梁雙鳳

(甘肅省測(cè)繪工程院，甘肅蘭州 730000 )

1 引 言

我國(guó)西部高原地區(qū)資源豐富，盛產(chǎn)天然氣、煤炭、石油以及地方特色土特產(chǎn)品。為了將這些資源和特產(chǎn)運(yùn)送出來(lái)，實(shí)現(xiàn)國(guó)家西部大開發(fā)戰(zhàn)略，國(guó)家做出了很多戰(zhàn)略部署。基礎(chǔ)測(cè)繪是戰(zhàn)略部署中非常重要的環(huán)節(jié)，但高原地區(qū)地形復(fù)雜，自然環(huán)境惡劣，人跡罕至，有廣袤的沙漠百里風(fēng)沙區(qū)，也有荒無(wú)人煙的戈壁灘“無(wú)人區(qū)”，給高原基礎(chǔ)測(cè)繪工作帶來(lái)了巨大的挑戰(zhàn)，特別是實(shí)施水準(zhǔn)測(cè)量非常困難。據(jù)有關(guān)報(bào)道，測(cè)繪生產(chǎn)單位開展的高原基礎(chǔ)測(cè)繪相關(guān)航空攝影測(cè)量?jī)?nèi)外業(yè)工作，包括構(gòu)建測(cè)量 C級(jí)GNSS網(wǎng)、布測(cè)像控點(diǎn)及像片外業(yè)調(diào)繪等，但水準(zhǔn)測(cè)量幾乎無(wú)法實(shí)施。GNSS測(cè)量得到的是大地高，而工程項(xiàng)目采用正常高程系統(tǒng)，因此，GNSS測(cè)量的大地高需要通過(guò)大地水準(zhǔn)模型轉(zhuǎn)換為正常高，才能用于工程項(xiàng)目?；A(chǔ)測(cè)繪航空攝影高程基準(zhǔn)的高程精度要求較高，為了解決高原基礎(chǔ)測(cè)繪高程基準(zhǔn)測(cè)量困難的問(wèn)題，需要一種更高效的獲取高精度像控點(diǎn)高程值的方法。

國(guó)內(nèi)外學(xué)者對(duì)GNSS高程轉(zhuǎn)換理論和實(shí)踐進(jìn)行了大量的探索和研究[1]。高偉[2]通過(guò)大量GNSS高程擬合理論與實(shí)踐研究，精化了似大地水準(zhǔn)面，其高程精度可達(dá)±2 cm。根據(jù)參考文獻(xiàn)[3]，深圳市似大地水準(zhǔn)面精化模型的高程異常的絕對(duì)精度達(dá)±1.4 cm。隨著數(shù)學(xué)與計(jì)算機(jī)技術(shù)的快速發(fā)展，地球重力場(chǎng)模型提出后，國(guó)內(nèi)外學(xué)者又提出了很多新方法、新模型來(lái)解決高程擬合問(wèn)題[4-5]。本文主要針對(duì)高原地區(qū)高程基準(zhǔn)難以解決的實(shí)際工程問(wèn)題，探索一種適用于高原地區(qū)的GNSS高程擬合方法，解決高原基礎(chǔ)測(cè)繪像控點(diǎn)高程基準(zhǔn)問(wèn)題。利用GNSS測(cè)量像控點(diǎn)平面坐標(biāo)得到精確大地高數(shù)據(jù)，結(jié)合測(cè)區(qū)少量已知公共點(diǎn)數(shù)據(jù)，通過(guò)高程擬合方法，解算出大量像控點(diǎn)的高精度高程異常值，進(jìn)而計(jì)算得到這些像控點(diǎn)的高程值。經(jīng)文獻(xiàn)研究發(fā)現(xiàn)，曲面擬合方法和“移去—恢復(fù)”法具有較高精度。因此，本文基于這兩種方法，利用MATLAB平臺(tái)自編程序，構(gòu)建多種擬合模型，并結(jié)合高原基礎(chǔ)測(cè)繪項(xiàng)目實(shí)例，得到適合高原基礎(chǔ)測(cè)繪的最佳高程擬合模型。本研究成果可以為高原基礎(chǔ)測(cè)繪以及其他相關(guān)測(cè)繪在高效解決高程基準(zhǔn)的問(wèn)題上提供參考。

2 GNSS高程擬合算法

2.1 多項(xiàng)式曲面擬合模型

當(dāng)測(cè)區(qū)范圍較大時(shí)，可采用多項(xiàng)式曲面擬合模型，通過(guò)已知點(diǎn)高程異常解算模型中的系數(shù)，然后計(jì)算待定點(diǎn)的高程異常，根據(jù)高程異常與正常高的關(guān)系計(jì)算待定點(diǎn)的高程值。該算法步驟如下[6]：

1)建立多項(xiàng)式曲面方程，如公式(1)所示。根據(jù)已知點(diǎn)的坐標(biāo)(x,y)和高程異常值ξ，求解多項(xiàng)式曲面函數(shù)的系數(shù)。

(1)

式中：ai為多項(xiàng)式曲面函數(shù)的系數(shù)，εi為高程異常殘差值。

2)根據(jù)步驟1)中求得的系數(shù)，將待求點(diǎn)坐標(biāo)代入多項(xiàng)式函數(shù)，解算出未知點(diǎn)的ξ，用該點(diǎn)的GNSS觀測(cè)的大地高，可求出該點(diǎn)的正常高的高程值。

多項(xiàng)式曲面擬合模型具有多種形式，根據(jù)公式(1)選用項(xiàng)的不同，可得不同的曲面擬合模型[7]：

1)當(dāng)選用式(1)中的前三項(xiàng)時(shí)，建立的擬合模型為平面擬合模型，該模型僅需3個(gè)已知點(diǎn)高程異常值，就可以解算方程系數(shù)。該模型適用于平緩地區(qū)，但在地形復(fù)雜地區(qū)精度不高。

2)當(dāng)選用式(1)中前三項(xiàng)和第五項(xiàng)進(jìn)行擬合時(shí)，可構(gòu)建相關(guān)平面擬合模型，解算模型系數(shù)需要至少4個(gè)已知點(diǎn)高程異常值。該模型仍是用平面代替附近高程異常值，適用于地勢(shì)平坦地區(qū)。

3)利用式(1)中前六項(xiàng)進(jìn)行擬合，所得到的數(shù)學(xué)模型是二次曲面擬合模型，模型系數(shù)求解需要6個(gè)已知點(diǎn)高程異常值，能夠較好地?cái)M合具有起伏變化的地形信息。

4)利用式(1)中所有項(xiàng)進(jìn)行擬合得到的數(shù)學(xué)模型為三次曲面擬合模型，模型較為復(fù)雜，需要至少10個(gè)已知點(diǎn)高程異常值，方可進(jìn)行求解系數(shù)。該模型對(duì)于較復(fù)雜的地形有很好的模擬效果。

已知高程點(diǎn)的個(gè)數(shù)，通常會(huì)多于模型所需要的數(shù)據(jù)個(gè)數(shù)，此時(shí)應(yīng)按最小二乘法平差方法，求解模型系數(shù)的最佳估值。

2.2 “移去—恢復(fù)”法高程擬合模型

李建成院士等諸多學(xué)者對(duì)我國(guó)似大地水準(zhǔn)面精化的理論和方法研究做出了巨大貢獻(xiàn)，研究發(fā)現(xiàn)重力場(chǎng)模型能夠提高GNSS的高程擬合精度，因此提出了“移去—恢復(fù)”法的高程異常擬合模型，在GNSS高程擬合中得到了廣泛的應(yīng)用[8]。

根據(jù)地球重力場(chǎng)的可疊加性，將高程異常分解為3個(gè)分量：

ξ=ξGM+ξΔG+ξT

(2)

式中：ξGM表示地球重力場(chǎng)中高程異常的長(zhǎng)波分量，ξΔG表示用斯托克斯積分得到的高程異常的中波分量，ξT表示地形起伏變化引起的高程異常的短波分量。3個(gè)分量的關(guān)系如圖1所示。

圖1 高程異常3個(gè)分量關(guān)系示意圖Fig.1 Schematic diagram of the relationship among three components of elevation anomaly

在一定范圍內(nèi)高程異常的中、長(zhǎng)波分量相對(duì)穩(wěn)定且變化??；而短波分量則受局部地形影響，小范圍內(nèi)也存在較大差異。因此，首先需要將高程異常中的中、長(zhǎng)波分量移去，然后對(duì)剩余的短波分量進(jìn)行擬合和插值計(jì)算，從待求點(diǎn)的短波分量中擬合高程異常值，最后再恢復(fù)待定點(diǎn)中、長(zhǎng)波分量。該方法可以有效地避免擬合和插值計(jì)算過(guò)程導(dǎo)致中、長(zhǎng)波分量的計(jì)算損失，進(jìn)而提高高程擬合的精度。模型計(jì)算基本步驟如下[9-10]：

1)移去：首先求出n個(gè)已知點(diǎn)的高程異常值，用重力場(chǎng)模型計(jì)算各點(diǎn)的中、長(zhǎng)波高程異常分量，然后計(jì)算高程異常的短波分量，即得到高程異常值與中、長(zhǎng)波分量的差值。

2)擬合：根據(jù)n個(gè)已知點(diǎn)的坐標(biāo)和短波分量值，利用MATLAB建立多種數(shù)學(xué)擬合模型，求解待求點(diǎn)的高程異常值短波分量的擬合和插值。

3)恢復(fù)：利用步驟2)中求解的高程異常短波分量，與步驟1)中求解的中、長(zhǎng)波分量進(jìn)行代數(shù)運(yùn)算，即可得到待求點(diǎn)的高程異常值的擬合值。

2.3 GNSS高程擬合精度分析

GNSS高程擬合精度通常需要根據(jù)內(nèi)符合精度、外符合精度和相關(guān)水準(zhǔn)等級(jí)精度要求進(jìn)行評(píng)定。

1)內(nèi)符合精度。為了檢驗(yàn)GNSS高程擬合的精度，利用已知高程異常值點(diǎn)與擬合高程異常值求殘差Vi，然后求其內(nèi)符合精度σ。

(3)

(4)

2)外符合精度。GNSS高程異常擬合精度還可用水準(zhǔn)實(shí)測(cè)方法來(lái)檢核，根據(jù)水準(zhǔn)測(cè)量結(jié)果求得高程異常值，再與擬合高程異常值計(jì)算殘差Vout，然后利用公式(4)計(jì)算得到外符合精度。

3)GNSS高程擬合精度評(píng)定。精度評(píng)定通常按照常規(guī)幾何水準(zhǔn)測(cè)量的精度評(píng)定方法，根據(jù)公共點(diǎn)、檢核點(diǎn)、待求點(diǎn)之間構(gòu)建的水準(zhǔn)網(wǎng)，設(shè)已知檢核點(diǎn)至公共點(diǎn)的距離L(km)，再根據(jù)國(guó)家水準(zhǔn)測(cè)量規(guī)范進(jìn)行公式計(jì)算和等級(jí)評(píng)定。

表1 水準(zhǔn)限差技術(shù)要求Tab.1 Levelingtolerancetechnicalrequirements等級(jí)允許殘差/mm三等水準(zhǔn)測(cè)量±12L四等水準(zhǔn)測(cè)量±20L五等水準(zhǔn)測(cè)量±30L

3 工程算例

3.1 工程概況

試驗(yàn)區(qū)為某高原基礎(chǔ)測(cè)繪項(xiàng)目，項(xiàng)目主要包括GNSS控制測(cè)量C級(jí)網(wǎng)及航空攝影測(cè)量外業(yè)工作。測(cè)區(qū)東西長(zhǎng)約165 km，南北寬約60 km，基礎(chǔ)測(cè)繪航空攝影測(cè)量設(shè)計(jì)為1∶10 000地形圖264幅(其中衛(wèi)星影像成圖157幅，航攝像片成圖107幅)。區(qū)域內(nèi)布設(shè)C級(jí)GNSS觀測(cè)點(diǎn)61個(gè)，其中檢核點(diǎn)25個(gè)(采用二、三等水準(zhǔn)測(cè)量)，像控點(diǎn)158個(gè)(待測(cè)高程點(diǎn))。C級(jí)GNSS控制測(cè)量，二、三等水準(zhǔn)測(cè)量，像控點(diǎn)平面坐標(biāo)測(cè)量均按國(guó)家相應(yīng)規(guī)范進(jìn)行測(cè)量。

3.2 多項(xiàng)式曲面擬合模型

根據(jù)多項(xiàng)式曲面擬合模型函數(shù)，自編MATLAB程序分別建立平面擬合、相關(guān)平面擬合、二次曲面擬合、三次曲面擬合4種擬合模型。利用已知點(diǎn)的數(shù)據(jù)得出各模型的擬合系數(shù)，再將像控點(diǎn)的大地坐標(biāo)代入該模型，解算出各像控點(diǎn)的高程異常。根據(jù)檢測(cè)點(diǎn)的高程異常值和擬合計(jì)算得到的高程異常值，計(jì)算殘差，再根據(jù)外符合精度進(jìn)行精度估算，并統(tǒng)計(jì)出各檢測(cè)點(diǎn)的殘差，結(jié)果如圖2所示。

圖2 多項(xiàng)式曲面法擬合殘差分布曲線Fig.2 Residual distribution curve of polynomial surface fitting

根據(jù)殘差分布曲線可得，平面擬合、相關(guān)平面擬合模型的殘差較大，而二次曲面擬合、三次曲面擬合模型殘差較小，三次曲面擬合方法總體趨勢(shì)優(yōu)于其他擬合模型，誤差分布在(-0.4 mm，0.8 mm)范圍。為了進(jìn)一步分析模型質(zhì)量，對(duì)檢核點(diǎn)高程異常的殘差進(jìn)行統(tǒng)計(jì)分析，結(jié)果如表2所示。

表2 高程異常的多項(xiàng)式曲面擬合殘差統(tǒng)計(jì)Tab.2 Statisticsofpolynomialsurfacefittingresidualsofelevationanomalies類別平面擬合模型/mm相關(guān)平面擬合模型/mm二次曲面擬合模型/mm三次曲面擬合模型/mm最大值0.7130.7720.3730.289最小值0.0160.0080.0330.006中誤差0.4810.5110.2280.181

根據(jù)表2分析，上述4種擬合方法中，平面擬合法與相關(guān)平面擬合法的中誤差較大，分別為0.481 m和0.511 m。而二次曲面擬合與三次曲面擬合方法較前兩種方法的精度有較大提升，分別為0.228 m和0.181 m。因此，三次曲面函數(shù)模型擬合效果最好。

3.3 “移去—恢復(fù)”擬合模型

根據(jù)EGM2008地球重力場(chǎng)模型，將所有公共點(diǎn)作為已知點(diǎn)，得到中、長(zhǎng)波高程異常分量。利用“移去—恢復(fù)”算法，解算像控點(diǎn)高程異常，并進(jìn)行精度分析。將公共點(diǎn)坐標(biāo)轉(zhuǎn)換為經(jīng)緯度，再利用AllTrans EGM2008 Calculator計(jì)算EGM2008地球重力場(chǎng)模型的高程異常。將已知點(diǎn)高程異常的中、長(zhǎng)波移除后，采用自編MATLAB程序構(gòu)建曲面擬合、多面函數(shù)擬合以及移動(dòng)曲面擬合3種模型，并分析模擬結(jié)果，最后利用同樣的原理計(jì)算外符合精度。4種模型的殘差分布結(jié)果如圖3所示。

圖3 “移去—恢復(fù)”法擬合殘差分布曲線Fig.3 Residual distribution curve fitted by “remove-restore”method

從圖3可以看出，曲面擬合方法與移動(dòng)曲面擬合方法稍差，多面函數(shù)法的總體外符合精度最好。但是總體來(lái)看，“移去—恢復(fù)”模型的精度要高于多項(xiàng)式曲面模型，誤差分布在(-0.3 mm，0.3 mm)范圍?！耙迫ァ謴?fù)”模型的殘差統(tǒng)計(jì)結(jié)果如表3所示。

表3 “移去—恢復(fù)”模型中誤差統(tǒng)計(jì)Tab.3 Theresidualstatisticsof“remove restore”model類 別曲面擬合/mm移動(dòng)曲面/mm多面函數(shù)/mm最大值0.2850.2920.141最小值0.0110.0060.000中誤差0.1660.1560.061

基于“移去—恢復(fù)”的三種擬合方法中，曲線擬合與移動(dòng)曲面方法的中誤差分別為0.166 mm和0.156 mm，而多面函數(shù)方法的精度最好，中誤差為0.061 mm。

通過(guò)本試驗(yàn)區(qū)的數(shù)據(jù)分析發(fā)現(xiàn)，高原基礎(chǔ)測(cè)繪高程基準(zhǔn)可采用高程擬合方法求解。通過(guò)對(duì)相對(duì)精度較高的集中模型進(jìn)行分析，發(fā)現(xiàn)“移去—恢復(fù)”擬合模型中多面函數(shù)方法精度最優(yōu)，且該方法的高程擬合精度能夠滿足高原基礎(chǔ)測(cè)繪像控點(diǎn)的高程精度要求。

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

利用MATLAB自編程序，結(jié)合高原基礎(chǔ)測(cè)繪項(xiàng)目，研究了GNSS高程擬合的多種擬合方法。本文利用少量公共點(diǎn)的GNSS高程值和高程異常值，通過(guò)高程擬合方法解算出像控點(diǎn)的高程值。結(jié)果表明，在高原地形復(fù)雜的地區(qū)，使用多面函數(shù)法的“移去—恢復(fù)”模型，可得到較高精度的GNSS擬合效果，精度最高可達(dá)到6 cm。高程擬合方法可為高原基礎(chǔ)測(cè)繪4D產(chǎn)品提供較高精度的高程基準(zhǔn)，并為高原測(cè)繪相關(guān)工作提供一種高程基準(zhǔn)測(cè)量的新方法。隨著數(shù)學(xué)和計(jì)算機(jī)技術(shù)的發(fā)展，高程擬合的算法會(huì)進(jìn)一步優(yōu)化，使其工作量更小，精度更高，更好地為高原測(cè)繪工作服務(wù)。