常用投影大地線的高效展繪及Mathematica實現(xiàn)

周東權(quán)，劉 敏，魏 沖，邊少鋒，張思遠(yuǎn)

常用投影大地線的高效展繪及Mathematica實現(xiàn)

周東權(quán)1，劉 敏2，魏 沖3，邊少鋒1，張思遠(yuǎn)1

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

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

第三扁率；Mathematica；大地線展繪；等間距內(nèi)插；自適應(yīng)控制模型

0 引言

圖1 大地主題解算示意圖

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

大地線由于本質(zhì)是一條空間曲面曲線，與平面曲線的繪制不同，函數(shù)方程復(fù)雜，需要綜合解算的方法來實現(xiàn)展繪[11-12]。“以直代曲”是大地線展繪的常見方法，即通過大地解算，得到大地線上一系列密集的點位，通過點與點的連接實現(xiàn)大地線的展繪。為此，相關(guān)學(xué)者在該思路的基礎(chǔ)上，提出了大地線的內(nèi)插方法，根據(jù)地球曲率計算最大插值間隔，建立絕對精度閾值約束的內(nèi)插算法，提高了大地線的絕對精度，卻由于過多的內(nèi)插點降低了大地線的展繪效率[13]。針對大地線展繪算法中內(nèi)插參數(shù)與精度自適應(yīng)匹配的問題，學(xué)者建立了大地線展繪長度誤差與拱高誤差的定量評估模型，提出了一種有效控制閾值的快速展繪算法。由于現(xiàn)有大地線展繪方法大多建立于大地線方程在投影中的精確表達(dá)，繼而建立的矢量數(shù)據(jù)模型[14]，形式單一，往往建立在墨卡托投影等角和固定比例尺的條件下，無法滿足許多場景的應(yīng)用。

Mathematica是基礎(chǔ)研究、應(yīng)用基礎(chǔ)研究以及工程技術(shù)領(lǐng)域流行的計算機代數(shù)系統(tǒng)[15-16]，在參考橢球數(shù)學(xué)計算分析、微分幾何分析、地圖投影和數(shù)據(jù)歸上有著系統(tǒng)的強大的數(shù)學(xué)分析能力、符號推導(dǎo)能力和可視化能力，具有便捷快速地建立任意投影的強大能力。基于此，本文在改化公式的基礎(chǔ)上，利用Mathematica建立大地線展繪算法，實現(xiàn)大地線在常用投影的高效展繪，挖掘可視化技術(shù)在大地測量學(xué)的應(yīng)用價值。

1 算法的基本原理

圖2 算法流程圖

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

2 大地主題解算

2.1 大地線長S計算

球面歸化緯度與球面大圓弧長的數(shù)學(xué)關(guān)系式如式（3）所示

則

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

式中：

2.2 大地方位角計算及修正

2.3 基于Mathematica的語法設(shè)計

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

2.4 精度自適應(yīng)控制模型

3 結(jié)果分析

3.1 大地線在空間直角坐標(biāo)下的展繪

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

表1 大地線展繪實驗數(shù)據(jù)

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

表2 大地線線展繪實驗結(jié)果

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

圖3 精度自適應(yīng)誤差計算結(jié)果

圖4 內(nèi)插點空間三維散點圖

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

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

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

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

4 結(jié)語

本文在利用Mathematica強大的數(shù)學(xué)分析功能探究大地問題正反解算法的基礎(chǔ)上，通過將等間距內(nèi)插展繪的應(yīng)用擴展至展繪精度及內(nèi)插間距自動調(diào)節(jié)的自適應(yīng)模型，同時利用Mathematica強大的繪圖能力，實現(xiàn)了常用投影的大地線高效展繪，得出了以下結(jié)論。

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

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

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

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

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

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

[4] 鄭天瑞. 利用Mathematica求解正軸圓錐投影的標(biāo)準(zhǔn)緯度[J]. 現(xiàn)代導(dǎo)航，2022，13（5）：383-386.

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

[6] 紀(jì)兵，邊少鋒. 大地主題問題的非迭代新解[J]. 測繪學(xué)報，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] 計宇陽，金立新，丁佳波. 復(fù)數(shù)高斯投影的大地線曲率與代曲直距[J]. 海洋測繪，2021，41（6）：31-35.

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

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

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

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

[16] 趙自豪，李鵬慧. 最小徑集求解方法及Mathematica實現(xiàn)[J]. 湖南科技大學(xué)學(xué)報（自然科學(xué)版），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。

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

國家自然科學(xué)基金項目（42074010）