DZW重力儀非線性零漂的貝葉斯估計方法

陳善勇 楊 江,2,3

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

地球重力場的精準測量對計量科學、防震減災、大地測量、地球物理等領域具有十分重要的科學意義。我國一直對重力相關領域投入巨大，也取得眾多研究成果。尤其是高精密重力儀器的開發，如超導重力儀、小型絕對重力儀和海空重力儀等均取得突破性進展[1]。高精密重力儀器的發展對重力數據處理方法提出更高的精度需求。

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

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

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

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

1 模型構建

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

1.1 矩陣定義

將經過潮汐校正的相對重力觀測數據表示為：

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

將每個樣本點對應的零漂表示為：

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

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

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

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

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

1.2 基本公式

對于經過潮汐校正后的相對重力觀測數據，其與零漂之間的關系式可表示為[7]：

x-v-g=δ

(1)

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

δ～N(0,Σ1)

(2)

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

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

x～N(v+g,Σ1)

(3)

1.3 先驗分布確定

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

(4)

假設零漂在連續時間中具有連續性[11]，則需要滿足：

Av=ξ

(5)

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

ξ～N(0，Σ2)

(6)

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

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

Av～N(0，Σ2)

(7)

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

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

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

(8)

先驗分布的概率密度函數為：

(9)

1.4 似然函數確定

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

(10)

1.5 邊緣密度函數確定

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

(11)

1.6 后驗分布計算

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

(12)

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

(13)

1.7 超參數確定

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

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

(14)

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

(15)

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

(16)

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

后驗分布的均值經過化簡可表示為：

(17)

2 算法數據驗證

本文算法在經過前文數學邏輯推導后，需要通過實際重力數據來進一步驗證。首先，對算法中超參數進行初步估計，通過數據進一步驗證零漂超參數的收斂性；然后，通過數據來確認均值的收斂極限性質。在算法本身邏輯性得到證明的基礎上，需要確定應用方案來模擬測量過程，對所得結果的可用性和準確性進行分析；最后，通過與傳統線性擬合去零漂方案進行對比分析，探討算法是否實現去非線性零漂的目標。

2.1 超參數初步估計

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

圖1 對1 d重力觀測數據進行線性擬合所得結果Fig.1 The results obtained by linear fitting of one-day gravity observation data

2.2 超參數收斂性驗證

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

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

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

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

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

2.3 均值收斂極限性質確認

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

圖4 觀測數據在去除不同計算次數所得零漂結果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，后者方差明顯小于前者，說明隨著計算次數的增加，均值不斷逼近，使去零漂結果為常數數值，均值的收斂極限即為真實的零漂值。

2.4 算法應用方案確定

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

圖5 算法應用方案Fig.5 Scheme of algorithm application

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

2.5 貝葉斯算法結果分析

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

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

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

2.6 干擾作用下的結果分析

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

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

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

2.7 線性掉格下的結果分析

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

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

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

2.8 經典算法結果對比

通過線性擬合對同樣2 d的連續重力觀測值進行零漂估計，擬合結果如圖9所示。

圖9 連續2 d重力觀測值線性擬合結果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)發現，貝葉斯算法可將重力觀測值中的非線性零漂去除，實現對非線性零漂的精準估計。

3 結 語

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

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