計及非穩態效應的地脈動功率譜密度包絡修正及應用

郝 雨， 陳學前， 范宣華， 王玉軍， 杜 強， 胡紹全

(中國工程物理研究院 總體工程研究所，四川 綿陽 621999)

地脈動作為一種隨機載荷，對很多大型精密科學裝置的穩定性會產生顯著影響[1-2]，為評估地脈動對科學裝置建筑地基及結構設計的影響，需要獲得相應的地脈動載荷譜曲線[3-4]，而地脈動的測試數據一般是時域信號，需要對信號進行甄別和相應處理，得到地脈動載荷的功率譜密度包絡，為建筑基礎及結構載荷輸入提供參考。

傳統的功率譜密度計算一般針對平穩信號進行[5]。但是，地脈動信號是一個典型的非平穩隨機過程，其功率譜密度是一個隨時間變化的函數，傳統的功率譜密度計算方法只能夠反映一段時間內信號的平均能量特性，如此會使結構偏于危險，而且，功率譜密度計算的時間長度選取會對計算結果產生較大影響。

學者們針對非平穩隨機信號的功率譜密度也展開了一定的研究。例如，王良曦等[6]應用Wiener過程模擬了不平路面的功率譜密度；孔凡等[7]應用局部平穩小波過程，提出了一種非平穩隨機過程的演變功率譜密度估計方法。這兩種方法雖然具有一定的通用性，并且能夠重復體現功率譜密度曲線的時間變化歷程，但是對于大量、長時非平穩時間序列而言，所需的時間和計算空間過大，難以在實踐中應用。史康等[8]通過將風速時間序列分解為時變確定性趨勢成分和平穩零均值隨機成分，分析了實測風場的非平穩特性。但是，這種分解方法是針對某一類問題得到的，難以推廣應用于地脈動信號中。林建生等[9]分析了泉州地震臺的地脈動數據，給出了若干時段的功率譜密度，但是其中沒有考慮功率譜密度的非平穩性。

本文以綿陽地震臺臺站的地脈動監測數據為例，分析了地脈動信號的數值特征,在此基礎上提出了一種考慮地脈動非平穩特性的功率譜密度修正算法，計算了綿陽地區地脈動信號的功率譜密度包絡。該算法消除了PSD(Power Spectral Denaity)計算中人為參數選取對計算結果的影響，同時避免了平均效應導致的地脈動PSD幅值估計偏低的現象。

1 研究對象

選取2015年—2016年綿陽地區花荄臺站的數字地震觀測資料作為研究對象，數據包括垂直(Z向)、東西(E向)、南北(N向)3個方向的地脈動數據。時間間隔為0.01 s，數據單位為counts/s，原始測量數據與加速度信號之間通過如下關系式聯系：

時域關系式為

(1)

根據式(1)，原始測量數據與加速度信號的功率譜之間的關系式為

(2)

式中：x(t)為原始速度測量信號；a(t)為轉換得到的加速度信號，g；Gx(f)和Ga(f)分別為它們的功率譜密度；η為測量靈敏度。

分析時一般選取無地震或其他干擾的時段[9]，去除烈度大于等于Ⅱ級的有感地震及非地震干擾引起的異常信號。具體方法如下：對第一類異常信號，根據文獻[10]，烈度為Ⅴ的地震峰值速度為0.02～0.04 m/s，且烈度每增加一級，地震峰值速度約增加1倍，選取峰值速度2.5 mm/s作為無感地震和有感地震的分界線。計算中,對數據按小時分段，如果某段數據的最大速度峰值大于2.5 mm/s，則去掉整段數據。

第二類異常信號的特點是振動幅值明顯高于常態地脈動，同時相比于無感地震信號，具有持續時間長的特點。根據地震局的測量數據，無感地震的持續時間一般在幾秒甚至更短的數量級[11]。以每10 s為一個窗，如果連續3個窗內有超過5%的數據點大于該小時或當天數據均方根值的5倍，則去掉整段數據。事實上，由于非地震干擾在量級上和持續時間上與其他信號有明顯差異，判斷閾值的選取具有較好的魯棒性。

2 地脈動信號的功率譜密度包絡

2.1 平穩信號的功率譜密度

平穩隨機信號的功率譜密度與時間無關，對零均值離散平穩時間序列xn=x(nΔt)，其單邊功率譜密度

(3)

式中：Δt為離散時間序列的時間步長；N為信號總長度。

為避免產生過多的毛刺，通常采用滑動平均的辦法，求解一段相對較長時間T內的平均功率譜密度。由于平穩隨機信號的功率譜密度與時間無關，時間長度T的選擇不影響功率譜密度計算的結果。

2.2 地脈動信號的非平穩特性

地脈動信號是一個典型的非平穩隨機過程。應用自相關函數(Augmented Dickey-Fuller,ADF)法[12]和游程法[13]對一段地脈動信號進行了平穩性檢驗，結果表明：地脈動信號在從幾十秒～1天的任意一段時間長度內都是非平穩的。經典的做法是將地脈動信號視作局部平穩的，在每時間長度T內視作平穩隨機過程，求出一條功率譜密度曲線，再將所有的曲線合并取包絡。這樣，時間長度T的選擇會對結果產生較大影響，時間長度T越大，求出的功率譜密度包絡幅值越低，原因在于振動量級較高的時間段會被振動量級較低的時間段平均。但是，時間長度T選擇過小時，又會影響功率譜密度的求解精度和頻率分辨率。

以花荄監測點E方向2016-03-01的地脈動信號為例，當時間長度取不同值時，得到的PSD包絡曲線如圖1所示。以不同時間長度計算出的功率譜密度譜形基本一致，但是幅值不同，T=30 s時的功率譜密度曲線比T=10 min時高約10 dB，并且沒有發現明顯的收斂趨勢。這提示我們必須考慮地脈動信號的非平穩特性，對經典的功率譜密度包絡求解方法進行修正。

圖1 時間長度T對功率譜密度包絡的影響Fig.1 Effect of time length T to the PSD envelope

2.3 地脈動信號的功率譜密度修正

在地震分析中，通常可以把非平穩地震信號近似為一個平穩隨機過程和一個反映振幅強度隨時間變化的函數的乘積[14]，即

x(t)=A(t)z(t)

(4)

式中：x(t)為非平穩地震信號；z(t)為平穩隨機過程；A(t)為反映振幅強度的量(見圖2)。

圖2 非平穩地脈動信號的分解Fig.2 Decomposition of non-stationary ground microtremor signal

由于A(t)相對于z(t)是緩變的，式(4)可以寫成頻域的形式

Gx(t,f)=A(t)·Gz(f)

(5)

式中：Gx和Gz分別為x(t)和z(t)的功率譜密度，由于z(t)是平穩隨機過程，Gz與時間無關。圖1中不同曲線的形狀相似性驗證了這一近似的合理性。

將一段時間T分成若干子區間，根據以上假設，在每一子區間內，其功率譜密度是成比例的，有

Gi(f)=AiGz(f)

(6)

式中：Gi為第i個子區間內的地脈動信號功率譜密度；Ai在每個子區間內是一個常數。

考慮到功率譜密度的積分與時域信號均方根值的關系，式(6)可以寫成

(7)

式中：RMS(xi)和RMS(x)分別為第i個子區間內和整段地脈動信號的均方根值。功率譜密度的包絡即可寫為

(8)

式中：Env為該段信號的功率譜密度包絡。

為防止高振幅信號被低振幅信號平均，導致遺漏危險載荷，子區間時間長度應當盡可能小。考慮地脈動信號相鄰兩次穿越0值之間的部分，稱之為一個游程，以每一個游程作為一個子區間，這可以作為子區間的最細劃分。在一個游程內，地脈動時域信號為一個半波，如圖3所示。假設半波的形狀為半正弦波，則該游程內信號的均方根值與最大值之間的關系為

(9)

為驗證這一假設，隨機選取若干段游程，統計其均方根值和峰值的比例，如圖4所示。結果表明，絕大多數游程的均方根值都位于0.5～0.7，式(9)偏于保守，且誤差小于3 dB。

綜合式(8)、式(9)，修正后的功率譜密度包絡可以寫成

(10)

圖3 地脈動信號游程示意圖Fig.3 Sketch of run-length of the ground microtremor signal

這一公式建立起功率譜密度包絡與平均功率譜密度的關系，公式中只需要計算整段信號的均方根值及最大峰值，不需要實際進行游程的劃分和穩態量/非穩態量的分解，對于大批量數據的處理顯著節約了計算量，同時避免了參數和游程選取的任意性。

圖5顯示了時間長度T不同時根據式(10)計算得到的典型地脈動修正功率譜密度包絡曲線，同時與未經修正的包絡曲線進行了對比。結果顯示，經過式(10)修正，功率譜密度包絡的計算與時間長度T的選取基本無關。包絡曲線與未修正的曲線形狀一致，但幅值更高，有效避免了局部高量級振動被平均化的風險。

3 分析結果

應用上述方法，計算了綿陽花荄監測站2015年—2016年地脈動數據的功率譜密度，分析頻率范圍為0～50 Hz，分析頻率間隔為0.048 8 Hz，對應的滑動平均窗長度為2 048個數據點，時間長度T選取為10 min。得到的功率譜密度包絡如圖6所示。在Intel Core i7-4770 CPU@3.4 GHz處理器上，使用MATLAB進行4線程并行計算，處理全部2年共計189億個數據點所需的時間約3.5 h。

從圖6中可見，地脈動信號的能量大部分集中在25～40 Hz的范圍內，最大峰值約為3.1×10-12g2/Hz，其中橫向(E向和N向)的振動幅度略大于垂直方向(Z向)。

圖6 花荄監測站2015年—2016年地脈動功率譜密度包絡Fig.6 PSD envelop of ground microtremor of Huagai monitoring station 2015—2016

4 結 論

本文研究了典型地脈動載荷的包絡計算方法，將非平穩地脈動信號分解為一個平穩隨機過程與時變函數的乘積，提出了考慮信號非平穩特性的功率譜密度包絡修正方法，為大型精密裝置結構穩定性設計中地脈動載荷的確定奠定了基礎。本文方法具有以下優點：

(1)本文方法避免了經典平穩或局部平穩假設下高振幅信號被平均化過程掩蓋的危險，所得結果可以充分反映地脈動信號的危險剖面，保證了結構設計中地脈動載荷的可靠性。

(2)經典平穩假設下的功率譜密度包絡幅值與時間長度的選取有關，時間長度過長時，高振幅更容易被掩蓋，時間長度過短時，數據量的缺乏會降低分析精度，本文通過修正方法，得到的地脈動PSD計算結果不依賴于計算時間長度的選擇。

(3)本文充分利用地脈動信號的時頻特征，在經典平穩功率譜密度的基礎上進行了改進，計算效率高，適于進行大量數據的分析。