基于絕對(duì)重力基準(zhǔn)控制的流動(dòng)重力觀測(cè)資料動(dòng)態(tài)平差方法研究

康開軒 李 輝 申重陽(yáng) 孫少安 郝洪濤，3

1 中國(guó)地震局地震研究所（地震大地測(cè)量重點(diǎn)實(shí)驗(yàn)室），武漢市洪山側(cè)路40號(hào)，430071

2 中國(guó)地震局地殼應(yīng)力研究所科技創(chuàng)新基地，武漢市洪山側(cè)路40號(hào)，430071

3 中國(guó)科學(xué)院大學(xué)，北京市玉泉路甲19號(hào)，100039

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

本文擬將動(dòng)態(tài)平差原理應(yīng)用于滇西重力測(cè)網(wǎng)1986～2009年共52 期流動(dòng)重力觀測(cè)數(shù)據(jù)的處理，采用線性速率模型擬合重力背景場(chǎng)長(zhǎng)期趨勢(shì)；優(yōu)化方程解算，對(duì)變化率參數(shù)作顯著性檢驗(yàn)，剔除變化不顯著的參數(shù)，使速率模型更接近真實(shí)物理場(chǎng)；采用Helmert方差分量估計(jì)方法實(shí)現(xiàn)基于驗(yàn)后信息調(diào)整不同期次測(cè)量中不同儀器觀測(cè)值的權(quán)重，保證精密重力網(wǎng)測(cè)量多期數(shù)據(jù)聯(lián)合解算結(jié)果的合理性。

1 原理與模型

相對(duì)重力聯(lián)測(cè)網(wǎng)平差的基本元素為相鄰測(cè)站i、j的重力段差觀測(cè)值Δgij，滿足基本觀測(cè)方程［5］：

式中，gi、gj分別為測(cè)站i、j的未知重力值?？紤]彈簧重力儀格值系統(tǒng)誤差及儀器漂移的影響，引入儀器格值函數(shù)ΔF（Lk，Pk）和儀器線性漂移率Dk，則：

式中，ΔF（Lk，Pk）為長(zhǎng)波項(xiàng)系數(shù)Lk和周期項(xiàng)系數(shù)Pk的函數(shù)，這里采用李輝等［5］提出的數(shù)學(xué)模型。

考慮到重力場(chǎng)在觀測(cè)過(guò)程中不是恒定不變的，即式（2）中未知數(shù)gi、gj均為時(shí)間的函數(shù)，本文擬引入動(dòng)態(tài)參數(shù)——重力變化率來(lái)描述時(shí)變重力場(chǎng)長(zhǎng)期變化趨勢(shì)。為此，選定基準(zhǔn)時(shí)間t0，時(shí)間t處的重力場(chǎng)模型為g（t）＝g（t0）＋（t－t0），其中g(shù)（t0）為基準(zhǔn)時(shí)間t0處的重力場(chǎng)，為重力變化率。由此得到測(cè)段Δgij在觀測(cè)時(shí)刻t處的觀測(cè)方程：

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

2 模型解算細(xì)節(jié)

2.1 數(shù)據(jù)描述

選取滇西重力網(wǎng)1986～2009年共52期相對(duì)重力聯(lián)測(cè)數(shù)據(jù)進(jìn)行整體動(dòng)態(tài)平差解算，并采用下關(guān)重力基準(zhǔn)點(diǎn)1986～2009年、麗江與昆明重力基準(zhǔn)點(diǎn)1986～2001年的絕對(duì)重力復(fù)測(cè)資料作基準(zhǔn)控制。該網(wǎng)以下關(guān)為中心布設(shè)，北到麗江、攀枝花一線，南抵云縣一帶，主要監(jiān)控紅河斷裂帶北段區(qū)域構(gòu)造活動(dòng)的重力效應(yīng)［6］。自1985年以來(lái)，滇西網(wǎng)重力測(cè)量采用LCR－G 型相對(duì)重力儀每年進(jìn)行2～3期流動(dòng)重力測(cè)量，前后共使用過(guò)的儀器達(dá)12臺(tái)，測(cè)量精度約為（7～10）×10－8ms－2。

2.2 觀測(cè)值協(xié)方差矩陣的構(gòu)建

2.2.1 重力段差觀測(cè)值協(xié)方差陣CΔg的構(gòu)建

多期相對(duì)重力聯(lián)測(cè)的段差觀測(cè)值協(xié)方差陣由重力儀性能、測(cè)網(wǎng)布設(shè)情況、儀器運(yùn)輸過(guò)程及野外觀測(cè)環(huán)境、觀測(cè)者個(gè)人技術(shù)等多種因素共同決定。各期觀測(cè)數(shù)據(jù)之間相互獨(dú)立，協(xié)方差陣可表示為：

其中，σi為第i期觀測(cè)數(shù)據(jù)的先驗(yàn)中誤差，Qi為第i期觀測(cè)數(shù)據(jù)的協(xié)方差陣。

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

另外，相對(duì)重力聯(lián)測(cè)的段差觀測(cè)值是一種非獨(dú)立觀測(cè)量。流動(dòng)重力野外觀測(cè)一般采用測(cè)線往返觀測(cè)的方案，相鄰段差觀測(cè)值均由公共重力測(cè)點(diǎn)相連，因此，相對(duì)重力聯(lián)測(cè)結(jié)果必須采用相關(guān)平差才是嚴(yán)密的。劉少府等［9］基于相關(guān)系數(shù)法確定相鄰段差觀測(cè)值理論相關(guān)系數(shù)為－0.5，實(shí)際觀測(cè)數(shù)據(jù)的相關(guān)系數(shù)在0.5左右，只有極個(gè)別的偏大或偏小。設(shè)相鄰段差觀測(cè)值相關(guān)系數(shù)為－0.5，第i期段差觀測(cè)數(shù)據(jù)的協(xié)方差陣為［9］：

2.2.2 絕對(duì)重力基準(zhǔn)觀測(cè)值協(xié)方差陣CG的構(gòu)建

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

圖1 絕對(duì)重力基準(zhǔn)點(diǎn)重力時(shí)變曲線Fig.1 The time series of absolute gravity measurements

2.3 模型初值選取及迭代計(jì)算

測(cè)網(wǎng)各相對(duì)聯(lián)測(cè)點(diǎn)重力值初值g（0）和重力變化率初值由多期重力相對(duì)聯(lián)測(cè)結(jié)果加權(quán)最小二乘回歸擬合獲得，權(quán)重采用各期數(shù)據(jù)實(shí)測(cè)精度，儀器格值系數(shù)和線性漂移率的初值均由儀器短基線標(biāo)定結(jié)果獲得。

為客觀分配觀測(cè)值協(xié)方差矩陣中段差觀測(cè)值和絕對(duì)重力觀測(cè)值的權(quán)重比例，采用Helmert方差分量估計(jì)原理［12］，根據(jù)驗(yàn)后信息使兩種獨(dú)立觀測(cè)量的權(quán)重趨于合理，迭代計(jì)算終止條件為（k1、k2為式（5）中的比例因子）。

3 平差結(jié)果及精度

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

圖2 重力年變率精度統(tǒng)計(jì)結(jié)果Fig.2 The statistic result of the gravity annual rate precision

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

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

與分期靜態(tài)平差模型比較，基于動(dòng)態(tài)平差模型的解算結(jié)果更為合理。這是因?yàn)椋?）分期靜態(tài)平差解算中沒(méi)有考慮觀測(cè)過(guò)程中重力點(diǎn)值隨時(shí)間的變化，該時(shí)間因素引起的系統(tǒng)誤差增加了其解算結(jié)果的不確定性，而動(dòng)態(tài)平差模型中考慮了這一時(shí)間因素的影響；2）基于分期靜態(tài)平差解算獲得的重力年變率為各期次觀測(cè)結(jié)果的等權(quán)擬合結(jié)果，沒(méi)有考慮不同儀器不同年代觀測(cè)資料的精度差異，而基于動(dòng)態(tài)平差模型的整體平差解算對(duì)參與計(jì)算的不同儀器不同年代觀測(cè)資料采用了不同的權(quán)重，并依據(jù)驗(yàn)后信息對(duì)各類觀測(cè)數(shù)據(jù)作了較為客觀的調(diào)整，同時(shí)采用相關(guān)穩(wěn)健估計(jì)有效減弱了粗差的影響，其解算結(jié)果較分期平差結(jié)果更為合理。

4 結(jié) 語(yǔ)

研究結(jié)果表明，隨著測(cè)區(qū)多年復(fù)測(cè)資料的積累，區(qū)域重力背景場(chǎng)長(zhǎng)期變化趨勢(shì)特征逐漸收斂，測(cè)點(diǎn)線性速率反映了測(cè)區(qū)重力多年累積變化趨勢(shì)。在經(jīng)典平差模型中引入初始時(shí)間歷元和測(cè)點(diǎn)的重力變化率，采用單點(diǎn)速率模型，聯(lián)合絕對(duì)觀測(cè)資料，建立段差觀測(cè)值的誤差方程，并依據(jù)最小二乘原理解算法方程，得到測(cè)網(wǎng)中每個(gè)測(cè)點(diǎn)的重力值及重力變化平均速率的平差值，實(shí)現(xiàn)多期觀測(cè)資料共同解算，同時(shí)也實(shí)現(xiàn)了單期網(wǎng)中各個(gè)測(cè)點(diǎn)時(shí)間基準(zhǔn)的統(tǒng)一歸算。

在處理較大面積的測(cè)網(wǎng)多期復(fù)測(cè)資料時(shí)，采用基于動(dòng)態(tài)平差原理的整體解算方法更為合理。本文采用線性速率模型可有效擬合區(qū)域重力場(chǎng)的長(zhǎng)期線性變化趨勢(shì)，為了能夠更準(zhǔn)確地模擬重力場(chǎng)的非線性時(shí)變信息，需要引入更精確的時(shí)變模型。以后的研究中將考慮基于高階項(xiàng)多項(xiàng)式擬合和周期擬合的動(dòng)態(tài)平差模型。

致謝：感謝中國(guó)地震局地震研究所和云南省地震局卓有成效的觀測(cè)工作和數(shù)據(jù)資料清理工作！

［1］李輝.中國(guó)大陸近期重力場(chǎng)動(dòng)態(tài)變化圖像［J］.大地測(cè)量與地 球 動(dòng) 力 學(xué)，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］孫少安，康開軒，黃邦武.關(guān)于區(qū)域重力場(chǎng)變化基準(zhǔn)的思考［J］.大地測(cè)量與地球動(dòng)力學(xué)，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］張祖勝，楊元喜，孫漢榮，等.地殼形變監(jiān)測(cè)中的水準(zhǔn)與重力資料聯(lián)合解算［J］.地震學(xué)報(bào)，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］孫和平，徐建橋，黎瓊.地球重力場(chǎng)的精細(xì)頻譜結(jié)構(gòu)及其應(yīng)用［J］.地球物理學(xué)進(jìn)展，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］李輝，劉冬至，劉紹府.地震重力監(jiān)測(cè)網(wǎng)統(tǒng)一平差模型的建立［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］孫少安，項(xiàng)愛(ài)民，李輝.滇西和北京區(qū)域重力場(chǎng)演化及其與地震的關(guān)系的探討［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］孫少安，項(xiàng)愛(ài)民，吳維日.LCR－G 型重力儀儀器參數(shù)的時(shí)變特征［J］.大地測(cè)量與地球 動(dòng)力 學(xué)，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］劉紹府，李輝，劉冬至.拉科斯特重力儀測(cè)量平差中的相關(guān)問(wèn)題［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］李輝，付廣裕，孫少安，等.滇西地區(qū)重力場(chǎng)動(dòng)態(tài)變化計(jì)算［J］.地殼形變與地震，2000（1）：60－66（Li Hui，F(xiàn)u 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］王勇，張為民，詹金剛，等.重復(fù)絕對(duì)重力測(cè)量觀測(cè)的滇西地區(qū)和拉薩點(diǎn)的重力變化及其意義［J］.地球物理學(xué)報(bào)，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ì)在重力網(wǎng)平差中的應(yīng)用［J］.大地測(cè)量與地球動(dòng)力學(xué)，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）