基于絕對重力基準控制的流動重力觀測資料動態平差方法研究

康開軒 李 輝 申重陽 孫少安 郝洪濤，3

1 中國地震局地震研究所（地震大地測量重點實驗室），武漢市洪山側路40號，430071

2 中國地震局地殼應力研究所科技創新基地，武漢市洪山側路40號，430071

3 中國科學院大學，北京市玉泉路甲19號，100039

我國流動重力觀測［1］自20世紀60年代起積累了大量寶貴資料，如何統一時間基準，整合不同空間尺度的重力網多期復測資料，得到重力空間變化趨勢背景場，是目前亟需解決的問題［2］。重力監測網是一種動態監測網，重力測量資料是一種時序觀測值［3］。目前，大多數重力測網解算采取的平差模型均為分期經典平差算法，該算法是基于靜態重力場背景下的，即假設在測網觀測期間各個觀測點的重力值是常數［2］。實際上，由于地球系統復雜的動力學過程，導致地球表面重力場隨時間發生變化［4］。在相對重力聯測數據處理中，時變重力場潮汐信號（如固體潮、海潮、極潮）以及氣壓、溫度等環境因素引起的區域重力非潮汐效應，可采用相對完善的數學物理模型加以改正［5－7］，但時變重力場中因地殼形變、構造運動等地球動力學因素引起的長期變化趨勢往往被忽略。動態平差法對靜態分期平差有良好的概括性，可以同時處理包括不完整觀測的多期觀測數據。建立基于絕對重力基準控制的多期流動重力觀測的動態平差方法與模型可有效解決多網多期資料整體解算問題［2］。

本文擬將動態平差原理應用于滇西重力測網1986～2009年共52 期流動重力觀測數據的處理，采用線性速率模型擬合重力背景場長期趨勢；優化方程解算，對變化率參數作顯著性檢驗，剔除變化不顯著的參數，使速率模型更接近真實物理場；采用Helmert方差分量估計方法實現基于驗后信息調整不同期次測量中不同儀器觀測值的權重，保證精密重力網測量多期數據聯合解算結果的合理性。

1 原理與模型

相對重力聯測網平差的基本元素為相鄰測站i、j的重力段差觀測值Δgij，滿足基本觀測方程［5］：

式中，gi、gj分別為測站i、j的未知重力值。考慮彈簧重力儀格值系統誤差及儀器漂移的影響，引入儀器格值函數ΔF（Lk，Pk）和儀器線性漂移率Dk，則：

式中，ΔF（Lk，Pk）為長波項系數Lk和周期項系數Pk的函數，這里采用李輝等［5］提出的數學模型。

考慮到重力場在觀測過程中不是恒定不變的，即式（2）中未知數gi、gj均為時間的函數，本文擬引入動態參數——重力變化率來描述時變重力場長期變化趨勢。為此，選定基準時間t0，時間t處的重力場模型為g（t）＝g（t0）＋（t－t0），其中g（t0）為基準時間t0處的重力場，為重力變化率。由此得到測段Δgij在觀測時刻t處的觀測方程：

式中，基準時間t0處的重力值為、重力變化率為，儀器線性漂移率為Dk，格值系數Lk、Pk為未知變量（若儀器格值系數由基線標定精密確定，則Lk、Pk為已知參數）。上述函數模型的系數陣列秩虧，為得到平差參數的最小二乘解，測網必須具有足夠的起算數據，即要采用合適的平差基準對測網進行約束［5］。約束條件采用絕對重力基準：

2 模型解算細節

2.1 數據描述

選取滇西重力網1986～2009年共52期相對重力聯測數據進行整體動態平差解算，并采用下關重力基準點1986～2009年、麗江與昆明重力基準點1986～2001年的絕對重力復測資料作基準控制。該網以下關為中心布設，北到麗江、攀枝花一線，南抵云縣一帶，主要監控紅河斷裂帶北段區域構造活動的重力效應［6］。自1985年以來，滇西網重力測量采用LCR－G 型相對重力儀每年進行2～3期流動重力測量，前后共使用過的儀器達12臺，測量精度約為（7～10）×10－8ms－2。

2.2 觀測值協方差矩陣的構建

2.2.1 重力段差觀測值協方差陣CΔg的構建

多期相對重力聯測的段差觀測值協方差陣由重力儀性能、測網布設情況、儀器運輸過程及野外觀測環境、觀測者個人技術等多種因素共同決定。各期觀測數據之間相互獨立，協方差陣可表示為：

其中，σi為第i期觀測數據的先驗中誤差，Qi為第i期觀測數據的協方差陣。

段差觀測值的協方差陣Qi的確定應考慮單程段差觀測值的精度以及相鄰段差觀測值之間的相關系數。Pagiatakis等［8］的研究結果表明，彈簧型重力儀觀測數據不可避免地受到儀器零漂的影響，區域重力觀測中重力儀在相鄰兩重力測點間漂移量與測段觀測時間相關，單程段差觀測數據的精度與段差觀測時間Δtij成反比，段差觀測值的權可依據其觀測時間的倒數確定。滇西測網測點平均間距約30km，段差觀測時間約為30 min～2h，均值約1h。為合理估算儀器零漂對各段差觀測值精度的影響，采取類似Pagiatakis［8］的方案：段差觀測時間Δtij在30min以內的段差觀測值比例系數factordrift為1；大于2h，factordrift為2；Δtij在上述區間的段差采用線性插值計算比例系數。

另外，相對重力聯測的段差觀測值是一種非獨立觀測量。流動重力野外觀測一般采用測線往返觀測的方案，相鄰段差觀測值均由公共重力測點相連，因此，相對重力聯測結果必須采用相關平差才是嚴密的。劉少府等［9］基于相關系數法確定相鄰段差觀測值理論相關系數為－0.5，實際觀測數據的相關系數在0.5左右，只有極個別的偏大或偏小。設相鄰段差觀測值相關系數為－0.5，第i期段差觀測數據的協方差陣為［9］：

2.2.2 絕對重力基準觀測值協方差陣CG的構建

不同類型的絕對重力儀的觀測精度不同。滇西重力網中昆明、麗江、下關絕對重力觀測在20世紀70年代末、80年代初主要采用意大利IMGC型和我國NIM 型絕對重力儀，測量精度約為（10～15）×10－8ms－2；90年代以來采用德國和芬蘭JILAG 型重力儀，測量精度約5×10－8ms－2；1995年之后引進由美國研制的FG5 絕對重力儀，測量精度可達（1～2）×10－8ms－2［10］。多期絕對重力觀測資料［11］按不同類型重力儀觀測精度組成權函數陣，采用加權最小二乘回歸擬合得到各基準點重力值、重力變化率及精度（圖1）。絕對重力基準觀測值協方差陣CG由上述最小二乘擬合參數誤差估值確定。

圖1 絕對重力基準點重力時變曲線Fig.1 The time series of absolute gravity measurements

2.3 模型初值選取及迭代計算

測網各相對聯測點重力值初值g（0）和重力變化率初值由多期重力相對聯測結果加權最小二乘回歸擬合獲得，權重采用各期數據實測精度，儀器格值系數和線性漂移率的初值均由儀器短基線標定結果獲得。

為客觀分配觀測值協方差矩陣中段差觀測值和絕對重力觀測值的權重比例，采用Helmert方差分量估計原理［12］，根據驗后信息使兩種獨立觀測量的權重趨于合理，迭代計算終止條件為（k1、k2為式（5）中的比例因子）。

3 平差結果及精度

對清理后的滇西重力測網1986～2009年共52期流動重力觀測資料進行整體動態平差解算，選定數據長度區間的中間年份（1998年）為基準時間，昆明、麗江、下關絕對重力觀測為空間基準控制，得到平差后各測點相對于基準時間的重力值和重力變化率。經統計，測區內各測點的重力年變率范圍為（－0.039 5±0.024 1）×10－5ms－2／a，精度范圍為（0.001±0.029 8）×10－5ms－2／a。由于流動重力多年復測的實測網形不可能完全一致，各測點的復測次數是不同的。由圖2的測點年變率解算精度與復測次數的統計結果可以看出，測點重力變化率精度與測點的復測次數相關，復測次數小于15的測點解算精度相對較差，最大可達±0.029×10－5ms－2／a；復測次數15～35的測點，精度約為±0.002×10－5ms－2／a；復測次數大于35 的測點，精度約為±0.001×10－5ms－2／a。這與平差理論是相符的，觀測次數多的測點其自由度也越大，解算結果精度越高。

圖2 重力年變率精度統計結果Fig.2 The statistic result of the gravity annual rate precision

圖3 測點重力年變率空間分布（全部測點）Fig.3 The spatial distribution of the gravity annual rate in the network（includes all the gravity sites）

圖4 測點重力年變率空間分布（剔除精度低于0.01×10－5 ms－2／a的測點）Fig.4 The spatial distribution of the gravity annual rate in the network（removes the gravity sites where the precision is less than 0.01×10－5 ms－2／a）

圖3為滇西測網中所有重力測點的重力年變率空間分布圖，圖4 為剔除了精度低于0.01×10－5ms－2／a的測點后的重力變化空間分布圖。通過比較分析，剔除精度較差的點值后，等值線部分畸點（尤其是麗江、鳳慶、彌渡等重力測點附近的畸點）消失，空間分布圖趨于平滑。由圖4可知，滇西地區南端的鳳慶一帶，中部的彌渡、姚安附近及雄楚地區，北部的麗江附近地區及攀枝花一帶，重力場變化呈減小趨勢，最大幅值可達－0.010×10－5ms－2／a；中部的賓川一帶，及南部的景東地區（馬鹿田附近）重力場變化呈小幅增加趨勢，變化幅值集中在（0.001～0.006）×10－5ms－2／a。

與分期靜態平差模型比較，基于動態平差模型的解算結果更為合理。這是因為：1）分期靜態平差解算中沒有考慮觀測過程中重力點值隨時間的變化，該時間因素引起的系統誤差增加了其解算結果的不確定性，而動態平差模型中考慮了這一時間因素的影響；2）基于分期靜態平差解算獲得的重力年變率為各期次觀測結果的等權擬合結果，沒有考慮不同儀器不同年代觀測資料的精度差異，而基于動態平差模型的整體平差解算對參與計算的不同儀器不同年代觀測資料采用了不同的權重，并依據驗后信息對各類觀測數據作了較為客觀的調整，同時采用相關穩健估計有效減弱了粗差的影響，其解算結果較分期平差結果更為合理。

4 結 語

研究結果表明，隨著測區多年復測資料的積累，區域重力背景場長期變化趨勢特征逐漸收斂，測點線性速率反映了測區重力多年累積變化趨勢。在經典平差模型中引入初始時間歷元和測點的重力變化率，采用單點速率模型，聯合絕對觀測資料，建立段差觀測值的誤差方程，并依據最小二乘原理解算法方程，得到測網中每個測點的重力值及重力變化平均速率的平差值，實現多期觀測資料共同解算，同時也實現了單期網中各個測點時間基準的統一歸算。

在處理較大面積的測網多期復測資料時，采用基于動態平差原理的整體解算方法更為合理。本文采用線性速率模型可有效擬合區域重力場的長期線性變化趨勢，為了能夠更準確地模擬重力場的非線性時變信息，需要引入更精確的時變模型。以后的研究中將考慮基于高階項多項式擬合和周期擬合的動態平差模型。

致謝：感謝中國地震局地震研究所和云南省地震局卓有成效的觀測工作和數據資料清理工作！

［1］李輝.中國大陸近期重力場動態變化圖像［J］.大地測量與地 球 動 力 學，2009（3）：1－10（Li Hui.Dynamic Gravity Change in Recent Years in China Continent［J］.Journal of Geodesy and Geodynamic，2009（3）：1－10）

［2］孫少安，康開軒，黃邦武.關于區域重力場變化基準的思考［J］.大地測量與地球動力學，2012（1）：17－20（Sun Shaoan，Kang Kaixuan，Huang Bangwu.Thinking on Datum of Regional Gravity Field Variation［J］.Journal of Geodesy and Geodynamic，2012（1）：17－20）

［3］張祖勝，楊元喜，孫漢榮，等.地殼形變監測中的水準與重力資料聯合解算［J］.地震學報，1998，20（1），76－85（Zhang Zusheng，Yang Yuanxi，Sun Hanrong，et al.Combination of the Leveling and Gravity Measurements in Crustal Deformation Detection［J］.Acta Seismologica Sinica，1998，20（1），76－85）

［4］孫和平，徐建橋，黎瓊.地球重力場的精細頻譜結構及其應用［J］.地球物理學進展，2006，21（2），345－352（Sun Heping，Xu Jianqiao，Li Qiong.Detail Spectral Structure of Earth’s Gravity Field and Its Application［J］.Progress in Geophysics，2006，21（2），345－352）

［5］李輝，劉冬至，劉紹府.地震重力監測網統一平差模型的建立［J］.地殼形變與地震，1991（增刊）：68－74（Li Hui，Liu Dongzhi，Liu Shaofu.Integtated Adjustment Models for the Seismic－Gravity Network［J］.Crustal Deformation and Earthquake，1991（Supp）：68－74）

［6］孫少安，項愛民，李輝.滇西和北京區域重力場演化及其與地震的關系的探討［J］.地震，1999，19（1）：97－106（Sun Shaoan，Xiang Aimin，Li Hui.Research on Evolution of Western Yunnan and Beijing Regional Gravity Fields and Their Relationship with Earthquake［J］.Earthquake，1999，19（1）：97－106）

［7］孫少安，項愛民，吳維日.LCR－G 型重力儀儀器參數的時變特征［J］.大地測量與地球 動力 學，2002（2）：101－105（Sun Shaoan，Xiang Aimin，Wu Weiri.Time－Varying Characteristics of LCR－G Gravimeter Parameters［J］.Journal of Geodesy and Geodynamic，2002（2）：101－105）

［8］Pagiatakis S D，Salib P.Historical Relative Gravity Observations and the Time Rate of Change of Gravity due to Postglacial Rebound and Other Tectonic Movements in Canada［J］.Journal of Geophysical Research：Solid Earth（1978－2012），2003，108（B9）

［9］劉紹府，李輝，劉冬至.拉科斯特重力儀測量平差中的相關問題［J］.地殼形變與地震，1990（2）：67－73（Liu Shaofu，Li Hui，Liu Dongzhi.Correlative Problem of Measurement Adjustment of LCR Gravimeter［J］.Crustal Deformation and Earthquake，1990（2）：67－73）

［10］李輝，付廣裕，孫少安，等.滇西地區重力場動態變化計算［J］.地殼形變與地震，2000（1）：60－66（Li Hui，Fu Guangyu，Sun Shaoan，et al.Computation on Dynamic Gravity Changes in the Western Area of Yunnan Province［J］.Crustal Deformation and Earthquake，2000（1）：60－66）

［11］王勇，張為民，詹金剛，等.重復絕對重力測量觀測的滇西地區和拉薩點的重力變化及其意義［J］.地球物理學報，2004，47（1）：95－100（Wang Yong，Zhang Weimin，Zhan Jingang，et al.Gravity Changes Detected by Repeated Absolute Gravity Measurements in the Western Yunnan and Lhasa，China and Its Implication［J］.Chinese Journal of Geophysics，2004，47（1）：95－100）

［12］康開軒，邢燦飛，李輝，等.抗差Helmert方差分量估計在重力網平差中的應用［J］.大地測量與地球動力學，2008（5）：115－119（Kang Kaixuan，Xing Canfei，Li Hui，et al.Application of Robust Helmert Method of Variance Component Estimation in Adjustment of Gravity Network［J］.Journal of Geodesy and Geodynamics，2008（5）：115－119）