大地线拱高限差约束的大地线内插方法

摘要

大地线拱高限差约束的大地线内插算法，包括以下步骤：读入大地线起点、终点的经纬度坐标；设定大地线拱高限差；计算未予标记的各折线段大地线长度、起始方位角及经纬度坐标中点；计算该点处的地球椭球面切点动径；计算未予标记的各折线段起点、终点的空间直角坐标；计算未予标记的各折线段大地线拱高；顺次判断未予标记的折线段大地线拱高与大地线拱高限差的关系，若折线段大地线拱高大于大地线拱高限差，则对折线段进行对分加密内插处理；若折线段大地线拱高小于等于大地线拱高限差，则该折线段不予内插处理。本发明通过大地线展绘精度与内插距离的自适应调整与匹配，实现了在指定精度阈值条件下任意长度大地线的精密展绘。

说明书

技术领域

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

背景技术

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

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

大地线的精确展绘，其最为关键的是构建精度阈值严密可控的大地线内插模型，实现大地线展绘精度与内插距离的自适应调整与匹配。如图4所示，P(B₁,L₁)、Q(B₂,L₂)表示地球椭球面上任意两点；w表示地球椭球面上由P、Q构成的大地线；w'表示P、Q及内插点相连构成的折线；h表示大地线w上任意点至直线w'的最大距离，称之为大地线拱高。受地球扁率f及P、Q处地球曲率的影响，大地线拱高h与大地线长度S(w)呈正相关关系。

理想情况下，当h→0，大地线w与折线w'重合。因此，对于大地线的精确展绘，可预先设定大地线展绘的精度阈值(大地线w拱高限差h_Ω)，判断折线w'中各折线段w_i'大地线拱高h_i与拱高限差h_Ω的数值大小关系，若折线段w_i'大地线拱高h_i大于大地线拱高限差h_Ω，则以折线段w_i'大地线长度S(w_i')的一半作为内插距离，对折线段w_i'进行对分加密内插处理；若折线段w_i'大地线拱高h_i小于等于长度大地线拱高限差h_Ω，则该折线段不予内插处理，并予以标记，避免重复运算。基于上述大地线的加密迭代内插处理过程，以折线段w_i'大地线拱高h_i为精度判断条件，可实现任意精度阈值(大地线w拱高限差h_Ω)条件下的大地线精密展绘。

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

发明内容

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

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

a、输入大地线起点P的经纬度坐标(B₁,L₁)和终点Q的经纬度坐标(B₂,L₂)；

b、确定起点P、终点Q在地球椭球面上构成的大地线w拱高限差h_Ω；

c、计算折线w'中未予标记的各折线段w_i'(i∈[0,n])(w'₀＝PT₁、w_i'＝T_iT_i+1、w'_n＝T_nQ)的大地线长度S(w_i')、起始方位角A_i及经纬度坐标中点起点P、终点Q及内插点T_i(i∈[1,n])相连构成折线w'，n表示大地线w上的内插点数，n＝0,1,2...；

d、依据折线w'中未予标记的各折线段w_i'(i∈[0,n])的经纬度坐标中点，计算折线w'中未予标记的折线段w_i'经纬度坐标中点处的地球椭球面切点动径r_i；

e、以球面代替地球椭球面，在替代球面上计算折线w'中未予标记的折线段w'大地线拱高h_i；

f、顺次判断折线w'中未予标记的折线段w_i'(i∈[0,n])大地线拱高h_i与大地线拱高限差h_Ω的关系，若h_i＞h_Ω，则内插出该折线段的大地线中点，遍历所有折线段后转至步骤c；若h_i≤h_Ω，则该折线段不予内插处理。

所述步骤b中，对于大地线展绘至纸质海图的情形，大地线w拱高限差h_Ω可采用纸质海图图上最小分辨距离d_chart对应的实地距离D_chart，海图比例尺分母为l，计算方法为：D_chart＝ld_chart/100。

所述步骤c中，已知折线段w_i'(i∈[0,n])起点T_i、终点T_i+1的经纬度坐标计算折线段w_i'的大地线长度S(w_i')、起始方位角A_i及经纬度坐标中点采用大地主题反解方法，即解算方程组：

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

所述步骤d中，依据折线段w_i'(i∈[0,n])的经纬度坐标中点，计算折线w'中未予标记的折线段w_i'经纬度坐标中点处的地球椭球面切点动径r_i：

计算折线段w_i'经纬度坐标中点令代入公式式中：a为椭球长半轴；e为椭球第一偏心率，计算折线段w_i'经纬度坐标中点处的地球椭球面切点动径r_i。

所述步骤e中，在替代球面上计算折线段w_i'大地线拱高h_i，大地坐标(B,L,H)至空间直角坐标(X,Y,Z)的转换计算方法为：

将折线段w_i'起点T_i、终点T_i+1的大地坐标分别代入上式计算各自空间直角坐标记折线段w_i'的长度则折线段w_i'大地线拱高h_i的计算方法为：

所述步骤f中，若h_i＞h_Ω，则内插出该折线段w_i'(i∈[0,n])大地线中点的计算方法为：以折线段w_i'起点T_i的经纬度坐标以及步骤c中计算所得的折线段w_i'起始方位角A_i、大地线长度S(w_i')的一半(作为内插距离)为已知条件，利用大地主题正解方法计算出折线段w_i'的大地线中点，大地主题正解方法为：已知折线段w_i'(i∈[0,n])起点T_i的经纬度坐标折线段w_i'大地线长度S(w_i')及其起始方位角A_i，计算折线段w_i'终点T_i+1的大地经纬度坐标采用正解方法，即解算方程组：

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

附图说明

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

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

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

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

图5是本发明局部范围内大地线拱高计算的简化图。

图6是本发明大地线拱高限差约束的大地线内插方法原理图。

图7是本发明大地线拱高限差约束的大地线内插方法流程图。

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

具体实施方式

本发明的地线拱高限差约束的大地线内插算法，采用计算机实现大地线的精密展绘，采用本发明方法实现大地线展绘精度与内插距离的自适应调整与匹配。假设大地线w起点、终点分别为点P和点Q；P、Q及内插点T_i(i∈[1,n])(n表示大地线w上的内插点数，n＝0,1,2...)相连构成的折线以w'表示；n初始值为0，表示P、Q相连的折线。

在给定大地线w拱高限差h_Ω的前提下，以球面代替地球椭球面计算折线w'中各折线段w_i'(i∈[0,n])(w'₀＝PT₁、w_i'＝T_iT_i+1、w'_n＝T_nQ)(n＝0时，w'₀＝PQ)的大地线拱高h_i，如图5所示，h_i可由直角ΔOT_iB间的三角关系解算获得。若h_i＞h_Ω，则以折线段w_i'大地线长度S(w_i')的一半对大地线进行对分内插，直至所有折线段w_i'的大地线拱高h_i≤h_Ω。包括以下步骤，如图7、图7所示：

步骤a，读入大地线w起点P的经纬度坐标(B₁,L₁)和终点Q的经纬度坐标(B₂,L₂)，P、Q在地球椭球面上构成的大地线以w表示；P、Q及内插点T_i(i∈[1,n])(n表示大地线w上的内插点数，n＝0,1,2...)相连构成的折线以w'表示；n初始值为0，表示P、Q相连的折线；

步骤b，设定大地线w拱高限差h_Ω(单位：m)。大地线w拱高限差h_Ω理论上可采用任意给定高度阈值(单位：m)。对于大地线展绘至纸质海图的情形，大地线w拱高限差h_Ω可采用纸质海图(海图比例尺分母为l)图上最小分辨距离d_chart(单位：cm)对应的实地距离D_chart(单位：m)，计算方法为：

D_chart＝ld_chart/100；

步骤c，利用大地测量学中的大地主题反解方法，计算折线w'中未予标记的折线段w_i'的大地线长度S(w_i')、起始方位角A_i及经纬度坐标中点已知折线段w_i'(i∈[0,n])(以w_i'＝T_iT_i+1为例)起点T_i、终点T_i+1的经纬度坐标计算折线段w_i'的大地线长度S(w_i')、起始方位角A_i及经纬度坐标中点的问题涉及大地测量学中的大地主题反解问题，即解算方程组：

式中：A_i'表示大地线反方位角。此外，已知折线段w_i'(i∈[0,n])(以w_i'＝T_iT_i+1为例)起点T_i的经纬度坐标折线段w_i'大地线长度S(w_i')及其起始方位角A_i，计算折线段w_i'终点T_i+1的大地经纬度坐标的问题归结为大地测量学中的大地主题正解问题，即解算方程组：

步骤d，令代入公式(式中：a为椭球长半轴；e为椭球第一偏心率)，计算折线w'中未予标记的折线段w_i'经纬度坐标中点处的地球椭球面切点动径r_i；以地球椭球面切点动径r_i作为替代地球椭球面的球面半径；如表1所示，若采用上述球面替代地球椭球面，则在200公里范围内，大地线长度S(w)与大圆弧长S_c(w)之间的差异δ_c小于0.5米；

表1球面代替地球椭球面产生的长度误差(基于克拉索夫斯基椭球参数)

注：表1及结论引自文献(华棠.海图数学基础[M].中国人民解放军海军司令部航海保证部,1985.)

步骤e，在替代球面上计算折线w'中未予标记的折线段w_i'大地线拱高h_i，涉及大地坐标(B,L,H)与空间直角坐标(X,Y,Z)的转换计算，大地坐标至空间直角坐标的转换(基于半径为地球椭球面切点动径r_i的球面)计算方法为：

将折线段w_i'(i∈[0,n])起点T_i、终点T_i+1的大地坐标分别代入上式计算各自空间直角坐标记折线段w_i'的长度

步骤f，半径为地球椭球面切点动径r_i的球面中，折线w'中未予标记的折线段w_i'大地线拱高h_i等价于图3中线段AB(AB＝OA-OB)的长度，从而线段OB长度的计算成为大地线拱高h_i计算的关键；在直角ΔOT_iB中，存在的三角关系；将已知条件OA＝OT_i＝r_i、代入上述公式，则大地线拱高

步骤g，顺次判断折线w'中未予标记的折线段w_i'大地线拱高h_i与大地线拱高限差h_Ω的关系，若h_i＞h_Ω，则内插出该折线段的大地线中点，遍历所有折线段后转至步骤c；若h_i≤h_Ω，则该折线段不予内插处理，并予以标记。

算例分析：

说明：以下算例中“拱高误差”参数为相邻内插点(包含大地线起点、终点)间小段大地线拱高的极大值统计。

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

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

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

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

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

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

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

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

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

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

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