向新桃，秦 堯，湯瑾璟

（上海船舶研究設(shè)計院，上海201203）

0 前 言

隨著高技術(shù)船舶不斷發(fā)展，智能船已成為各界關(guān)注的重點。 《全球港航信息化發(fā)展報告（2018版）》［1］指出，中國、日本、韓國以及不少歐洲國家和地區(qū)正在加快智能船舶的研發(fā)，全球各大船級社也積極參與其中并紛紛頒布了相關(guān)規(guī)范。 中國船級社將智能船功能分為智能航行、智能船體、智能機艙、智能能效管理、智能貨物管理和智能集成平臺六大類［2］。 其中智能航行功能要求對航線和航行策略進行優(yōu)化，達到減少油耗等目的。 這就要求智能系統(tǒng)能夠精確計算地球表面各幾何要素，實現(xiàn)恒向線、大圓弧、大地線等主題計算等。

在地球橢球模型中，大地線可以表征地球表面兩點間最短距離，其主要幾何特征為曲線各點主法線與橢球面法線重合［3］。顯然，子午線和赤道上的弧段都是大地線，此時大地線各點位于同一平面內(nèi)，為二維平面曲線。 一般情況下，大地線各點主法線并不在同一平面，此時大地線為三維空間曲線。 由于大地線方向、距離等幾何要素的計算過程較復(fù)雜，在經(jīng)典航海學(xué)中應(yīng)用不多。 近年來隨著船舶智能化趨勢的確立，智能船對精確感知和精確計算提出了更高要求，加之智能系統(tǒng)的算力也大幅提高，大地線主題計算在智能船航線管理和數(shù)據(jù)分析等功能模塊中發(fā)揮了日益重要的作用。

大地線主題計算包括正解和反解兩類。 正解是已知起點坐標、起點大地方位角和大地線距離，求終點坐標和終點大地方位角；反解是已知起始點坐標，求兩個點的大地方位角和大地線距離。 目前，大地線主題計算的方法有70 多種，求解的出發(fā)點也各不相同。 勒讓德級數(shù)公式、高斯平均引數(shù)公式［4］等是通過對大地線微分方程進行離散近似之后得到大地線長，其計算精度與距離有關(guān)，距離越長精度越低［3］。 對大地線微分方程直接進行數(shù)值積分的方法雖然可以保證計算精度， 但是其計算時間較長，且有不收斂的風(fēng)險。 為此，采用貝塞爾大地主題解算方法來計算大地線長。 該方法適用于長距離計算，具有精度高、計算速度快等優(yōu)點。 將計算模型應(yīng)用于氣象預(yù)報數(shù)據(jù)插值匹配，為智能船數(shù)據(jù)分析打下堅實的基礎(chǔ)。

1 大地線主題計算

1.1 計算目標

大地線長度可表征兩點間最短球面距離，可用于最短航線設(shè)計、偏航監(jiān)測、數(shù)據(jù)插值等多個功能模塊。 大地線長度的精確計算也是算法中最具應(yīng)用價值的技術(shù)點，因而本文著重介紹大地線長度的計算方法和應(yīng)用案例。 大地線正反解問題完整算法可參閱文獻［3］和［5］。

1.2 地球橢球模型

地球橢球模型采用WGS84 標準。 以B、L 表示地理緯度和經(jīng)度，約定東經(jīng)與北緯坐標取正值，西經(jīng)與南緯坐標取負值，暫不將極點納入計算范圍。 因而緯度與經(jīng)度的取值范圍為：

B∈（-π/2，π/2），L∈（-π，π）

結(jié)合橢球幾何公式，需引用的橢球幾何要素列于表1。

表1 地球橢球模型幾何要素

1.3 貝塞爾投影變換參數(shù)關(guān)系

采用貝塞爾投影法求大地線長。 其基本思路為：將橢球面上的幾何元素按貝塞爾投影條件映射到輔助球面上， 繼而在球面上進行大地主題解算，最后再將球面計算結(jié)果變換回橢球面。 貝塞爾投影需滿足以下條件：（a）橢球面大地線投影到球面為大圓?。唬╞）大地線和大圓弧上對應(yīng)點的方位角相等；（c） 球面上任意點的緯度等于橢球面上對應(yīng)點的歸化緯度。 如圖1 所示，設(shè)PN為橢球北極點，大地線P1P2交赤道于P0，設(shè)pn為輔助球面北極點。 將大地線P0P1P2按貝塞爾投影條件變換為輔助球面上的大圓弧p0p1p2。 貝塞爾投影變換下，大地元素滿足以下關(guān)系：

圖1 貝塞爾投影幾何要素

1）P1、P2的地理坐標（B1，L1），（B2，L2）對應(yīng)p1、p2的球面坐標（φ1，λ1），（φ2，λ2）。其中：緯度坐標滿足式（2），進一步還可求得正弦和余弦函數(shù)，見式（6）；經(jīng)度坐標滿足貝塞爾微分方程，見式（7）。

2）P0、P1的大地方位角A0、A1對應(yīng)p0、p1的大圓弧方位角A0、A1。 投影前后方位角保持不變。

3）P1和P2間大地線長S（長度）對應(yīng)p1和p2間大圓弧長σ（弧度）。 S 和σ 滿足貝塞爾微分方程：

1.4 貝塞爾投影法求大地線長

1.4.1 輔助計算

為了求得大地線長S， 需完成以下兩步輔助計算：

第一步：在輔助球面上計算sinA0和σ。

令Δλ=λ2-λ1，由球面三角形p1p2pn可得：

tan A1=q1/q2（9）

其中：

假定：

q1≠0，q2≠0 （11）

那么可求得A1：

進一步求得sin A0和σ：

由球面三角形p0p1pn可得：

sin A0=cos φ1sin A1（13）

由球面三角形p1p2pn得：

聯(lián)立式（9）～（14），獲得σ 與Δλ 以及sin A0與dλ 之間的關(guān)系：

sin A0=F1Δλ （15）

σ=F2Δλ （16）

第二步：通過貝塞爾微分方程求經(jīng)差改正量。

令ΔL=L2-L1，對公式（7）積分并展開為級數(shù)形式，那么經(jīng)差改正量可表示為：

考慮大圓弧p0p1p2上任意點p（φ，λ），設(shè)該點大圓弧方位角為A，p0p 大圓弧長為σ， 子午線pnp 交赤道于p′，如圖2 所示。

圖2 大圓弧上任意一點p 及相關(guān)幾何要素

由球面三角形p0ppn可得：

再由球面三角形p0pp′可得：

sinφ=cosA0sinσ?cos2φ=1-cos2A0sin2σ （19）

進而將式（17）的被積函數(shù)替換為：

e 為小量，略去e8項及e8的高次項，再應(yīng)用三角恒等式求式（20）的積分形式，可得：

δ=sin A0［ασ-2β（2sin φ1sin φ2-cos2A0cos σ）sin σ］（21）

其中：

1.4.2 求大地線長

輔助計算公式建立后，聯(lián)立（15）、（16）和（21）式，可采用迭代法同時求得δ、sinA0和σ。再將σ 代入貝塞爾微分方程公式（8）可得S 的級數(shù)形式：

其中k=e′cos A0為小量， 略去k8項及k8的高次項，仿照1.4.1 節(jié)經(jīng)差改正量的積分推導(dǎo)過程，可得S 與δ 之間的關(guān)系：

S=C1σ+2C2（2sin φ1sin φ2/cos2A0-cos σ）sin σ+2C3［1-2（2sin φ1sin φ2/cos2A0-cos σ）2］cos σ sin σ（23）

其中：

1.5 特殊工況大地線長度計算方法

1.4 節(jié)推導(dǎo)的大地線算法并不能適用于任意工況，還需對以下特殊工況單獨處理。

1.5.1 經(jīng)差等于180°

根據(jù)1.2 節(jié)的約定，ΔL=L2-L1取值范圍為-2π＜ΔL＜2π。 當=π 時，等同于大地線經(jīng)過極點，不在本文討論范疇，計算程序應(yīng)當報錯退出。

1.5.2 經(jīng)差大于180°

1.5.3 經(jīng)差為0

若ΔL=0，式（11）不能滿足，1.4 節(jié)的算法失效。此時大地線長就是子午線上的弧段，那么可調(diào)用式（5）求得大地線長

1.5.4 起始點同位于赤道

若B1=B2=0，式（11）不能滿足，1.4 節(jié)的算法失效。 此時大地線就是赤道上的弧段，那么可調(diào)用公式（3）求得大地線長：S=。

經(jīng)過以上修正之后，算法適用于跨越180°經(jīng)線計算。 同時，式（14）和式（23）表明算法也適用于跨赤道計算。

1.6 算法的適用條件

貝塞爾投影并不是同胞映射。 具體來說，由公式（7）可知：

那么，貝塞爾投影算法在赤道附近的適用范圍為：

需要指出的是，航?；顒邮艿匦蜗拗?，一般不會求解超長距離大地線，式（25）作為約束條件是合適的。 但是在宇航、測繪等領(lǐng)域，常需要繪制繞地球半圈的大地線，當貝塞爾投影算法不再適用時，還需要進行特別處理［6］。

2 計算結(jié)果驗證

選用文獻［6］的計算案例，與本文算法的計算結(jié)果進行對比，如表2 所示。 計算結(jié)果表明，本文提供的算法可以準確計算短距離和長距離大地線，也可以跨赤道計算。 計算結(jié)果與文獻［6］的計算結(jié)果偏差約-0.002%。

表2 大地線主題計算模型驗證

3 大地線主題計算在智能船數(shù)據(jù)分析中的應(yīng)用

獲得大地線算法之后，可應(yīng)用于智能船數(shù)據(jù)分析。 以氣象數(shù)據(jù)后處理為例，介紹大地線算法的應(yīng)用。

圖3 為“明遠”輪某段航跡，記錄時間為2019年3 月8 日9 時至2019 年3 月13 日9 時。“明遠”輪是招商系列第3 艘400 000 DWT 級超大型智能礦砂船，也是智能船1.0 專項首批示范船。 通過實時獲取全船運營數(shù)據(jù)， 該船實現(xiàn)了輔助自動駕駛、能效管理、設(shè)備運維、船岸一體、礦物液化監(jiān)測等五大智能功能。

圖3 “明遠”輪航行軌跡

以T、B 和L 分別表示時間坐標、緯度坐標以及經(jīng)度坐標。 基于航測數(shù)據(jù)庫，航跡線可以表示為：

B=fB（T），L=fL（T） （26）

另一方面， 船舶航測數(shù)據(jù)庫記錄了隨船風(fēng)速w和風(fēng)向θ（對地風(fēng)），它們是時空坐標［T，B，L］的函數(shù)，分布于航跡線上，航測數(shù)據(jù)的時間間隔為1 s。同時，氣象預(yù)報數(shù)據(jù)庫也提供了相同時間段內(nèi)絕對風(fēng)速W 和絕對風(fēng)向Θ 的預(yù)測值，它們是時空坐標［T，B，L］的函數(shù)，分布于全平面，氣象預(yù)報數(shù)據(jù)的時間間隔為3 h，空間間隔為1°。

顯然，基于［w，θ］和［W，Θ］，可評估氣象預(yù)報的準確度。 在此之前，需將兩份數(shù)據(jù)進行坐標匹配，即對數(shù)據(jù)庫［W，Θ］做矢量插值，求得航跡線［式（26）］上的氣象預(yù)報值。 基于大地線主題計算，可設(shè)計匹配算法如下：

第一步：將矢量數(shù)據(jù)庫［W，Θ］變換為標量數(shù)據(jù)庫，即

U=WcosΘ，V=WsinΘ （27）

式中：U 和V 分別表示沿正北和正東方向的風(fēng)速。

第二步：時間坐標匹配。 氣象預(yù)報數(shù)據(jù)的時間分辨率為3 h，提供了時間序列［T0，T1，T2，…，T40］上的全球氣象預(yù)報信息。 具體來說：

T0=′2019 年3 月8 日9 時′，T1=′2019 年3 月8日12 時′，…，T40=′2019 年3 月13 日9 時′

為了對比預(yù)報數(shù)據(jù)和實測數(shù)據(jù)，在航測數(shù)據(jù)庫中篩選出時間序列［T0，T1，T2，…，T40］上的船位和風(fēng)速風(fēng)向：

從而獲得41 個船位點：（Bi，Li），i=0，1，2， …，40 和41 組 實 測 數(shù) 據(jù)：wi（Ti，Bi，Li），θi（Ti，Bi，Li），i=0，1，2，…，40。

第三步：空間坐標匹配。 氣象預(yù)報數(shù)據(jù)的空間分辨率為1°， 提供了整數(shù)經(jīng)緯坐標上的氣象預(yù)報值，如圖4 所示。 為了與實測數(shù)據(jù)進行對比，需求得時空坐標（［Ti，Bi，Li］，i=0，1，2，…，40）上的氣象預(yù)報值。以T0時刻為例，此時船位為［B0，L0］，對B0，L0向上和向下取整，找到距離船位最近的4 個氣象點，它們與船之間的距離即大地線長S1，S2，S3，S4。 在4 個氣象點上，風(fēng)信息［U，V］是已知的，故而可采用線性插值估算船位點的風(fēng)信息：

圖4 氣象預(yù)報數(shù)據(jù)空間匹配

圖5 實測數(shù)據(jù)與氣象預(yù)報數(shù)據(jù)對比

從而獲得41 組預(yù)報數(shù)據(jù)Ui（Ti，Bi，Li），Vi（Ti，Bi，Li），i=0，1，2，…，40

第 四步：將Ui（Ti，Bi，Li），Vi（Ti，Bi，Li），i=0，1，2，…，40 變換回矢量形式，得：

Wi（Ti，Bi，Li），Θi（Ti，Bi，Li），i=0，1，2，…，40

實現(xiàn)匹配計算之后，對氣象預(yù)報數(shù)據(jù)和實測值進行比較，如圖5 所示。 分析表明：對23 號之前的船位點（3 天內(nèi)），氣象預(yù)報數(shù)據(jù)和實測數(shù)據(jù)較吻合，氣象預(yù)報能夠合理地預(yù)測未來3 天的風(fēng)速風(fēng)向信息。24 號船位點之后的預(yù)報數(shù)據(jù)與實測數(shù)據(jù)差異較大。 鑒于氣象預(yù)報工作本身的復(fù)雜性，建議謹慎對待長期和超長期預(yù)報數(shù)據(jù)，盡可能只用3 天內(nèi)短期預(yù)報數(shù)據(jù)服務(wù)于智能船。 在確保預(yù)測準確度的前提下，提出運營建議。

4 結(jié) 語

將地球視為橢圓體，利用貝塞爾投影理論建立大地線長度的計算模型。 對原始算法進行補充和修正，使算法適用于跨赤道和跨180°經(jīng)線計算。 通過理論分析提出算法成立的充分條件。 通過與文獻資料對比驗證了算法的準確性。 利用程序設(shè)計語言將計算模型應(yīng)用于智能船數(shù)據(jù)分析，討論了氣象預(yù)報數(shù)據(jù)的準確度并提出建議。