基于加權(quán)最小二乘譜迭積算法探測超導(dǎo)重力觀測弱信號＊

尹 暉 王艷艷

(武漢大學(xué)測繪學(xué)院,武漢 430079)

基于加權(quán)最小二乘譜迭積算法探測超導(dǎo)重力觀測弱信號＊

尹 暉 王艷艷

(武漢大學(xué)測繪學(xué)院,武漢 430079)

提出采用加權(quán)最小二乘譜迭積算法探測超導(dǎo)重力觀測數(shù)據(jù)隱藏的弱周期信號,該算法基于不等權(quán)的觀測序列和序列的自相關(guān)性和互相關(guān)性,以高分辨率提取序列潛在的共同信號,削弱或抵消噪聲信號。并以實例進(jìn)行了分析和驗證。

加權(quán)最小二乘譜分析;譜迭積;超導(dǎo)重力數(shù)據(jù);弱信號;噪聲

1 引言

超導(dǎo)重力儀具有穩(wěn)定性強、靈敏度高、噪音低、漂移小等特點,可以提供精度為 10-11ms-2的重力觀測數(shù)據(jù),為地球深內(nèi)部構(gòu)造和核模弱信號檢測等研究提供了有效途徑。但由于超導(dǎo)重力儀觀測數(shù)據(jù)受大氣環(huán)境、地震和儀器干擾如斷電,維修、調(diào)試、系統(tǒng)升級、氦注入等多種因素的影響,觀測數(shù)據(jù)會出現(xiàn)不連續(xù)和不平穩(wěn)等現(xiàn)象,對存在斷鏈、尖峰、基準(zhǔn)偏移、不等權(quán)、不等間隔等情況下的重力觀測數(shù)據(jù),采用最小二乘譜分析進(jìn)行處理已顯示出較大的優(yōu)越性[1]。

但由于核模等信號很弱,對掩蓋在噪聲背景下的弱信號僅采用譜分析仍很難檢測到。本文提出采用基于加權(quán)最小二乘譜迭積算法來探測超導(dǎo)重力觀測數(shù)據(jù)隱藏的弱周期信號。該算法能考慮觀測序列和序列的自相關(guān)性和互相關(guān)性,以高分辨率提取序列潛在的共同信號,削弱或抵消噪聲信號。并以加拿大超導(dǎo)重力觀測站實測數(shù)據(jù)進(jìn)行了分析驗證。

2 加權(quán)最小二乘譜迭積算法

2.1 加權(quán)最小二乘譜分析

給定一組觀測值向量 f=f(t){fi}(i=1,2,…, n),該觀測值向量可能存在斷鏈、尖峰、基準(zhǔn)偏移、不等間隔、不等權(quán)等情況,設(shè)觀測值向量具有協(xié)方差陣 Cf,建立如下模型：

根據(jù)最小二乘譜分析原理[2],得加權(quán)最小二乘譜的計算公式為

若設(shè)定基函數(shù)形式為一組頻率為ωi,i=1,2,…,m的三角函數(shù),則對于不同的頻率,可以得到不同的頻譜值:

2.2 譜迭積算法

諧波譜積算法是將交叉譜的思想推廣應(yīng)用到多個不同時間序列的聯(lián)合譜分析上,以反映多個相同時段的不同時間序列在頻域變化上的相互關(guān)系,將不同時間序列在相同頻率段對應(yīng)的頻譜值相乘,能夠?qū)Σ煌瑫r間序列計算得到的同一頻率處的信號起放大作用,從而將掩藏在噪音中的信號迭積檢測出來。

Smylie D E[3]1993年應(yīng)用快速傅立葉變換和譜迭積方法檢測歐洲 3個臺站的超導(dǎo)重力觀測數(shù)據(jù)中的弱重力共振信號。類似地,Vanicek P[4,5]提出將最小二乘譜分析應(yīng)用到多個時間序列的聯(lián)合譜分析上,進(jìn)行譜迭積,由公式 (4),給出在指定頻率區(qū)間[ω1,ω2]內(nèi),譜迭積的計算公式為:

式 (5)是N個時間序列的譜積,si(ω)是第 i個時間序列的最小二乘譜,從式 (5)可以看出,對于序列中潛在的共同信號,經(jīng)過譜迭積計算后,譜值以高分辨率得到增強,而對于序列中的噪聲信號,由于其出現(xiàn)的隨機(jī)性和不確定性,其迭積后的譜值將得以削弱或抵消。

3 實例分析

3.1 數(shù)據(jù)處理與分解

加拿大超導(dǎo)重力觀測站實測數(shù)據(jù)為 1 s的采樣率,對實測數(shù)據(jù)采取如下的處理步驟[4]：1)數(shù)據(jù)獲取。對以天為 1個記錄文件、采樣率為 1 s的原始數(shù)據(jù)(單位為V)進(jìn)行轉(zhuǎn)化 (乘以 -78.48×10-5m/ s2V轉(zhuǎn)化為重力單位 10-5ms-2),刪除重疊數(shù)據(jù),使數(shù)據(jù)序列的時間單位呈單調(diào)遞增,合并成以年為單位的數(shù)據(jù)文件;2)重力潮汐影響剔除。根據(jù)重力潮汐預(yù)測程序 G-wave[6]計算得到重力潮汐理論結(jié)果來剔除超導(dǎo)重力原始觀測數(shù)據(jù)中的重力潮汐信號的影響;3)數(shù)據(jù)過濾。采用滯后窗為 50 s的 Parzen窗作為權(quán)函數(shù),最大步長為 100 s進(jìn)行隨機(jī)采樣,為避免數(shù)據(jù)間的相關(guān)性,要求兩個連續(xù)窗不能疊置,這樣每小時大約得到 24個不等間隔的過濾數(shù)據(jù),由窗函數(shù)按誤差傳播定律計算得到其方差,因而得到不等權(quán),不等間隔的隨機(jī)觀測序列;4)數(shù)據(jù)分解。針對本項目研究感興趣的弱信號位于周期 2～10小時之間,按照 Nyquist采樣定理,我們將全年數(shù)據(jù)分解為以月為單位時段(12個月),其目的是利用時間序列的隨機(jī)性、非平穩(wěn)性和自相關(guān)性,運用譜迭積多尺度地檢測出隱藏在所有時間序列中的共同信號。

3.2 計算與分析

對 1998年 12個月不等間隔和不等權(quán)值的超導(dǎo)重力觀測數(shù)據(jù)序列,運用加權(quán)最小二乘譜分析模型進(jìn)行計算,得到最小二乘譜圖,圖 1、2給出計算的部分譜圖,從圖中可見,由于超導(dǎo)重力觀測受噪聲信號的影響,信號具有明顯的隨機(jī)性和不確定性,要分辨和檢測到超導(dǎo)重力觀測數(shù)據(jù)中的弱信號很難,不同月份譜分析計算的結(jié)果沒有明顯的對應(yīng)比較關(guān)系,一致性較差,根本無法判斷檢測信號的結(jié)果。

對此,在譜分析計算的基礎(chǔ)上,運用譜迭積算法對超導(dǎo)重力觀測數(shù)據(jù)的最小二乘譜分析結(jié)果進(jìn)行譜迭積,圖 3給出了 1998年超導(dǎo)重力觀測數(shù)據(jù)的最小二乘譜迭積后的結(jié)果,從圖中可見,在經(jīng)過譜迭積計算后,譜值的分辨率得到大大增強,而序列中的噪聲信號,由于其出現(xiàn)的隨機(jī)性和不確定性,其迭積后的譜峰得到削弱或抵消。應(yīng)用最小二乘譜迭積分析方法分析加拿大 6年 (1997—2002)實測超導(dǎo)重力觀測數(shù)據(jù)的弱信號,進(jìn)行了地球內(nèi)核 Slichter模譜峰信號的檢測[1],表 1給出了 90%置信水平下的檢測結(jié)果,其結(jié)果與 Smylie的結(jié)果[7]較為接近。

圖 2 1998年 2月的最小二乘譜圖Fig.2 Least square spectrum in Feb.1998

圖 3 1998年加權(quán)最小二乘譜迭積圖Fig.3 Product ofweighted least square spectrum in 1998

表 1 Slichter模譜峰檢測結(jié)果(單位:小時)Tab.1 Test results for spectrum peak of Slichter mode(un it:hour)

4 結(jié)語

1)基于最小二乘譜分析技術(shù),無須內(nèi)插缺失的數(shù)據(jù),無須平滑尖鋒和漂移,直接由加權(quán)最小二乘譜分析得到譜值,從而避免了不必要的數(shù)據(jù)處理,而使隱藏在序列中的有效信息受到損失。

2)對超大量數(shù)據(jù)進(jìn)行過濾,在整個數(shù)據(jù)段上對所有數(shù)據(jù)進(jìn)行融和,并以加權(quán)處理提取時間序列的有效信息,是數(shù)據(jù)融合和信息挖掘的一種有效方法,過濾數(shù)據(jù)為不等間隔,其方差由窗函數(shù)按誤差傳播定律得到。

3)基于加權(quán)最小二乘譜迭積算法探測超導(dǎo)重力觀測數(shù)據(jù)隱藏的弱周期信號,可以顧及觀測序列和序列的自相關(guān)性和互相關(guān)性,從而削弱或抵消了超導(dǎo)重力觀測數(shù)據(jù)中的噪聲信號,增強譜值的分辨率。

4)Slichter模信號檢測是目前地學(xué)界研究的熱點之一,本文利用超導(dǎo)重力觀測數(shù)據(jù),運用最小二乘譜迭積方法檢測 Slichter模信號的研究結(jié)果表明,該方法是地球深內(nèi)部弱信號檢測的有效方法,為進(jìn)一步研究地球深內(nèi)部結(jié)構(gòu)及其動力學(xué)現(xiàn)象提供了新的途徑。

1 SpirosD,et al.Least-squares self-coherency analysis of superconducting gravi meter records in search for the Slichter triplet[J].Phys Earth Planet Int.,2007,160：108-123.

2 Pagiatakis SD.Stochastic significance of peaks in the leastsquares spectrum[J].Journal of Geodesy,1999 73：67-78.

3 Smylie D E,et al.The product spectra of gravity and barometric pressure in Europe[J].Physics of the Earth and Planetary Interiors,1993,80：135-157.

4 Vanicek P.Approximate spectral analysis by least-squares fit [J].Astrophysics and Space Science,1969,4：387-391.

5 Vanicek P.Further development and properties of the spectral analysis by least-squares[J].Astrophysics and Space Science,1971,12：10-33.

6 Yin Hui and SpirosD Pagiatakis.Least squares spectral analysis and its application to superconducting gravimeter data analysis[J].Geo-Spatial Information Science,2004,7(4)：279-283.

7 Smylie D E,et al.The Inner core translational triplet and the desity near earth’s center[J].Science,1992,255：1 678 -1 682.

DETECTI ON OF W EAK SUPERCONDUCTING GRAVI M ETRIC SIGNALSW ITH PRODUCT SPECTRUM OF W EIGHTED LSSA

Yin Hui andWang Yanyan
(School of Geodesy and Geom atics,W uhan University,W uhan 430079)

Least Squares SpectralAnalysis(LSSA)has been shown more suitable than Fourier analysis for analyzing time serieswith spikes,gaps,datum shrifts(offsets)and other disturbances in a variety of applications in geodesy and geophysics.To detect the weak superconducting gravimetric signals,a methodology of product spectrum based on weighted LSSA is presented,which can take consideration the correlations among the signals and highlight the potential spectra power of high resolution in the common figures as well as weaken or eliminate the power of noise.A data set of superconducting gravi meter observations at the Canada Superconducting Gravi meter Installation(CSGI)has been assembled for analysis to show its efficacy.

weighted LSSA;product spectrum;superconducting gravimetric data;weak signals;noise

1671-5942(2011)04-0063-03

2011-02-19

國家自然科學(xué)基金(51077105);國家電網(wǎng)公司科技項目(SGKJ2009-0609-51)

尹暉,女,1962年生,教授,博士,博導(dǎo),主要從事變形監(jiān)測分析與預(yù)報、超導(dǎo)重力數(shù)據(jù)分析及地球動力學(xué)研究.E-mail：hyin @sgg.whu.edu.cn

P207

A