丁阿鹿，田 婕，聶建亮，劉曉云，郭鑫偉，趙大江，張海平

(1. 北京市勘察設(shè)計研究院有限公司，北京 100038； 2. 自然資源部大地測量數(shù)據(jù)處理中心，陜西 西安 710054； 3. 山東省國土測繪院，山東 濟南 250102)

GNSS與水準等大地測量觀測技術(shù)是國內(nèi)外學者研究地殼垂直運動的重要技術(shù)手段[1-6]。精密水準測量是地殼垂直運動監(jiān)測的經(jīng)典方法，是定量獲得垂直形變信息精確的技術(shù)方法，20世紀初該技術(shù)就應用于東京及周邊區(qū)域監(jiān)測，但其觀測成本高、觀測周期長、觀測條件要求高等因素決定了大區(qū)域精密水準測量的復測頻率不高。GNSS可以進行實時觀測，獲取區(qū)域垂直形變信息，GNSS連續(xù)運行基準站系統(tǒng)可以提供連續(xù)形變信息，提高垂直形變分析控制點密度，目前該技術(shù)已成為垂直形變監(jiān)測工作重要的觀測技術(shù)。然而上述兩種技術(shù)獲得的垂直運動之間存在一個垂線偏差，文獻[1]分析了小范圍內(nèi)監(jiān)測點的大地高變化速度近似替代正常高變化速度；兩類觀測數(shù)據(jù)所采用不同的觀測手段、噪聲模型、起算基準、數(shù)據(jù)處理方法，導致兩類垂直運動之間可能存在一定的系統(tǒng)偏差。為了充分發(fā)揮GNSS和水準對垂直運動的貢獻，需要根據(jù)GNSS、水準精度實現(xiàn)多源數(shù)據(jù)融合。文獻[6]提出了大地測量數(shù)據(jù)融合的基于觀測信息和基于平差結(jié)果的融合模式，證明兩種模式理論上是等價的。考慮到從原始觀測值層面融合難度較大，可以采用水準高差觀測值和GNSS垂直運動速度結(jié)果進行融合。這樣如何有效確定水準高差觀測值權(quán)矩陣與GNSS垂直運動速度權(quán)矩陣之間的比例關(guān)系成為實現(xiàn)數(shù)據(jù)融合的關(guān)鍵。目前，國內(nèi)外學者常用兩種方法確定權(quán)比：①方差分量估計方法，該方法能夠有效協(xié)調(diào)不同觀測之間的權(quán)重，已應用于地震反演[8]、地球重力場[9]、GPS網(wǎng)平差[10]、自適應濾波[11]、地面沉降[12]等方面；②利用反演結(jié)果分辨率與反演殘差的折中曲線確定權(quán)比值的方法，已應用于地震反演[13]、地殼垂直運動[14]。而方差分量估計方法是國內(nèi)外學者青睞的權(quán)比確定方法。在以上研究的基礎(chǔ)上，將水準原始高差觀測值和GNSS垂直運動速度結(jié)果作為數(shù)據(jù)融合觀測值，將一定數(shù)量GNSS水準點的垂直運動速度作為限制條件，利用方差分量估計方法自適應調(diào)整GNSS、水準觀測值權(quán)比，采用相關(guān)抗差估計方法抑制異常誤差影響，從而實現(xiàn)GNSS與水準的垂直運動融合，以提高區(qū)域垂直運動模型的精度。

1 GNSS水準聯(lián)合平差

水準觀測可以獲取高精度正常高方向形變信息，GNSS觀測能夠獲得大地高方向變化量，兩類形變信息之間存在較小差異，可用大地高變化代替正常高變化[1]。文獻[1，8]聯(lián)合GNSS連續(xù)站的高程分量變化速率與復測水準資料進行動態(tài)平差，解決了水準高程變化基準問題。

GNSS與水準聯(lián)合平差的誤差方程由水準高差觀測值誤差方程與GNSS高程方向速度值誤差方程組成。選取一定數(shù)量的GNSS水準點作為約束點，則聯(lián)合平差誤差方程為

(1)

聯(lián)合平差觀測值的權(quán)矩陣為

(2)

式中，n0為水準高差觀測值個數(shù)；m為GNSS大地高速度個數(shù)；α為兩類觀測值間的權(quán)比；水準高差觀測值權(quán)矩陣PL按照測站數(shù)定權(quán)；GNSS垂直運動速度值的權(quán)矩陣PG為垂直運動速度平差精度協(xié)因數(shù)矩陣的逆矩陣。為了控制觀測數(shù)據(jù)中異常誤差影響，采用相關(guān)觀測抗差估計最小二乘方法[15]降低異常觀測值的權(quán)值，構(gòu)造雙因子相關(guān)等價權(quán)矩陣[15]。

2 自適應聯(lián)合平差

理論上基于觀測信息、中間結(jié)果和最終平差結(jié)果的數(shù)據(jù)融合是一致的，為了得到可靠的計算結(jié)果，確定不同觀測數(shù)據(jù)之間合理的權(quán)比α顯得尤為重要。折中曲線確定權(quán)比方法依靠經(jīng)驗確定兩類觀測值之間的權(quán)比[12-13]，無法實現(xiàn)自適應糾正，而方差分量估計方法依據(jù)驗后估計精度逐步調(diào)整不同觀測值之間的權(quán)比，使權(quán)比更接近真實比例關(guān)系。

Helmert方差分量估計方法得到的第i0類觀測值的殘差平方和為

(3)

Helmert方差分量估計的Baumker簡化公式為

(4)

式中，ni0為第i0類觀測值的個數(shù)。不同類觀測值之間的權(quán)比α為

(5)

經(jīng)過多次迭代，直到α趨于穩(wěn)定。

3 算例分析

算例采用2006、2015年山東省似大地水準面精化GNSS、水準觀測數(shù)據(jù)。水準網(wǎng)觀測按照二等水準網(wǎng)觀測綱要于2004—2006、2015年進行觀測，GNSS控制網(wǎng)按照C級網(wǎng)觀測綱要于2004—2006、2015年同步進行數(shù)據(jù)采集(水準路線及C級點示意圖如圖1所示)，其中兩期分別有1337、1105個水準高差觀測值，兩期GNSS觀測共有153個C級重合點。利用兩期GNSS控制網(wǎng)、水準網(wǎng)觀測數(shù)據(jù)，采用聯(lián)合平差方法實現(xiàn)GNSS水準融合，從而獲得山東省區(qū)域垂直形變信息。

為了實現(xiàn)水準觀測值與GNSS垂直運動速度聯(lián)合平差，采用如下步驟處理：①使用GAMIT/GLOBK軟件處理GNSS數(shù)據(jù)，獲得GNSS垂直運動速度及協(xié)方差矩陣；②整理兩期水準網(wǎng)重合點數(shù)據(jù)，水準網(wǎng)平差起算點采用青島國家水準原點，其速度為+2.2 mm/a[4]；③組建聯(lián)合平差觀測誤差方程，水準觀測值根據(jù)測站數(shù)定權(quán)，GNSS根據(jù)協(xié)方差矩陣定權(quán)，兩類觀測值之間的權(quán)比α采用Helmert方差分量估計方法確定，從而獲取垂直運動速度。

為了比較GNSS水準融合效果，采用4種方案進行比較分析：

(1) 利用最小二乘方法進行水準網(wǎng)整體平差。

(2) 利用抗差最小二乘進行水準網(wǎng)整體平差。

(3) 選取兩類高程變化一致性較好的60個GNSS水準并置點作為約束點，利用聯(lián)合平差方法進行融合，權(quán)比利用Helmert方差分量估計方法確定。

(4) 同方案(3)，構(gòu)造雙因子相關(guān)等價權(quán)函數(shù)削弱異常觀測值的權(quán)值。

4種方案計算結(jié)果見表1，誤差曲線如圖2所示。由計算結(jié)果可以得到如下初步結(jié)論:

表1 各方案誤差統(tǒng)計信息

(1) GNSS水準聯(lián)合平差與水準網(wǎng)平差相比，聯(lián)合平差精度有一定程度提高。聯(lián)合平差充分利用了水準高差觀測值和GNSS平差結(jié)果，通過約束GNSS水準點垂直運動速度建立水準垂直運動速度與GNSS垂直運動速度之間的關(guān)聯(lián)，實現(xiàn)水準高差觀測值與GNSS垂直運動速度的有效融合，從而提高水準環(huán)包圍區(qū)域的垂直形變信息精度。

(2) 方案2、4與方案1、3相比，對應平差精度都有一定程度提高。方案2、4采用抗差估計方法降低了異常觀測值的權(quán)值，削弱了異常誤差影響，提高了垂直運動估計精度；方案3、4與方案1、2相比，利用Helmert方差分量估計方法根據(jù)驗后殘差自動調(diào)整兩類觀測值之間的權(quán)比，實現(xiàn)了GNSS與水準的有效數(shù)據(jù)融合，提高了垂直形變信息的可靠性；由于誤差轉(zhuǎn)移、處理策略的不同，聯(lián)合平差與水準網(wǎng)動態(tài)平差殘差最大值、最小值的位置可能會存在一定不同。

(3) 利用GNSS水準數(shù)據(jù)獲取了山東省西降東升、北降南升的整體垂直運動趨勢。西北、西南部地區(qū)沉降嚴重，整體呈下降趨勢，垂直運動速度平均為20.0 mm/a；德州、菏澤、廣饒地區(qū)出現(xiàn)3個沉降中心，主要由于地下水過度開采導致地面沉降，其中廣饒地區(qū)沉降最嚴重，10年平均沉降速度70.0 mm/a左右；而在膠東半島呈5.0 mm/a的上升趨勢。

致謝：感謝山東國土測繪院與陸態(tài)網(wǎng)絡中心提供的GNSS試驗觀測數(shù)據(jù)。