顧及地球重力場模型的高程轉換方法

王式太，肖守圣，殷 敏*，華 磊，姜新偉，程 波

(1. 桂林理工大學 測繪地理信息學院，廣西 桂林 541006； 2. 廣西空間信息與測繪重點實驗室，廣西 桂林 541006； 3. 中國地質大學(武漢)公共管理學院，湖北 武漢 430074; 4. 山東省地質測繪院，山東 濟南 250013)

0 引言

無人機攝影測量具有分辨率高、成本低和任務周期短等優點[1-2]，通過無人機攝影測量，能夠快速得到數字正射影像圖、大比例尺地形圖和數字高程模型(Digital Elevation Model,DEM) 等數字化成果[3-5]，由于定位定向系統(Position and Orientation System，POS)獲得的坐標數據是通過全球導航衛星系統(Global Navigation Satellite System, GNSS)測量得到的，所獲得的大地高與正常高之間需要通過高程異常值進行轉換[6-8]。通常，求取高程異常的方法有多項式函數擬合、多面函數擬合和神經網絡擬合等[9-12]，但幾何擬合方法要求擬合區域內要有一定數量且分布合理的水準聯測GNSS點，該方法在水準點稀少的測區(尤其是山區)實施困難且山區高程異常變化較大時精度會受到較大影響[13-14]。任超等[15]使用不同組合方法對高程異常擬合進行討論，發現WLSSVM組合和非負變權組合的擬合效果較好，轉換精度較高。劉斌等[16]提出了一種基于EGM2008地球重力場模型和地形變化的影響，并結合二次曲面函數來進行高程擬合，該方法能提高GPS高程擬合精度。馬雷等[17]使用EGM2008，EIGEN6C4，GOCO05S三種地球重力場模型與克里金插值和一些函數擬合方法來進行精度分析，證明了地球重力場模型與數學函數相結合的“移去—恢復”法的優勢。然而，對于不同幾何模型與不同地球重力場模型的結合效果未進一步研究，也未將其應用到無人機攝影測量領域。本文采用平面函數、二次曲面函數和多面函數3種幾何模型以及EGM2008，EIGEN6C2，EIGEN6C4三種重力場模型。利用已知數據，通過設計不同組合方案進行分析和對比，得到了精度較高的組合模型，并將該模型應用于無人機攝影測量，驗證模型得到的高程異常精度。結果表明，EIGEN6C4地球重力場模型與二次曲面函數相結合的擬合模型能較好地提高無人機攝影測量成果的高程轉換精度，給工程應用提供了一定的參考價值。

1 重力場模型與幾何模型結合的高程異常擬合算法

1.1 高程異常與重力場模型的關系

通過Bruns公式來計算地球重力場與高程異常二者之間的關系：

Sn,msin(mλ))Pn,m(sinθ) ，

(1)

式中，ζG為該點用地球重力場模型計算得到的高程異常值；(ρ,θ,λ)分別為地心向徑、地心緯度和地心經度；GM為引力常數和地球質量的乘積；γ為待定點的正常重力值；a為參考橢球的長半軸；Cn,m,Sn,m為完全規格化得到的位系數；Pn,m(sinθ)為完全規格化的Laggardness函數；N為模型的最高展開階數。

1.2 重力場模型

隨著測量技術的不斷發展，不同地球重力場模型的分辨率和精度持續提高。本文選用了廣泛使用的3種典型地球重力場模型進行實驗分析。

① EGM2008地球重力場模型：由美國國家地理空間情報局2008年研制的全球超高階重力場模型，采用WGS-84坐標系，階次達到2 190。相比于以前的低階重力場模型，EGM2008重力場模型精度有著大幅度提高[16]。在一些精度要求不高的場合，可以使用它來獲得相應的數據。

② EIGEN地球重力場模型：EIGEN6C2重力場模型是德國地學中心2012年發布的1 949階的重力場模型；EIGEN6C4重力場模型是德國地學中心在2014年發布的2 190階的超高階重力場模型[13]。該模型的中長波精度有較大提升，是目前精度最高的重力場模型。

1.3 “移去—恢復”法原理

將高程異常分為3部分：長波項、中波項和短波項。短波項是地形起伏引起的，需要結合測區的數字高程模型來進行地形改正。但在實際應用過程中數字高程模型通常難以獲得，因此一般不進行地形改正。可以將高程異常分為2部分：

ζ=ζG+Δζ，

(2)

式中，ζG為用地球重力場模型計算得到的高程異常值；Δζ為實際測量得到的高程異常與地球重力場模型計算得到的高程異常殘差值。

“移去—恢復”法的具體實現過程與技術路線如下：

第1步“移去”：測區中選擇n個GNSS點，其中m個GNSS點含有大地高H與正常高Hγ，可以求解得到m個水準點的高程異常值ζ=H-Hγ，通過相應的地球重力場模型計算得到m個GNSS點的近似高程異常值ζG，m個GNSS點的高程異常殘差項Δζ(Δζ=ζ-ζG)；

第2步“擬合”：對第1步求得m個GNSS點的Δζ，利用常規的幾何函數對此殘差進行擬合，求解得到相應模型的擬合系數，將剩余點的Δζ利用上述的幾何函數進行擬合求解；

第3步“恢復”：將第1步利用地球重力場模型解算得到待定點的ζG與第2步中擬合求得待定點的Δζ相加，求得待求點的最終高程異常值ζ=ζG+Δζ，得到各待求點的正常髙。

2 實例分析

2.1 高程轉換實例方案

以文獻[18]中的丘陵地區測量數據作為本文的實驗數據。該測區西北方向地勢較高，東南方向地勢較低。測區內部以海拔650～850 m的丘陵為主，還有一些淺丘和平壩，其中海拔850 m以上的高山也占了很大的比例。實驗選取了測區中17個GNSS點作為實驗點，其中有3個國家一等水準點，其余點以一等水準點為基準點，按照國家四等水準測量規范進行水準測量，得到質量可靠的正常高數據。名山測區控制點分布示意如圖1所示。實驗結果采用中誤差、內符合精度和外符合精度進行轉換精度評定。

圖1 名山測區控制點分布示意Fig.1 Distribution of control points in Mingshan

2.2 方案1

選取測區內均勻分布的6個控制點作為已知點進行實驗，其余11個點作為檢核點。采用不同的擬合函數與不同的重力場模型進行擬合，包括幾何函數擬合、平面函數擬合(Planar)、二次曲面擬合(Quadratic)和多面函數擬合(Polyhedral)；顧及EGM2008地球重力場模型進行高程轉換；顧及EIGEN6C2地球重力場模型進行高程轉換；顧及EIGEN6C4地球重力場模型進行高程轉換。精度分析匯總如表1所示。點位分布如圖2所示。殘差如圖3～圖6所示。

表1 方案1精度分析匯總Tab.1 Summary of precision analysis of scheme 1 單位：m

圖2 方案1點位分布Fig.2 Point distribution of scheme 1

圖3 幾何函數擬合殘差Fig.3 Residual graph of geometric function fitting

圖4 EGM2008地球重力場的殘差Fig.4 Residual graph of EGM2008 earth gravity field

圖5 EIGEN6C2地球重力場的殘差Fig.5 Residual graph of EIGEN6C2 earth gravity field

圖6 EIGEN6C4地球重力場的殘差Fig.6 Residual graph of EIGEN6C4 earth gravity field

由表1可以看出，基于二次曲面函數擬合法的精度較其余2種擬合方案高，中誤差為0.025 m，外符合精度為0.032 m，平均值為0.005 m。通過比較基于不同地球重力場模型的“移去—恢復”法的高程轉換精度可以看出，基于3種不同地球重力場模型的高程轉換結果中殘差的最大值、最小值、中誤差、平均值以及外符合精度都比單一的幾何函數擬合得到的結果小。通過殘差圖同樣可以看出，得到的精度較高。在3種不同的重力場模型之間進行比較發現，使用EIGEN6C2重力場模型的精度較差，使用EGM2008重力場模型的精度與EIGEN6C4重力場模型的擬合精度相當，但基于EIGEN6C4地球重力場模型進行高程轉換結果的擬合殘差值中平均誤差較小，得到的精度更加可靠，更適合用作高程異常擬合計算。顧及EIGEN6C4重力場的二次曲面函數擬合的中誤差為0.023 m，外符合精度0.030 m，平均值0.003 m，約為單一的二次曲面函數擬合的平均值的1/2，得到的結果更加可靠。

2.3 方案2

選取測區內均勻分布的9個已知高程異常的點進行擬合函數系數計算，其余8個點作為檢核點。采用的擬合函數和重力場模型同方案1。點位分布如圖7所示。將上述4個方法計算得到的高程異常殘差值進行精度分析，不同方法的二次曲面擬合殘差如圖8所示，精度分析如表2所示。

圖7 方案2點位分布Fig.7 Point distribution of scheme 2

圖8 幾種不同方法的二次曲面擬合殘差Fig.8 Residual graph of quadric surface fitting with several different methods

表2 方案2精度分析Tab.2 Summary of precision analysis of scheme 2 單位：m

當選取測區內數量較多、分布均勻的已知點時，使用幾何函數擬合中的二次曲面擬合的中誤差能夠達到0.022 4 m，能較好地進行高程異常擬合。顧及不同地球重力場模型和二次曲面擬合相結合時的精度對比其他2種方法的擬合精度仍有較大提升，其他2種函數擬合的精度可以達到與二次曲面函數相當的水平，中誤差均縮小到了3 cm之內，相對幾何函數中的7 cm誤差，是一個較大的提升。由圖8可以看出，顧及地球重力場模型的二次曲面擬合比單一的二次曲面擬合波動和殘差平均值都較小，這也說明了考慮地球重力場模型的高程轉換精度有一定的提升。

綜上所述，通過以上2個方案的結果可以看出，在地形起伏較大的山區進行高程異常擬合時，選用二次曲面函數和EIGEN6C4地球重力場模型相結合的擬合方法來進行高程轉換計算是最佳的選擇。

3 高程轉換在無人機攝影測量中應用

實驗使用大疆御Mavic Pro無人機進行影像數據采集，測區位于廣西壯族自治區欽州市石甌山水庫；覆蓋面積約為1.4 km2，測區以小山丘為主，最大高差約為100 m，無人機可視性良好。在內業數據處理中，利用Pix4Dmapper軟件進行影像處理，通過ArcGIS進行點云數據的處理，得到數字高程模型。選取的點位位置分布如圖9所示。

圖9 選取的點位位置分布Fig.9 Distribution of selected points

為分析高程轉換精度，通過正射影像選取具有代表性的特征點進行處理。本文選用測區中6個像控點作為已知點，選擇8個特征點進行精度評定。分別采用3種不同的方法進行高程異常計算和精度評定：① 通過Pix4Dmapper軟件與GNSS測得的像控點坐標數據經軟件直接處理，得到特征點的正常高；② 利用二次曲面函數根據像控點的坐標數據來求解特征點的高程異常并獲得正常高；③ 利用上述顧及EIGEN6C4地球重力場模型與二次曲面函數擬合的方法，求解特征點的高程異常并獲得正常高。點位與不同方案得到的高程異常如10所示。

由圖10可以看出，直接使用Pix4Dmapper軟件進行影像處理，得到的特征點正常高與真實值偏差較大、精度較低且波動較大。使用函數擬合的方法進行特征點高程轉換時，8個特征點的高程轉換精度有明顯提升，誤差較小，誤差波動較平緩。

圖10 不同方案得到的高程異常Fig.10 Height anomalies obtained by different schemes

無人機特征點高程轉換分析如表3所示。

表3 無人機特征點高程轉換分析Tab.3 Height conversion analysis of UAV feature points 單位：m

由表3可以看出，直接使用軟件處理得到的高程異常中誤差為22 cm，使用二次曲面函數進行高程轉換的中誤差為3.73 cm，顧及EIGEN6C4地球重力場模型時，中誤差進一步縮小為1.154 cm。由此可知，顧及EIGEN6C4重力場模型進行擬合的結果更貼近真實值，能夠提高高程異常內插的精度。因此，在地形起伏較大的測區進行無人機攝影測量和內業數據處理時，通過應用上述相關的高程轉換方法，能夠提高無人機攝影測量計算得到正常高的精度。

4 結束語

綜合上述實驗分析，可以得到以下結論：

① 進行高程異常擬合時，應該結合實際測區情況采用多種方案進行比較分析，得到最佳的擬合模型。綜合來說，當測區的已知點分布較均勻而且已知點個數較多時，通常選用二次曲面函數擬合。

② 在顧及地球重力場模型的轉換方法中，相比于幾何函數擬合法，精度和可靠性均有一定提升，充分利用GNSS的高程信息和地球重力場模型提供的中長波部分，可以獲得較大的精度提升，較大地減少水準測量的外業工作，提高作業效率。

③ 在無人機攝影測量數據處理中，相對于軟件直接處理輸出得到的正常高，使用高程轉換模型得到的正常高有明顯的精度提升。說明了顧及地球重力場模型的高程異常內插在無人機攝影測量數據處理中的可行性，對提高無人機攝影測量正常高的精度提供了一定的思路。

本文采用平面函數、二次曲面函數和多面函數3種幾何模型和3種重力場模型進行高程異常的計算，未來可在此基礎上結合機器學習等算法進一步實驗。