隧道突涌水防突結構微震響應初至拾取方法研究

李利平，陳彥好，靳昊，成帥，孫子正，胡慧江，陳迪楊

（山東大學 齊魯交通學院，山東 濟南 250002）

0 引言

隧道突涌水災害是嚴重威脅我國交通基礎設施建設安全的地質災害之一，具有強破壞性、突發性等特點[1]，利用先進技術實現對地質災害的時空監測預警是當前研究的熱門方向。微震監測技術因其遠程無損監測的優勢得以在隧道工程、油田服務、礦山安全監測等領域廣泛應用[2]。在隧道工程領域，我國巖溶發育的隧址地區，特別是高山險峻的西部，突涌水災害由于高水壓、高地應力因素使破壞力激增，災變機理更加復雜。利用微震監測技術對防突結構的微震活動進行監測和分析，實現防突結構整體穩定性評價，是突涌水災害監測的有效方法。因此，有必要開展一系列微震監測技術在突涌水災害監測方面的應用研究，建立防突結構穩定性評價與微震結果之間的聯系，為災害監測預警奠定基礎。

在微震監測技術應用過程中，震源定位和震源機制分析的關鍵是P 波初至時間拾取。突涌水過程存在水力壓裂、通道流體撞擊等微小振動信號，需要從中提取清晰的表征微震事件并實現精準定位。理論方法在理想信號或高信噪比信號中效果較好，但在低信噪比信號中具有一定難度，誤差較大；算法選擇、拾取閾值等人為因素也會導致誤差加大。

對此，各國學者開展了大量研究。針對瞬態非平臺的復雜微震信號，常規識別檢測方法有時域法、頻域法、時頻域法及綜合法[3]。基于頻域、時頻域的算法也被相繼提出[4]，如時間-頻率域能量比[5]、時頻分析[6]等。此外，運用分形理論[7-8]、互相關[9]、神經網絡[10]、數字圖像[11]等技術亦可進行地震波初至時間拾取。通過對波形數據的統計分析，獲取信號的偏斜度和峰度，進而拾取初至時間的PAI-S/K 法，可獲得較準確的初至時間。在時域法——數值P 波初至拾取中，長短時均值比拾取技術是應用最廣泛的技術之一[12]，該技術利用信號特征函數STA 和LTA 之比識別P 波初至相位，當STA/LTA 超過預定或動態指定的閾值時觸發信號識別。隨著現代計算機技術的迅猛發展，人工神經網絡[13]、小波變換[14-15]、波偏振分析[16]等也逐漸應用到P 波初至時間拾取中，一定程度上給P波初至檢測帶來新的突破口，但這些方法需要1 個或多個預設標準，如檢測間隔、閾值設置，不同的研究人員對此難以形成共識標準，因此分析結果存在一定差異。

在單自由度振動體系結構抗震研究領域中，已有研究以結構的總地震輸入能及滯回耗能占總耗能的比例評估結構滯回耗能水平，相關理論研究在該領域已趨于成熟[17]，并且用于抗震設計的彈塑性單自由度體系的總能量譜研究也有較大進展[18]，許多學者采用以“單”代“多”的方法，研究單自由度體系和多自由度體系輸入能的相關性，通過單自由度體系等效疊加得到多自由度體系的相關結論[19]。以上研究成果在微震信號分析方面尚無成熟應用，采用單自由度體系這一相對穩定的結構處理微震信號，利用體系能量譜的形式估計振動輸入能，嘗試實現對微震信號P 波初至特征的捕捉，提供了一種能量分析拾取算法的新角度。

1 單自由度體系振動規律分析

單自由度體系分為無阻尼和有阻尼2 種，前者振動總是以動能和勢能交換為特征，并未考慮體系能量的耗散，即結構體系在振動過程中總能量保持不變，因而與能量大小密切相關的振幅始終不變[20]。實際上，單個微震事件信號振幅逐漸衰減，最終消失于背景噪聲中，在單自由度體系中表現為質量塊m 逐漸靜止在靜力平衡位置，即質量塊m 從微震波動能量輸入后進行的是有阻尼振動。

在單自由度體系中，阻尼主要由2 部分因素構成：一是外部介質的摩擦阻尼力；二是結構內部變形的內耗。阻尼力的性質比較復雜，為了進行阻尼結構的振動計算，采用線性粘滯阻尼理論求解。最簡單的單自由度阻尼結構模型見圖1。

圖1 單自由度阻尼結構模型

假定阻尼力與變形速度成正比，但是方向與速度方向相反，公式為：

式中：D(t)為阻尼力；y為質量塊的運行速度；c為結構系統的阻尼系數。結構的阻尼力由多種因素引起，機理復雜，多數情況將阻尼系數理解為實際結構中各種因素的綜合阻尼系數，不影響最終計算結論。

設ξ為阻尼比，定義為：

式中：ω為結構自振角頻率。

將式（3）代入式（2），設方程解為y=ert，則：

式（4）方程特征根為：

式（6）表達了單自由度含阻尼體系質量塊在外界振動能量輸入條件下，位移y隨時間t的變化規律。當阻尼比設置為0.05，自振頻率為1 時，含阻尼自由振動位移-時間曲線見圖2。

圖2 典型含阻尼自由振動位移-時間曲線

為了理解微震信號能量的輸入過程，并理解體系內能量轉化，對以下關鍵參量進行說明：

（1）振幅衰減。

由式（6）可知，含阻尼自由振動的振幅為Ae-ξωt，且振幅值按照指數e-ξωt的規律迅速衰減。如圖2 所示，對于每1 個周期振幅的衰減情況，采用衰減率η表示振幅的衰減：

由此可見，阻尼造成振幅不斷衰減，從初始時刻開始外界振動輸入能逐漸耗散。

（2）頻率減小。

為進一步表明結構體固有阻尼參數和自振周期對輸入能衰減過程不同的反應，基于Matlab 完成相同初始位移、初始速度及時間步長條件下不同阻尼比的單自由度體系位移-時間曲線計算（見圖3）。

圖3 不同阻尼比的單自由度體系位移-時間曲線

如圖3 所示，在輸入體系能量一致的情況下，含有不同阻尼比的體系對輸入能的振動反應有所不同。阻尼比越大，原始振動振幅和頻率衰減越明顯，最終趨于靜止，阻尼比對輸入能的耗散作用明顯，同時也反映了單自由度模型對復雜振動過程的表征能力。本能量分析拾取法即利用單自由度體系的固有阻尼對微震信號輸入能的耗散作用來表征有效微震信號的能量變化過程，從而拾取P 波初至時間。

2 單自由度體系振動能量反應

研究單自由度體系的微震能量輸入，通過結構彈塑性地震反應的時程分析法，計算結構在微震作用下的振動輸入能、阻尼耗能和滯回耗能時程曲線[21]，為微震事件信號中P 波初至時間拾取提供依據。單自由度體系在微震作用下的動力方程為：

將式（7）進行從y(0)到y(t0)的積分：

上式簡化為：

式中：EK為單自由度體系的相對動能；ED為單自由度體系的阻尼耗能；EA為單自由度體系的變形能，由2 部分組成，一是可恢復的彈性應變能，另一部分是不可恢復的塑性累計滯回耗能；EI為微震事件輸入能。事實上，如將包含微震事件信號的加速度振動譜作為輸入能，其因阻尼逐漸發生能量耗散，直至結構停止振動，即體系EK=0。

3 微震信號能量辨識與初至拾取

計算截取2 種典型微震監測的現場數據作為輸入能，基于Matlab 計算結構體系相對于地面的振動輸入能、動能、阻尼耗能和變形能等各種能量時程曲線，解析微震信號時間序列上的能量變化過程。

將微震彈性波能量作為單自由度體系輸入的總能量，受體系內阻尼力和變形作用，總能量以阻尼能和變形能的形式耗散，轉化過程可由Matlab 計算實現。單自由度體系能量響應計算流程見圖4。

圖4 單自由度體系能量響應計算流程

采集某隧道微震監測現場數據，微震信號加速度時間序列見圖5，輸入能時程曲線見圖6。

圖5 微震信號加速度時間序列

圖6 輸入能時程曲線

由圖5 可見，高信噪比微震信號P 波振幅突跳明顯，波形輪廓清晰；低信噪比微震信號P 波突跳振幅不明顯，P 波波至前存在噪聲干擾。針對圖5 的2 種現場典型微震信號進行分析，可得微震信號計算處理結果，輸入能在體系能轉換后的動能響應時程曲線見圖7，阻尼耗能時程曲線見圖8。

圖7 動能響應時程曲線

圖8 阻尼耗能時程曲線

在圖6 所示各項能量中，輸入能和阻尼耗能的表征曲線整體呈遞增狀。如圖5（a）所示，動能曲線在200 ms 前遞增明顯，但維持時間較短，到達頂峰值后迅速衰減，在400 ms 后微震信號持續時間內在零線之上微小波動。從微震信號原始加速度波形分析，體系質量塊受微震輸入能影響，質點從靜止開始振動，阻尼能量耗散作用開始，相應動能迅速增加，與微震信號轉化為單自由度能量響應域的理論預期相符。如圖6（a）所示，較高信噪比微震信號可清晰識別P 波初至振幅突跳位置，在輸入能時程曲線相應時刻出現陡增。阻尼耗能隨體系輸入能一開始即發生能量耗散作用，所以阻尼耗能與輸入能時程曲線變化趨勢相對一致。動能是體系質量塊運動狀態的直接反映，體系動能曲線相對阻尼耗能變化復雜。在變化趨勢上，動能時程曲線與微震信號加速度譜幾乎同步達到峰值，然后迅速衰減。

低信噪比的微震信號處理是研究難點，微震事件的初至受噪聲影響容易被掩蓋。如圖7 所示，微震信號來自同一次微震事件的加速度波形，但圖7（b）相對于圖7（a），有效信號在抵達之前存在一個噪聲段，幅值相當于微震事件信號段峰值的1/2，容易造成算法識別錯誤。

如計算結果所見，體系輸入能和阻尼耗能隨時間累積，整體保持遞增趨勢。阻尼耗能激增點位置見圖9。動能時程曲線隨微震信號時間序列出現噪聲信號段平穩、微震信號段內激增后迅速衰減的規律。如圖7（b）所示，廣泛應用的STA/LTA 算法拾取結果受噪聲干擾，初至識別位置有所提前。這是因為算法觸發機制設定在長短時均值比值超過預設閾值時發生拾取，而一般噪聲在加速度時間序列內有著和微震事件相似的地震屬性，存在誤判的可能。當微震加速度譜轉入單自由度體系作為能量變化考慮時，由于噪聲信號具有隨機特性，頻散嚴重，計算的能量累積處于穩定水平，如圖7（b）噪聲段的動能時程曲線部分，但微震有效信號具有特定主頻段和震源極性，所以噪聲和微震事件信號在體系各部分能量時程曲線（見圖8（b）、圖9）變化上有清晰的辨識特征。

圖9 阻尼耗能激增點位置

目前，常用的STA/LTA、AIC 等算法對微震初至的拾取，要求設定檢測間隔或觸發門檻，一定程度上增加了研究人員的主觀因素[22-23]。將微震引起的質點運動加速度記錄轉化為體系能量域，使得在時間序列中的波至特征在體系輸入能、阻尼耗能及彈性動能序列上得到凸顯，以能量出現特征反應的時刻作為震波初至時間，進行定位計算。

通過上述計算結果分析可知，阻尼耗能在P 波相位到達之前為0 或接近0，并在P 波相位之后迅速積累。在阻尼作用下，隨時間推移產生光滑的能量耗散累積曲線。因為阻尼能函數在低信噪比信號中反應仍然敏感，所以阻尼能變化過程可作為追蹤和檢測P 波初至時間的指標。根據阻尼耗能變化曲線，取耗能激增點位置0.16 s 為微震初至時刻，比采用長短時均值比檢測延后了0.03 s（見圖10（a））。針對低信噪比的微震信號，微震信號波至在阻尼耗能曲線上會出現第2 次增幅更大的跳躍，由此將微震初至定為1.85 s，與長短時均值比檢測結果相差較大（見圖10（b））。這是因為長短時均值比的檢測方式一旦檢測閾值確定，對于低信噪比的微震信號來說，其容易受噪聲干擾，這就導致誤將噪聲識別為有效波至。

4 拾取誤差分析

為了驗證微震初至拾取方法的精度，通過建立震源定位雙曲線控制的定位誤差計算模型[24-25]，基于核密度估計函數得到能量分析拾取法與常用STA/LTA 在同一微震定位模型下的定位誤差分布密度函數，從而判斷2 種方法的拾取精度。以震波從震源到2 個傳感器i和j為例，波至到時差曲線模型見圖11。

圖10 拾取結果對比

圖11 波至到時差曲線模型

微震彈性波從震源到達傳感器i和j之間的到時差（二維平面）表示為：

式中：vP表示彈性波P 波波速。式（10）是震源定位雙曲線（平面范圍）控制方程，形成的曲線為ti-tj到時差軌跡線，線上每點到2 個傳感器的距離為vP（ti-tj）。在概率論方法中，概率密度函數是一個描述此變量的輸出值，在某個確定取值點附近的可能性的函數。對模型臺網內所有臺站添加了滿足相同的正態分布ξ～N（0，σ2I）的旅行時隨機誤差，I為單位矩陣，σ為隨機誤差的方差。利用能量分析拾取法與STA/LTA 計算同一觀測系統下的定位雙曲線，檢驗2 種算法在含隨機誤差情況下的拾取結果，計算定位誤差概率分布。通過獲得的誤差概率密度分布函數判斷不同拾取結果下的定位精度，從而實現對2 種方法P 波初至拾取精度的比較。經Matlab 計算，進行2 種方法的反演定位結果與精確值對比，可得P 波初至拾取結果誤差分析（見圖12）及各自定位誤差概率密度分布（見圖13）。

圖12 拾取結果誤差分析

圖13 水平向定位誤差概率密度分布曲線

如圖13 所示，能量分析拾取法的結果在水平向上誤差控制更好，誤差分布帶寬更窄，定位誤差在7.8 m范圍內波動，STA/LTA 的定位誤差則在13.3 m 范圍內波動。較小的誤差分布帶寬，表示在零誤差范圍內具有較高集中度，說明基于能量分析的拾取方法在反演計算微震事件位置時，對于定位誤差具有更好的控制力。

5 結束語

基于單自由度振動體系的能量分析過程實現有效波動初至的拾取，并采用雙曲線定位模型進行驗證，研究結果表明：

（1）不同信噪比的微震信號在單自由度振動體系能量分析中可以捕獲明顯特征點。高信噪比信號段在阻尼耗能、動能曲線上激增點凸顯，輸入能曲線激增幅度大。低信噪比信號在阻尼耗能曲線上出現階段性遞增，為微震事件波至拾取提供參考依據。

（2）由于能量特征點的檢測不需要預設閾值，避免了其他預設閾值方法中的人為誤差，具有更好的誤差控制，通過數值計算驗證，能量分析拾取法可獲得微震反演更窄的定位誤差分布。

隧道工程監測領域所研究的微震監測對事件定位精度比天然地震監測的精度要求更高，眾多研究旨在解決這項技術在尺度變化下的適用性問題。本研究的微震事件初至拾取作為復雜系統的一環，力求提供一種新方法和思路。要做好定位精度的控制仍需考慮前后環節的銜接處理和相互影響，如數字濾波技術和定位計算，未來很多相關研究工作必定繼續開展。