王 宇，肖 烜，劉績寧，鄧志紅，王 博

（北京理工大學，北京 100089）

2021年3月發(fā)布的《中華人民共和國國民經(jīng)濟和社會發(fā)展第十四個五年規(guī)劃和2035年遠景目標綱要》中提出，要將深海探測作為科技前沿領域攻關，將精密重力測量研究設施作為國家重大科技基礎設施。重力場是地球的固有物理場[1]，重力儀等測量儀器可以實時獲得所在位置的重力異常信息，由于重力異常信息具有良好的時空位置特征，且不破壞水下航行器的隱蔽性，因此使用重力輔助慣性導航系統(tǒng)，能夠安全、有效地提高慣性導航系統(tǒng)的導航定位精度。

適配區(qū)選取是水下重力輔助慣性導航系統(tǒng)的關鍵技術之一[2,3]。傳統(tǒng)適配區(qū)選取算法主要有閾值法和單一特征值選取法等，主要對重力場統(tǒng)計特征參數(shù)進行比較分析[4]，而忽略了重力場的空間特點和幾何特征，導致漏掉潛在的適配區(qū)域。

本文提出了一種基于分割嵌套三角剖分的重力場三維模型，首先利用墨卡托投影和重力異?？臻g校正，將傳統(tǒng)重力場柵格信息變換為三維高程信息。再利用分割嵌套的思想，不斷從重力場最小環(huán)形域中分割出最優(yōu)三角形，從而形成重力場三角網(wǎng)絡。對網(wǎng)絡中每一個小三角形，選取了坡度、坡向、起伏度和粗糙度作為重力場局部幾何特征，采用聚類算法，將重力場分為強適配區(qū)、弱適配區(qū)和非適配區(qū)，據(jù)此進行重力場適配區(qū)選取，使得適配區(qū)的選取更為有效。

1 重力場柵格信息預處理

1.1 經(jīng)緯度數(shù)據(jù)投影變換

墨卡托變換可以使得投影坐標系更合適地描述不規(guī)則的重力場特征，對于重力場特征變化較為明顯、起伏程度更大的區(qū)域可以表示得更加精確。將赤道作為標準緯線，本初子午線作為中央經(jīng)線，赤道與本初子午線的交點作為坐標原點，取東向北向為正方向。根據(jù)等角條件，將地球表面的信息投影到一個平面上，地球表面的每一部分都對應相應平面上的微積線段，求得對應微積線段的關系式。然后先將地球近似視為一個標準的球體：在地球表面取一點A，假設A點在經(jīng)緯度坐標系下的坐標為(λ,)φ，在墨卡托投影坐標系下對應的投影坐標為(x,y)，將赤道設置為基準緯線，用R表示地球半徑?？傻媚ㄍ型队胺匠淌綖椋?/p>

最后推至更一般的表達式，將地球視為旋轉中的橢球體，對應的(X,Y)=f(B,L)的墨卡托正解變換由式(2)表示：

式中，a為橢球體長半軸長度，取6378136.49 m，b為橢球體短半軸長度，取6356755.00 m，e為第一偏心率，取0.0818，e′為第二偏心率，取0.0821，N為卯酉圈曲率半徑長度。

1.2 重力異常值z 因子轉換

重力場數(shù)據(jù)在投影坐標系中的橫、縱坐標單位為米，而重力異常的單位為毫伽（mGal），由于單位不同，因此將z因子作為重力異常值的轉換因子，通過映射關系將重力異常值信息與高度信息之間建立聯(lián)系，這樣可以保證重力場三角剖分得到的三維曲面單位一致，都用長度單位米來表示。

根據(jù)重力異?？臻g校正的泰勒展開式，可以將空間校正看作將大地水準面上的正常重力值γ0()φ進行延拓至高度h上的正常重力值與原來高度的差異。將高程信息h作為地球半徑R的增量，可得：

其中誤差項即為高度校正值：(?γ/?r)0h+(1/2)(?2γ/?h2)0h+…。

若將地球看作扁率為f的橢球，可得與平靜海平面處于一致狀態(tài)的大地水準面上的一點正常重力值g0與該點高程為h處的正常重力異常值gh之間滿足以下方程：

其中ge為地球赤道上的重力值；R為地球的平均半徑；q為赤道上慣性離心力與引力之比；φ為地理緯度。

地球參數(shù)采用ge=9780300mGal，R=6371.2km，f=1/298.3，q=0.003467。得到空間校正公式近似為：

若進行進一步的簡化，將地球看作半徑為R的均質圓球，則距離地面高度h處的重力值為：

對應的大地水準面上點的引力為：

將式(6)(7)作差可得：

將地球半徑和地球平均重力代入可得一階近似公式：

由此可得，每1 m 的高程位移產(chǎn)生約0.3086 mGal的重力異常變化。由式(2)(10)可以將經(jīng)緯度信息、重力異常值信息進行數(shù)據(jù)預處理，轉化為以米為單位的坐標。如將圖1的重力場傳統(tǒng)柵格信息經(jīng)過墨卡托投影變換后得到圖2所示區(qū)域。

圖1 傳統(tǒng)重力場等值線平面結構Fig.1 Plane structure of isoline of traditional gravity field

圖2 墨卡托投影變換后的等值線圖Fig.2 Contour map after Mercator transformation

2 基于分割嵌套三角剖分的重力場三角網(wǎng)絡構建

重力場三角網(wǎng)是根據(jù)已知的傳統(tǒng)重力場重力異常柵格點及其中蘊含的位置等信息，通過三角剖分等處理方式，將這些點按照一定規(guī)則、一定順序串聯(lián)起來，構成由互不相交、互不重疊的三角形形成的三角網(wǎng)絡[6]。這種三角網(wǎng)絡對于處理復雜重力曲面信息，分析其中的拓撲結構具有重要作用，因此可以用構建重力場的三角網(wǎng)絡來對重力場幾何拓撲特點進行分析和提取[7]。

在構建重力場三角網(wǎng)絡時，需要利用三角剖分算法。常用的三角剖分算法例如Delaunay 三角剖分[8,9]，這種算法更側重于對三維立體結構進行剖分，更適用于對某個場立體特征的分析。而本文提出的基于分割嵌套環(huán)形域三角剖分算法能夠較為準確地對重力場表面的信息進行描述，可以將重力場背景圖中提供的重力場點集的幾何特征表現(xiàn)出來，且算法的可行性較高，復雜度較小，還可以解決許多較常見三角剖分算法無法解決的最小權三角剖分的實現(xiàn)問題，所以在構造重力場三角網(wǎng)絡時采用這種基于分割嵌套的三角剖分算法。

將重力場背景圖中提供的重力場點集經(jīng)坐標變換投影到經(jīng)度-緯度坐標平面上，由此得到一組平面重力場點集，再將平面重力場點集中全部點組成一個凸殼，通過不斷篩選凸殼中的頂點，找到位于最內層的嵌套凸殼，從此開始可以分為三種情況：最內層凸殼不包含、包含1 個，或者包含2 個平面重力場點集中的點。包含3 個或3 個以上平面重力場點集中的點的情況是不可能發(fā)生的，因為如果平面重力場點集中有3 個點都位于篩選得到的最內層凸殼中，則該凸殼中還可以繼續(xù)劃分出一個位于其內側的凸殼，也就意味著它不能作為嵌套區(qū)域中最內層的凸殼。

接下來對三種情況分別處理。若最內層凸殼中不包含平面重力場點集中的點，則求出該凸殼的直徑，從中取點作為判斷時的頂點，判斷凸殼中點的共線情況，將不共線的點與其他點連成線段，比較所得線段之間的長度，刪去長度較長的線段，保留長度較短的線段，對每個不共線的點都進行這樣的步驟，直到將整個凸殼完全分割成三角形。若最內層凸殼中包含1個平面重力場點集中的點，則將該點分別與凸殼的各頂點相連接，將這些連成的線段按照本身的長度由大到小的順序排序，得到對應的凸殼頂點的序列，再將這些凸殼頂點按順序作為判斷時的頂點，判斷凸殼內部的共線情況，將不共線的點與其他點連成線段，通過比較線段長度的大小關系，對線段進行添加和刪除，直到將整個凸殼完全劃分為三角形。若最內層凸殼中包含平面重力場點集中的2 個點，則將這2 個點連接成一條線段，取凸殼中位于該線段右側的點，將該點分別與最內層凸殼包含的平面重力場中的2 個點相連接，這樣就形成了一個完整的凸殼，在新得到的凸殼中判斷點與線的位置關系，再通過線段長度的比較進行篩選，將該凸殼進行劃分，直至使之成為全部由三角形構成的凸殼。

在對平面重力場點集中最內層嵌套凸殼處理后，再對中間的環(huán)形區(qū)域進行劃分，可以取出一條固定邊，分別選出環(huán)形區(qū)域中位于固定邊同一側的點鏈，將固定邊的兩個端點分別與該點鏈相連接，由此形成一個新的凸殼，再用與處理最內層相同的方法對其進行劃分，即可得到完整的三角剖分結果。分割嵌套環(huán)形域的三角剖分算法流程圖如圖3所示。

圖3 分割嵌套環(huán)形域的三角剖分算法流程圖Fig.3 Flow chart of triangulation algorithm for segmenting nested ring domain

3 重力場三角網(wǎng)絡幾何特征選取

從經(jīng)過三角剖分后的重力場三角網(wǎng)絡中，可以提取重力場幾何空間特征參數(shù)來對重力場特征，尤其是局部特征進行分析，進而更加合理、更加精確地描述重力場三維表面特性。本文提取的重力場幾何空間特征包括：坡度、坡向、起伏度和粗糙度[10]。其中，坡度和坡向能夠表現(xiàn)重力場局部三角形法向量在俯仰角及旋轉角上的變化。起伏度和粗糙度能夠表征重力場局部曲面變化特征及起伏程度。

3.1 坡 度

坡度為重力場三角面的法線矢量與z軸矢量的夾角，可以表示重力場曲面局部的最陡坡傾斜程度，是每個三角面中的最大重力異常變化率。坡度值的范圍為0°～90°。坡度值越小，重力場曲面變化越平緩，坡度值越大，重力場曲面變化越陡峭。坡度slope的表達式如下：

3.2 坡 向

坡向為重力場三角面的法線矢量在水平面的投

影矢量與正北方向矢量的夾角，可以表示曲面的水平方向，是坡度所面對的方向。坡向按照順時針方向進行測量，坡向角度范圍介于0°（正北）到360°（仍是正北）之間，即完整的圓。坡向可以用于識別重力場曲面某一位置處的最陡下坡方向。坡向aspect的單位為°，表達式如下：

其中，nx表示標準法線矢量在x軸的分量，ny表示標準法線矢量在y軸的分量。

3.3 起伏度

起伏度是指重力場三角面內最大的重力異常值zmax與最小重力異常值zmin之差，是描述重力場特征的一個宏觀性指標，起伏度反映的是區(qū)域內重力場起伏變化的特征，起伏度越大，則說明重力場特征變化越豐富。起伏度relief的單位為mGal，表達式如下：

3.4 粗糙度

粗糙度是重力場空間三角面積S與其在水平面上的投影面的面積S' 之比，它是描述重力場曲面的侵蝕程度的重要指標，粗糙度越大，則說明重力場空間表面越曲折，特征變化越豐富。粗糙度rough表達式如下：

4 基于k-means 的適配區(qū)聚類分析

為了利用重力場三角網(wǎng)絡的幾何特征，需要對其進行數(shù)據(jù)挖掘。k-means方法是一種經(jīng)典的迭代求解聚類分析算法，被廣泛應用在機器學習、圖形分析領域。本文利用k-means方法對4 個重力場幾何特征值進行聚類分析，從而將重力場按照適配性劃分為3 種區(qū)域：強適配區(qū)、弱適配區(qū)和非適配區(qū)。

從重力場三角網(wǎng)絡中隨機選取3 個三角面的4 個標準化處理后的幾何特征值作為初始聚類中心，計算其余三角面的標準化幾何特征值到這些聚類中心的相似性度量。本文采用的歐氏距離度量如下：

其中，i,j為重力場三角網(wǎng)絡中局部三角形編號，xij為重力場三角面的標準化幾何特征值。

根據(jù)相似性度量的計算結果，將不同的重力場三角形分配給相似度最高的聚類中心，更新聚類中心并重新分配幾何特征值，在聚類過程中使用的聚類準則函數(shù)為：

其中，mij為聚類中心特征向量，E為均方差之和。重復迭代更新聚類中心和特征向量分配環(huán)節(jié)，當E的值不變時，代表算法收斂，最終將重力場按照適配性劃分為3 種適配區(qū)域。

圖4 基于三角剖分的適配區(qū)選取流程圖Fig.4 Flow chart of area selection based on triangulation

5 仿真實驗結果及分析

選取北緯24°～28°，東經(jīng)126°～130°范圍的重力場背景圖進行實驗，分辨率為1′×1 ′，經(jīng)改進克里金算法將其插值為0.5′×0.5′，其中重力異常值最大為91.09 mGal，最小為-75.97 mGal。在仿真中，x軸、y軸坐標為墨卡托投影變換后的經(jīng)度、緯度坐標，z軸為經(jīng)過重力異?？臻g校正變換后的重力異常值。采用本文提出的基于分割嵌套三角剖分算法，將傳統(tǒng)的重力場等值線平面格網(wǎng)結構，變換為重力場立體三角網(wǎng)結構，如圖5所示。形成由115198 個小三角形組成的連續(xù)三角面，來描述重力場的幾何拓撲特征。

圖5 嵌套分割三角剖分后的局部重力場立體三角網(wǎng)結構Fig.5 Three dimensional triangulation structure of local gravity field after nested segmentation triangulation

在經(jīng)過三角剖分后的重力場三角網(wǎng)中，選擇每一個局部三角形的坡度、坡向、起伏度和粗糙度作為重力場幾何空間特征參數(shù)。4 個特征參數(shù)在仿真區(qū)域內的空間分布如圖6所示。

圖6 重力場坡度、坡向、起伏度和粗糙度空間分布圖Fig.6 Spatial distribution of gradient,aspect,fluctuation and roughness of gravity field

經(jīng)過對坡度、坡向、起伏度和粗糙度的k-means聚類后，得到的該區(qū)域內重力場適配性如圖7所示，其中深藍色區(qū)域代表強適配區(qū)，黃色區(qū)域代表弱適配區(qū)，綠色區(qū)域代表非適配區(qū)。同時采用傳統(tǒng)的閾值法多指標融合方法，選取重力場標準差和經(jīng)緯度方向相關系數(shù)作為重力特征參數(shù)[14]，并通過經(jīng)驗確定閾值k=10，應用閾值法得到最終選取標準，并在同一片區(qū)域內進行適配區(qū)選取，得到的結果如圖8所示。

圖7 3 種重力場適配區(qū)空間分布圖Fig.7 Distribution map of three gravity adaptation areas

圖8 傳統(tǒng)閾值法適配區(qū)選取結果Fig.8 Selection results of adaptation region of traditional threshold method

結合圖7-8 可以看出，同一片重力場內，傳統(tǒng)方法與本文所提方法在適配區(qū)的選擇上有較大不同，主要體現(xiàn)在圖8的三片區(qū)域內。傳統(tǒng)閾值法將區(qū)域1 和區(qū)域2 劃分為非適配區(qū)，區(qū)域3 作為適配區(qū)，而本文所提出的適配區(qū)選取方法，將區(qū)域1 作為強適配區(qū)，區(qū)域2 作為弱適配區(qū)，而區(qū)域3 作為非適配區(qū)。傳統(tǒng)閾值法過于強調參數(shù)閾值的作用，而缺少多指標融合的手段，忽略了重力場的空間特點和幾何特征，導致漏掉潛在的適配區(qū)域。

本文算法選出來的強適配區(qū)（如區(qū)域1），重力場的坡度較大，坡向變化明顯，同時起伏度和粗糙度較大，表明該區(qū)域內，重力場的幾何空間特征明顯，重力異常信息豐富，適合進行重力匹配定位。根據(jù)仿真結果，選擇圖7中紅圈內區(qū)域（即圖8中的區(qū)域1）作為適配區(qū)，規(guī)劃水下航行器的仿真軌跡。設航行器在水下進行勻速直線運動，航行速度為20 節(jié)（約10.289 m/s），航向為西北方向，設定航行時間為150 min。陀螺常值偏移為0.01°/h，陀螺隨機誤差為0.001°/h，加速度計零偏為100 μg，隨機誤差為50 μg。

如圖9所示，黑色的曲線是水下航行器的真實軌跡，藍色的曲線是慣性導航系統(tǒng)輸出軌跡，可以看出經(jīng)過一段時間的航行后，INS 輸出軌跡逐漸發(fā)散，偏離真實軌跡。而紅色的曲線代表基于ICCP 算法的重力匹配航跡，與真實軌跡幾乎重合，匹配定位平均位置誤差與重力場背景圖網(wǎng)格分辨率相當。表1為重力匹配定位解算結果。仿真結果表明，本文所提算法能夠選擇被傳統(tǒng)適配區(qū)選取算法漏掉的強適配區(qū)，在此適配區(qū)內，重力匹配定位效果好，精度高。

圖9 強適配區(qū)內重力匹配定位航跡圖Fig.9 Gravity matching positioning track map in strong adaptation area

表1 強適配區(qū)內重力匹配定位誤差Tab.1 Gravity matching positioning error in strong adaptation area

6 結 論

本文針對傳統(tǒng)重力場適配區(qū)選取算法缺乏對重力場三維信息、局部信息充分利用的問題，提出一種基于分割嵌套三角剖分的重力場適配區(qū)選取算法。首先利用墨卡托投影和重力異?？臻g校正，將傳統(tǒng)重力場柵格信息變換為三維高程信息。再利用分割嵌套的思想，不斷從重力場最小環(huán)形域中分割出最優(yōu)三角形，從而形成重力場三角網(wǎng)絡。對網(wǎng)絡中每一個小三角形，本文選取了坡度、坡向、起伏度和粗糙度作為重力場局部幾何特征，采用k-means聚類算法，將重力場分為強適配區(qū)、弱適配區(qū)和非適配區(qū)。在強適配區(qū)內，重力場的坡度較大，坡向變化明顯，同時起伏度和粗糙度較大。仿真實驗結果表明，基于分割嵌套三角剖分的重力場適配區(qū)選取算法彌補了傳統(tǒng)方法描述重力場特征時單一、片面的缺陷，能夠挖掘出重力場更多的局部信息減少了人工干預和人工經(jīng)驗判斷，更加全面客觀地劃分適配區(qū)。未來可以繼續(xù)針對三角剖分后的重力場進行幾何特征挖掘，同時將軌跡信息也進行相似處理，利用重力場三角形和軌跡三角形進行基于計算幾何的匹配算法研究。