超導重力觀測噪聲水平的極大似然估計*

韋進 李輝 劉子維 康開軒

1)武漢大學測繪學院，武漢4300792)中國地震局地震研究所，武漢4300713)地殼運動與地球觀測實驗室，武漢430071

超導重力觀測噪聲水平的極大似然估計*

韋進1，2，3)李輝2，3)劉子維2，3)康開軒2，3)

1)武漢大學測繪學院，武漢430079

2)中國地震局地震研究所，武漢430071
3)地殼運動與地球觀測實驗室，武漢430071

利用極大似然估計和3種互不相關的噪聲(白噪聲，閃爍噪聲，游走噪聲)模型及其組合，對SG-053超導重力儀產(chǎn)出的重力固體潮整時值殘差時間序列進行分析，得出殘差時間序列中明顯存在彩色噪聲(閃爍噪聲、游走噪聲):白噪聲1.40×10－8ms－2，閃爍噪聲1.85×10－8ms－2，游走噪聲為2.40×10－8ms－2。該方法估計的白噪聲水平與潮汐分析方法得出的結(jié)果一致。

超導重力儀;殘差時間序列;極大似然估計;噪聲水平;白噪聲

AbstractThe instrument noise level is one of the most important factor for the evaluation of SG-053.The tidal gravity residual time series of the SG-053 has been analyised with maximum likelihood estimation(MLE)and three noise(white noise，flicker noise，walk noise)models and their combination.It is shown that there is color noise in the tidal gravity residual，the white noise level is 1.40×10－8ms－2，the flicker noise 1.85×10－8ms－2and the random walk noise 2.40×10－8ms－2.This analysis result is consistent with the tidal analysis.

Key words:superconducting gravimeter;residual time series;maximum likelihood estimation;noise level;white noise

1 引言

超導重力儀的噪聲水平的量級是儀器性能的重要指標，它是選擇臺址［1］，影響潮汐分析的重要因素［2，3］。文獻［3，4］1997年利用武漢基準站超導重力儀確定了一種所謂的“單邊”噪聲，據(jù)經(jīng)驗估計量級在(0.2～0.8)×10－8ms－2。目前，分析超導重力儀噪聲水平的主要方法是潮汐分析［1］。用兩套潮汐分析軟件(VAV與ETERNA)可以估計儀器的白噪聲。文獻［2］在利用超導重力儀精密測定地球潮汐常數(shù)時，認為武漢臺站背景噪聲較大，約為1.5× 10－8ms－2。而在很多地球物理現(xiàn)象中，時間序列噪聲可以被描述成為一個冪率過程［5］。在連續(xù)GPS觀測中采用冪率噪聲模型來估計測站位移速率，能夠很好地減少由于冪率噪聲影響的低頻振幅的影響［7－11］。實際上利用超導重力儀觀測數(shù)據(jù)進行分析時，同樣需要進行復雜的數(shù)據(jù)處理(濾波、潮汐分析、潮汐改正、氣壓改正等)。而這樣的過程可能會引入和連續(xù)GPS分析類似的噪聲。因此，如果對殘差時間序列依然采用白噪聲模型進行估計，也可能引入冪率(彩色)噪聲，進而可能影響儀器性能的判定和臺址選擇。

本文利用潮汐模型、氣壓觀測時間序列改正后的重力固體潮后的殘差時間序列研究超導重力儀白噪聲和冪率(彩色)噪聲之間的規(guī)律，并利用極大似然估計方法估計各噪聲分量的噪聲水平。

2 噪聲數(shù)據(jù)的極大似然估計

2.1 極大似然估計

MLE(極大似然估計)［10］可以用于計算時間序列中白噪聲和冪率(彩色)噪聲的量級。在考慮冪率過程時，MLE可根據(jù)不同冪率(彩色)噪聲模型來估計對應模型的噪聲量級。用MLE估計噪聲項時，平差數(shù)據(jù)的協(xié)因數(shù)陣具最大的概率密度。概率密度方程為:

式中，like是極大似然值，det是矩陣的行列式。兩邊取自然對數(shù):

這里ln是自然對數(shù)，N是采樣數(shù)，C是數(shù)據(jù)的協(xié)因數(shù)，v是利用C協(xié)因數(shù)和加權(quán)最小二乘對原始數(shù)據(jù)估計出來的殘差。

利用極大似然估計的方法對扣除實際模型的理論固體潮值和氣壓影響后的殘差進行白噪聲分析的模型如下:

式中，x(ti)是經(jīng)過潮汐和氣壓改正后的殘差，x0是重力的常數(shù)項，r是重力儀的漂移?？紤]到殘差中也存在周期信號，則可寫為:

當ε(ti)=aαj(t)+bβj(t)時，利用極大似然估計，殘差的協(xié)因數(shù)為Cx=a2I+b2Jκ。

2.2 模型定義

白噪聲模型認為觀測數(shù)據(jù)之間沒有相關性。而其他冪率(彩色)噪聲模型通過不同的表達方式反映了觀測數(shù)據(jù)之間存在相互關系［11］。各模型協(xié)因數(shù)陣C、噪聲水平a、b和中誤差σ如下:

1)白噪聲模型:b=0，Cx=a2I，C－1x=(1/a2)I

2)閃爍噪聲模型:a=0，Cx=b2Jκ，C－1x=(1/

這里，N為采樣數(shù)，T為時間跨度［13］。

4)白噪聲加閃爍噪聲:a≠0，b≠0，Cx=a2I+ b2Jκ，κ=1。

5)白噪聲、游走噪聲和閃爍噪聲:a≠0，b≠0，Cx=a2I+b2J1+c2J2。

復合噪聲模型用Brent［14］算法實現(xiàn)參數(shù)估計。

3 數(shù)據(jù)處理計算與結(jié)果分析

3.1 殘差時間序列

用國際固體潮中心［15］提供的VAV(V03.11)和ETERNA34潮汐分析軟件對SG-053超導重力儀2009-03-01—09-01預處理后的整時值和氣壓觀測數(shù)據(jù)聯(lián)合進行潮汐分析。除了分析各潮波的潮汐因子和相位外，還分析了殘差的振幅中誤差和氣壓導納值(表1)。

表1 VAV和ETERNA的潮汐分析白噪聲分析結(jié)果Tab.1Tidal analysis results of the VAV and ETERNA

VAV與ETERNA軟件均采用最小二乘方法估計潮汐參數(shù)。但在漂移模型、氣壓回歸因子和擾動數(shù)據(jù)剔除等方面各自采用了不同的方法消除系統(tǒng)誤差和粗差對參數(shù)估計的影響［1］。從上述分析結(jié)果來看兩套軟件計算的氣壓導納值之間和白噪聲水平之間趨于同一量級。氣壓導納值在(－0.322 29～－0.302 361)×10－8ms－2/hPa之間，白噪聲水平在(1.1～1.4)×10－8ms－2水平范圍內(nèi)。而且從圖1可以看到，利用兩個模型的殘差序列的差不超過0.004×10－8ms－2的量級，表明利用兩個潮汐模型和氣壓導納值模型改正后的序列幾乎是相同的。

對重力固體潮進行潮汐改正和氣壓改正后的殘差時間序列以及它們之間的差值見圖1。

3.2 不同模型的極大似然估計

和潮汐分析方法不同，極大似然估計在已知估計參數(shù)的數(shù)學期望的情況下，能夠無偏地估計出系統(tǒng)誤差中和時間相關的誤差影響［6－8］。然而，儀器在觀測過程中不可避免地受到外界因素的干擾，出現(xiàn)直接影響估計結(jié)果的粗差。因此，在已知時間序列變化速率、各波群振幅相位的數(shù)學期望值的前提下，利用極大似然估計選擇最優(yōu)的估計參數(shù)，從而無偏地估計出不同噪聲模型的噪聲水平。

圖1 SG-053超導重力儀2009-03-01—09-01殘差時間序列(a)和兩個模型殘差的差結(jié)果(b)Fig.1Residual time serial of SG-053 between 2009-03-01 and 2009-09-01(a)and difference between the results of two model residuals(b)

3.2.1 模型分波頻率的關系

利用殘差時間序列和5種噪聲模型進行極大似然估計。計算出從起始的2 000小時觀測數(shù)據(jù)不同周期下各模型的噪聲水平、極大似然值之間的關系(表2)。

進行極大似然估計，在利用6個周期時，各模型均出現(xiàn)了極大值拐點現(xiàn)象;超過7個周期后，模型全部出現(xiàn)系統(tǒng)粗差現(xiàn)象;但采用不同的周期對極大值及其噪聲分量的估計影響不大。同時，組合模型的估計明顯優(yōu)于單一模型的估計(極大值普遍較大)，組合模型估計出的各噪聲分量在數(shù)量級上基本一致(白噪聲水平在(1.6～1.8)×10－8ms－2，閃爍噪聲在(2.1～2.8)×10－8ms－2)。而且利用3種模型的估計極值優(yōu)于其他各模型。分析表明，利用6周期和復合噪聲模型的估計能更加有效地估計出觀測時間序列的各噪聲水平。

3.2.2 模型采樣數(shù)和采樣時段的關系

根據(jù)殘差時間序列，利用3種噪聲的組合模型進行極大似然估計。計算周期數(shù)為6的情況下，兩個組合模型、采樣數(shù)和采樣時段的關系如圖2所示。

1)模型和采樣數(shù)的關系

圖2(a，b)是利用兩種組合噪聲模型，以殘差時間序列的起點作為起點，估計出的各模型的噪聲分量和觀測數(shù)據(jù)長度之間的關系。在采樣數(shù)為1 500～3 000的各噪聲水平均趨于平穩(wěn)(變化不超過0.5 ×10－8ms－2)。而當采樣數(shù)在800～1 000時范圍內(nèi)時，3類噪聲水平普遍不穩(wěn)定，其主要原因是采樣數(shù)不夠。觀測采樣數(shù)超過該范圍到3 500時，各模型的噪聲水平開始同步上揚，其主要原因是后期的觀測數(shù)據(jù)出現(xiàn)了不明原因的噪聲(圖1)。

2)模型和采樣時段的關系

圖2(c，d)是利用兩種組合噪聲模型，以100為滑動步長取不同的滑動起點，估計出的各組合模型的各噪聲分量和采樣時段之間的關系。而每個點是不同觀測數(shù)據(jù)長度下極大估值的個噪聲分量的均值和方差。該時間序列表明:兩個組合模型的各噪聲分量都表現(xiàn)出先減小后增大的過程，各噪聲分量的變化也一致。在1 300小時作為起點估計的各噪聲分量達到了極小值。

3)模型、采樣數(shù)和采樣時段之間的關系

從模型、采樣數(shù)和采樣時段之間的關系表明，無論是WN+FN+RWN還是WN+FN模型，利用極大似然估計后，它們都表現(xiàn)出了在1 300小時為起始點，1 000～1 500小時窗口達到各噪聲分量的極小值。但是兩個組合模型估計出的圖3描述了兩個組合模型的白噪聲分量、采樣數(shù)量和采樣時段之間的關系。

3.2.3 組合模型的極大似然估計

分析表明，利用組合模型和最優(yōu)估計參數(shù)(分波周期、組合模型、采樣數(shù)量和采樣時段)采用Brent algorithm算法［12］實現(xiàn)了兩種組合模型的各噪聲水平分量的估計(圖4)。分析結(jié)果表明在橫坐標為53°，縱坐標為46°時達到極大值:白噪聲水平為1.40×10－8ms－2，閃爍噪聲水平為1.85×10－8ms－2，游走噪聲為2.40×10－8ms－2。該分析結(jié)果中白噪聲的分量與潮汐分析的結(jié)果一致。

4 結(jié)論與討論

1)利用SG-053預處理后的4 422個采樣數(shù)據(jù)進行潮汐分析和進行了潮汐改正和氣壓改正后的殘差序列進行極大似然估計結(jié)果表明，兩個軟件分析的殘差序列相差不超過0.004×10－8ms－2;白噪聲量級趨于一致，在(1.1～1.4)×10－8ms－2范圍內(nèi)。

2)在最優(yōu)的分析參數(shù)下，利用Brent算法計算兩組合模型的各噪聲分量，結(jié)果都趨于一致。但是用WN+FN+RWN估計出的當前重力固體潮殘差時間序列的噪聲量級為:白噪聲水平為1.40×10－8ms－2，閃爍噪聲水平為1.85×10－8ms－2，游走噪聲為2.40×10－8ms－2。這一分析結(jié)果中白噪聲的分量和潮汐分析的分析結(jié)果趨于一致。

表2 周期、數(shù)噪聲水平和極大值之間的關系Tab.2Relations among period，noise level and maximum value

圖2 兩種組合模型采樣數(shù)和各噪聲水平關系序列Fig.2Relation between sample numbers and the noise level with the two types of combined model

圖3 組合模型采樣數(shù)和采樣時段之間的關系Fig.3Relations between sample numbers and sample period with the two types of combined model

圖4 Brent算法的極大似然估計Fig.4Brent maximum likelihood estimation algorithm

1田桂娥，等.VAV和ETERNA潮汐分析方法的比較和研究［J］.大地測量與地球動力學，2009，(2):96－99.(Tian Guie，et al.Comparison and investigation of VAV and ETERNA tidal analysis methods［J］.Journal of Geodesy and Geodynamics，2009，(2):96－99)

2孫和平，等.用超導重力儀觀測數(shù)據(jù)精度測定地球潮汐常數(shù)［J］.地殼形變與地震，1997，(4):17－25.(Sun Heping，et al.The constant precision determination of the earth tides with superconducting gravimeter datums［J］.Crustal Deformation and Earthquake，1997，(4):17－25)

3徐建橋.重力固體潮汐理論及分析方法——武昌臺超導重力儀觀測資料的分析處理［D］.中國科學院測量與地球物理研究所，1997.(Xu Jianqiao.The theory and analysis of gravity tidal——superconducting gravimeter data analysis and processing in Wuchang station［D］.Institute of Geodesy and Geophysics Chinese Academy of Sciences，1997)

4Kroner C，et al.Long-term stability and tidal parameters of the SCG-record at Wuhan［A］.Proceedings 12th International Symposium on Earth Tides［C］.1995:265－276.

5Mandelbrot B and Van Ness J.Fractional brownian motions，fractional noises，and applications［R］.SIAM Rev.1968，10，422－439.

6Simon D and Williams P.Error analysis of continuous GPS position time series［J］.Journal of Geophysical Research，2004，109(B03412):19.

7Simon D and Williams P.The effect of coloured noise on the uncertainties of rates estimated from geodetic time series［J］.Journal of Geodesy，2003，76:483－494.

8John Langbein.Correlated errors in geodetic time series:Implications for time－dependent deformation［J］.Journal of Geophysical Research，1997，102(B1):591－603.

9Ailin Mao，et al.Noise in GPS corrdinate time series［J］.Journal of Geophysical Research，1999，104(B2):2 797－2 816.

10Simon D and Williams P.CATS:GPS coordinate time series analysis software［J］.GPS solut，2008，12:147－153.

11Jie Zhang，et al.Southern California permanent GPS geodetic array:error analysis of daily position estimates and site velocities［J］.Journal of Geophysical Research，2007，102: 18 035－18 055.

12Teferle F N，et al.Crustal motions in Great Britain:evidence from continuous GPS，absolute gravity and Holocene sea level data［J］.Geophy.J.Int.，2009，178:23－46.

13Johnson H O and Wyatt FK.Geodetic network design for fault-mechanics studies［J］.Manuscr.Geod.，1994，(10): 309－323.

14Press WH，et al.Numerical recipes［M］.New York:Cambridge University Press，2007.

15Jean-Pierre.ETERNA［EB/OL］.http://www.observatoire.be/ICET/soft/index.html， 2010-07-01/2010-12-15.

NOISE LEVEL MEASUREMENT OF SG-053 WITH MLE

Wei Jin1，2，3)，Li Hui2，3)，Liu Ziwei2，3)and Kang Kaixuan2，3)
1)School of Geodesy and Geometics，Wuhan University，Wuhan430079
2)Institute of Seismology，CEA，Wuhan430071
3)Crustal Movement Laboratory，Wuhan430071

P207;P203

A

1671-5942(2011)03-0069-06

2010-12-25

國家自然科學基金(41004030);中國地震局地震研究所所長基金(IS200951041)

韋進，男，1981年生，博士研究生，助理研究員，主要從事重力臺網(wǎng)管理和重力固體潮分析研究.E－mail:pierce212@163.com