顧及地球重力場(chǎng)模型的高程轉(zhuǎn)換方法

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

(1. 桂林理工大學(xué) 測(cè)繪地理信息學(xué)院，廣西 桂林 541006； 2. 廣西空間信息與測(cè)繪重點(diǎn)實(shí)驗(yàn)室，廣西 桂林 541006； 3. 中國(guó)地質(zhì)大學(xué)(武漢)公共管理學(xué)院，湖北 武漢 430074; 4. 山東省地質(zhì)測(cè)繪院，山東 濟(jì)南 250013)

0 引言

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

1 重力場(chǎng)模型與幾何模型結(jié)合的高程異常擬合算法

1.1 高程異常與重力場(chǎng)模型的關(guān)系

通過(guò)Bruns公式來(lái)計(jì)算地球重力場(chǎng)與高程異常二者之間的關(guān)系：

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

(1)

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

1.2 重力場(chǎng)模型

隨著測(cè)量技術(shù)的不斷發(fā)展，不同地球重力場(chǎng)模型的分辨率和精度持續(xù)提高。本文選用了廣泛使用的3種典型地球重力場(chǎng)模型進(jìn)行實(shí)驗(yàn)分析。

① EGM2008地球重力場(chǎng)模型：由美國(guó)國(guó)家地理空間情報(bào)局2008年研制的全球超高階重力場(chǎng)模型，采用WGS-84坐標(biāo)系，階次達(dá)到2 190。相比于以前的低階重力場(chǎng)模型，EGM2008重力場(chǎng)模型精度有著大幅度提高[16]。在一些精度要求不高的場(chǎng)合，可以使用它來(lái)獲得相應(yīng)的數(shù)據(jù)。

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

1.3 “移去—恢復(fù)”法原理

將高程異常分為3部分：長(zhǎng)波項(xiàng)、中波項(xiàng)和短波項(xiàng)。短波項(xiàng)是地形起伏引起的，需要結(jié)合測(cè)區(qū)的數(shù)字高程模型來(lái)進(jìn)行地形改正。但在實(shí)際應(yīng)用過(guò)程中數(shù)字高程模型通常難以獲得，因此一般不進(jìn)行地形改正?？梢詫⒏叱坍惓７譃?部分：

ζ=ζG+Δζ，

(2)

式中，ζG為用地球重力場(chǎng)模型計(jì)算得到的高程異常值；Δζ為實(shí)際測(cè)量得到的高程異常與地球重力場(chǎng)模型計(jì)算得到的高程異常殘差值。

“移去—恢復(fù)”法的具體實(shí)現(xiàn)過(guò)程與技術(shù)路線如下：

第1步“移去”：測(cè)區(qū)中選擇n個(gè)GNSS點(diǎn)，其中m個(gè)GNSS點(diǎn)含有大地高H與正常高Hγ，可以求解得到m個(gè)水準(zhǔn)點(diǎn)的高程異常值ζ=H-Hγ，通過(guò)相應(yīng)的地球重力場(chǎng)模型計(jì)算得到m個(gè)GNSS點(diǎn)的近似高程異常值ζG，m個(gè)GNSS點(diǎn)的高程異常殘差項(xiàng)Δζ(Δζ=ζ-ζG)；

第2步“擬合”：對(duì)第1步求得m個(gè)GNSS點(diǎn)的Δζ，利用常規(guī)的幾何函數(shù)對(duì)此殘差進(jìn)行擬合，求解得到相應(yīng)模型的擬合系數(shù)，將剩余點(diǎn)的Δζ利用上述的幾何函數(shù)進(jìn)行擬合求解；

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

2 實(shí)例分析

2.1 高程轉(zhuǎn)換實(shí)例方案

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

圖1 名山測(cè)區(qū)控制點(diǎn)分布示意Fig.1 Distribution of control points in Mingshan

2.2 方案1

選取測(cè)區(qū)內(nèi)均勻分布的6個(gè)控制點(diǎn)作為已知點(diǎn)進(jìn)行實(shí)驗(yàn)，其余11個(gè)點(diǎn)作為檢核點(diǎn)。采用不同的擬合函數(shù)與不同的重力場(chǎng)模型進(jìn)行擬合，包括幾何函數(shù)擬合、平面函數(shù)擬合(Planar)、二次曲面擬合(Quadratic)和多面函數(shù)擬合(Polyhedral)；顧及EGM2008地球重力場(chǎng)模型進(jìn)行高程轉(zhuǎn)換；顧及EIGEN6C2地球重力場(chǎng)模型進(jìn)行高程轉(zhuǎn)換；顧及EIGEN6C4地球重力場(chǎng)模型進(jìn)行高程轉(zhuǎn)換。精度分析匯總?cè)绫?所示。點(diǎn)位分布如圖2所示。殘差如圖3～圖6所示。

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

圖2 方案1點(diǎn)位分布Fig.2 Point distribution of scheme 1

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

圖4 EGM2008地球重力場(chǎng)的殘差Fig.4 Residual graph of EGM2008 earth gravity field

圖5 EIGEN6C2地球重力場(chǎng)的殘差Fig.5 Residual graph of EIGEN6C2 earth gravity field

圖6 EIGEN6C4地球重力場(chǎng)的殘差Fig.6 Residual graph of EIGEN6C4 earth gravity field

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

2.3 方案2

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

圖7 方案2點(diǎn)位分布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

當(dāng)選取測(cè)區(qū)內(nèi)數(shù)量較多、分布均勻的已知點(diǎn)時(shí)，使用幾何函數(shù)擬合中的二次曲面擬合的中誤差能夠達(dá)到0.022 4 m，能較好地進(jìn)行高程異常擬合。顧及不同地球重力場(chǎng)模型和二次曲面擬合相結(jié)合時(shí)的精度對(duì)比其他2種方法的擬合精度仍有較大提升，其他2種函數(shù)擬合的精度可以達(dá)到與二次曲面函數(shù)相當(dāng)?shù)乃?，中誤差均縮小到了3 cm之內(nèi)，相對(duì)幾何函數(shù)中的7 cm誤差，是一個(gè)較大的提升。由圖8可以看出，顧及地球重力場(chǎng)模型的二次曲面擬合比單一的二次曲面擬合波動(dòng)和殘差平均值都較小，這也說(shuō)明了考慮地球重力場(chǎng)模型的高程轉(zhuǎn)換精度有一定的提升。

綜上所述，通過(guò)以上2個(gè)方案的結(jié)果可以看出，在地形起伏較大的山區(qū)進(jìn)行高程異常擬合時(shí)，選用二次曲面函數(shù)和EIGEN6C4地球重力場(chǎng)模型相結(jié)合的擬合方法來(lái)進(jìn)行高程轉(zhuǎn)換計(jì)算是最佳的選擇。

3 高程轉(zhuǎn)換在無(wú)人機(jī)攝影測(cè)量中應(yīng)用

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

圖9 選取的點(diǎn)位位置分布Fig.9 Distribution of selected points

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

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

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

無(wú)人機(jī)特征點(diǎn)高程轉(zhuǎn)換分析如表3所示。

表3 無(wú)人機(jī)特征點(diǎn)高程轉(zhuǎn)換分析Tab.3 Height conversion analysis of UAV feature points 單位：m

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

4 結(jié)束語(yǔ)

綜合上述實(shí)驗(yàn)分析，可以得到以下結(jié)論：

① 進(jìn)行高程異常擬合時(shí)，應(yīng)該結(jié)合實(shí)際測(cè)區(qū)情況采用多種方案進(jìn)行比較分析，得到最佳的擬合模型。綜合來(lái)說(shuō)，當(dāng)測(cè)區(qū)的已知點(diǎn)分布較均勻而且已知點(diǎn)個(gè)數(shù)較多時(shí)，通常選用二次曲面函數(shù)擬合。

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

③ 在無(wú)人機(jī)攝影測(cè)量數(shù)據(jù)處理中，相對(duì)于軟件直接處理輸出得到的正常高，使用高程轉(zhuǎn)換模型得到的正常高有明顯的精度提升。說(shuō)明了顧及地球重力場(chǎng)模型的高程異常內(nèi)插在無(wú)人機(jī)攝影測(cè)量數(shù)據(jù)處理中的可行性，對(duì)提高無(wú)人機(jī)攝影測(cè)量正常高的精度提供了一定的思路。

本文采用平面函數(shù)、二次曲面函數(shù)和多面函數(shù)3種幾何模型和3種重力場(chǎng)模型進(jìn)行高程異常的計(jì)算，未來(lái)可在此基礎(chǔ)上結(jié)合機(jī)器學(xué)習(xí)等算法進(jìn)一步實(shí)驗(yàn)。