概率密度函數在超導重力臺站背景噪聲中的應用

張凌云,薛秀秀,田桂娥,程海港

(華北理工大學 礦業工程學院,河北 唐山 063210)

引言

超導重力儀(Superconducting Gravimeter,簡稱SG)是一種高度靈敏且穩定的重力測量儀器,具有低噪聲水平、穩定性和廣泛的頻率響應范圍。靈敏度可達1 ngal (0.01 nm/s2),年零漂小于10×10-8m/s2,可以準確監測重力隨時間的微小變化,有助于深入了解地球內部結構和各種地球動力學過程[1,2]。SG臺站的背景噪聲水平是評估觀測數據質量和臺站觀測環境的關鍵指標,對于臺站選址、儀器校準和地球動力學研究至關重要。全球SG背景噪聲水平分析為全球動力學信號的探測提供了數據選擇依據[3]。通過選擇噪聲水平較低的SG觀測數據,并應用有效的疊加方法,可以進一步降低噪聲水平,并有利于提取微弱信號[4-10],Peterson通過分析全球75個地震儀臺站的觀測記錄,得出了低噪聲模型NLNM(New Low Noise Model)。NLNM是全球地震儀觀測臺站噪聲水平的下包絡線,代表了儀器性能和臺站綜合條件最佳的地震儀臺站的噪聲水平[11,12]。背景噪聲水平可以從功率譜密度(Power Spectral Density) 中提取,其中簡正模頻段為200～600 s的平均功率譜密度被稱為地震噪聲水平(seismologicalnoisemagnitude,SNM)[13]。SNM是描述地震頻段噪聲大小的無量綱參數,可用來評估臺站地震頻段噪聲水平的高低,SNM值越大,表示地震頻段的噪聲水平越高[14]。

在研究噪聲隨功率密度的變化時,通常需要選擇沒有地震干擾且外部環境干擾較小的時間段進行功率譜密度 (PSD)分析。然而,這種方法的局限性在于無法全面評估臺站噪聲水平,且工作量巨大。為了解決這個問題,Mcnamara和Raymonnd發展了Peterson的臺站噪聲模型估算方法,引入了概率密度函數統計法(probability density function,PDF)來分析地震觀測臺站的噪聲水平[15]。使用概率密度函數方法可以保持數據的連續性,并直觀地將影響背景噪聲的因素體現在概率密度函數 (PDF) 的概率值中。這種方法能夠充分反映臺站噪聲水平的實時變化情況。目前,國際地震臺網(IRIS)等機構已將這種方法應用于日常地震儀數據質量跟蹤,但在超導重力臺站的背景噪聲分析方面應用較少。因此,本項目將概率密度函數方法應用于超導重力臺站的背景噪聲分析,并與傳統的背景噪聲分析方法進行比較。這將作為評估超導重力臺站背景噪聲的輔助方法,以提供更全面的評價和分析。通過比較2種方法的結果,可以更好地了解超導重力臺站背景噪聲特性,并為進一步的研究和數據分析提供重要依據。

1 傳統方法分析拉薩臺站地震頻段背景噪聲

選擇2012年拉薩超導重力臺站采樣間隔為1 s和1 min的觀測數據,進行預處理并計算SNM:

(1)重力數據標定,將2012年采樣間隔為1 s和1 min的連續重力數據分割成單天數據段;

(2)格值轉換,將電壓值轉換成重力值(μgal)和氣壓值(hPa),圖1(a)為拉薩超導重力臺站2012年1月份原始重力記錄,單位為μgal;

圖1 拉薩超導重力臺站單日觀測數據的功率譜密度

圖1 拉薩超導重力臺站2012年1月份觀測數據預處理

圖3 數據長度分別為1h、3h和24h的平均PSD

(3)去除均值和線性趨勢;

(4)扣除重力潮汐理論值和改正大氣負荷效應,大氣導納值采用的是-0.3 μgal/hPa;

(5)扣除擬合得到的9階多項式,去除殘余潮汐和儀器漂移等長周期項,圖1(b)為重力殘差;

(6)人工剔除單天異常數據如間斷、尖峰、階躍等;

(7)計算單天重力殘差數據的均方根RMS;

(8)選取該月中RMS值最小的5天數據作為該月的平靜期數據,對其進行傅里葉分析,獲得這5天的平均振幅譜;

(9)對選取的5天平靜期數據200～600 s頻率范圍內的平均功率譜 (mean PSD) 計算作為該月的噪聲水平;

(10)根據平均功率譜得到SNM(見式1),其中,PSD單位為μgal2/Hz。新低背景噪音模型(NLNM)為利用超導重力儀計算地震頻段背景噪聲提供參考和約束,其SNM接近于0。

SNM=log10(meanPSD)+2.5

(1)

由式(1)可得,采樣率為1 s的超導重力數據的SNM為1.226,2012年全年最平靜5天為323、209、031、100和305年積日(如表1所示),平靜期5天的平均功率譜密度為0.053 23 μgal2/Hz;采樣率為1 min的SNM為0.909,全年最平靜5天為049、006、013、050和322年積日,平靜期5天的平均功率譜密度為0.025 66 μgal2/Hz。需要注意的是,采樣率為1 min的數據是通過對采樣率為1 s的數據進行低通濾波得到的。該濾波抑制了高于Nyquist頻率(即1/120 s=8.3 mHz)的能量,導致得到的功率譜衰減更快,因此得到的SNM(0.909)也較低。

表1 基于功率譜概率密度函數方法得到的平靜期與傳統方法對比

2 拉薩臺站背景噪聲概率密度函數分析

2.1 計算功率譜密度

為了最大程度地減少重疊后的頻譜泄露效應并增加頻峰的寬度,對各臺各分量的數據進行以下步驟處理:將數據劃分為1 h的數據段,在數據段之間進行50%的重疊,這意味著每個數據段將與相鄰數據段共享一部分數據,這樣做有助于減少重疊效應,并提高數據的連續性和穩定性;將每小時的連續數據分為13個時間段,各段之間進行75%的重疊。這樣可以進一步減小重疊效應,并確保頻譜數據的完整性;利用FFT(快速傅里葉變換)計算每個時間段的功率譜密度值;對于1 h的PSD值,選取13個時間段的PSD平均值,這樣可以得到每小時更加穩定和可靠的PSD值,以用于后續分析和比較;將PSD換算成分貝。

通過以上處理,能夠盡可能減少重疊后的頻譜泄露效應,增加頻峰的寬度,并獲得更準確和可靠的功率譜密度數據。這種方法在數據處理中具有重要意義,可以提高數據的質量和可解釋性。通過上述方法可產生大量的PSD光滑曲線,進而統計得到臺站地震噪聲功率譜密度函數。圖2所示為拉薩超導重力臺站單日觀測數據的功率譜密度。以2012年1月2日的拉薩超導重力臺站觀測數據為例,得到功率譜密度如圖 2(a)所示,圖中包含了47條PSD最大值、最小值、均值、中值和10百分位統計值,90百分位統計值,其中(PSD)平均值為(橘色直線)、中位數(綠色實線)以及第10(黑色直線)和第90百分位統計數據(黑色虛線)。圖2(b)為2012年100年積日拉薩超導重力臺站功率譜密度數據,圖中包含(PSD)平均值、中位數、第20和第95百分位統計值,灰色區域為平均值±95百分位統計值,由圖可知譜密度(PSD)平均值和中位數相差較小,PSD數據的分布比較對稱、集中,所以采用平均PSD作為每天的PSD值。

傳統評估超導重力儀背景噪聲的方法是對一整天的觀測數據進行加窗傅里葉變換,然后計算其PSD。將數據長度為1 h、3 h和24 h的平均PSD統計值進行比較(見圖 3),可以觀察到小時平均值PSD更為平滑,能夠反映整天PSD的變化趨勢。由于超導重力儀的采樣率是1 s,1 h只包含3 600個采樣點,因此頻率間隔較大,頻率分辨率較低?？梢哉J為1 h的數據長度PSD是一天數據PSD的平滑結果。數據長度越長,頻率分辨率越高,PSD細節越豐富。地震頻段背景噪聲關注的是200～600 s的頻段,可選擇使用1 h的數據長度來計算概率密度函數(PDF),計算頻段范圍內的PSD平均值。

2.2 計算概率密度函數PDF

對得到的每小時PSD數據,以1/8倍頻為單位間隔,滑動計算全頻段內的平均功率譜。平均功率值取PSD在Tc處的值

(2)

其中,TS為起始點周期,TL=2×Ts,Ts=Ts×20.125,每次均以20.125倍數增加,滑動計算每一間隔的平均功率值,當TL達到原始數據的最長有效周期停止計算。通過計算可以得到大量的PSD曲線,之后對每條曲線以1 dB為間隔進行頻數統計,這樣對于某一周期Tc,它的概率密度函數可表示為:

P(Tc)=NPTc/NTc

(3)

其中NPTc是在某個1 dB區間內PSD值總數。NTc是在周期Tc上,所有PSD值的總數。對于給定的周期,利用(3)式可得出某一PSD值的概率密度P(Tc)。即可在周期-功率密度圖上將其概率以不同顏色表示,得到功率譜密度概率密度函數 (見圖4),色度條表征PSD的概率密度值,范圍為0%至30%。通過這種方法能直接反映出噪聲的時間變化分布特征。

圖4 2012年100年積日拉薩超導重力臺站功率譜密度(PSD)

為了確認該種方法的準確性,仿照式(1),定義

DNM=log10(PSD)+3.5

(4)

DNM (Daily noise magnitude) 方法是一種用于評估超導重力臺站背景噪聲的方法。與傳統方法不同的是,DNM方法利用每天的所有數據,而不僅僅選擇最平靜期的天數。這種方法更全面地反映了每天的噪聲情況。DNM方法得到的SNM值為0.945(如表2所示),平均功率譜密度為0.027 85 μgal2/Hz,其介于1 s和1 min 2種采樣的平均功率譜密度之間且更傾向于1 min。最平靜5天為052、050、029、031和014年積日,與1 s和1 min的最平靜期分布有交叉。因此,DNM方法可以作為傳統評估超導重力臺站背景噪聲方法的有效補充。通過使用DNM指標,可以更全面地了解每日噪聲幅度信息,從而更準確地評估臺站的背景噪聲水平。

表2 基于功率譜概率密度函數方法得到的平靜期與傳統方法對比

3 結論

(1)1 s采樣得到的SNM為1.226,1 min采樣得到的SNM為0.909,與前人的研究結果一致,以此作為評判新方法的基準。將連續的超導重力數據分割成長度為1 h、3 h和24 h的數據段,比較了不同數據長度的PSD均值,據此計算得到了臺站地震噪聲功率譜密度函數。

(2)引入DNM方法,該方法得到的SNM值位于1 s和1 min采樣結果的中間。因此,該方法可以作為傳統評估超導重力臺站背景噪聲方法的有效補充。通過該方法可以更全面地了解超導重力臺站的背景噪聲特征,并準確評估其背景噪聲水平。