加權(quán)遞推最小二乘在絕對(duì)重力儀數(shù)據(jù)處理中的應(yīng)用

羅棋，龐聰，李查瑋

（中國(guó)地震局地震研究所（地震大地測(cè)量重點(diǎn)實(shí)驗(yàn)室），湖北武漢430071）

加權(quán)遞推最小二乘在絕對(duì)重力儀數(shù)據(jù)處理中的應(yīng)用

羅棋，龐聰，李查瑋

（中國(guó)地震局地震研究所（地震大地測(cè)量重點(diǎn)實(shí)驗(yàn)室），湖北武漢430071）

針對(duì)絕對(duì)重力儀高精度的需求，在抗差參數(shù)估計(jì)的基礎(chǔ)上提出了加權(quán)遞推算法及其改進(jìn)算法，加權(quán)遞推和改進(jìn)算法改良抗差權(quán)重因子并且避免求復(fù)雜多維的逆矩陣。通過(guò)仿真表明，該算法對(duì)含有白噪聲或者異常數(shù)據(jù)都起到很好的仰制作用并提高了重力測(cè)量精度。

絕對(duì)重力儀；抗差估計(jì)；加權(quán)遞推；最小二乘

絕對(duì)重力測(cè)量在計(jì)量學(xué)，地球物理和大地測(cè)量有著重要的作用，其精度的提高對(duì)地震活動(dòng)、地質(zhì)勘測(cè)都有不可估量的作用[1]。一般測(cè)量絕對(duì)重力加速度g的方法是基于落體每下落半個(gè)波長(zhǎng)距離產(chǎn)生一個(gè)干涉條紋的規(guī)律，獲得物體下落的位移信息，對(duì)比銣原子鐘給出的時(shí)間，得到時(shí)間距離對(duì)，擬合得到重力加速度值[2]。

目前普遍采用最小二乘和非線性擬合時(shí)間距離對(duì)[3]，與非線性相比，最小二乘原理簡(jiǎn)單，易于理解和掌握，且在編程中易于實(shí)現(xiàn)，但是最小二乘算法適用于近似服從正態(tài)分布的測(cè)量數(shù)據(jù)，參數(shù)擬合結(jié)果對(duì)非正態(tài)分布的觀測(cè)粗差極為敏感[4]，沒(méi)有考慮觀測(cè)誤差的不同對(duì)算法解的影響[5]。用雙采樣過(guò)零檢測(cè)，時(shí)刻數(shù)據(jù)將越來(lái)越密集，時(shí)間間隔越來(lái)越短，測(cè)量的誤差影響將會(huì)放大，因此需要對(duì)最小二乘進(jìn)行改進(jìn)。文獻(xiàn)[4]提出了抗差估計(jì)，但是涉及到求解復(fù)雜多維的逆矩陣。基于上述問(wèn)題，為了進(jìn)一步提供算法的可行性和精度要求，提出加權(quán)遞推形式，避免了復(fù)雜逆矩陣的求解。文中基于文獻(xiàn)[4]中的迭代算法，采用文獻(xiàn)[5]的加權(quán)陣，對(duì)擬合的時(shí)間距離對(duì)進(jìn)行加權(quán)處理，擬合誤差方差越大的賦予的權(quán)重越小，誤差方差小的賦予的權(quán)重大，這也充分體現(xiàn)了加權(quán)遞推最小二乘算法的優(yōu)越性[5]。

1 加權(quán)遞推最小二乘算法

1.1 加權(quán)遞推最小二乘

絕對(duì)重力儀研究的對(duì)象是做自由落體運(yùn)動(dòng)的落體，針對(duì)落體運(yùn)動(dòng)，絕對(duì)重力儀數(shù)據(jù)符合二次多項(xiàng)式曲線（包含重力梯度和光速有限性請(qǐng)參考文獻(xiàn)[6]中的算法，在擬合公式中增加正弦項(xiàng)和余弦項(xiàng)以仰制激光器調(diào)制頻率對(duì)輸出的影響參考文獻(xiàn)[4]，將其整理后依舊可以帶入最小二乘擬合參數(shù)，并不影響算法的通用性），通常假設(shè)時(shí)間距離對(duì)數(shù)據(jù)為(ti,si)(i=1,2,…m)，最小二乘的基本思想即選擇一個(gè)n次多項(xiàng)式s(t)用一種使均方誤差極小的方式來(lái)逼近(ti,si)，其中n遠(yuǎn)小于m。

寫成向量模型，模型偏差為：

其中ε為偏差，。

有m個(gè)等式，將其合成矩陣向量形式：

選取適當(dāng)?shù)膶?duì)角加權(quán)矩陣W，得到加權(quán)非遞推最小二乘解：

初始化：P(0)=aIm，a為很大的實(shí)數(shù)，Im為m×m單位矩陣，（0）=0 。

圖1給出了對(duì)應(yīng)的加權(quán)遞推算法計(jì)算框圖，絕對(duì)重力儀數(shù)據(jù)處理系統(tǒng)每讀取一個(gè)時(shí)間距離對(duì)(ti,si)，根據(jù)式（2）算出偏差ε(ti)，取權(quán)值帶入加權(quán)遞推最小二乘中，通過(guò)迭代，擬合出高精度的參數(shù)估計(jì)。對(duì)式（4）取加權(quán)矩陣W(m)=Im，則上述遞推公式則為遞推最小二乘。

圖1 加權(quán)遞推算法計(jì)算框圖

1.2 加權(quán)遞推算法改進(jìn)

在上述遞推公式中，每次讀取到新的時(shí)間距離對(duì)數(shù)據(jù)，立刻參與計(jì)算得到參數(shù)，更新的參數(shù)取代舊參數(shù)參與下一步的擬合，即可以立刻判斷參數(shù)的優(yōu)異程度，但是數(shù)據(jù)精度與初始化參數(shù)有關(guān)，并且需要大量擬合數(shù)據(jù)來(lái)消去初始化參數(shù)帶來(lái)的影響。實(shí)際運(yùn)用中，先存取所有的時(shí)間距離對(duì)，再取出合適數(shù)據(jù)來(lái)參與計(jì)算，很明顯這種分段式處理數(shù)據(jù)，其精度將得到提高。文中將結(jié)合這種方式，對(duì)加權(quán)遞推方式做出改進(jìn)。

將公式（5）做變形處理

上式中s0,v0,g0分別為初始位移、初始速度、當(dāng)?shù)刂亓铀俣?，在?shí)際中只需要求出重力加速度參數(shù)g值，因此對(duì)式（8）做變形寫成矩陣形式：

式中δi=w(i)為權(quán)值 ，當(dāng)δi=1 時(shí)，上式（9）（10）為改進(jìn)最小二乘遞推。明顯看出改進(jìn)之后的算法避免了求逆矩陣，由于時(shí)間距離對(duì)的選取可以是非連續(xù)的，在數(shù)據(jù)選取和處理上是非常容易實(shí)現(xiàn)的。

2 仿真實(shí)驗(yàn)

為驗(yàn)證算法在絕對(duì)重力儀數(shù)據(jù)處理中的可行性，文中基于MATLAB軟件分別針對(duì)最小二乘，加權(quán)最小二乘遞推，改進(jìn)最小二乘遞推和改進(jìn)加權(quán)最小二乘遞推算法這4種情況進(jìn)行了模擬仿真。對(duì)文獻(xiàn)[7]中的掃頻正弦信號(hào)模擬干涉條紋進(jìn)行雙采樣提取時(shí)間距離對(duì)，往信號(hào)中添加18%的幅值噪聲。

圖2給出了這4種算法的g值仿真結(jié)果，其中直線代表真實(shí)值g=9.8 m/s2，曲線代表估計(jì)值，用上述4種算法及其抗差最小二乘得到的數(shù)據(jù)處理結(jié)果見(jiàn)表1。

表1 4種算法及抗差加權(quán)仿真結(jié)果

如圖2和表1可以看出，加權(quán)后的算法精度和準(zhǔn)度比未加權(quán)的高出兩個(gè)量級(jí)，同時(shí)改進(jìn)后的算法比未改進(jìn)的算法更快速的收斂于真實(shí)值，即改進(jìn)算法只需要相對(duì)少量數(shù)據(jù)便可得到高精度的估計(jì)值。在干擾嚴(yán)重時(shí)，無(wú)論是抗差還是加權(quán)遞推都有效地提高了精度，當(dāng)只加入幅值白噪聲時(shí)，抗差效果比最小二乘提高不大，但是加權(quán)遞推依舊提高了精度，這充分的說(shuō)明了本文提出的算法能夠滿足準(zhǔn)確性。

圖2 4種算法參數(shù)估計(jì)收斂圖

3 結(jié)束語(yǔ)

文中借鑒文獻(xiàn)[4]中抗差的思想和迭代流程，就算法的改進(jìn)提出了加權(quán)遞推最小二乘算法，采用平差準(zhǔn)則作為加權(quán)陣，經(jīng)過(guò)改進(jìn)算得到的估計(jì)g值不僅具有抗差參數(shù)估計(jì)的優(yōu)越性，還具有運(yùn)算簡(jiǎn)單方便，易于實(shí)現(xiàn)、移植的特征。結(jié)合計(jì)算機(jī)運(yùn)算特點(diǎn)，改進(jìn)后的算法通過(guò)仿真實(shí)驗(yàn)驗(yàn)證，具有快速收斂特性，不依賴初始值設(shè)定的特點(diǎn)，經(jīng)過(guò)多次實(shí)驗(yàn)，改進(jìn)算法非常穩(wěn)定，具有工程應(yīng)用的可行性。該算法還可以在參數(shù)細(xì)節(jié)上進(jìn)行修改，以達(dá)到實(shí)際絕對(duì)重力儀數(shù)據(jù)處理中最佳的估計(jì)值。

[1]Svitlov S.Frequency-domain analysis of absolute gravimeters[J].Metrologia，2012，49（6）:706-726.

[2]吳瓊.高精度絕對(duì)重力儀關(guān)鍵技術(shù)研究[D].北京：中國(guó)地震局地球物理研究所，2011.

[3]Svitlov S，Maslyk P，Rothleithner C，et al.Comparison of Three Digital Fringe Signal Processing Methods in a Ballistic Free-Fall Absolute Gravimeter[J].Metrologia，2010，47（6）:677-689.

[4]胡明，張為民.抗差估計(jì)在自由落體式絕對(duì)重力儀中的應(yīng)用[J].大地測(cè)量與地球動(dòng)力學(xué)，2016，9（9）:822-825.

[5]劉謝進(jìn).遞推加權(quán)最小二乘算法的研究[J].系統(tǒng)仿真學(xué)報(bào)，2007，21（14）:4248-4250.

[6]吳書清，吉望西，劉達(dá)倫，等.絕對(duì)重力儀數(shù)值計(jì)算及不確定度評(píng)定方法[J].計(jì)量學(xué)報(bào)，2009，30（3）：212—215.

[7]楊厚麗，鄒彤，郭唐永，等.復(fù)外差法處理絕對(duì)重力儀數(shù)字干涉信號(hào)[J].大地測(cè)量與地球動(dòng)力學(xué)，2015，4（2）:346-348.

[8]周姍姍，柴金廣，李丹.實(shí)時(shí)遞推的最小二乘預(yù)測(cè)跟蹤算法[J].電子設(shè)計(jì)工程，2011，19（10）：53-55.

[9]強(qiáng)明輝，張京娥.基于MATLAB的遞推最小二乘法識(shí)別與仿真[J].自動(dòng)化與儀器儀表，2008（6）:4-5，39.

[10]賈沛璋.誤差分析與數(shù)據(jù)處理[M].北京：國(guó)防工業(yè)出版社，1992.

[11]秦延，陳宗海，李衍杰.遞推最小二乘算法的補(bǔ)充證明[J].系統(tǒng)仿真學(xué)報(bào)，2004，16（10）:2159-2160.

[12]Orlob M，Braun A.Impact Estimation and Filtering ofDisturbancesin FG5 AbsoluteGravimeter Observations[J].International Journal of Geosciences，2013，4（2）:302-308.

[13]鄧自立，王欣，高媛.建模與估計(jì)[M].北京：科學(xué)出版社，2007.

[14]邵建新.最小二乘法線性擬合中參數(shù)的確定問(wèn)題[J].大學(xué)物理，2003，22（1）.

[15]卓金武.MATLAB在數(shù)學(xué)建模中的應(yīng)用[M].北京：北京航空航天大學(xué)出版社，2014.

Application of weighted recursive least square to data processing of absolute gravimeter

LUO Qi，PANG Cong，LI Cha?wei
（Institute of Seismology，CEA，Wuhan430071，China）

In view of the high precision demand of absolute gravimeter，the weighted recursive algorithm and its improved algorithm are proposed based on the robust parameter estimation，weighted recursive algorithm and improved robust weighting factor and inverse matrix for complex multidimensional.Simulation results show that the proposed algorithm has a good effect on the white noise or abnormal data and improves the accuracy of gravity measurement.

absolute gravimeter；robust estimation；weighted recursive；least square

P315

A

1674-6236（2017）22-0065-04

2016-09-21稿件編號(hào)：201609188

羅棋（1990—），男，湖北天門人，碩士。研究方向：絕對(duì)重力儀落體控制。