大地线长度限差约束的大地线内插方法与流程

本发明涉及海洋测绘技术领域。

背景技术：

由大地测量学知识可知，任意给定的地球椭球面上的两点，必然存在唯一一条大地线。利用这一特性可以在给定的地球椭球面上实现对直线领海基线等各种“直线”的准确表达，因而在海洋划界中将各种“直线”视为大地线已被越来越多的海洋国家所接受。如图1-3所示，由于大多数大地线通常在各种投影平面上都表现为十分复杂的曲线，且不容易按严格的数学模型实现其在海图平面上的表示，因此必须设法实现其精确展绘的问题。

为解决上述问题，国内外大地测量学家都进行了深入而详尽的研究并提出了各自的解决方案。当前普遍采用的是绘制任意曲线时采用的“以直代曲”内插方法，如图4所示，其关键步骤在于获取最大内插距离，对于该参数的求解可依据地球曲率所确定出的由直线(折线段)替代大地线所对应的地球椭球面上的长度获得。然而大地线各位置处的地球曲率变化不一，且长距离大地线各位置处的地球曲率相差较大，最大内插距离的获取涉及较多的大地线微分方程的解算，方法相对复杂且不易实现。此外，该类方法尽管可以保证任意小段直线(折线段)逼近大地线的绝对精度，但对于大地线展绘的整体累积误差却未予探讨，且无法实现任意精度阈值条件下的大地线精确展绘。

大地线的精确展绘，其最为关键的是构建精度阈值严密可控的大地线内插模型，实现大地线展绘精度与内插距离的自适应调整与匹配。如图4所示，p(b1,l1)、q(b2,l2)表示地球椭球面上任意两点；w表示地球椭球面上由p、q构成的大地线；w'表示p、q及内插点相连构成的折线。受地球扁率f及p、q处地球曲率的影响，大地线w的长度s(w)(大地线长度)恒大于折线w'的长度d(w')(d(w')＝∑d(wi')，wi'表示p、q及内插点相连所构成的任意单根折线段)，且两者之间的长度差值δ(s(w)-d(w'))与大地线w长度s(w)呈正相关关系。

理想情况下，当δ→0，大地线w与折线w'重合。因此，对于大地线的精确展绘，可预先设定大地线展绘的精度阈值(大地线w长度限差δω)，判断长度差值δ与长度限差δω的数值大小关系，若长度差值δ大于长度限差δω，则分别计算各折线段wi'的大地线距离s(wi')并以作为局部内插距离，对各折线段wi'进行对分加密内插处理；若长度差值δ小于等于长度限差δω，则停止内插。基于上述大地线的加密迭代内插处理过程，以大地线长度与大地线起点、终点及内插点所构折线的长度差值为精度判断条件，可实现任意精度阈值(大地线w长度限差δω)条件下的大地线精密展绘。

在长期的实践中，对大地线的精确展绘，主要依赖于大地线最大内插距离的计算及小段直线(折线段)逼近大地线的局部精度评价，尚未有公开的文献资料能够定量地评估大地线展绘的整体精度指标，并实现精度阈值严密可控的大地线加密内插。

技术实现要素：

为了克服传统定性分析方法存在的上述问题，本发明提供了大地线长度限差约束的大地线内插算法。

本发明为实现上述目的所采用的技术方案是：大地线长度限差约束的大地线内插方法，包括以下步骤：

a、读入大地线w起点p的经纬度坐标(b1,l1)和终点q的经纬度坐标(b2,l2)；

b、设定大地线w长度限差δω；

c、根据p、q的经纬度坐标(b1,l1)和(b2,l2)，计算大地线起点为p(b1,l1)、终点为q(b2,l2)的大地线长度s(w)及其起始方位角a12；

d、计算内插距离以内插距离对大地线w进行无约束条件的等间距内插处理，p、q及内插点ti(i∈[1,2ⁿ-1])相连构成的折线以w'表示，n表示大地线w的内插次数，n＝1,2,3...；

f、计算大地线w长度s(w)与p、q及内插点ti(i∈[1,2ⁿ-1])所构成的折线w'长度d(w')的差值δi；

f、判断差值δi与限差δω的关系，若δi＞δω，则转至步骤d；若δi≤δω，则停止内插。

所述步骤b中，对于大地线展绘至纸质海图的情形，长度限差δω可采用纸质海图图上最小分辨距离dchart对应的实地距离代替，即：δω＝ldchart/100，海图比例尺分母为l。

所述步骤c中，利用大地测量学中的大地主题反解方法，计算大地线起点为p(b1,l1)、终点为q(b2,l2)的大地线长度s(w)及其起始方位角a12，计算方法如下：解算方程组：

式中：a21表示大地线反方位角。

所述步骤d中，以内插距离对大地线w进行无约束条件的等间距内插处理的方法为：

①对于大地线w未进行内插处理的情形(w'＝{p,q})：以内插距离起始方位角a12、p点的经纬度坐标(b1,l1)为已知条件，利用大地主题正解方法计算出该小段大地线末端的经纬度坐标(b'2,l'2)，计算大地线w终点q的大地经纬度坐标(b'2,l'2)的正解方法如下：解算方程组：

以(b'2,l'2)作为新的起始经纬度坐标(b1',l'1)与q点的经纬度坐标(b2,l2)配对，利用大地主题反解方法计算新的起始方位角a1'2，以新的起始方位角a1'2、经纬度坐标(b1',l'1)及内插距离δs(n)为已知条件，利用大地主题正解方法计算出下一小段大地线末端的经纬度坐标(b”2,l”2)；依此类推，即每次以前一次求得的小段大地线末端经纬度坐标及起始方位角计算下一小段大地线的起始方位角及末端经纬度坐标，直至点位内插完毕；此时，大地线w可由p、q及内插点ti所构成折线w'的有序点集表示，即：

②对于大地线w已进行内插处理的情形遍历折线w'的有序点集，以w'上任意连续节点ti、ti+1(i∈[1,2ⁿ-2])的经纬度坐标和为已知条件，利用大地主题反解方法计算各折线段titi+1(wi')的大地线距离s(wi')及起始方位角ai；以起始方位角ai、节点ti的经纬度坐标及内插距离为已知条件，利用大地主题正解方法计算出各折线段titi+1(wi')的大地线中点(对分内插点)；折线段pt1、大地线中点的计算与折线段titi+1(wi')大地线中点的计算原理一致。

所述步骤e中，大地坐标至空间直角坐标的转换计算方法为：

式中：为卯酉圈曲率半径；a为椭球长半轴；e为椭球第一偏心率，大地坐标(b,l,h)、空间直角坐标(x,y,z)；将p、q及内插点ti(i∈[1,2ⁿ-1])的大地坐标(b1,l1,0)、(b2,l2,0)和分别代入上式计算各自空间直角坐标(x1,y1,z1)、(x2,y2,z2)和则p、q及内插点ti(i∈[1,2ⁿ-1])所构成的折线w'长度d(w')中各项的计算方法为：

大地线w长度s(w)与p、q及内插点ti(i∈[1,2ⁿ-1])所构成的折线w'长度d(w')的差值δi为：

本发明的大地线长度限差约束的大地线内插方法，通过大地线展绘精度与内插距离的自适应调整与匹配，实现了在指定精度阈值条件下任意长度大地线的精密展绘。

附图说明

图1是大地线在墨卡托投影平面上表现形式图。

图2是大地线在墨卡托投影平面上表现形式图。

图3是大地线在墨卡托投影平面上表现形式图。

图4是大地线“以直代曲”的展绘原理图。

图5是本发明大地线长度限差约束的大地线内插算法的主流程图。

图6是本发明无约束条件的等间距大地线(对分)内插算法原理图。

具体实施方式

本发明的大地线长度限差约束的大地线内插方法，采用计算机实现大地线的精密展绘，采用本方法实现大地线展绘精度与内插距离的自适应调整与匹配。

假设大地线w起点、终点分别为点p和点q，在给定大地线w长度限差δω的前提下，每次以局部折线段wi'大地线距离s(wi')的一半对折线段wi'进行对分加密内插处理；计算大地线w长度s(w)与p、q及内插点所构成的折线长度∑d(wi')的差值δi；判断差值δi与限差δω的关系，当δi≤δω时停止大地线对分内插。如图5所示，包括以下步骤：

步骤a，读入大地线w起点p的经纬度坐标(b1,l1)和终点q的经纬度坐标(b2,l2)；

步骤b，设定大地线w长度限差δω(单位：m)。大地线w长度限差δω理论上可采用任意给定长度阈值(单位：m)，对于大地线展绘至纸质海图的情形，δω可采用纸质海图(海图比例尺分母为l)图上最小分辨距离dchart(单位：cm)对应的实地距离代替，即：δω＝ldchart/100；

步骤c，利用大地测量学中的大地主题反解方法，计算大地线起点为p(b1,l1)、终点为q(b2,l2)的大地线长度s(w)及其起始方位角a12，已知大地线w起点p、终点q的经纬度坐标(b1,l1)、(b2,l2)，计算大地线w长度s(w)及其起始方位角a12的问题归结为大地测量学中的大地主题反解问题，即解算方程组：

式中：a21表示大地线反方位角。此外，已知大地线w起点p的经纬度坐标(b1,l1)、大地线w长度s(w)及其起始方位角a12，计算大地线w终点q的大地经纬度坐标(b2,l2)的问题归结为大地测量学中的大地主题正解问题，即解算方程组：

步骤d，对于大地线w未进行内插处理的情形(w'＝{p,q})：以内插距离(n表示大地线w的内插次数，n＝1,2,3...)、起始方位角a12、p点的经纬度坐标(b1,l1)为已知条件，利用大地主题正解方法计算出该小段大地线末端的经纬度坐标(b'2,l'2)，并转至步骤f；

步骤e，对于大地线w已进行内插处理的情形遍历折线w'的有序点集，以w'上任意连续节点ti、ti+1(i∈[1,2ⁿ-2])的经纬度坐标和为已知条件，利用大地主题反解方法计算各折线段titi+1(wi')的大地线距离s(wi')及起始方位角ai；以起始方位角ai、节点ti的经纬度坐标及内插距离为已知条件，利用大地主题正解方法计算出各折线段titi+1(wi')的大地线中点(对分内插点)；折线段pt1、大地线中点的计算与折线段titi+1(wi')大地线中点的计算原理一致，并转至步骤i。

步骤f，以步骤d中计算所得的(b'2,l'2)作为新的起始经纬度坐标(b1',l'1)与q点的经纬度坐标(b2,l2)配对，利用大地测量学中的大地主题反解方法计算新的起始方位角a1'2；

步骤g，以新的起始方位角a1'2、经纬度坐标(b1',l'1)及内插距离δs(n)为已知条件，利用大地测量学中的大地主题正解方法计算出下一小段大地线末端的经纬度坐标(b”2,l”2)；

步骤h，依此类推，重复步骤d至步骤f，即每次以前一次求得的小段大地线末端经纬度坐标及起始方位角计算下一小段大地线的起始方位角及末端经纬度坐标，直至点位内插完毕(如图6所示)；此时，大地线w可由p、q及内插点ti所构成折线w'的有序点集表示，即：

步骤i，将p、q及内插点ti(i∈[1,2ⁿ-1])的大地坐标(b1,l1,0)、(b2,l2,0)和转换为空间直角坐标(x1,y1,z1)、(x2,y2,z2)和转换公式如下：

式中：为卯酉圈曲率半径；a为椭球长半轴；e为椭球第一偏心率。

步骤j，计算p、q及内插点ti(i∈[1,2ⁿ-1])所构成的折线w'长度d(w')各项的计算方法为：

步骤k，计算大地线w长度s(w)与p、q及内插点ti(i∈[1,2ⁿ-1])所构成的折线w'长度d(w')的差值

步骤l，判断差值δi与限差δω的关系，若δi＞δω，则转至步骤e；若δi≤δω，则停止内插。

采用该方法的内插效果及内插误差如下表所示：

表1大地线长度限差约束的大地线内插算法原理

算例分析：

大地线展绘算例(一)引自文献(黄继文.椭球大地测量学[m].郑州:解放军测绘学院,1991.)，文献中算例的计算结论如表2所示，本发明中算例计算结论及与文献的计算偏差如表3所示。

表2大地线展绘算例(一)计算结论

表3大地线展绘算例(一)计算偏差统计

大地线展绘算例(二)引自文献(华棠.海图数学基础[m].中国人民解放军海军司令部航海保证部,1985.)，文献中算例的计算结论如表4所示，本发明中算例计算结论及与文献的计算偏差如表5所示。

表4大地线展绘算例(二)计算结论

表5大地线展绘算例(二)计算偏差统计

大地线展绘算例(三)引自文献(黄继文.椭球大地测量学[m].郑州:解放军测绘学院,1991.)，文献中算例的计算结论如表6所示，本发明中算例计算结论及与文献的计算偏差如表7所示。

表6大地线展绘算例(三)计算结论

表7大地线展绘算例(三)计算偏差统计

大地线展绘算例(一)、(二)、(三)中的大地线长度分别为28230.936m、225310.09m和15000000.1m。对比分析不难看出：未进行内插处理的大地线展绘误差随着大地线长度的增大而急剧增大；传统的“以直代曲”的固定距离内插方法(无约束条件的等间距(对分)内插方法)，在大地线长度较小(算例(一)、(二))时能保证较高的展绘精度，但对于长距离大地线(算例(三))的展绘，则会出现大地线展绘整体累积误差无法有效控制的情况；本发明中的大地线长度限差约束的大地线内插算法通过大地线展绘精度与内插距离的自适应调整与匹配，实现了在指定精度阈值条件下任意长度大地线的精密展绘。

本发明是通过实施例进行描述的，本领域技术人员知悉，在不脱离本发明的精神和范围的情况下，可以对这些特征和实施例进行各种改变或等效替换。另外，在本发明的教导下，可以对这些特征和实施例进行修改以适应具体的情况及材料而不会脱离本发明的精神和范围。因此，本发明不受此处所公开的具体实施例的限制，所有落入本申请的权利要求范围内的实施例都属于本发明的保护范围。