利用EGM2008模型與加權(quán)組合模型進(jìn)行高程異常擬合

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

(1. 中國(guó)礦業(yè)大學(xué)環(huán)境與測(cè)繪學(xué)院， 江蘇 徐州 221116； 2. 紹興市城市規(guī)劃測(cè)繪院，浙江 紹興312000)

高程異常擬合是GPS高程測(cè)量應(yīng)用于實(shí)際工作的重要技術(shù)手段，高精度地求解高程異常是技術(shù)關(guān)鍵。對(duì)于這一熱點(diǎn)問題，國(guó)內(nèi)外學(xué)者進(jìn)行了大量的研究探索，目前主要的方法有數(shù)學(xué)模型擬合法、利用地球重力場(chǎng)模型求高程異常和神經(jīng)網(wǎng)絡(luò)方法。這些方法各有優(yōu)點(diǎn)，但是對(duì)應(yīng)用的地形和分布情況有一定要求，不能滿足復(fù)雜的環(huán)境。

針對(duì)這一問題，國(guó)內(nèi)外學(xué)者進(jìn)行了大量的研究探索。邢志斌等[1]提出了快速構(gòu)建利用垂線偏差計(jì)算高程異常差法方程的方法，使高程異常差計(jì)算更加高效；劉斌等[2]提出了利用EGM2008重力場(chǎng)模型和地形改正進(jìn)行高程異常計(jì)算的方法，該方法能夠達(dá)到較高的精度但是不能適應(yīng)較為復(fù)雜的應(yīng)用環(huán)境；黎劍[3]提出了一種加權(quán)綜合模型高程異常擬合方法，該方法兼顧多種方法的優(yōu)勢(shì)具有良好的適應(yīng)性，但是沒有引入EGM2008重力場(chǎng)模型進(jìn)行改正導(dǎo)致高程異常擬合精度不高。

本文提出一種顧及EGM2008重力場(chǎng)模型的加權(quán)組合模型高程異常擬合方法，該方法具有較高的擬合精度并且能夠適應(yīng)較為復(fù)雜的應(yīng)用環(huán)境。通過實(shí)際數(shù)據(jù)驗(yàn)證可知，該方法的解算精度能夠達(dá)到厘米級(jí)，使程異常擬合的效果提升10%～30%。

1 原理與方法

1.1 原理與流程

本文通過“移去—恢復(fù)”法求解高程異常，該方法可以根據(jù)EGM2008重力場(chǎng)模型和聯(lián)測(cè)少量水準(zhǔn)點(diǎn)對(duì)GPS高程異常進(jìn)行擬合而不需要測(cè)量額外的重力參數(shù)，能夠適應(yīng)地形起伏較大且水準(zhǔn)聯(lián)測(cè)點(diǎn)較少的區(qū)域，具有高效率和較好的適用性。根據(jù)物理大地測(cè)量學(xué)理論，可以得出GPS高程擬合的表達(dá)式[4]為

ζ=ζGM+ζCT+ζRES

(1)

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

ζ=ζGM+ζ0

(2)

式中，ζ0表示高程異常的短波項(xiàng)和中波項(xiàng)，統(tǒng)稱殘余高程異常值。“移去—恢復(fù)”法的核心思想為，將已知點(diǎn)高程異常中通過EGM2008重力場(chǎng)模型求解的長(zhǎng)波項(xiàng)ζGM移去，剩下殘余高程異常值ζ0運(yùn)用加權(quán)組合模型對(duì)其進(jìn)行擬合；通過加權(quán)組合模型對(duì)未知點(diǎn)的殘余高程異常ζ0進(jìn)行求解，再加上高程異常長(zhǎng)波項(xiàng)ζGM就可以得到未知點(diǎn)的高程異常ζ。具體的解算流程如圖1所示。

圖1 “移去—恢復(fù)”法

1.2 EGM2008重力場(chǎng)模型求解高程異常長(zhǎng)波項(xiàng)

EGM2008全球重力場(chǎng)模型是由美國(guó)國(guó)家地理空間情報(bào)局，通過先進(jìn)建模技術(shù)與算法，結(jié)合測(cè)高數(shù)據(jù)和重力數(shù)據(jù)，得到拓展至2190次球諧系數(shù)的全球超高階地球重力場(chǎng)模型。EGM2008全球重力場(chǎng)模型的覆蓋率高達(dá)83.8%，模型分辨率達(dá)到5′，全球5′×5′網(wǎng)格的估算精度最小誤差為0.304 5 dm，平均誤差為1.114 dm。

利用EGM2008重力場(chǎng)模型求解地球上任意一點(diǎn)的GPS高程異常長(zhǎng)波項(xiàng)ζGM，首先要計(jì)算重力擾動(dòng)位，再根據(jù)Burns公式計(jì)算ζGM。Burns[4]公式如下

(3)

1.3 加權(quán)組合模型

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

圖2 加權(quán)組合模型算法

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

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

(4)

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

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

(5)

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

(6)

式中，δ為光滑系數(shù)，系數(shù)的取值根據(jù)實(shí)際情況而定，也就是說若想獲得一個(gè)測(cè)區(qū)的最佳多面函數(shù)模型就需要對(duì)δ的取值作出準(zhǔn)確判斷，通常是取一系列數(shù)據(jù)分析模型精度選取最佳光滑系數(shù)值。本文取多面函數(shù)模型的光滑系數(shù)δ=5。

加權(quán)組合模型的公式為[3]

(7)

2 算例分析

2.1 試驗(yàn)區(qū)及試驗(yàn)分組介紹

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

試驗(yàn)區(qū)內(nèi)的高程點(diǎn)的坐標(biāo)值、高程異常值(ζ)、高程異常長(zhǎng)波項(xiàng)(ζGM)及殘余高程異常值(ζRES)見表1。

圖3 點(diǎn)位分布

圖4 高程異常曲面

點(diǎn)號(hào)緯度(B)/(° ′ ″)經(jīng)度(L)/(° ′ ″)高程異常(ζ)/m高程異常長(zhǎng)波項(xiàng)(ζ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重力場(chǎng)模型的二次曲面擬合方法、多面函數(shù)擬合方法及其加權(quán)組合模型對(duì)試驗(yàn)區(qū)內(nèi)的數(shù)據(jù)進(jìn)行GPS高程擬合，然后采用顧及EGM2008重力場(chǎng)模型的二次曲面擬合方法、多面函數(shù)擬合方法及其加權(quán)組合模型對(duì)試驗(yàn)區(qū)內(nèi)的數(shù)據(jù)進(jìn)行GPS高程擬合。具體試驗(yàn)?zāi)Ｐ秃头椒ㄟx擇及其分組編號(hào)見表2。

表2 試驗(yàn)?zāi)Ｐ团c方法組合及其編號(hào)

2.2 試驗(yàn)結(jié)果及分析

本文根據(jù)間接平差原理進(jìn)行算法設(shè)計(jì)并實(shí)現(xiàn)，通過算法對(duì)6種組合試驗(yàn)進(jìn)行驗(yàn)證，計(jì)算得到高程殘差值見表3。

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

表3 驗(yàn)證點(diǎn)高程異常擬合殘差 m

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

圖5 EGM2008重力場(chǎng)模型對(duì)高程擬合的影響

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

圖6 加權(quán)組合模型對(duì)高程擬合的影響

m

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

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

2.3 工程實(shí)例

將本方法應(yīng)用于廣東某市城市測(cè)量工程中，通過已知高程點(diǎn)的水準(zhǔn)高程和GPS高程，對(duì)測(cè)區(qū)內(nèi)的高程點(diǎn)進(jìn)行高程擬合計(jì)算。計(jì)算結(jié)果見表5。分析可知，加權(quán)組合模型方法達(dá)到GPS高程測(cè)量精度水準(zhǔn)，滿足GPS高程的要求并且具有較好的實(shí)用性。

表5 應(yīng)用實(shí)例數(shù)據(jù) m

3 結(jié) 語

通過引入EGM2008重力場(chǎng)模型能夠計(jì)算高程異常長(zhǎng)波項(xiàng)的值，在高程異常擬合計(jì)算中可以移出高程異常長(zhǎng)波項(xiàng)的值，從而提高高程異常擬合的精度。加權(quán)組合模型可以兼顧二次曲面擬合方法和多面函數(shù)擬合方法的優(yōu)點(diǎn)，既能保證高程異常整體的平滑又能顧及高程異常極小誤差，保證高程異常結(jié)果精度。顧及EGM2008重力場(chǎng)模型的加權(quán)組合模型擬合方法能夠顯著提升高程異常擬合的精度，擬合效果可以提升10%～30%。工程實(shí)例的應(yīng)用說明，加權(quán)組合模型能夠滿足GPS高程測(cè)量精度要求，具有較好的實(shí)用性。

[1] 邢志斌，李?yuàn)檴?，王偉，?利用垂線偏差計(jì)算高程異常差法方程的快速構(gòu)建方法[J].武漢大學(xué)學(xué)報(bào)(信息科學(xué)版)，2016，41(6)：778-783.

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

[3] 黎劍.區(qū)域GPS高程異常擬合及建模方法研究[D].昆明：昆明理工大學(xué)，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高程分區(qū)擬合轉(zhuǎn)換正常高的研究[J].武漢大學(xué)學(xué)報(bào)(信息科學(xué)版)，2004，29(10)：908-911.

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

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

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

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

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

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

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