丁阿鹿,田 婕,聶建亮,劉曉云,郭鑫偉,趙大江,張海平
(1. 北京市勘察設計研究院有限公司,北京 100038; 2. 自然資源部大地測量數據處理中心,陜西 西安 710054; 3. 山東省國土測繪院,山東 濟南 250102)
GNSS與水準等大地測量觀測技術是國內外學者研究地殼垂直運動的重要技術手段[1-6]。精密水準測量是地殼垂直運動監測的經典方法,是定量獲得垂直形變信息精確的技術方法,20世紀初該技術就應用于東京及周邊區域監測,但其觀測成本高、觀測周期長、觀測條件要求高等因素決定了大區域精密水準測量的復測頻率不高。GNSS可以進行實時觀測,獲取區域垂直形變信息,GNSS連續運行基準站系統可以提供連續形變信息,提高垂直形變分析控制點密度,目前該技術已成為垂直形變監測工作重要的觀測技術。然而上述兩種技術獲得的垂直運動之間存在一個垂線偏差,文獻[1]分析了小范圍內監測點的大地高變化速度近似替代正常高變化速度;兩類觀測數據所采用不同的觀測手段、噪聲模型、起算基準、數據處理方法,導致兩類垂直運動之間可能存在一定的系統偏差。為了充分發揮GNSS和水準對垂直運動的貢獻,需要根據GNSS、水準精度實現多源數據融合。文獻[6]提出了大地測量數據融合的基于觀測信息和基于平差結果的融合模式,證明兩種模式理論上是等價的。考慮到從原始觀測值層面融合難度較大,可以采用水準高差觀測值和GNSS垂直運動速度結果進行融合。這樣如何有效確定水準高差觀測值權矩陣與GNSS垂直運動速度權矩陣之間的比例關系成為實現數據融合的關鍵。目前,國內外學者常用兩種方法確定權比:①方差分量估計方法,該方法能夠有效協調不同觀測之間的權重,已應用于地震反演[8]、地球重力場[9]、GPS網平差[10]、自適應濾波[11]、地面沉降[12]等方面;②利用反演結果分辨率與反演殘差的折中曲線確定權比值的方法,已應用于地震反演[13]、地殼垂直運動[14]。而方差分量估計方法是國內外學者青睞的權比確定方法。在以上研究的基礎上,將水準原始高差觀測值和GNSS垂直運動速度結果作為數據融合觀測值,將一定數量GNSS水準點的垂直運動速度作為限制條件,利用方差分量估計方法自適應調整GNSS、水準觀測值權比,采用相關抗差估計方法抑制異常誤差影響,從而實現GNSS與水準的垂直運動融合,以提高區域垂直運動模型的精度。
水準觀測可以獲取高精度正常高方向形變信息,GNSS觀測能夠獲得大地高方向變化量,兩類形變信息之間存在較小差異,可用大地高變化代替正常高變化[1]。文獻[1,8]聯合GNSS連續站的高程分量變化速率與復測水準資料進行動態平差,解決了水準高程變化基準問題。
GNSS與水準聯合平差的誤差方程由水準高差觀測值誤差方程與GNSS高程方向速度值誤差方程組成。選取一定數量的GNSS水準點作為約束點,則聯合平差誤差方程為
(1)
聯合平差觀測值的權矩陣為
(2)
式中,n0為水準高差觀測值個數;m為GNSS大地高速度個數;α為兩類觀測值間的權比;水準高差觀測值權矩陣PL按照測站數定權;GNSS垂直運動速度值的權矩陣PG為垂直運動速度平差精度協因數矩陣的逆矩陣。為了控制觀測數據中異常誤差影響,采用相關觀測抗差估計最小二乘方法[15]降低異常觀測值的權值,構造雙因子相關等價權矩陣[15]。
理論上基于觀測信息、中間結果和最終平差結果的數據融合是一致的,為了得到可靠的計算結果,確定不同觀測數據之間合理的權比α顯得尤為重要。折中曲線確定權比方法依靠經驗確定兩類觀測值之間的權比[12-13],無法實現自適應糾正,而方差分量估計方法依據驗后估計精度逐步調整不同觀測值之間的權比,使權比更接近真實比例關系。
Helmert方差分量估計方法得到的第i0類觀測值的殘差平方和為
(3)

Helmert方差分量估計的Baumker簡化公式為
(4)
式中,ni0為第i0類觀測值的個數。不同類觀測值之間的權比α為
(5)
經過多次迭代,直到α趨于穩定。
算例采用2006、2015年山東省似大地水準面精化GNSS、水準觀測數據。水準網觀測按照二等水準網觀測綱要于2004—2006、2015年進行觀測,GNSS控制網按照C級網觀測綱要于2004—2006、2015年同步進行數據采集(水準路線及C級點示意圖如圖1所示),其中兩期分別有1337、1105個水準高差觀測值,兩期GNSS觀測共有153個C級重合點。利用兩期GNSS控制網、水準網觀測數據,采用聯合平差方法實現GNSS水準融合,從而獲得山東省區域垂直形變信息。
為了實現水準觀測值與GNSS垂直運動速度聯合平差,采用如下步驟處理:①使用GAMIT/GLOBK軟件處理GNSS數據,獲得GNSS垂直運動速度及協方差矩陣;②整理兩期水準網重合點數據,水準網平差起算點采用青島國家水準原點,其速度為+2.2 mm/a[4];③組建聯合平差觀測誤差方程,水準觀測值根據測站數定權,GNSS根據協方差矩陣定權,兩類觀測值之間的權比α采用Helmert方差分量估計方法確定,從而獲取垂直運動速度。
為了比較GNSS水準融合效果,采用4種方案進行比較分析:
(1) 利用最小二乘方法進行水準網整體平差。
(2) 利用抗差最小二乘進行水準網整體平差。
(3) 選取兩類高程變化一致性較好的60個GNSS水準并置點作為約束點,利用聯合平差方法進行融合,權比利用Helmert方差分量估計方法確定。
(4) 同方案(3),構造雙因子相關等價權函數削弱異常觀測值的權值。
4種方案計算結果見表1,誤差曲線如圖2所示。由計算結果可以得到如下初步結論:

表1 各方案誤差統計信息
(1) GNSS水準聯合平差與水準網平差相比,聯合平差精度有一定程度提高。聯合平差充分利用了水準高差觀測值和GNSS平差結果,通過約束GNSS水準點垂直運動速度建立水準垂直運動速度與GNSS垂直運動速度之間的關聯,實現水準高差觀測值與GNSS垂直運動速度的有效融合,從而提高水準環包圍區域的垂直形變信息精度。
(2) 方案2、4與方案1、3相比,對應平差精度都有一定程度提高。方案2、4采用抗差估計方法降低了異常觀測值的權值,削弱了異常誤差影響,提高了垂直運動估計精度;方案3、4與方案1、2相比,利用Helmert方差分量估計方法根據驗后殘差自動調整兩類觀測值之間的權比,實現了GNSS與水準的有效數據融合,提高了垂直形變信息的可靠性;由于誤差轉移、處理策略的不同,聯合平差與水準網動態平差殘差最大值、最小值的位置可能會存在一定不同。
(3) 利用GNSS水準數據獲取了山東省西降東升、北降南升的整體垂直運動趨勢。西北、西南部地區沉降嚴重,整體呈下降趨勢,垂直運動速度平均為20.0 mm/a;德州、菏澤、廣饒地區出現3個沉降中心,主要由于地下水過度開采導致地面沉降,其中廣饒地區沉降最嚴重,10年平均沉降速度70.0 mm/a左右;而在膠東半島呈5.0 mm/a的上升趨勢。
致謝:感謝山東國土測繪院與陸態網絡中心提供的GNSS試驗觀測數據。