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

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

關(guān)鍵詞：漂移率；貝葉斯估計；流動重力觀測；平差結(jié)果；云南

中圖分類號：P315.62文獻(xiàn)標(biāo)識碼：A文章編號：1000-0666（2023）04-0511-10

doi：10.20015/j.cnki.ISSN1000-0666.2023.0060

0引言

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

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

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

1數(shù)學(xué)模型

1.1漂移率

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

1.2最小二乘平差方法

1.3貝葉斯平差方法

2數(shù)據(jù)分析

2.1漂移率及點值誤差

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

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

2.2漂移率變化分析

由表1可以看出，一個觀測周期內(nèi)，使用BAY計算的非線性漂移率的數(shù)學(xué)期望就是使用ADJ計算的線性漂移率，說明BAY計算結(jié)果是可靠的。本文以2020-03期為例，進(jìn)一步分析1個測期內(nèi)漂移率的變化情況，圖1為2020-03期觀測網(wǎng)示意圖。

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

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

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

2.3點值精度對比

從表1能夠發(fā)現(xiàn)，使用CG-5型相對重力儀后，BAY全網(wǎng)計算結(jié)果的平均點值精度均優(yōu)于ADJ計算結(jié)果，本文采用2020-03測期數(shù)據(jù)，分析兩種方法計算結(jié)果的具體點值精度情況。

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

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

2.4實例分析

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

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

3討論

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

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

4結(jié)論

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

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

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

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

（4）在測網(wǎng)邊緣和遠(yuǎn)離控制點、支線等控制較弱的位置，BAY處理結(jié)果的點值精度優(yōu)于DAJ。

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

參考文獻(xiàn)：

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

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

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

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

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

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

李小孟，曾華霖.1996.高精度重力資料在勝利油區(qū)油氣藏探測中的應(yīng)用［J］.現(xiàn)代地質(zhì)，10（2）：250-259.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

祝意青，劉芳，李鐵明，等.2015.川滇地區(qū)重力場動態(tài)變化及其強(qiáng)震危險含義［J］.地球物理學(xué)報，58（11）：4187-4196.

祝意青，劉芳，徐云馬，等.2018.重力監(jiān)測在地震預(yù)報中的應(yīng)用與展望［J］.國際地震動態(tài)，（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