云南流動重力觀測中相對重力儀漂移估計方法

黃江培　杜家云　曹穎　鄭秋月　劉東　吳宇琴　王青華

摘要：針對在陸地百千米以上的時變流動重力觀測網絡中，相對重力儀非線性漂移率和儀器性能偏差是影響觀測結果精度的主要誤差來源這一問題，以1990—2022年云南流動重力觀測數據為例，對比傳統最小二乘平差方法和貝葉斯平差方法得到的非線性漂移率結果，系統測評了貝葉斯估計方法對于提高重力平差結果精度的有效性。結果表明，可以通過貝葉斯優化方法給出合適的儀器權重和非線性漂移率，減少多臺相對重力儀性能差異造成的誤差和傳統線性漂移計算方法造成的有效漂移遺失，提高平差結果的精度。

關鍵詞：漂移率；貝葉斯估計；流動重力觀測；平差結果；云南

中圖分類號：P315.62文獻標識碼：A文章編號：1000-0666（2023）04-0511-10

doi：10.20015/j.cnki.ISSN1000-0666.2023.0060

0引言

由于地球的不可入性，高精度時變重力被視為探測地球內部密度變化的主要手段，被廣泛地應用于各個科學研究領域，包括礦物探測（李小孟，曾華霖，1996）、地殼構造特征（陳石，王溪身，2015，汪健等，2015）、地下水變化特征（佘雅文等，2015）和地震研究（盧造勛等，1978；陳運泰等，1980；吳國華等，1995，1997；祝意青等，2015，2018；孫少安等，2015；胡敏章等，2019；Xing et al，2021；劉洪良等，2021；黃江培等，2022）等方面。但在高精度重力獲取過程中，受外部環境、儀器設備性能等影響，其數據質量成為眾多學者最關心的問題之一。外部環境中的溫度、氣壓、固體潮等影響因素目前已經有被國際學會接受的經驗改正公式進行解算；由于個體差異極大，儀器性能方面的一次項、格值系數、漂移率等特征，無法進行統一規定，只能通過儀器廠商的出廠設置、觀測前的儀器標定、觀測數據的似然估計綜合獲取。

2010年以前，中國大陸的流動重力觀測網絡主要采用以Lacoste（拉科斯特）和Burries（貝爾雷斯）為主的金屬彈簧重力儀進行流動重力觀測工作，金屬彈簧型的漂移率基本為（1～3）×10^-8（m·s^-2）/h（邢樂林等，2010；趙云峰等，2018），所以這個階段主要是針對儀器一次項及格值系數的改正研究。2010年以后，CG-5型石英彈簧相對重力儀被引進并廣泛應用，其采用自動傾斜補償、高精度溫度控制、電子讀數等先進技術，出廠設置將格值進行了處理，參數估計中只需要考慮一次項及漂移率即可，而且在一個觀測周期內，可以采用相同的一次項系數解算（郝洪濤等，2016，黃江培等，2020）。

2009年6 月，中國地震局在武漢及江西九江廬山基線場對10 臺CG-5 相對重力儀進行測試，顯示其靜態漂移率呈較好的線性，平均靜態漂移率與平均動態漂移率符合性較好，個別儀器線性漂移率較大，最大超過100×10^-8（m·s^-2）/h（邢樂林，李輝，2010）。通過本次標定，發現平均動態漂移率與靜態線性漂移率吻合較好，認為將動態漂移進行線性化計算是可行的。但在實際工作中，CG-5相對重力儀的零點漂移并不是線性的（汪健等，2016；楊雅慧等，2021），并且漂移率變化可達到10×10^-8（m·s^-2）/h（隗壽春等，2016，2017），因此在一個觀測周期內將漂移率當做一個定值進行平差的方法是不恰當的。針對這種情況，學者們提出了分段平差法（隗壽春等，2016）、基于貝葉斯準則的似然估計法（Chen et al，2019）等，其中貝葉斯平差方法的有效性得到了部分學者的驗證（王林海等，2020；楊錦玲等，2021；Zheng et al，2022）。本文基于云南地區近30年的流動重力觀測數據，分別采用基于傳統最小二乘原理的平差方法計算線性漂移及基于貝葉斯原理的平差方法計算非線性漂移率，提出合理的相對重力儀零漂參數估計方案，以期為高精度時變重力數據處理提供參考。

1數學模型

1.1漂移率

相對重力儀零點漂移是指儀器零點隨著彈簧的老化、彈性疲勞、測值段的突變或在運輸過程中發生顛簸等因素出現的偏移，在不考慮觀測誤差和測點重力實際變化的情況下，可以理解為同一測點在不同時間段重復觀測的差值。漂移率則是在單位時間內的零點漂移量，一般取1 h為時間單位。設觀測值為g，觀測時間間隔為t，漂移率v計算公式為（隗壽春等，2017）：

1.2最小二乘平差方法

1.3貝葉斯平差方法

2數據分析

2.1漂移率及點值誤差

本文收集1990—2022年云南地區的流動重力觀測數據，采用最小二乘平差方法（ADJ）及貝葉斯平差方法（BAY）分別計算其漂移率，統計結果見表1，表中ADJ即為整期漂移率，BAY取整期觀測中每日漂移率的平均值。其中1990—2013年所用儀器為Lacoste金屬彈簧相對重力儀，2014年以后采用CG-5石英彈簧相對重力儀，2014年所用數據為在儀器中設置漂移率后的觀測數據，2015年以后所用數據為未進行漂移率修正的數據。從表1可見，不管使用哪種型號的儀器，通過BAY和ADJ計算的每日漂移率的平均值都非常接近。

從表1中可以看出，2014年以前，儀器漂移率基本都在4×10^-8（m·s^-2）/h以內，兩種方法計算的點值平均誤差也沒有明顯差距；2014年以后，CG-5相對重力儀的漂移率較大，且漂移率明顯是變化的，在一個觀測周期內，將漂移率當做一個定值處理是欠妥的，ADJ計算結果的點值誤差基本大于BAY計算結果。

2.2漂移率變化分析

由表1可以看出，一個觀測周期內，使用BAY計算的非線性漂移率的數學期望就是使用ADJ計算的線性漂移率，說明BAY計算結果是可靠的。本文以2020-03期為例，進一步分析1個測期內漂移率的變化情況，圖1為2020-03期觀測網示意圖。

野外數據采集過程中，采用A→B→C…N…C→B→A的測線往返觀測模式，一條測線往返觀測時間不超過72 h，其中A→B與B→A構成一個往返測段，2020-03測期共計完成82條測線、275個測段，歷時91 d。如果不考慮偶然誤差，根據式（1）可以計算出每個測段及每條測線的儀器漂移率，然后根據大量的漂移率樣本量，擬合出漂移率的變化趨勢，如圖2所示。為了與BAY計算結果對應，測段及測線編號均按照時間先后排序。

根據式（1）計算的每個測段漂移率，是將每個往返兩次觀測的不等差均當做漂移處理，即誤差全帶入，此時如果觀測間隔較短，實際漂移不大，求取的漂移率中大部分是由誤差引起的。2020-03期全網觀測點值平均精度為8.8×10^-8m/s²，全網漂移率為26.3×10^-8（m·s^-2）/h。根據誤差理論，一般取大于2倍誤差的變化值作為有效數據，即變化大于17.6×10^-8m/s²才能視為計算的漂移率是有效的，通過與漂移率取商，得出測段觀測時間間隔需達到40 min以上，才能計算有效漂移率。在實際觀測中有部分測段觀測間隔低于40 min，所以圖2a中，存在部分離散幅度非常大的測點，可理解為觀測誤差大于實際漂移。圖2b中采用往返測線作為計算基礎，每條測線均包含2段以上測段，部分高達10多個，每個測段均是一次獨立觀測，偶然誤差能夠疊加消除一部分，計算采用的有效間隔時間也較長，所以計算的漂移率離散度明顯降低。雖然測段及測線計算的結果均是將誤差全帶入的結果，不過當樣本量較大時，通過擬合得出的趨勢是準確的。從圖2可以看出，CG-1170儀器的漂移率近似于線性，變化不大，但CG-1169儀器的漂移率有明顯的變化趨勢。圖3為整網計算的使用ADJ和BAY得到的漂移率。

從圖3可以看出，BAY結果與所有測段和測線計算的漂移率變化趨勢（圖2）一致，其中使用BAY計算的CG5-1170儀器的漂移率從4月7日以后一直處于ADJ計算結果的上方，且差值在不斷增大；在5月17日之前，使用BAY計算的CG5-1169儀器的漂移率處于ADJ計算結果附近上下波動的狀態，但是從5月18日之后，一直處于ADJ計算結果之上，并且差值也在不斷增大。在ADJ平差計算時，偏離部分的有效漂移率被當做誤差處理了，而BAY算法能夠極大的擬合這些有效漂移率。

2.3點值精度對比

從表1能夠發現，使用CG-5型相對重力儀后，BAY全網計算結果的平均點值精度均優于ADJ計算結果，本文采用2020-03測期數據，分析兩種方法計算結果的具體點值精度情況。

從表2可以看出，在σ＜5區間，BAY計算結果的點值精度明顯大于ADJ計算結果；在σ＜7區間，BAY結果測點占比達到86％，ADJ測點占比為60.2％；在σ≥9區間，BAY只有4個測點，ADJ有19個測點，所以，BAY計算方法對點值精度有優化效果。本文在空間分布上分析其優化能力如圖4所示。

結合圖1和圖4可以看出，ADJ計算結果點值誤差主要分布在測網邊緣和遠離控制點的地方，特別是支線上，點值精度較差，這是符合最小二乘平差中邊緣效應原理的；BAY計算結果點值誤差的空間分布趨勢與ADJ計算結果是一致的，但是在邊緣效應上有明顯的優化結果，特別是在支線的處理上，能夠明顯降低點值誤差。

2.4實例分析

2021年云南漾濞MS6.4地震后，劉東等（2021）和黃江培等（2022）采用相同的相對重力歷史數據，分別采用ADJ和BAY方法進行解算，并分別對震前相對重力變化特征進行了分析，本文對兩份解算結果進行對比分析，如圖5所示。由圖5可以看出，兩種方法計算結果采用相同的繪圖范圍和歷史數據，并使用相同的色標?？傮w來看，兩種方法解算的重力變化趨勢是一致的，漾濞地震震中西南側的局部區域異常不一致是兩位學者對重力變化高頻噪聲的取舍不同所致。

從圖5b可以看出，在震中附近，NW-SE方向上形成了一個弱四象限趨勢，沿紅河斷裂帶是以負變化為主。震中附近，紅河斷裂兩側，有正變化趨勢，地震發生于零值線附近，符合學者們的部分經驗結論（祝意青等，2018；胡敏章等，2019），所以認為本次地震前相對重力變化特征分析時，基于BAY所得結果是能夠有效捕捉震前異常的。

3討論

陸地高精度相對重力觀測中，觀測段差往返不符值中包含零點漂移及觀測誤差，如何恰當地區分兩者一直是一個難題。在使用Lacoste儀器的年代，漂移率較小且變化也小，經典最小二乘平差方法由于計算方便，計算機語言也比較簡潔，無疑是最優的平差方法。但是CG-5型相對重力儀投入使用后，其漂移率較大，且非線性漂移明顯，個體差異也大。從表1可以看出，2015—2022年CG5-1169和CG5-1170儀器的漂移率有明顯的下降，上下半年也有明顯的差距，但是隨著時間的推移，下降速率有收斂的跡象，上下半年差距也在減小。所以，從長時間來看，漂移率是變化的，且漂移率變化也是非線性的，這種情況下，把漂移率當做一個定值處理是不恰當的。而BAY能夠較好地反映漂移率的變化趨勢（圖3、4），并且受到光滑矩陣的約束，加強了其穩健性，部分離散度較大的漂移率對其影響并不明顯。

地震重力的研究是基于高精度重力觀測，研究10×10^-8m/s級的重力變化，對于云南測區這種數百千米以上的大尺度構造環境的監測，閉合時間長、儀器多、測網復雜的情況下，漂移率的變化對計算結果影響較大，更需要BAY這種能夠計算非線性漂移率，并且根據儀器性能自動分配儀器權重的算法。

4結論

本文利用傳統最小二乘平差方法和貝葉斯平差方法計算1990—2022年云南流動重力觀測數據的漂移率，并以2021年云南漾濞MS6.4地震為例進行分析，主要得出以下結論：

（1）流動重力觀測數據的測段時間間隔需要達到40 min以上，才能夠計算有效漂移率。

（2）對于漂移較小且非線性變化不明顯的儀器觀測結果，ADJ及BAY計算結果是一樣的。

（3）對于漂移較大且非線性變化的儀器觀測結果，ADJ把漂移率當做一個定值進行平差，會把部分有效漂移率當做誤差處理，造成計算結果點值誤差增大；BAY能夠計算出每日漂移率的極大似然值，較好地擬合漂移率變化趨勢，具有更好的數據自洽性，計算結果點值精度能夠得到明顯優化。

（4）在測網邊緣和遠離控制點、支線等控制較弱的位置，BAY處理結果的點值精度優于DAJ。

綜上所述，BAY適用于云南測區這種空間跨度大、時間周期長的測網，由于BAY需要解算的參數較多，計算機解算過程較長，對于儀器數量少、省網級的解算是值得推薦的，但是對于儀器數量多且型號不一的全國聯測網的解算，由于對計算機配置要求較高，并且解算耗時較長，目前還沒有學者進行過嘗試，可在將來進一步研究。

參考文獻：

陳石，王謙身.2015.蒙古及周邊地區重力異常和地殼不均勻體分布［J］.地球物理學報，58（1）：79-91.

陳運泰，顧浩鼎，盧造勛，等.1980.1975年海城地震與1976年唐山地震前后的重力變化［J］.地震學報，2（1）：21-31.

郝洪濤，李輝，孫和平，等.2016.CG-5重力儀零漂改正及格值系數檢測應用研究［J］.武漢大學學報（信息科學版），41（9）：1265-1271.

胡敏章，郝洪濤，李輝，等.2019.地震分析預報的重力變化異常指標分析［J］.中國地震，35（3）：417-430.

黃江培，曹穎，劉東，等.2022.漾濞MS6.4地震前后的重力變化特征及其孕震含義分析［J］.地震地質，44（6）：1557-1573.

黃江培，王青華，徐聲鑫，等.2020.CG-5重力儀一次項系數變化特性分析及其對觀測數據的影響研究［J］.地震研究，43（1）：101-108.

李小孟，曾華霖.1996.高精度重力資料在勝利油區油氣藏探測中的應用［J］.現代地質，10（2）：250-259.

劉東，郝洪濤，王青華，等.2021.2021年云南漾濞MS6.4地震前重力變化［J］.地震地質，43（5）：1157-1170.

劉洪良，王青華，張展偉，等.2021.河北唐山古冶5.1級地震前的重力變化［J］.華南地震，41（2）：71-75.

盧造勛，方昌流，石作亭，等.1978.重力變化與海城地震［J］.地球物理學報，21（1）：1-8.

佘雅文，付廣裕，韋進，等.2015.十三陵地震臺gPhone重力儀的儀器性能與水文響應分析［J］.大地測量與地球動力學，35（5）：901-905.

孫少安，郝洪濤，韋進，等.2015.云南景谷M6.6地震前重力場變化的區域性特征［J］.大地測量與地球動力學，35（4）：613-615.

汪健，孫少安，邢樂林，等.2016.CG-5重力儀的漂移特征［J］.大地測量與地球動力學，36（6）：556-560.

汪健，王安怡，申重陽，等.2015.南北地震帶南段莫霍面重力反演研究［J］.大地測量與地球動力學，35（6）：931-935.

王林海，陳石，莊建倉，等.2020.精密重力測量中相對重力儀格值系數的貝葉斯估計方法［J］.測繪學報，49（12）：1543-1553.

隗壽春，徐建橋，郝洪濤，等.2017.零漂改正對中國地殼運動觀測網絡重力數據處理的影響［J］.大地測量與地球動力學，37（4）：403-406.

隗壽春，徐建橋，周江存，等.2016.重力網的分段線性動態平差［J］.測繪學報，45（5）：511-520

吳國華，羅增雄，賴群，等.1995.1988年瀾滄─耿馬地震與滇西實驗場的重力變化［J］.地殼形變與地震，（2）：66-73.

吳國華，羅增雄，賴群，等.1997.麗江7.0級地震前后滇西實驗場的重力異常變化特征［J］.地震研究，（1）：103-109.

邢樂林，李輝，夏正超，等.2010.CG-5重力儀零漂特性研究［J］.地震學報，32（3）：369-373.

邢樂林，李輝.2010.Burris重力儀性能研究［C］//中國地球物理2010——中國地球物理學會第二十六屆年會、中國地震學會第十三次學術大會論文集，北京，851.

楊錦玲，陳石，王林海，等.2021.華南陸地時變重力觀測數據質量評估［J］.測繪學報，50（3）：333-342.

楊雅慧，劉洪良，郝洪濤，等.2021.CG-5重力儀格值系數修正與應用研究［J］.華南地震，41（4）：63-68.

趙云峰，祝意青，梁偉鋒，等.2018.Burris型重力儀性能分析［J］.地震地磁觀測與研究，39（2）：178-185.

祝意青，劉芳，李鐵明，等.2015.川滇地區重力場動態變化及其強震危險含義［J］.地球物理學報，58（11）：4187-4196.

祝意青，劉芳，徐云馬，等.2018.重力監測在地震預報中的應用與展望［J］.國際地震動態，（8）：9-10.

Chen S，Zhuang J C，Li X Y，et al.2019.Bayesian approach for network adjustment for gravity survey campaign：Methodology and model test［J］.Journal of Geodesy，93（5）：681-700.

Nelder J，Mead R.1965.A simplex method for function minimization［J］.The Computer Journal，8（1）：308-313.

Xing L L，Liu Z W，Jia J G，et al.2021.Far－filed coseismic gravity changes related to the 2015 MW7.8 Nepal（Gorkha）earthquake observed by superconducting gravimeters in China continent［J］.Earth and Planetary Physics，5（2）：141-148.

Zheng Q，Yao X，Chen S，et al.2022.Data quality assessment of time-variable surface microgravity surveys in the Southeastern Tibetan Plateau［J］.Appl Sci，12：3310.

The Estimation Method of the Zero Drift of the Relative Gravimeter

in High－precision Time－varying Gravity Observation in Yunnan

HUANG Jiangpei DU Jiayun CAO Ying ZHENG Qiuyue LIU Dong WU Yuqing WANG Qinghua

（1.Yunnan Earthquake Agency，Kunming 650224，Yunnan，China；2.Mile Earthquake Agency，Mile 652399，Yunnan，China）

Abstract

In the time－varying gravity observation network covering hundreds kilometers of land，the nonlinear drift rate of the relative gravimeter and the instrument performance deviation are the main error sources that affect the accuracy of the observational results.In this paper，the drift rates of the mobile－gravity observational data from 1990 to 2022 in Yunnan are calculated respectively with the Bayesian gravity adjustment method and the traditional least squares adjustment method.Then the non－linear dirft rate results from these two methods are compared，and the effectiveness of the Bayesian estimation method in improving the accuracy of gravity adjustment results is systematically evaluated.The results show that the appropriate instrument weight and nonlinear drift rate can be given by the Bayesian optimization method，in order to reduce the error caused by the performance difference of multiple relative gravimeters，and reduce the effective－drift loss caused by the traditional linear drift calculation method.In this way the accuracy of the adjustment results is improved.

Keywords：drift rate；the Bayesian estimation；mobile－gravity observation；adjustment results；Yunnan