韋進(jìn) 李輝 劉子維 康開軒

1)武漢大學(xué)測(cè)繪學(xué)院，武漢4300792)中國(guó)地震局地震研究所，武漢4300713)地殼運(yùn)動(dòng)與地球觀測(cè)實(shí)驗(yàn)室，武漢430071

超導(dǎo)重力觀測(cè)噪聲水平的極大似然估計(jì)*

韋進(jìn)1，2，3)李輝2，3)劉子維2，3)康開軒2，3)

1)武漢大學(xué)測(cè)繪學(xué)院，武漢430079

2)中國(guó)地震局地震研究所，武漢430071
3)地殼運(yùn)動(dòng)與地球觀測(cè)實(shí)驗(yàn)室，武漢430071

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

超導(dǎo)重力儀;殘差時(shí)間序列;極大似然估計(jì);噪聲水平;白噪聲

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 引言

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

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

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

2.1 極大似然估計(jì)

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

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

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

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

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

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

2.2 模型定義

白噪聲模型認(rèn)為觀測(cè)數(shù)據(jù)之間沒有相關(guān)性。而其他冪率(彩色)噪聲模型通過不同的表達(dá)方式反映了觀測(cè)數(shù)據(jù)之間存在相互關(guān)系［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為時(shí)間跨度［13］。

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

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

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

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

3.1 殘差時(shí)間序列

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

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

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

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

3.2 不同模型的極大似然估計(jì)

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

圖1 SG-053超導(dǎo)重力儀2009-03-01—09-01殘差時(shí)間序列(a)和兩個(gè)模型殘差的差結(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 模型分波頻率的關(guān)系

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

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

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

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

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

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

2)模型和采樣時(shí)段的關(guān)系

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

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

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

3.2.3 組合模型的極大似然估計(jì)

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

4 結(jié)論與討論

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

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

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

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

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

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

1田桂娥，等.VAV和ETERNA潮汐分析方法的比較和研究［J］.大地測(cè)量與地球動(dòng)力學(xué)，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孫和平，等.用超導(dǎo)重力儀觀測(cè)數(shù)據(jù)精度測(cè)定地球潮汐常數(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徐建橋.重力固體潮汐理論及分析方法——武昌臺(tái)超導(dǎo)重力儀觀測(cè)資料的分析處理［D］.中國(guó)科學(xué)院測(cè)量與地球物理研究所，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

國(guó)家自然科學(xué)基金(41004030);中國(guó)地震局地震研究所所長(zhǎng)基金(IS200951041)

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