陳善勇 楊 江,2,3

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

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

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

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

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

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

1 模型構(gòu)建

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

1.1 矩陣定義

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

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

將每個(gè)樣本點(diǎn)對(duì)應(yīng)的零漂表示為：

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

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

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

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

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

1.2 基本公式

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

x-v-g=δ

(1)

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

δ～N(0,Σ1)

(2)

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

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

x～N(v+g,Σ1)

(3)

1.3 先驗(yàn)分布確定

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

(4)

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

Av=ξ

(5)

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

ξ～N(0，Σ2)

(6)

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

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

Av～N(0，Σ2)

(7)

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

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

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

(8)

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

(9)

1.4 似然函數(shù)確定

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

(10)

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

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

(11)

1.6 后驗(yàn)分布計(jì)算

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

(12)

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

(13)

1.7 超參數(shù)確定

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

對(duì)于零漂后驗(yàn)分布均值和方差的確定，方法如下：

(14)

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

(15)

第n次后驗(yàn)分布的方差為：

(16)

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

后驗(yàn)分布的均值經(jīng)過(guò)化簡(jiǎn)可表示為：

(17)

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

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

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

將九峰臺(tái)站DZW重力儀在2019年的連續(xù)觀測(cè)數(shù)據(jù)作為算法的驗(yàn)證數(shù)據(jù)源。該數(shù)據(jù)已經(jīng)過(guò)精度為10-5μGal的固體潮校正，可不考慮溫度、氣壓等因素的影響。通過(guò)線性擬合，得到圖1中表示線性零漂的直線，與潮汐校正后的數(shù)據(jù)點(diǎn)進(jìn)行對(duì)比。擬合直線斜率即為估計(jì)的零漂率，為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的對(duì)角矩陣。由于零漂的先驗(yàn)分布的協(xié)方差矩陣Σ2為未知參數(shù)，但在多層運(yùn)算后將會(huì)趨近于0，因此不妨設(shè)其為與Σ1同數(shù)量級(jí)的矩陣，令其為0.01的對(duì)角矩陣。

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

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

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

圖2 后驗(yàn)分布方差變化Fig.2 Variance variation of posteriori distribution

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

圖3 樣本點(diǎn)均值變化Fig.3 Meanvalues change of sample points

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

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

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

圖4 觀測(cè)數(shù)據(jù)在去除不同計(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，后者方差明顯小于前者，說(shuō)明隨著計(jì)算次數(shù)的增加，均值不斷逼近，使去零漂結(jié)果為常數(shù)數(shù)值，均值的收斂極限即為真實(shí)的零漂值。

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

第1天的1 440個(gè)樣本點(diǎn)采用貝葉斯模型進(jìn)行預(yù)處理，通過(guò)對(duì)其進(jìn)行1 000次重復(fù)計(jì)算，得到精準(zhǔn)的1 440個(gè)零漂值。對(duì)于新樣本數(shù)據(jù)，通過(guò)左移操作實(shí)現(xiàn)樣本數(shù)組的更新，為確保樣本點(diǎn)與零漂值一一對(duì)應(yīng)，零漂的均值矩陣與協(xié)方差矩陣中的對(duì)角元素?cái)?shù)組同時(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代表與新樣本值對(duì)應(yīng)的均值，是需要預(yù)測(cè)的參數(shù)。Sn+1代表與新樣本值對(duì)應(yīng)的方差，等于先驗(yàn)分布的方差。該方案的優(yōu)點(diǎn)是對(duì)于新移入的樣本值，其均值由移入至最后輸出會(huì)經(jīng)歷1 440次貝葉斯算法計(jì)算，可保證結(jié)果的高精度。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 結(jié) 語(yǔ)

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

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