陳善勇 楊 江,2,3

1 中國地震局地震研究所，武漢市洪山側(cè)路40號，430071 2 武漢地震科學(xué)儀器研究院有限公司，武漢市洪山側(cè)路40號，430071 3 湖北省重大工程地震監(jiān)測與預(yù)警處置工程技術(shù)研究中心，湖北省咸寧市青龍路11號，437000

地球重力場的精準(zhǔn)測量對計量科學(xué)、防震減災(zāi)、大地測量、地球物理等領(lǐng)域具有十分重要的科學(xué)意義。我國一直對重力相關(guān)領(lǐng)域投入巨大，也取得眾多研究成果。尤其是高精密重力儀器的開發(fā)，如超導(dǎo)重力儀、小型絕對重力儀和?？罩亓x等均取得突破性進(jìn)展[1]。高精密重力儀器的發(fā)展對重力數(shù)據(jù)處理方法提出更高的精度需求。

DZW重力儀是我國自行研發(fā)的第一臺高精度重力儀，分辨率可達(dá)到μGal級。由于其采用的是垂直懸掛彈性系統(tǒng)，彈簧張力的衰減以及未被補償或屏蔽的外界作用會產(chǎn)生零漂[2]。DZW重力儀1 d的零漂有時可高達(dá)40 μGal，直接影響到重力儀的觀測精度與準(zhǔn)確性。

經(jīng)典的線性零漂處理方法認(rèn)為，零漂率在觀測過程中為恒值，在1個月以上的長期實際外業(yè)測量中呈非線性[2]。因此，線性零漂處理方法在重力觀測過程中的應(yīng)用方案是以24 h作為時間節(jié)點進(jìn)行分段線性擬合，得到各段的零漂率[3]。然而對于處理幾個月以上的長期重力觀測數(shù)據(jù)，該方法會耗費較多時間和物力。如果1 d數(shù)據(jù)中有溫度、氣壓等其他非線性因素，也會導(dǎo)致經(jīng)過線性零漂校正后觀測值的精度難以保證。

文獻(xiàn)[4-5]通過建立貝葉斯平差模型來解決非線性漂移觀測誤差，從而為提高數(shù)據(jù)精度提供新的解決方案。本文在其基礎(chǔ)上進(jìn)行優(yōu)化嘗試，貝葉斯平差模型中的研究對象為零漂率，本文貝葉斯模型將其變?yōu)榱闫?。這是由于零漂與重力觀測值之間只涉及到加減運算，不涉及與時間的乘法運算，所需要的參量只有重力觀測值以及理想的重力儀讀數(shù)，因此可以簡化運算并增加算法的計算效率。

為將本文算法運用到實際重力觀測中，我們設(shè)計一套關(guān)于該算法的應(yīng)用方案，使得該算法可以在每次獲得新觀測數(shù)據(jù)后自動進(jìn)行計算，增加算法的實用性。該方案將算法運用到重力觀測數(shù)據(jù)的實時處理中，在采樣重力觀測值的同時進(jìn)行零漂計算，可以減少后期數(shù)據(jù)處理中消除零漂所消耗的人力和物力[6]。

1 模型構(gòu)建

零漂的貝葉斯估計模型需要先定義與算法相關(guān)的矩陣以及基本公式，模型的核心公式為貝葉斯公式，需要分別確定先驗分布形式、似然函數(shù)以及邊緣密度函數(shù)。

1.1 矩陣定義

將經(jīng)過潮汐校正的相對重力觀測數(shù)據(jù)表示為：

x=[x1,x2,…,xn]T

將每個樣本點對應(yīng)的零漂表示為：

v=[v1,v2,…,vn]T

將每個樣本點的重力儀讀數(shù)(理想值)表示為：

g=[g1,g2,…,gn]T

將用于得到一組數(shù)據(jù)的梯度變化量的矩陣A(A∈Rn×n)表示為：

式中，ζ為方陣的補充值，取值趨近于0。

1.2 基本公式

對于經(jīng)過潮汐校正后的相對重力觀測數(shù)據(jù)，其與零漂之間的關(guān)系式可表示為[7]：

x-v-g=δ

(1)

式中，δ為公式計算產(chǎn)生的誤差，滿足以下分布：

δ～N(0,Σ1)

(2)

式中，Σ1為高斯分布的協(xié)方差矩陣。

由式(1)可推出以下分布：

x～N(v+g,Σ1)

(3)

1.3 先驗分布確定

對零漂的先驗分布形式進(jìn)行確認(rèn)，每個時刻的零漂是未知的，因此存在隨機性[8-9]，可用隨機變量表示：

(4)

假設(shè)零漂在連續(xù)時間中具有連續(xù)性[11]，則需要滿足：

Av=ξ

(5)

式中，ξ為趨向于0的數(shù)組，滿足以下分布：

ξ～N(0，Σ2)

(6)

式中，Σ2為高斯分布的協(xié)方差矩陣。

由式(5)可推出以下分布：

Av～N(0，Σ2)

(7)

已知v為高斯過程，根據(jù)高斯過程的線性性質(zhì)，由式(7)可推算出v的具體分布為：

v～N(A-1·0,A-1Σ2(AT)-1)?

v～N(0,A-1Σ2(AT)-1)

(8)

先驗分布的概率密度函數(shù)為：

(9)

1.4 似然函數(shù)確定

似然函數(shù)是樣本信息的集中體現(xiàn)，由式(3)可推出其表達(dá)式為：

(10)

1.5 邊緣密度函數(shù)確定

邊緣密度函數(shù)本身并不會影響后驗分布的分布形式，但會影響后驗分布的幅值變化[12]，可通過式(9)、(10)推出邊緣密度函數(shù)的表達(dá)式為：

(11)

1.6 后驗分布計算

將先驗分布的概率密度函數(shù)、似然函數(shù)和邊緣密度函數(shù)代入貝葉斯公式，推出后驗分布的概率密度函數(shù)：

(12)

正態(tài)均值(方差已知)的共軛先驗分布為正態(tài)分布。后驗分布滿足以下分布形式：

(13)

1.7 超參數(shù)確定

由式(13)可知，要得到確定的后驗分布，需要確定兩組超參數(shù)值，分別是Σ1和Σ2。其中，Σ1為樣本所產(chǎn)生的方差，代表重力儀的屬性，視為恒定值[13]，可以取其中1 d的數(shù)據(jù)進(jìn)行線性擬合得到零漂率，從而得到零漂校正后的數(shù)據(jù)方差；Σ2為零漂先驗分布的超參數(shù)，可通過增加計算次數(shù)得到收斂極限進(jìn)行確定，初始值可以取與Σ1同一量級的數(shù)值。

對于零漂后驗分布均值和方差的確定，方法如下：

(14)

經(jīng)過貝葉斯公式計算得到的零漂后驗分布方差是兩者的綜合體現(xiàn)。將該后驗分布作為下一次貝葉斯公式的先驗分布，則得到的第二次后驗分布的方差可表示為：

(15)

第n次后驗分布的方差為：

(16)

由式(16)可知，后驗分布方差會隨著計算層數(shù)的增加而不斷減小并趨向于零，表明后驗分布具有收斂性，并且收斂極限為0。

后驗分布的均值經(jīng)過化簡可表示為：

(17)

2 算法數(shù)據(jù)驗證

本文算法在經(jīng)過前文數(shù)學(xué)邏輯推導(dǎo)后，需要通過實際重力數(shù)據(jù)來進(jìn)一步驗證。首先，對算法中超參數(shù)進(jìn)行初步估計，通過數(shù)據(jù)進(jìn)一步驗證零漂超參數(shù)的收斂性；然后，通過數(shù)據(jù)來確認(rèn)均值的收斂極限性質(zhì)。在算法本身邏輯性得到證明的基礎(chǔ)上，需要確定應(yīng)用方案來模擬測量過程，對所得結(jié)果的可用性和準(zhǔn)確性進(jìn)行分析；最后，通過與傳統(tǒng)線性擬合去零漂方案進(jìn)行對比分析，探討算法是否實現(xiàn)去非線性零漂的目標(biāo)。

2.1 超參數(shù)初步估計

將九峰臺站DZW重力儀在2019年的連續(xù)觀測數(shù)據(jù)作為算法的驗證數(shù)據(jù)源。該數(shù)據(jù)已經(jīng)過精度為10-5μGal的固體潮校正，可不考慮溫度、氣壓等因素的影響。通過線性擬合，得到圖1中表示線性零漂的直線，與潮汐校正后的數(shù)據(jù)點進(jìn)行對比。擬合直線斜率即為估計的零漂率，為0.014 25 μGal/min=0.855 μGal/h，由給出的均方根誤差RMSE(root mean squared error)推出樣本的數(shù)據(jù)方差Σ1=RMSE2=0.184 92μGal2=0.034 19 μGal2，Σ1即為0.034 19的對角矩陣。由于零漂的先驗分布的協(xié)方差矩陣Σ2為未知參數(shù)，但在多層運算后將會趨近于0，因此不妨設(shè)其為與Σ1同數(shù)量級的矩陣，令其為0.01的對角矩陣。

圖1 對1 d重力觀測數(shù)據(jù)進(jìn)行線性擬合所得結(jié)果Fig.1 The results obtained by linear fitting of one-day gravity observation data

2.2 超參數(shù)收斂性驗證

為確保計算次數(shù)足夠多，將計算次數(shù)設(shè)置為1 000。1 000次計算中后驗分布的方差變化如圖2所示。

圖2 后驗分布方差變化Fig.2 Variance variation of posteriori distribution

在前100次計算中，方差下降速度較快；在后900次計算中，下降速度趨于平緩，總體呈收斂態(tài)勢。為進(jìn)一步驗證均值的收斂性，需要研究均值的收斂方向和收斂極限。為遍歷所有點的均值收斂情況，將計算次數(shù)減少至100。經(jīng)過研究發(fā)現(xiàn)，每個樣本點的收斂方向不同，收斂極限也存在差異。以第11個和第1 003個樣本點的均值變化為例，結(jié)果見圖3(a)和3(b)。

圖3 樣本點均值變化Fig.3 Meanvalues change of sample points

通過上述分析可知，隨著計算次數(shù)不斷增加，后驗分布的方差趨近于0，各樣本點的均值沿著各自方向趨近于不同的收斂極限。但收斂極限是否為零漂的真實值需進(jìn)一步確認(rèn)。

2.3 均值收斂極限性質(zhì)確認(rèn)

將重力觀測值與后驗分布均值作差，理論結(jié)果為真實重力觀測值，是一個恒值，結(jié)果如圖4所示。

圖4 觀測數(shù)據(jù)在去除不同計算次數(shù)所得零漂結(jié)果Fig.4 Drift results of observed data after removing different calculation times

圖4(a)中方差為0.000 343 μGal2，圖4(b)中方差為1.803 78×10-8μGal2，后者方差明顯小于前者，說明隨著計算次數(shù)的增加，均值不斷逼近，使去零漂結(jié)果為常數(shù)數(shù)值，均值的收斂極限即為真實的零漂值。

2.4 算法應(yīng)用方案確定

第1天的1 440個樣本點采用貝葉斯模型進(jìn)行預(yù)處理，通過對其進(jìn)行1 000次重復(fù)計算，得到精準(zhǔn)的1 440個零漂值。對于新樣本數(shù)據(jù)，通過左移操作實現(xiàn)樣本數(shù)組的更新，為確保樣本點與零漂值一一對應(yīng)，零漂的均值矩陣與協(xié)方差矩陣中的對角元素數(shù)組同時進(jìn)行左移操作(圖5)。

圖5 算法應(yīng)用方案Fig.5 Scheme of algorithm application

圖5中v0作為零漂數(shù)組的輸出數(shù)據(jù)按序存放，形成的數(shù)據(jù)組即為精準(zhǔn)的零漂數(shù)組。xn+1為新樣本值，vn+1代表與新樣本值對應(yīng)的均值，是需要預(yù)測的參數(shù)。Sn+1代表與新樣本值對應(yīng)的方差，等于先驗分布的方差。該方案的優(yōu)點是對于新移入的樣本值，其均值由移入至最后輸出會經(jīng)歷1 440次貝葉斯算法計算，可保證結(jié)果的高精度。

2.5 貝葉斯算法結(jié)果分析

貝葉斯算法得到連續(xù)2 d的零漂以及零漂校正后的重力觀測值如圖6所示。

圖6 貝葉斯算法所得結(jié)果Fig.6 Bayesian algorithm results

由圖6(b)可知，樣本矩陣、零漂矩陣、協(xié)方差矩陣的不斷變化，會使第2天的重力觀測值方差更大，第1天重力觀測值方差為8.458 01×10-12μGal2，而第2天為1.595 53×10-10μGal2，方差數(shù)量級已符合高精度的要求。

2.6 干擾作用下的結(jié)果分析

人為添加白噪聲作為干擾信號加入到重力觀測值中，模擬在觀測過程中遇到外部干擾的情況。通過算法得到連續(xù)2 d的零漂結(jié)果以及零漂校正后的重力觀測值如圖7所示。

圖7 加入干擾后的算法結(jié)果Fig.7 Algorithm results after adding interference

由圖7(a)可知，干擾信號會對零漂結(jié)果造成一定影響，影響在可控范圍內(nèi)的原因在于先驗分布的協(xié)方差矩陣會約束零漂的變動幅度。由圖7(b)可知，經(jīng)過零漂校正后的重力觀測值仍會將干擾信號特征體現(xiàn)出來，同時也說明算法具有一定的抗干擾能力，在干擾結(jié)束后仍然可以精準(zhǔn)地計算零漂值。

2.7 線性掉格下的結(jié)果分析

人為在重力觀測數(shù)據(jù)中添加線性參數(shù)，模擬重力儀由于人員操作或溫度驟變等外界因素而使儀器出現(xiàn)線性掉格的情況，結(jié)果如圖8(a)所示。

由圖8(b)可知，重力儀出現(xiàn)掉格后，該算法的去零漂結(jié)果會出現(xiàn)微小跳變，經(jīng)過300個數(shù)據(jù)點后零漂校正結(jié)果回歸正常，說明該算法需要300個數(shù)據(jù)點進(jìn)行恢復(fù)零漂計算，但由于跳動較小，因此可以視為在出現(xiàn)線性掉格時，該算法仍然可以正常工作并保持準(zhǔn)確度。

圖8 線性掉格情況下算法結(jié)果Fig.8 Algorithm results in case of linear drifting

2.8 經(jīng)典算法結(jié)果對比

通過線性擬合對同樣2 d的連續(xù)重力觀測值進(jìn)行零漂估計，擬合結(jié)果如圖9所示。

圖9 連續(xù)2 d重力觀測值線性擬合結(jié)果Fig.9 Linear fitting results of gravity observations for two consecutive days

線性擬合的零漂率為0.014 02 μGal/min ，由此得到去除零漂的重力觀測值如圖10所示。

圖10 線性擬合算法零漂校正Fig.10 Drift correction of linear fitting algorithm

對比圖10與圖6(b)發(fā)現(xiàn)，貝葉斯算法可將重力觀測值中的非線性零漂去除，實現(xiàn)對非線性零漂的精準(zhǔn)估計。

3 結(jié) 語

本文提出一種基于貝葉斯公式，以零漂的連續(xù)性作為先驗約束條件，考慮重力觀測數(shù)據(jù)浮動產(chǎn)生的偏差，將每個采樣點零漂的均值和方差作為超參數(shù)進(jìn)行求解，從而獲得零漂估計值的算法。通過零漂與重力觀測數(shù)據(jù)的關(guān)系式以及零漂的連續(xù)性得到零漂的分布規(guī)律，采用邏輯推導(dǎo)和數(shù)據(jù)驗證證明，貝葉斯算法得到零漂后驗分布的均值會隨著計算次數(shù)的增加而向零漂真實值收斂。在模擬的數(shù)據(jù)采樣過程中，通過貝葉斯算法得到零漂校正結(jié)果的方差為1.595 53×10-10μGal2，符合DZW重力儀的精度需求。在添加人為干擾后得到的零漂校正結(jié)果表明，該算法可以在保留重力信號特征的前提下精準(zhǔn)去除零漂值。

本文提出的非線性零漂貝葉斯估計方法能有效地去除重力觀測值中的非線性零漂，可以在重力分采樣過程中實時對重力觀測值進(jìn)行預(yù)處理，從而為重力觀測精度提供保障。