三種非線性時變識別方法的對比研究

陳長龍 劉啟方

摘要：

利用2021年日本福島地震（MJMA 7.3）中6個KiK-net臺站的觀測記錄，分別基于移動窗解卷積法、移動窗譜比法和移動窗自相關函數法對場地非線性時變過程進行研究，分析非線性發生的閾值和程度。研究表明：三種方法識別場地非線性時變特征的穩定性不同。移動窗解卷積法更容易獲得較為穩定的土體非線性時變過程，但其反映的是地表至井下平均波速的變化，識別的非線性程度偏低。對于存在強波阻抗比的淺表層土體，移動窗譜比法和移動窗自相關函數法識別的非線性程度更強;對于較低波阻抗比的淺表層土體，則不容易識別出穩定的非線性時變結果。三種方法識別的非線性最強時刻基本一致，皆位于整條記錄加速度峰值附近。基于三種方法識別的6個臺站的非線性閾值約在40～140 gal，剪切波速下降比在16%～53%之間，場地非線性程度較明顯。

關鍵詞：

移動窗解卷積; 移動窗譜比; 自相關函數; 場地非線性; 時變

中圖分類號： P631????? 文獻標志碼：A?? 文章編號： 1000-0844（2023）02-0372-10

DOI：10.20000/j.1000-0844.20211110002

Comparative study on three identification methods of the nonlinear time-

varyin? process： a case study of the 2021 Fukushima earthquake in Japan

CHEN Changlong， LIU Qifang

（Suzhou University of Science and Technology， Suzhou 215011， Jiangsu， China）

Abstract：

Based on the records of six KiK-net stations in the 2021 Fukushima earthquake （MJMA7.3）， the nonlinear time-varying process of the site was studied using the moving-time-window deconvolution method， moving-time-window spectral ratio method， and moving-time-window autocorrelation function method. Then， the threshold and degree of nonlinearity were analyzed. Results show that the ability of the three methods to identify the nonlinear time-varying characteristics of the site differs. The moving-time-window deconvolution method can easily obtain a more stable nonlinear time-varying response for soil， but the degree of nonlinearity is low. The method reflects the change of average wave velocity from the surface to the underground. For shallow soils with a strong wave impedance ratio， the nonlinear degree of identification results derived from the moving-time-window spectral ratio method and moving-time-window autocorrelation function method are strong. These two methods cannot provide stable nonlinear time-varying results of shallow soils with a low wave impedance ratio. The strongest nonlinear moments identified by the three methods are basically the same and are located near the peak acceleration of the whole record. The nonlinear threshold values of the six stations identified by the three methods are about 40-140 gal， and the shear wave velocity drops between 16% and 53%， indicating that the degree of site nonlinearity is high.

Keywords：

moving-time-window deconvolution; moving-time-window spectral ratio; autocorrelation function; site nonlinearity; time-varying

0 引言

場地非線性是地震工程領域的研究熱點，土體在強震動作用下會產生損傷，從而導致應變幅值增加，剪切模量降低及阻尼比增強［1-7］，非線性發生的閾值、程度及恢復過程均是場地非線性研究的主要問題。目前，基于強震動觀測數據直接分析場地非線性的方法主要有譜比法［2-3］、廣義反演法［8］、S變換法［9］、地震干涉測量法［10］和自相關函數法［11］等。這些傳統分析方法一般是利用主震之前的小震記錄確定場地處于線性狀態時的水平，利用主震記錄確定非線性狀態時的水平，利用余震記錄分析非線性恢復情況。基于以上方法，研究人員判斷的不同場地非線性閾值差別較大［1-3］，在100～300 gal之間，非線性的恢復過程也有很大差異，既有震后幾分鐘甚至更短時間內恢復［10］，也有需幾年恢復［2，12］。另外，研究人員還提出了一些非線性程度指標［5，13-14］（DNL、ADNL、PNL等）。

場地經歷的強地震動一般會從較低幅值逐漸增加到較高的震動水平，然后再降低到弱震動水平。因此可根據一條強震記錄作用下不同時段場地的反應，獲得場地土的非線性時變特性。近十幾年來，研究人員將移動時間窗技術與傳統分析方法結合，提出了移動窗譜比法、移動窗解卷積法和移動窗自相關函數法等，針對一條記錄分析不同時間窗內非線性程度的變化，從而獲得非線性時變特性。如Wu等［15］ 利用移動窗解卷積法，針對2011 Tohoku 地震（MW9.0）KiK-net臺網中6個臺站主震記錄的分析表明，場地共振頻率明顯下降，等效剪切波速最大下降到主震前的70%。Sawazaki等［16］利用移動窗譜比法和2000年日本Western-Tottori地震（MW6.7）的主震記錄，對KiK-net臺網的SMNH01臺站場地土的分析表明，在峰值加速度最大值（PGA）附近，剪切模量下降到主震前的52%，等效剪切波速下降30%。Wu等［17］利用移動窗譜比法對KiK-net 臺網中6個臺站1999年1月至2008年5月2 204個地震的觀測記錄分析表明，PGA在60～100 gal時場地發生非線性反應，明顯低于傳統方法的結果。Federico等 ［18］ 利用移動窗解卷積法和2011 Tohoku 地震（MW9.0）主震記錄，分析了IBRH16臺站土體剪切模量的時變特性，Bonilla等［19］則進一步利用該臺站同一記錄，分析了自相關函數法和移動窗解卷積法在識別非線性時變特性的差異，結果表明，兩種方法均可有效識別出場地非線性的發展過程，但自相關函數法獲得的非線性程度大于移動窗解卷積法，且非線性恢復速度更快。

為了分析移動窗解卷積法、譜比法和自相關函數法在識別場地非線性時變特性的差異，本文從KiK-net臺網中選擇6個臺站，其中4個臺站的表層與下臥基巖波阻抗比很大，另外2個臺站土層波速則由表層向下逐漸增大，無強波阻抗比分界面。利用2021年日本福島地震（MJMA7.3）中接收到的主震記錄，分析三種方法對于這兩類不同場地的非線性時變特性識別的差異和可行性。

1 數據選取

2021年2月13日23時8分日本福島縣東部近海地區發生7.3級地震，震源深度約為55 km，震中位于（37.728°N，141.698°E），日本防災技術研究所（NIED）所屬的KiK-net臺網有超過384個臺站記錄到了完整的主震記錄。此次地震與2011年東日本大地震（MW9.0）的震中距離約75 km，是2011年東日本大地震的余震［20］。本文選擇獲得此次地震主震記錄的6個臺站識別場地的非線性時變特性，圖1為臺站和震中位置。表1給出了這6個臺站的土層地質資料［21］，其中1個臺站為C類場地，其余臺站均為D類場地，這些臺站記錄到的地表峰值加速度在300～1 100 gal之間，易發生非線性。其中FKSH17、FKSH18、FKSH19和FKSH11臺站存在強波阻抗比的分界面，而FKSH20和MYGH10臺站剪切波速變化較平緩。

2 研究方法

本文采用移動窗解卷積法、譜比法和自相關函數法分析場地土的非線性時變特征。移動窗解卷積法是對固定寬度時窗內的地表和井下記錄解卷積，獲得一個時間窗內的等效剪切波速，然后沿時間軸按一定的窗口距離移動時間窗，獲得下一個時間窗內結果［如圖2（a）中紅色框表示一個時間窗，黑色框為下一個時間窗］。通過計算每個窗口的PGA（亦稱為窗內PGA），獲得不同時段的土體波速隨時間或窗內PGA的變化。根據等效剪切波速的變化，判斷土體是否發生非線性及其程度。該方法首先將井下和地表處的地震動記錄進行解卷積計算，解卷積結果G（ω）為：

G（ω）=U1（ω）U2（ω）=U1（ω）U*2（ω）U′2（ω）2 （1）

式中：地表和井下加速度記錄的傅里葉譜分別為U1（ω）和U2（ω）;* 表示復共軛;分子為地表與井下加速度在頻域內的互功率譜，分母為井下加速度的功率譜。為穩定計算，一般采用U′2（ω）對U2（ω）的進行校正：

U′2（ω）2=U2（ω）2;當U2（ω）≥0.1maxU2（ω）

0.1maxU2（ω）2;當U2（ω）<0.1maxU2（ω）（2）

然后對解卷積結果進行傅里葉逆變換，獲得井下脈沖輸入時地表的響應。通過脈沖響應獲得井下臺站到地表地震波的傳播時間（一般取正或負時間軸上最大值所對應的時刻），利用式（3）計算土體等效剪切波速。

vS=hΔt （3）

式中：vS代表等效剪切波速;h代表井下臺站到地表的距離;Δt代表井下臺站到地表脈沖波的傳播時間。

移動窗譜比法中窗口的選擇與移動窗解卷積類似，不同的是通過計算每一時間窗內地表和井下記錄傅氏譜比，獲得時間窗內的峰值頻率，得到峰值頻率隨時間或窗內PGA的變化。在進行譜比前，一般需要對傅里葉譜進行平滑，本文采用0.5 Hz的帕曾窗進行平滑。通過峰值頻率的變化，判斷土體是否發生非線性及其程度。設v0和f0是線性狀態土體波速和固有頻率，某一時間窗內識別的波速和峰值頻率為v和f，由于每個時間窗內的峰值頻率為其下覆土層的固有頻率，正比與土層的波速，則土體波速的變化Δv正比峰值頻率的變化Δf［式（4）］，因此，可以通過固有頻率的變化判斷剪切波速的變化。

Δvv0=Δff0 （4）

移動窗自相關函數法［11］是利用一定時間窗內地表記錄進行自相關，從地震事件波形中提取土層某一界面反射波信號，觀察時間延遲，由于反射波延時正比于界面之上土層波速，因此可以通過延時的變化計算波速的變化，如式（5）

Δtt=t0-tt0=v0-vv0=-Δvv0 （5）

式中：t0和v0是土層線性狀態下延時和波速;t和v是某一時間窗自相關函數法計算的延時和波速;Δt是延時變化;Δv波速變化。

采用這三種方法進行分析時，時間窗寬和窗距的選取是移動窗技術的關鍵，其直接決定了結果的穩定性和分辨率。時間窗口寬度過窄，會導致結果的不穩定。時間窗口寬度過長，時間分辨率太低，無法很好地反映頻譜特性隨時間的變化特征。為此，我們參考了相關研究［16-19］，并進行了大量的窗長和窗距測試，最終確定10 s窗長、2 s窗距可獲得穩定結果，上述三種方法采用相同的窗長和窗距。為消除截斷誤差的影響，本文對截取窗口兩端采用10%窗口長度的邊瓣余弦函數進行處理。由于所有觀測記錄未濾波的加速度傅氏譜對應的能量主要集中在1～13 Hz范圍，本文選擇此頻段進行濾波。

3 分析與討論

圖2～圖7給出了利用移動窗解卷積法、譜比法和自相關函數法識別的4個存在強波阻抗比分界面臺站的非線性時變過程。

FKSH17臺站的鉆孔深度100 m，vS30為544 m/s，屬于美國NEHRP規范中的C類場地，鉆孔資料表明（表1），表層6 m是剪切波速為180 m/s 的砂礫，其下是剪切波速為1 100 m/s的風化花崗巖和花崗巖，兩層之間的波阻抗比為9.9（表1），存在強波阻抗比分界面，強波阻抗比分界面之上的土層總厚度為6 m，等效剪切波速為180 m/s。由于風化花崗巖波速較高，非線性很可能主要發生在淺表層6 m的砂礫中，地表的峰值加速度（PGA）為333 gal。三種方法在該臺站識別非線性時變過程大致相同，即非線性程度隨著窗內PGA的增大而增強，當PGA達到整條記錄的峰值時，非線性特征最明顯，隨著PGA的逐漸減小，非線性逐步恢復。從非線性發生的閾值看，三種方法識別出結果差別不大，在60～70 gal之間。非線性下降程度也基本一致，識別的波速最大下降分別為22%、24%和22%（表2），均位于PGA最大值附近。

FKSH18臺站鉆孔深度100 m，vS30為307 m/s，屬于美國NEHRP規范中的D類場地，表層4 m是剪切波速為140 m/s的風化花崗巖，第二層8 m是剪切波速為220 m/s的風化花崗巖，第三層18 m是剪切波速為550 m/s的風化花崗巖。二三層之間的波阻抗比為3.8。強波阻抗比分界面之上的土層總厚度為12 m，等效剪切波速為185 m/s。地表的PGA為424.8 gal，三種方法識別出的非線性閾值相差不大，在56～66 gal之間。三種方法識別的非線性下降程度不同，識別的波速最大下降為29%、40%和43%（表2），譜比法和自相關法的結果較為接近，且明顯大于解卷積法的結果。

FKSH19臺站鉆孔深度100 m，vS30為338 m/s，屬于美國NEHRP規范中的D類場地。表層2 m是剪切波速為170 m/s的表土，第二層18 m是剪切波速為270 m/s的花崗巖，第三層20 m是剪切波速為970 m/s的花崗巖，二三層之間的阻抗比為5.0，強波阻抗比分界面之上的土層總厚度為20 m，等效剪切波速為255 m/s。地表的PGA為513.2 gal，解卷積法和譜比法識別的非線性閾值較接近，分別為141 gal和133 gal，從自相關函數法得到的延時結果看，突變點不明顯，不易識別準確結果，本文識別為55 gal，三種方法識別的波速最大下降為16%、28%和23%（表2），仍表現為譜比法和自相關函數法的結果較為接近，且大于解卷積法的結果。

FKSH11臺站鉆孔深度115 m，vS30為240 m/s，屬于美國NEHRP規范中的D類場地，表層1 m是剪切波速為110 m/s 的礫質土，第二層33 m是剪切波速為250 m/s的粉土混砂礫，第三層22 m是剪切波速為1 200 m/s的溶結凝灰巖，二三層之間的波阻抗比為5.2。強波阻抗比分界面之上的土層總厚度為34 m，等效剪切波速為241 m/s。地表的PGV為353.9 gal，三種方法識別非線性時變特性的能力不同，解卷積法和自相關函數法識別的非線性時變特性較明顯，識別的非線性閾值分別為59 gal和63 gal，波速最大下降為18%、32%（表2），隨著PGA的逐漸減小，波速逐步恢復。但自相關函數法識別的恢復過程振蕩較大。而譜比法則幾乎識別不出非線性的時變過程，該臺站譜比結果在1 Hz、5 Hz、10 Hz處顯示3條峰值曲線，但各峰值曲線隨時間變化均不明顯。

圖8給出了利用移動窗解卷積法、譜比法和自相關函數法識別的2個無強波阻抗比分界面臺站的非線性時變過程。MYGH10臺站鉆孔深度205 m，vS30為348 m/s，屬于美國NEHRP規范中的D類場地。表層1m是剪切波速為110 m/s 的黏性土，第二層2 m是剪切波速為250 m/s的砂巖，第三層31 m是剪切波速為390 m/s的礫巖和砂巖，第四層80 m是剪切波速為590 m/s的凝灰質砂巖，第五層91 m是剪切波速為770 m/s的砂巖，土層剪切波速逐漸增大，無強分界面。該臺站地表PGA很高，達1041.5 gal，值很大。從圖8可以看出，三種方法識別出的非線性閾值相差不大，分別識別的閾值為103、93和98 gal。三種方法識別的非線性下降程度不同，識別的波速最大下降為18%、53%和42%，均位于PGA最大值附近。譜比法和自相關函數法的結果較為接近，且顯著大于解卷積法的結果，約為2.5～3倍。

FKSH20臺站鉆孔深度109 m，vS30為350 m/s，屬于美國NEHRP規范中的D類場地。表層32 m是剪切波速為350 m/s的腐殖土和砂礫，第二層28 m是剪切波速為500 m/s的砂礫，第三層49 m是剪切波速為610 m/s的砂巖，一二層和二三層的波阻抗比為1.4和1.3，阻抗比較小，剪切波速逐漸增大，無強波阻抗比分界面。地表的PGA為345 gal。從圖8可以看出，解卷積法、譜比法分別識別的閾值為39 gal和35 gal，識別的波速最大下降為32%和35%，均很接近。自相關函數法則無法識別非線性時變過程。

從上述分析可見，對于存在強阻抗比分層的4個臺站，三種方法識別的閾值總體上差別不大，但在識別非線性程度上，譜比法和自相關函數法的結果較為接近，除FKSH17臺站，后兩種方法識別的非線性程度均且明顯大于解卷積法的結果，這很可能因為譜比法和自相關函數法更容易反映強阻抗比界面之上的淺層波速變化，而解卷積法則反映地表至井下平均波速的變化，由于強阻抗比界面之下的波速較高，非線性程度可能很低（或者不發生），所以平均波速的變化會低于強阻抗比界面之上土層波速的變化。Bonilla等［19］利用獲得2011年東日本大地震主震記錄的IBRH16臺站，分析了自相關函數法和移動窗解卷積法在識別非線性時變特性的差異。結果表明，自相關函數法獲得的非線性程度大于移動窗解卷積法，本文結果與其基本一致。

對于無強阻抗比界面的土層，波速隨深度變化不明顯，因而通過譜比法和自相關函數法不易獲得穩定的峰值頻率和土層界面的反射波延時，不易獲得非線性時變過程。此時解卷積法仍可得到穩定結果，主要由于該方法給出的是井下臺站到地表的平均波速，強震動作用下土層平均剪切波速隨震動水平的增加仍出現了明顯的下降，因此仍可以較好地反映出土層非線性的平均變化，如FKSH20臺站。同時我們發現，盡管土層中不存在明顯強阻抗比界面，但當場地PGA很高，如MYGH10臺站達到1 041.5 gal，三種方法也均能識別出較穩定非線性時變過程。

4 結論

本文利用2021年日本福島地震（MJMA 7.3）中6個KiK-net 臺站的觀測記錄，分別用移動窗解卷積法、譜比法和自相關函數法分析了場地非線性時變過程，研究表明：

（1） 三種方法識別場地非線性時變特征的穩定性不同。解卷積法更容易獲得較為穩定的土體非線性時變過程，但其反映的是地表至井下平均波速的變化，識別的非線性程度偏低。但對于淺表層存在強阻抗比的土體，譜比法和自相關函數法可以獲得更強的非線性程度，對于淺表層不存在強阻抗比的土體則不容易識別出穩定的非線性時變過程。三種方法識別的非線性閾值基本上相差不大，對于非線性最強時刻的判定也基本相同，皆位于整條記錄加速度峰值附近。

（2） 6個臺站地表峰值加速度在300～1 100 gal之間，三種方法識別的非線性閾值在40～140 gal之間，剪切波速下降比在16%～53%之間，表明場地產生較嚴重的非線性反應。對于PGA最大（1 041.5 gal）的MYGH10臺站，剪切波速最大下降達到53%。

本文僅用6個臺站初步分析了三種方法的優缺點，由于場地非線性的復雜性，擬在今后的研究中，我們將通過更多的鉆井臺陣記錄，進一步分析上述方法的可靠性和有效性。

參考文獻（References）

［1］ BERESNEV I A，WEN U K L.Nonlinear soil response：a reality？［J］.Bulletin of the Seismological Society of America，1996，86（6）：1964-1978.

［2］ WEN K L，BERESNEV I A，YEH Y T.Nonlinear soil amplification inferred from downhole strong seismic motion data［J］.Geophysical Research Letters，1994，21（24）：2625-2628.

［3］ HARTZELL S.Variability in nonlinear sediment response during the 1994 Northridge，California，earthquake［J］.Bulletin of the Seismological Society of America，1998，88（6）：1426-1437.

［4］ 王海云.基于強震觀測數據的土層場地反應的研究現狀［J］.地震工程與工程振動，2014，34（4）：42-47.

WANG Haiyun.A review of study on soil site response estimating from strong motion data［J］.Earthquake Engineering and Engineering Dynamics，2014，34（4）：42-47.

［5］ 任葉飛.基于強震動記錄的汶川地震場地效應研究［D］.哈爾濱：中國地震局工程力學研究所，2014.

REN Yefei.Study on site effect in the Wenchuan earthquake using strong-motion recordings［D］.Harbin：Institute of Engineering Mechanics，CEA，2014.

［6］ 李小軍.非線性場地地震反應分析方法的研究［D］.哈爾濱：中國地震局工程力學研究所，1993.

LI Xiaojun.Study on seismic response analysis method of nonlinear site［D］.Harbin：Institute of Engineering Mechanics，CEA，1993.

［7］ 陳學良.土體動力特性、復雜場地非線性地震反應及其方法研究［D］.哈爾濱：中國地震局工程力學研究所，2006.

CHEN Xueliang.Study on soil dynamic characteristics，nonlinear seismic response of complex site and its methods［D］.Harbin：Institute of Engineering Mechanics，CEA，2006.

［8］ REN Y F，WEN R Z，YAMANAKA H，et al.Site effects by generalized inversion technique using strong motion recordings of the 2008 Wenchuan earthquake［J］.Earthquake Engineering and Engineering Vibration，2013，12（2）：165-184.

［9］ BONILLA L F，TSUDA K，PULIDO N，et al.Nonlinear site response evidence of K-NET and KiK-net records from the 2011 off the Pacific coast of Tohoku earthquake［J］.Earth，Planets and Space，2011，63（7）：785-789.

［10］ YAMADA M，MORI J，OHMI S.Temporal changes of subsurface velocities during strong shaking as seen from seismic interferometry［J］.Journal of Geophysical Research：Solid Earth，2010，115（B3）：B03302.

［11］ SCHERBAUM F.Seismic imaging of the site response using microearthquake recordings.Part Ⅱ.application to the Swabian Jura，southwest Germany，seismic network［J］.Bulletin of the Seismological Society of America，1987，77（6）：1924-1944.

［12］ 王蘇陽.基于日本KiK-net地震動數據的場地反應研究［D］.哈爾濱：中國地震局工程力學研究所，2017.

WANG Suyang.Study on site effect using ground motion data from KiK-net in Japan［D］.Harbin：Institute of Engineering Mechanics，CEA，2017.

［13］ NOGUCHI S，SASATANI T.Quantification of degree of nonlinear site response［C］// Proceedings of the 14th World Conference on Earthquake Engineering.Beijing：Seismological Press，2008.

［14］ REGNIER J，CADET H，BONILLA L F，et al.Assessing nonlinear behavior of soils in seismic site response：statistical analysis on KiK-net strong-motion data［J］.Bulletin of the Seismological Society of America，2013，103（3）：1750-1770.

［15］ WU C Q，PENG Z G.Temporal changes of site response during the 2011 MW9.0 off the Pacific coast of Tohoku earthquake［J］.Earth，Planets and Space，2011，63（7）：791-795.

［16］ SAWAZAKI K，SATO H，NAKAHARA H，et al.Time-lapse changes of seismic velocity in the shallow ground caused by strong ground motion shock of the 2000 western-Tottori earthquake，Japan，as revealed from coda deconvolution analysis［J］.Bulletin of the Seismological Society of America，2009，99（1）：352-366.

［17］ WU C Q，PENG Z G，BEN-ZION Y.Refined thresholds for non-linear ground motion and temporal changes of site response associated with medium-size earthquakes［J］.Geophysical Journal International，2010，182（3）：1567-1576.

［18］ FEDERICO D，BONILLA L F，FOTI S.In-situ shear modulus reduction computation using seismic interferometry by deconvolution from borehole and surface data：theory and examples［C］//Proceedings of the 5th IASPEI/IAEE International Symposium：Effects of Surface Geology in Seismic Motion.Taiwan：［s.n.］，2016：15-17.

［19］ BONILLA L F，GUGUEN P，BEN-ZION Y.Monitoring coseismic temporal changes of shallow material during strong ground motion with interferometry and autocorrelation［J］.Bulletin of the Seismological Society of America，2019，109（1）：187-198.

［20］ 張喆，許力生.2021年2月13日日本本州東海岸MW7.2地震矩心矩張量解［J］.地震學報，2021，43（2）：255-259.

ZHANG Zhe，XU Lisheng.The centroid moment tensor solution of the 13 February 2021 MW7.2 earthquake in the east coast of Honshu，Japan［J］.Acta Seismologica Sinica，2021，43（2）：255-259.

［21］ GARDNER G H F，GARDNER L W，GREGORY A R.Formation velocity and density：the diagnostic basics for stratigraphic traps［J］.Geophysics，1974，39（6）：770-780.

［22］ NIED.NIED K-NET，KiK-net，National Research Institute for Earth Science and Disaster Resilience［Z］.2019.doi：10.17598/NIED.0004.