虛擬球諧方法逼近市級大地水準面計算分析

王建強 孫云龍

1 東華理工大學測繪工程學院，南昌市廣蘭大道418號，330013

高程基準的統一是大地測量的長期目標[1]，對工業、經濟及大型工程建設都具有重要的支撐作用。隨著GNSS技術的廣泛應用及我國北斗系統的全面建成，構建高精度高程基準成為高程測量實時化建設需要突破的關鍵技術。發達國家近幾十年不斷精化本國大地水準面，美國最新大地水準面模型Geoid18[2]的分辨率優于2 km，內符合精度優于2 cm；澳大利亞發布的AUSGeoid2020內符合精度達到亞cm級[3]；加拿大同美國正合作構建Geoid2022，屆時可全面覆蓋加拿大區域[4]。我國的大地水準面研究技術不論是在理論還是工程應用方面均已走在世界前列[5-6]，李建成[1]在2020-10中國測繪學會年會上指出，中國似大地水準面的擬合精度已達2～3 cm，但構建高精度(似)大地水準面仍有提升空間，還需開展更多研究分析。本文研究構建市級區域似大地水準面的虛擬球諧方法，并以南昌市為例進行計算分析。

1 位模型計算大地水準面

地球外部引力位的球諧表達式[7]為：

(1)

(2)

式中，J2n為系數值，當J2已知時，其他數值可以計算得出；P2n(θ)為勒讓德函數。將式(1)減去式(2)可得擾動位重力異常位T：

T(r,θ,λ)=V(r,θ,λ)-U(r,θ,λ)

(3)

(4)

(5)

(6)

式中，γ為正常重力值，可通過公式直接計算獲取。

利用式(6)計算出的大地水準面有2個誤差源，分別為位系數誤差和截斷到Nmax階所引起的誤差，后者被稱為模型截斷誤差，可采用重力異常進行估算。大地水準面上的重力異常表達式為:

(7)

(8)

由此可得：

(9)

將式(9)代入式(6)可得：

(10)

擾動引力分量的截斷誤差表達式為：

(11)

式中，R為地球平均半徑；Cn可利用Moritz兩分量模型和Lapp參數進行計算：

(12)

其中，C2=7.5 mGal2，計算階次的上限為100 000，R=6 371 km，結果如圖1所示。由圖可知，要達到cm級精度的大地水準面，截斷階次要達到2 000以上。

圖1 大地水準面模型截斷誤差Fig.1 Geoid model truncation error

2 虛擬球諧理論與方法

球冠諧模型是構建區域大地水準面模型的典型方法[8]。在球冠坐標系下，球冠半徑為θ0，任一點的坐標為(r,θ,λ)，θ為余緯，λ為經度。在球冠諧分析中，地球外部重力場位滿足Laplace方程，同樣可用分離變量法獲得方程的解[9]。由于球冠諧分析中，余緯的邊界條件[9]為：

T(r,θ0,λ)=f(r,λ)

(13)

(14)

可以看出，兩式等號右端的函數均與θ無關，而在球諧分析中，僅有締合勒讓德函數與θ有關。式(13)和式(14)的基函數可以通過以下2個方程分別滿足：

(15)

(16)

由于其經度范圍和球諧分析中的經度范圍相同，因此對應的本征值m為整數變量。當給定式(15)和式(16)的θ0時，可單獨確定一組對應m的n值序列，此時的n為非整數，非整階締合勒讓德函數的特征值可通過Muller方法計算得到[10]。通過球冠來確定非整階勒讓德函數序列[11]時, 需要2個正交基數，令k為n值序列的下標，并定義k-m為奇數時，采用式(15)獲取的本征值序列；k-m為偶數時采用式(16)獲取的本征值序列。可以看出，球冠諧函數描述的地球重力場是局部區域的，與全球區域的球諧函數存在差異：在求解關于余緯的偏微分方程時，球諧分析中方程的本征值是整數，而球冠諧分析中的本征值是非整數。擾動位球冠諧展開式為：

(17)

由于受非整階締合勒讓德函數的限制，球冠諧模型很難達到高分辨率，球冠半徑較小(如5°以下)的函數變化更為迅速。圖2為球冠半徑為5°，m=0、1、2、4時的函數值，可以看出，能計算出的零根值數量在10個左右。Haines[13]利用逼近方法改進算法，但仍有較大誤差。利用球冠諧模型逼近區域大地水準面可達到cm級的逼近精度，但隨著區域的擴大及地形的復雜化，球冠諧模型仍受階次擴展的限制。為克服這一限制，借用球冠諧映射方法，將球冠坐標系變換后采用整階次締合勒讓德函數進行計算[14]。

圖2 締合勒讓德函數值Fig.2 Legendre values

確定逼近區域后，根據球冠諧理論設計球冠半徑和球冠中心，在保持相當精度的前提下，將球冠諧函數改進為虛擬球諧函數。將余緯θ的定義域(0，θ0)映射到(0，π)上，用整階次締合勒讓德函數代替非整階締合勒讓德函數。虛擬球諧方法需要將原坐標系轉換到新坐標系中：

(18)

3 數值實驗

以南昌市為例進行實驗(圖3)，數據點由主要交通路線隨機采樣獲得。通過EGM2008模型計算該區域內階次分別為2～60、61～360、361～1 000、1 001～2 100的大地水準面，計算結果如圖4所示。可以看出，60階次內的重力場模型計算的大地水準面起伏范圍達到2 m，61～360階次范圍內的大地水準面起伏范圍達到40 cm，360階次以上的大地水準面起伏范圍大約為20 cm，此時的起伏范圍已經較小。對圖3中的采樣點進行數值計算，得到剩余大地水準面數據范圍為-42～-30 cm，具體見圖5。在此基礎上，增加1 cm 的白噪聲誤差。

圖3 數據采集范圍Fig.3 Data collection range

圖4 不同階次大地水準面起伏情況Fig.4 The geoid undulates in different orders

圖5 去除360階次重力場模型后的大地水準面數據Fig.5 The geoid after removing the 360-ordergravitational field model

通過重力場模型濾波后得到新的觀測值，然后將地理坐標轉換為球冠坐標。實驗數據范圍大約在1°左右，其中半徑為0.5°，添加0.5°的過渡帶。根據虛擬球諧理論設置虛擬半徑[15]為1°，模型階次設定為13，利用式(18)將球冠坐標轉換為新的虛擬坐標，并在新坐標系下利用球諧函數構建模型。為突出虛擬球諧理論的優越性，加入多項式擬合和神經網絡模型進行對比，結果如圖6所示。

圖6 誤差分布Fig.6 Error distribution

實驗結果表明，虛擬球諧方法的擬合精度最高，達到1 cm，RMS值和STD值是4種方法中最小的，都只有0.308 cm；其次是神經網絡模型，擬合精度基本優于2 cm；4次多項式的擬合誤差較大，接近5 cm。具體實驗結果如表1所示。

表1 精度表

4 結 語

本文研究構建市級大地水準面的虛擬球諧理論及方法，并以南昌市為例，對比分析了虛擬球諧方法、多項式擬合方法及神經網絡模型。從實驗結果可以看出，虛擬球諧方法的RMS值只有0.308 cm，低于多項式擬合和神經網絡模型，表明虛擬球諧方法偏差更小、效果更好。本文僅采用重力場位模型數據對觀測值進行簡單的移去和恢復操作，沒有綜合利用重力數據和地形數據，因此未來仍需作進一步探索研究。