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

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

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

（1 中國(guó)地質(zhì)大學(xué)（武漢）地質(zhì)探測(cè)與評(píng)估教育部重點(diǎn)實(shí)驗(yàn)室，武漢 430074；2 中國(guó)人民解放軍91001部隊(duì)，北京 100830；3 中國(guó)人民解放軍92823部隊(duì)，三亞 572000）

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

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

0 引言

圖1 大地主題解算示意圖

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

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

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

1 算法的基本原理

圖2 算法流程圖

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

2 大地主題解算

2.1 大地線(xiàn)長(zhǎng)S計(jì)算

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

則

根據(jù)牛頓二項(xiàng)式展開(kāi)定理的推廣，可進(jìn)一步展開(kāi)，如式（7）所示

式中：

2.2 大地方位角計(jì)算及修正

2.3 基于Mathematica的語(yǔ)法設(shè)計(jì)

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

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

3 結(jié)果分析

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

根據(jù)上述算法描述、內(nèi)插方法和繪制要求，通過(guò)Mathematica12.0作為測(cè)試環(huán)境，分別對(duì)不同實(shí)驗(yàn)數(shù)據(jù)進(jìn)行了大地線(xiàn)展繪實(shí)驗(yàn)，實(shí)驗(yàn)中所使用的參考橢球?yàn)榭死鞣蛩够鶛E球，實(shí)驗(yàn)數(shù)據(jù)如表1所示。

表1 大地線(xiàn)展繪實(shí)驗(yàn)數(shù)據(jù)

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

表2 大地線(xiàn)線(xiàn)展繪實(shí)驗(yàn)結(jié)果

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

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

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

3.2 大地線(xiàn)在常用投影下的展繪

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

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

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

4 結(jié)語(yǔ)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[16] 趙自豪，李鵬慧. 最小徑集求解方法及Mathematica實(shí)現(xiàn)[J]. 湖南科技大學(xué)學(xué)報(bào)（自然科學(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—），廣東化州人，碩士研究生，主要研究方向?yàn)闄E球大地測(cè)量計(jì)算機(jī)代數(shù)分析。

國(guó)家自然科學(xué)基金項(xiàng)目（42074010）