利用EGM2008模型與加權組合模型進行高程異常擬合

羅陶榮，王中元，梁 寧，馬效申

(1. 中國礦業大學環境與測繪學院， 江蘇 徐州 221116； 2. 紹興市城市規劃測繪院，浙江 紹興312000)

高程異常擬合是GPS高程測量應用于實際工作的重要技術手段，高精度地求解高程異常是技術關鍵。對于這一熱點問題，國內外學者進行了大量的研究探索，目前主要的方法有數學模型擬合法、利用地球重力場模型求高程異常和神經網絡方法。這些方法各有優點，但是對應用的地形和分布情況有一定要求，不能滿足復雜的環境。

針對這一問題，國內外學者進行了大量的研究探索。邢志斌等[1]提出了快速構建利用垂線偏差計算高程異常差法方程的方法，使高程異常差計算更加高效；劉斌等[2]提出了利用EGM2008重力場模型和地形改正進行高程異常計算的方法，該方法能夠達到較高的精度但是不能適應較為復雜的應用環境；黎劍[3]提出了一種加權綜合模型高程異常擬合方法，該方法兼顧多種方法的優勢具有良好的適應性，但是沒有引入EGM2008重力場模型進行改正導致高程異常擬合精度不高。

本文提出一種顧及EGM2008重力場模型的加權組合模型高程異常擬合方法，該方法具有較高的擬合精度并且能夠適應較為復雜的應用環境。通過實際數據驗證可知，該方法的解算精度能夠達到厘米級，使程異常擬合的效果提升10%～30%。

1 原理與方法

1.1 原理與流程

本文通過“移去—恢復”法求解高程異常，該方法可以根據EGM2008重力場模型和聯測少量水準點對GPS高程異常進行擬合而不需要測量額外的重力參數，能夠適應地形起伏較大且水準聯測點較少的區域，具有高效率和較好的適用性。根據物理大地測量學理論，可以得出GPS高程擬合的表達式[4]為

ζ=ζGM+ζCT+ζRES

(1)

式中，ζGM是由EGM2008重力場模型求得的高程異常長波項；ζRES是殘余高程異常，為高程異常的中波項；ζCT是通過地形改正解算獲得的高程異常短波項。本文引入EGM2008重力場模型來求解高程異常中的長波項，對于其中的短波項和中波項通過加權組合模型進行求解。故將式(1)變換成[5]

ζ=ζGM+ζ0

(2)

式中，ζ0表示高程異常的短波項和中波項，統稱殘余高程異常值?！耙迫ァ謴汀狈ǖ暮诵乃枷霝?，將已知點高程異常中通過EGM2008重力場模型求解的長波項ζGM移去，剩下殘余高程異常值ζ0運用加權組合模型對其進行擬合；通過加權組合模型對未知點的殘余高程異常ζ0進行求解，再加上高程異常長波項ζGM就可以得到未知點的高程異常ζ。具體的解算流程如圖1所示。

圖1 “移去—恢復”法

1.2 EGM2008重力場模型求解高程異常長波項

EGM2008全球重力場模型是由美國國家地理空間情報局，通過先進建模技術與算法，結合測高數據和重力數據，得到拓展至2190次球諧系數的全球超高階地球重力場模型。EGM2008全球重力場模型的覆蓋率高達83.8%，模型分辨率達到5′，全球5′×5′網格的估算精度最小誤差為0.304 5 dm，平均誤差為1.114 dm。

利用EGM2008重力場模型求解地球上任意一點的GPS高程異常長波項ζGM，首先要計算重力擾動位，再根據Burns公式計算ζGM。Burns[4]公式如下

(3)

1.3 加權組合模型

本文根據黎劍的加權綜合模型提出了一種權重動態變化的加權組合模型，該模型能夠根據兩種高程擬合方法的擬合效果動態確定兩種方法的權重，具有更好的擬合精度和計算效率。本文選用二次曲面擬合方法與多面函數擬合方法作為組合模型的基本擬合單元，動態分配兩種擬合單元的權重，最后得到加權后的GPS高程異常模型，如圖2所示。

圖2 加權組合模型算法

二次曲面擬合方法的公式為[6]

ξ(x,y)=a0+a1Δx+a2Δy+a3Δx2+a4Δy2+a5ΔxΔy

(4)

式中，ai(i=0,1,2,…,5)為二次曲面擬合模型的系數；Δx=x-x0；Δy=y-y0；(x0，y0)一般為重心點坐標。

多面函數擬合方法的公式為[7-8]

(5)

式中，aj(j=1,2,3,…,n)為各個核函數的比例系數，是需要求取的關鍵參數；Q(x,y,xj,yj)為核函數。本文選用了對稱性較好的正雙曲面函數作為核函數[8]

(6)

式中，δ為光滑系數，系數的取值根據實際情況而定，也就是說若想獲得一個測區的最佳多面函數模型就需要對δ的取值作出準確判斷，通常是取一系列數據分析模型精度選取最佳光滑系數值。本文取多面函數模型的光滑系數δ=5。

加權組合模型的公式為[3]

(7)

2 算例分析

2.1 試驗區及試驗分組介紹

本文以江蘇省某市C級GPS測量成果和一等水準測量成果為研究數據對提出的高程擬合新方法進行實際數據的驗證。該試驗區長21 km，寬20 km，點位均勻分布，具體分布如圖3所示。由圖3可知，共有30個已知點，其中25個作為已知點(三角形表示)，最大高程為481.810 0 m，最小高程為260.915 8 m，平均高程為372.785 0 m；5個作為檢核點(圓形表示)，最大高程為448.320 3 m，最小高程為75.327 7 m，平均高程為222.848 4 m。高程異常曲面如圖4所示。

試驗區內的高程點的坐標值、高程異常值(ζ)、高程異常長波項(ζGM)及殘余高程異常值(ζRES)見表1。

圖3 點位分布

圖4 高程異常曲面

點號緯度(B)/(° ′ ″)經度(L)/(° ′ ″)高程異常(ζ)/m高程異常長波項(ζGM)/m殘余高程異常(ζRES)/m1334．280473683116．2126554108．9989．257-0．2592034．243935403116．2539383018．6248．939-0．3152234．261696960116．2346416808．7879．105-0．3182434．262193619116．1907406089．1489．386-0．2382634．282952237116．1843342529．2089．43-0．222

本文采用未顧及EGM2008重力場模型的二次曲面擬合方法、多面函數擬合方法及其加權組合模型對試驗區內的數據進行GPS高程擬合，然后采用顧及EGM2008重力場模型的二次曲面擬合方法、多面函數擬合方法及其加權組合模型對試驗區內的數據進行GPS高程擬合。具體試驗模型和方法選擇及其分組編號見表2。

表2 試驗模型與方法組合及其編號

2.2 試驗結果及分析

本文根據間接平差原理進行算法設計并實現，通過算法對6種組合試驗進行驗證，計算得到高程殘差值見表3。

由表3可知，高程異常擬合的精度達到了厘米級，滿足GPS定位對于高程的要求。未顧及EGM2008重力場模型時3種方法擬合的精度：二次曲面擬合方法擬合結果中誤差1.39 cm，平均殘差0.94 cm；多面函數擬合方法擬合結果中誤差9.16 cm，平均殘差6.77 cm；加權組合模型擬合方法擬合結果中誤差1.1 cm，平均殘差0.83 cm。顧及EGM2008重力場模型時3種方法擬合的精度：二次曲面擬合方法擬合結果中誤差0.62 cm，平均殘差0.5 cm；多面函數擬合方法擬合結果中誤差2.61 cm，平均殘差2.4 cm；加權組合模型擬合方法擬合結果中誤差0.51 cm，平均殘差0.38 cm。

表3 驗證點高程異常擬合殘差 m

顧及EGM2008重力場模型時3種擬合方法的擬合精度明顯好于未顧及EGM2008重力場模型時3種擬合方法的擬合精度，通過引入EGM2008重力場模型能夠有效提高高程異常的擬合精度。顧及EGM2008重力場模型和未顧及EGM2008重力場模型的擬合效果如圖5所示。

圖5 EGM2008重力場模型對高程擬合的影響

由圖5可知，在顧及EGM2008重力場模型時二次曲面擬合法和加權組合模型擬合法的高程異常擬合殘差整體要小于未顧及EGM2008重力場模型時高程異常擬合殘差，且高程異常殘差的離散程度也較小；多面函數擬合方法的高程異常擬合殘差整體要大于未顧及EGM2008重力場模型時高程異常擬合殘差，且高程異常殘差相近。由圖6可知，加權組合模型擬合效果在顧及EGM2008重力場模型和未顧及EGM2008重力場模型時都優于二次曲面擬合方法和多面函數擬合方法。

圖6 加權組合模型對高程擬合的影響

m

注：提升率ρ=(α-β)/α，α為二次曲面法與多面函數法中最高精度，β為加權組合模型精度。

由表4可知，加權組合模型擬合效果整體可以提升10%～30%，具有很好的擬合效果和實際意義。

2.3 工程實例

將本方法應用于廣東某市城市測量工程中，通過已知高程點的水準高程和GPS高程，對測區內的高程點進行高程擬合計算。計算結果見表5。分析可知，加權組合模型方法達到GPS高程測量精度水準，滿足GPS高程的要求并且具有較好的實用性。

表5 應用實例數據 m

3 結 語

通過引入EGM2008重力場模型能夠計算高程異常長波項的值，在高程異常擬合計算中可以移出高程異常長波項的值，從而提高高程異常擬合的精度。加權組合模型可以兼顧二次曲面擬合方法和多面函數擬合方法的優點，既能保證高程異常整體的平滑又能顧及高程異常極小誤差，保證高程異常結果精度。顧及EGM2008重力場模型的加權組合模型擬合方法能夠顯著提升高程異常擬合的精度，擬合效果可以提升10%～30%。工程實例的應用說明，加權組合模型能夠滿足GPS高程測量精度要求，具有較好的實用性。

[1] 邢志斌，李姍姍，王偉，等.利用垂線偏差計算高程異常差法方程的快速構建方法[J].武漢大學學報(信息科學版)，2016，41(6)：778-783.

[2] 劉斌，郭際明，史俊波，等.利用EGM2008模型與地形改正進行GPS高程擬合[J].武漢大學學報(信息科學版)，2016，41(4)：554-558.

[3] 黎劍.區域GPS高程異常擬合及建模方法研究[D].昆明：昆明理工大學，2013.

[4] ITOH S，YASUDA I.Characteristics of Mesoscale Eddies in the Kuroshio-Oyashio Extension Region Detected from the Distribution of the Sea Surface Height Anomaly[J].Journal of Physical Oceanography，2010，40(5):1018-1034.

[5] SREEJITH K M，RAJESH S，MAJUMDAR T J，et al.High-resolution Residual Geoid and Gravity Anomaly Data of the Northern Indian Ocean——An Input to Geological Understanding[J].Journal of Asian Earth Sciences，2013，62(62):616-626.

[6] LI Feng，ZHOU X，LI Xin，et al.The Application of Ridge Estimation in GPS Height Fitting with Polynomial Surface[J].Coastal Engineering，2012.

[7] XU N，YAN Z，MA F.Multi-surface Function Model of GPS Height Fitting Based on Genetic Algorithm[J].Geotechnical Investigation & Surveying，2013.

[8] WU Y，ZHU J，ZUO T.The Application of Vondrak Filter Combined with Multisurface Function in Height Anomaly Fitting[J].Gohnal Nvgaon & Rvyng，2012，40(7):58-61.

[9] 高偉，徐紹銓.GPS高程分區擬合轉換正常高的研究[J].武漢大學學報(信息科學版)，2004，29(10)：908-911.

[10] 何海清，黃聲享，陳婷.一種顧及粗差的徑向神經網絡高程曲面擬合法[J].武漢大學學報(信息科學版)，2015，40(4):547-551.

[11] 唐健林，徐陶，龍盈.小波函數應用于GPS高程擬合的分析[J].測繪通報，2016(3)：58-60.DOI:10.13474/j.cnki.11 -2246.2016.0086.

[12] 田曉，鄭洪艷，許明元，等.一種改進的適用于不同地形的GPS高程擬合模型[J].測繪通報，2017(1):35-38.DOI:10.13474/j.cnki.11-2246.2017.0008.

[13] 周長志，孫佳龍，郭淑艷.基于移動-多面函數高程異常擬合方法[J].測繪通報，2016 (12)：25-27.DOI:10.13474/j.cnki.11-2246.2016.0394.

[14] 陳艷紅，馬全明，劉瑞敏，等.基于相對高程異常的GPS高程擬合方法探討[J].測繪通報,2014(12):67-69.

[15] 肖杰，張錦，鄧增兵，等.礦區似大地水準面精化方法研究[J].測繪通報,2015(2):14-18.DOI:10.13474/j.cnki.11-2246.2015.0032.

[16] 赫林，李建成.褚永海.1985國家高程基準與全球高程基準之間的垂直偏差[J].測繪學報，2016，45(7):768-774.DOI:10.11947/j.AGCS.2016.20160029.