改進的Prony算法提取舒曼諧振參數

周洪娟，于海雁，喬曉林

(哈爾濱工業大學(威海)信息工程研究所，264209 山東威海)

改進的Prony算法提取舒曼諧振參數

周洪娟，于海雁，喬曉林

(哈爾濱工業大學(威海)信息工程研究所，264209 山東威海)

為在更短的時間內提取舒曼諧振參數，引入對時域采樣數據進行直接處理的Prony算法.針對該算法受噪聲影響大的缺點，對算法的核心過程——超定方程的求解方法進行了改進，采用Martin等人提出的最小殘數法，并與頻域處理結果進行對比分析，驗證該算法的有效性.

舒曼諧振;Prony算法;非線性最小平方擬合;諧振頻率;最小殘數法

舒曼諧振(SR)是由全球閃電活動激勵在地-電離層腔中固定存在的一種諧振現象，各階諧振波長與地球周長成比例，前四階諧振頻率近似為8、14、20和26 Hz，舒曼諧振各階參數的研究對全球閃電活動、氣溫的微弱變化以及大氣層和電離層參數的變化具有重要意義［1－5］，另外，最近也有關于 SR與地震關系的研究［6－7］.但傳感器在接收低頻環境磁場數據過程中，通常會受到其他干擾和噪聲的影響，如人類在周圍空間的活動、近處的雷暴活動以及強烈的地殼運動等，都會使得測量的背景噪聲增強，從而淹沒極其微弱的舒曼諧振信號.當然，這些干擾因素是不可避免的，我們能做的就是選擇遠離工業區的偏遠地方作為觀測站，以及采用抗干擾和噪聲比較好的算法來提取舒曼諧振特征參數.

對舒曼諧振參數的提取，最常規的方法就是在頻域上處理，但對于信噪比很低的觀測數據往往要進行長時間的積累平滑以得到比較清晰的頻譜圖，這對于一些持續時間比較短的突發自然現象的研究，如太陽耀斑、磁暴等對SR參數的影響的研究是不利的.本文引入基于時域采樣數據的諧波估計算法——Prony算法，該算法在電力系統諧波和間諧波的檢測方面應用比較多［8－10］，可以在較短的時間內提取諧波參數，但是經典Prony算法的諧波估計性能受噪聲影響較大.本文首先對Prony算法在低信噪比數據中的應用進行討論，并提出改進算法，以位于云南某觀測站的低頻磁場數據為例，分別在頻域和時域估計舒曼諧振的參數，再進行對比分析，得到比較可靠的舒曼諧振參數時域估計方法.

1 數據說明

觀測站位于比較偏遠的云南山區，主要觀測ELF頻段環境磁場數據，觀測帶寬40 Hz左右，采樣率為100 Hz.觀測數據中的干擾和噪聲很多，舒曼諧振信號是淹沒在這些干擾和噪聲之下的，圖1(a)為10 s的采樣數據，數據中噪聲污染嚴重，對應的頻譜如圖1(b)藍線所示，幾乎看不出SR的輪廓，這通常需要一段時間的積累才能得到比較清晰的頻譜，圖1(b)紅線為10 min的數據做傅里葉變換后的平均頻譜，頻率分辨率為0.1 Hz，可明顯看到前四階舒曼諧振信號的輪廓，諧振頻率分別在8、14、20和26 Hz左右，但頻譜數據極不平滑，不利于對舒曼諧振參數的進一步提取與分析.

2 Prony算法分析

2.1 算法簡介

Prony算法是用復指數函數的線性組合來描述等間距采樣數據的一種數學方法，從而根據采樣數據直接估計出信號的頻率、衰減、幅值和初相位.對實數采樣數據，算法的主要思想如下:設時域信號為一系列存在衰減的正弦信號之和，由三角函數的變換，可以寫成一系列復指數函數之和，如下式所示:

進一步，令

其中:p=2l為估計的復指數函數階數;bm=Amejθm或 bm= Ame-jθm， 互 為 共 軛 對;Zm=e-?mts+j2πfmts或 Zm=e-?mts-j2πfmts，也是互為共軛對.

2.2 算法改進

對任意超定方程如下所示:

其中:G∈Rm×n，且m ＞ n;a ∈Rn，y∈Rn.式(5)的過程相當于對式(6)求最小平方和估計，即有

為了避免ρ0項的影響，Martin等人采用如下的方法，稱為最小殘數法［12］:即構造矩陣G的子集G1，同時取對應向量y的子集y1，有G1a=y1，兩邊同乘以1個矩陣，有

其中G2也為矩陣G的子集，同時要確保矩陣G1中不包含ρ0項.

按照這一思路，任取1個值n，使其滿足p+1≤n≤N-p，將式(3)寫成式(6)的標準矩陣形式，取

按照式(8)可以構造如下的關系式:

可見式(9)與式(5)相比，避免了部分噪聲的影響.

2.3 算法仿真

采用3個頻率分別為8、14和20 Hz的衰減信號之和為例，再加入一定的噪聲，分析噪聲對基于式(5)和式(9)的Prony算法的影響.采樣率設為100 Hz，分析時長為1 s.

圖2中實線、帶點實線和帶星號實線分別表示未加噪聲σ=0、加入噪聲σ=0.2和σ=0.6三種情況下，按照式(5)和式(7)的傳統Prony算法估計出的信號波形以及頻譜圖，估計出的參數見表1，其中估計階數均定為p=10.可見傳統的Prony方法對理想的信號估計精度很高，但對噪聲非常敏感，更不適用于SR信號的處理.而表2為同樣的測試數據采用式(9)估計出的參數，明顯發現該算法性能優于傳統的Prony方法，尤其是對衰減因子估計性能的改善.

圖2 傳統Prony算法在不同的噪聲污染情況下恢復的信號波形和頻譜圖

表1 基于傳統的Prony算法估計結果

表2 基于改進算法的Prony算法估計結果

3 SR特征參數提取

SR參數包含各階的中心頻率、幅度、以及帶寬等.本節首先采用頻域方法對SR參數進行提取，以此驗證Prony算法的有效性.

3.1 頻域提取

對圖1(b)的頻譜采用非線性最小平方擬合方法進行擬合，擬合后的曲線如圖3紅色粗線所示，藍色曲線為擬合前的曲線，前四階舒曼諧振的中心頻率和幅度如表3所示.這種方法估計結果比較準確，但時間分辨率比較低，在一些特殊的分析中，如太陽耀斑時引發的SR頻率的短暫變化等應用中會有局限性.

圖3 擬合后的SR頻譜圖

表3 基于頻域方法提取的SR參數

3.2 時域提取

Prony算法的估計階數采用基于奇異值分解的階數確定方法［13］，對噪聲的有效抑制采用式(9)的方法.通常情況下，若對波形進行估計，Prony算法一次參與計算的數據點數不能太多，因為隨著時間的增大，衰減因子和頻率估計的誤差就會被逐漸放大，而使得估計波形的誤差越來越大，尤其對干擾比較嚴重的數據，這種情況更明顯［14］.因為波形本身并不被關心，而是希望得到比較準確的參數估計，因此本文采用20 s數據為一處理時長，對圖1(b)的10 min平滑頻譜對應的數據進行分段Prony諧波估計，得到前兩階諧振的頻率變化如圖4所示.

其中紅線代表30個估計值的平均值，對應前兩階SR的平均值分別為7.359 1 Hz和13.991 6 Hz，與頻域估計值基本吻合，但存在差別，主要因為觀測數據的噪聲污染很嚴重.Prony算法雖然估計性能會受到噪聲的影響，但估計的時間分辨率可以很高，要得到比較穩定的估計值，下一步還需要對估計階數采用更靈活的確定方法，采用抗噪聲能力更強的算法.

圖4 采用Prony算法估計的SR前兩階頻率

4 結論

本文利用擴展的Prony算法對信噪比接近0 dB的SR觀測數據進行參數估計，采用了Martin等人對含有噪聲的超定方程的求解思路，避免了超定方程系數矩陣中的自相關函數項，提高了Prony算法的抗噪聲性能.實際測試表明，這種方法估計出的SR頻率與頻域方法的結論基本吻合.但是，為了得到Prony算法對SR觀測數據的穩定估計，還需要做進一步的研究，包括對數據的預處理、估計階數的自適應選擇以及更魯棒的抗噪聲算法等.

［1］SHVETS A V，HPBARA Y，HAYALAWA M.Variations of the global lightning distribution revealed from three-station schumann resonance measurements［J］.Journal of Geophysical Research，2010，115(A12)，doi:10.1029/2010JA015851.

［2］DE S S，DE B K，BANDYOPADHYAY B，et al.Studies on the shift in the frequency of the first Schumann resonance mode during a solar proton event［J］.Journal of Atmospheric and Solar-Terrestrial Physics，2010，72(11/12):829－836.

［3］WILLIAM E R，SATORI G.Solar radiation-induced changes in ionospheric height and the Schumann resonance waveguide on different timescales［J］.Radio Science，2007，42(RS2S11)，doi:10.1029/2006RS003494.

［4］ROLDUGIN V C，MALTSEV Y P，VASILJEV A N，et al.Changes of schumann resonance parameters during the solar proton event of 14 July 2000［J］.Journal of Geophysical Research，2003，108(A3)，doi:10.1029/2002JA009495.

［5］SIMOES F，GRARD R，HAMELIN M，et al.The schumann resonance:a tool for exploring the atmospheric environment and the subsurface of the planets and their satellites［J］.Icarus，2008，194(1):30 － 41.

［6］OHTA K，WATANABE N，HAYAKAWA M.Survey of anomalous schumann resonance phenomena observed in Japan，in possible association with earthquakes in Taiwan［J］.Physics and Chemistry of the Earth，2006，31(4/5/6/7/8/9):397－402.

［7］OHTA K，IZUTSU J，HAYAKAWA M.Anomalous excitation of schumann resonances and additional anomalous resonances before the 2004 mid-niigata prefecture earthquake and the 2007 Noto Hantou earthquake［J］.Physics and Chemistry of the Earth，2008，34(6/7):441 －448.

［8］沈楊，戴本祁，張惠東.基于小波和擴展Prony算法的電壓閃變檢測新方法［J］.電力系統保護與控制，2010，38(10):43 －47.

［9］ZYGARLICKI J，ZYGARLICKA M，MROCZKA J，et al.A reduced Prony’s method in power-quality analysis-parameters selection［J］.IEEE Transactions on Power Delivery，2010，25(2):979 －986.

［10］CHEN C I，CHANG G W.An effective time-domain approach based on Prony’s method for time-varying power system harmonics estimation［J］.Power and Energy Society General Meeting，2009(7):1 －6.

［11］張賢達.現代信號處理［M］.北京:清華大學出版社，2002:119－126.

［12］FULLEKRUG M.On the minimization of correlated residuals［J］.Geophysical Journal International，1996，126(1):63－68.

［13］熊俊杰，邢衛榮，萬秋蘭.Prony算法的低頻振蕩主導模式辨識［J］.東南大學學報(自然科學版)，2008，38(1):64－68.

［14］董航，劉滌塵，鄒江峰.基于Prony算法的電力系統低頻振蕩分析［J］.高電壓技術，2006，32(6):97－100.

Improved prony algorithm for abstracting schumann resonance parameters

ZHOU Hong-juan，YU Hai-yan，QIAO Xiao-lin

(Institute of Information Engineering，Harbin Institute of Technology at Weihai，264209 Weihai Shandong，China)

In order to abstract Schumann Resonance parameters in shorter time，Prony algorithm is adopted and analyzed here to process the sampling data directly without FFT.Due to the fact that the algorithm is noise-sensitive，the minimization of correlated residuals proposed by Martin is adopted for solving the overdetermined equations and verified with simulated signals.It is proved to be feasible for SR signal processing after compared with the results from frequency domain.

schumann resonance;prony algorithm;non-linear least-square fit;resonance frequency;minimization of correlated residuals

TN911.72

A

0367－6234(2012)07－0083－04

2011－04－12.

山東省自然科學基金資助項目(ZR2010DM013);科技部中日合作基金資助項目(2012DFG20510).

周洪娟(1980—)，女，講師;

喬曉林(1948—)，男，教授，博士生導師.

周洪娟，hongjue_zhou@sina.com.

(編輯 張 宏)