常用投影大地線的高效展繪及Mathematica實現

周東權，劉 敏，魏 沖，邊少鋒，張思遠

常用投影大地線的高效展繪及Mathematica實現

周東權1，劉 敏2，魏 沖3，邊少鋒1，張思遠1

（1 中國地質大學（武漢）地質探測與評估教育部重點實驗室，武漢 430074；2 中國人民解放軍91001部隊，北京 100830；3 中國人民解放軍92823部隊，三亞 572000）

針對大地線方程的復雜性及高效展繪難以實現的問題，以第三扁率改化大地問題正反解的計算公式，根據法截線方位角與大地方位角的數學關系提高大地方位角計算精度，并提出了等間距內插的大地線展繪方法；利用Mathematica建立了精度自適應控制模型，提高了等間距內插的效率和精度；利用Mathematica強大的繪圖功能實現了不同距離和應用場景下不同地圖投影的大地線展繪，挖掘了大地線展繪的應用價值。實驗結果表明：該方法在精度可控的條件下能實現對任意長度大地線的高效率展繪，不同地圖投影的大地線展繪在不同的應用場景下具有獨特的意義。

第三扁率；Mathematica；大地線展繪；等間距內插；自適應控制模型

0 引言

圖1 大地主題解算示意圖

隨著計算機代數系統的不斷發展，大地主題解算的精度和效率也得到了提升[3-4]。與傳統的研究不同的是，計算機代數系統對大地解算公式進行了系統的革新，并能創新性地討論貝塞爾大地解算等傳統思路的新方法，簡化了大地主題正反解的公式，并提高效率[5-6]。同時符號計算的出現也能夠針對過往大地問題研究中缺乏考慮的問題，在建立新的參考橢球的基礎上，提出新的大地問題解算方法，如顧及高程的大地問題解算方法[7-8]。大地問題解算算法的不斷完善，為大地線的高精度展繪提供了算法基礎[9-10]。

大地線由于本質是一條空間曲面曲線，與平面曲線的繪制不同，函數方程復雜，需要綜合解算的方法來實現展繪[11-12]?！耙灾贝笔谴蟮鼐€展繪的常見方法，即通過大地解算，得到大地線上一系列密集的點位，通過點與點的連接實現大地線的展繪。為此，相關學者在該思路的基礎上，提出了大地線的內插方法，根據地球曲率計算最大插值間隔，建立絕對精度閾值約束的內插算法，提高了大地線的絕對精度，卻由于過多的內插點降低了大地線的展繪效率[13]。針對大地線展繪算法中內插參數與精度自適應匹配的問題，學者建立了大地線展繪長度誤差與拱高誤差的定量評估模型，提出了一種有效控制閾值的快速展繪算法。由于現有大地線展繪方法大多建立于大地線方程在投影中的精確表達，繼而建立的矢量數據模型[14]，形式單一，往往建立在墨卡托投影等角和固定比例尺的條件下，無法滿足許多場景的應用。

Mathematica是基礎研究、應用基礎研究以及工程技術領域流行的計算機代數系統[15-16]，在參考橢球數學計算分析、微分幾何分析、地圖投影和數據歸上有著系統的強大的數學分析能力、符號推導能力和可視化能力，具有便捷快速地建立任意投影的強大能力?；诖?，本文在改化公式的基礎上，利用Mathematica建立大地線展繪算法，實現大地線在常用投影的高效展繪，挖掘可視化技術在大地測量學的應用價值。

1 算法的基本原理

圖2 算法流程圖

在計算得出各個點的大地坐標后，同時可以得出各點的空間坐標，利用Mathematica的強大內核，完成在不同投影的展繪。

2 大地主題解算

2.1 大地線長S計算

球面歸化緯度與球面大圓弧長的數學關系式如式（3）所示

則

根據牛頓二項式展開定理的推廣，可進一步展開，如式（7）所示

式中：

2.2 大地方位角計算及修正

2.3 基于Mathematica的語法設計

Mathematica有強大的符號計算功能和制圖能力，可以利用簡潔的語言設計大地主題正反解計算的代數分析，如在大地方位角計算及修正等復雜的公式中，Mathematica能利用簡單的語句設計，實現兩點間大地坐標相關函數的計算。在計算過程中，與Matlab、Python等數值計算分析編程語言相比，Mathematica只需要最簡單的函數自定義，即可將代數計算中復雜代數關系進行關聯，從而實現復雜的代數推導及結果輸出，涉及多元矩陣等計算時，代數系統能在操作頁面建立并識別矩陣，與數值計算工具相比，更加簡單便捷。其計算精度之高對于大地線的推導、改化和模型構建都提供了重大幫助。同時，Mathematica還能利用簡單的語言建立復雜的地圖投影，并實現自定義設計地圖投影。對于Matlab、Python等數值分析語言，自定義地圖投影需要復雜的嵌套函數以及程序設計，而在代數系統上，能利用內置的強大的數學工具及投影函數，對常用投影如Mercator投影，能夠一句代碼實現地圖投影的建立，自定義參數的任意設定，同時能過通過代數推導，并利用Mathematica強大的繪圖能力，實現自定義地圖投影的建立，是大地線代數推導及展繪實現的堅實基礎。

2.4 精度自適應控制模型

3 結果分析

3.1 大地線在空間直角坐標下的展繪

根據上述算法描述、內插方法和繪制要求，通過Mathematica12.0作為測試環境，分別對不同實驗數據進行了大地線展繪實驗，實驗中所使用的參考橢球為克拉索夫斯基橢球，實驗數據如表1所示。

表1 大地線展繪實驗數據

大地方位角可從上述大地方位角計算及修正方法得出，即大地經緯度的變化會帶動平行圈的變化，從而帶動大地方位角的變化，因此大地線的展繪精度與兩點大地經緯度以及大地方位角有關。三個算例的實驗結果如表2所示，從中可以獲知，算例2中大地方位角大，同時由于在大地線的展繪過程中，隨著緯度的增加，大地方位角的減少幅度也在不斷增加，會對展繪精度造成一定的影響。 算例3中，由于大地線的展繪從南半球到北半球展繪的過程中經過，大地方位角也由逐步變小到逐步變大，因此在展繪的過程中精度也在波動變化。

表2 大地線線展繪實驗結果

此次展繪的實驗結果還表明精度自適應控制模型根據大地線的展繪情況，選擇了滿足精度的內插間距，且滿足了大地線繪制高效率的要求，同時可以得知，大地線展繪過程中，內插點數量與大地線長無關?？偟膩碚f，大地線的內插過程中大地方位角變化幅度越小，其展繪精度越高，三個算例的誤差計算結果如圖3所示，可以看出，算例1兩點位于中低緯度，因此在51次的自適應調整中，精度都能保持較好，而算例2由于是從低緯度向高緯度展繪，而且終點緯度較高，精度較低，且精度變化較大。算例3由于是從南半球向北半球進行展繪，自適應模型隨著調整次數的增多，即內插距離的逐步增大，精度得到了較好地調整，基本能滿足展繪要求。理論條件下，算例2和算例3由終點向起點進行展繪能夠得到精度更高的展繪結果，盡管如此，在精度自適應控制模型的調整下，都在滿足條件的情況下實現了大地線高效率的展繪，這在未來大地問題的解算中，以及大地問題可視化算法進一步優化上，都具有一定的作用。

圖3 精度自適應誤差計算結果

圖4 內插點空間三維散點圖

3.2 大地線在常用投影下的展繪

大地線是地球橢球體表面一條空間曲面曲線，在辨識和使用時往往要將其投影在地圖上，其是一個曲面向平面投影的過程，不同的數學映射關系會構成性質不同的投影，按變形性質主要分為等距投影、等角投影、等積投影和任意投影。由于應用場景不同，往往采用不同的地圖投影，地圖投影的不同會影響大地線的展繪及使用。本文在不同場景的地圖投影下對大地線進行了展繪，探究其特性。

Mercator投影是正軸等角圓柱投影，其沒有角度變形，恒向線在其投影上表現為一條直線。三個算例在Mercator投影上的展繪如圖5所示，藍色為大地線，綠色為恒向線。通過三個算例在Mercator投影上的繪制可以得知：Mercator投影在兩極變形大，在低緯度地區變形小，算例1大地線距離較短，大地方位角變化也不大，同時位于較低緯度，經差也較小，其投影變形小，因此在地圖上近似為恒向線。而算例2和算例3距離遠，大地方位角變化大，算例3大地線還跨越南北半球，投影變形大，可以看出與恒向線具有較大差別?？偟膩碚f，大地線在地圖上展繪的表現形式受到地圖投影變形以及大地方位角的影響，大地線與恒向線的關系和性質在船舶大洋航行、精確制導等場景應用中具有一定的意義。

大地線是地球橢球體上的最短曲線，其在地圖平面上的視覺表現會因為不同的地圖投影而發生變化。常用的地圖投影除了Mercator投影還有Gauss投影和Lambert投影。Gauss投影雖然較大程度地降低了地圖投影產生的變形，但是并沒有保持真實的方向，因此圖6可以看出恒向線在該投影上表現為一條曲線，大地線則因為投影變形較小更加滿足其為最短路徑的視覺定義。Lambert投影是擬定的正形圓錐投影，常用的包括：Lambert Conic Conformal投影和Lambert Azimuthal EqualArea投影，前者適合中緯度東西方向的地圖繪制，變形較小，因此可以看出圖7中算例2的大地線繪制中，能夠較好地滿足大地線為最短曲線的視覺效果?？偠灾?，在短距離的大地線展繪中，不同投影大地線與恒向線的差別不大，而長距離的大地線展繪，不同投影性質下的地圖投影會影響其展繪的視覺效果，同時地圖投影參數如Central Scale Factor等的設置也會進一步影響其視覺效果，Mathematica強大的繪圖能力進一步提高了大地線的應用價值。

4 結語

本文在利用Mathematica強大的數學分析功能探究大地問題正反解算法的基礎上，通過將等間距內插展繪的應用擴展至展繪精度及內插間距自動調節的自適應模型，同時利用Mathematica強大的繪圖能力，實現了常用投影的大地線高效展繪，得出了以下結論。

1）本文利用第三扁率對大地問題正反解的算法進行改化，提高了大地問題正反解的計算精度和計算效率，同時利用法截線方位角與大地方位角的數學關系，通過修正項的引入降低了大地方位角的計算誤差。

2）針對大地線是一條空間曲面曲線的特殊性，提出了基于等間距的內插方法，通過內插點線段連接的方法提出了大地線展繪的基本思路，同時利用Mathematica設計了內插間距的自適應模型，在提高精度的情況下同時提高了展繪效率。

3）結合大地線可視化的具體應用場景，利用Mathematica強大的制圖功能展示了大地線與恒向線的關系，并根據不同距離和不同應用場景對不同投影下的大地線展繪進行了初步探究，提高了其應用價值，對后續具體科學問題的探究打下了基礎。其使得展繪算法的應用不再局限于墨卡托投影平面，具有一定的科學性。

[1] 邊少鋒，李厚樸. 大地測量計算機代數分析[M]. 北京：科學出版社，2018：9-23+38-42.

[2] 董箭，李彬彬，彭認燦，等. 顧及拱高誤差的墨卡托大地線快速展繪算法[J]. 測繪科學，2020，45（9）：43-51.

[3] 張曉東，涂玲，劉寶. 基于Mathematica的控制系統根軌跡探究性實驗[J]. 實驗技術與管理，2022，39（12）：163-168.

[4] 鄭天瑞. 利用Mathematica求解正軸圓錐投影的標準緯度[J]. 現代導航，2022，13（5）：383-386.

[5] 吳祖新，鄭中義. 大地線航法在智能船舶上的應用[J]. 艦船科學技術，2023，45（1）：180-185.

[6] 紀兵，邊少鋒. 大地主題問題的非迭代新解[J]. 測繪學報，2007（3）：269-273.

[7] CHANG S，JI B，WU M，et al. Evaluation of Height Correction on Loran Signal's Groundwave Transmission Delay Model[J]. IEEE Antennas and Wireless Propagation Letters，2022：1-5.

[8] SEIF M R，GHALENOEI E. Numerical determination of the geodesic curves: the solution of a two-point boundary value problem[J]. Earth Observation and Geomatics Engineering，2018，2（1）：26-35.

[9] PE?KA P，MIKE? J，RYPAROVá L. Almost Geodesic Curves as Intersections of n-Dimensional Spheres[J]. Lobachevskii Journal of Mathematics，2022，43（3）：687-690.

[10] BASELGA S，MARTíNEZ-LLARIO J C. Intersection and point-to-line solutions for geodesics on the ellipsoid[J]. Studia Geophysica Et Geodaetica，2017：1-11.

[11] 計宇陽，金立新，丁佳波. 復數高斯投影的大地線曲率與代曲直距[J]. 海洋測繪，2021，41（6）：31-35.

[12] 周江華，苗育紅，成文生，等. 貝塞爾大地反解問題的高效率算法[J]. 測繪學報，2002（2）：108-111.

[13] 唐紅濤，王微，楊永崇，等. 橢球面上繪制大地線的算法[J]. 測繪科學，2015，40（4）：7-10.

[14] 劉文超，卞鴻巍，王榮穎，等. 大橢圓航線設計的空間矢量方法[J]. 測繪學報，2015，44（7）：741-746.

[15] 王海鷹，閆浩文，梁曉鶴. 經典范疇認知下地圖的局限性[J]. 海洋測繪，2023，43（1）：78-82.

[16] 趙自豪，李鵬慧. 最小徑集求解方法及Mathematica實現[J]. 湖南科技大學學報（自然科學版），2022，37（2）：27-32.

Efficient Mapping of Common Projected Geodetic Lines and Mathematica Implementation

ZHOU Dongquan, LIU Min, WEI Chong, BIAN Shaofeng, ZHANG Siyuan

In view of the complexity of geodetic line equation and the difficulty of efficient plotting, the formula of the forward and backward solutions of geodetic problem with third flatteningis changed, the accuracy of geodetic azimuth calculation according to the mathematical relationship between normal transversal azimuth and geodetic azimuth is improved, and an equidistant interpolation geodetic line plotting method is proposed. A precision adaptive control model was established by using Mathematica to improve the efficiency and precision of equal spacing interpolation. Using the powerful mapping function of Mathematica, the geodetic line plotting of different map projections under different distances and application scenarios is realized, and the application value of geodetic line plotting is explored. Experimental results show that the method can achieve high efficiency geodetic line drawing of arbitrary length under the condition of controlled accuracy. Geodetic line plotting of different map projections has unique significance in different application scenarios.

Third Flattening; Mathematica; Geodetic Line Plotting; Equidistant Interpolation; Adaptive Control Model

P226

A

1674-7976-(2023)-06-422-09

2023-06-20。

周東權（1999.08—），廣東化州人，碩士研究生，主要研究方向為橢球大地測量計算機代數分析。

國家自然科學基金項目（42074010）