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

（山東大學(xué) 齊魯交通學(xué)院，山東 濟(jì)南 250002）

0 引言

隧道突涌水災(zāi)害是嚴(yán)重威脅我國(guó)交通基礎(chǔ)設(shè)施建設(shè)安全的地質(zhì)災(zāi)害之一，具有強(qiáng)破壞性、突發(fā)性等特點(diǎn)[1]，利用先進(jìn)技術(shù)實(shí)現(xiàn)對(duì)地質(zhì)災(zāi)害的時(shí)空監(jiān)測(cè)預(yù)警是當(dāng)前研究的熱門方向。微震監(jiān)測(cè)技術(shù)因其遠(yuǎn)程無損監(jiān)測(cè)的優(yōu)勢(shì)得以在隧道工程、油田服務(wù)、礦山安全監(jiān)測(cè)等領(lǐng)域廣泛應(yīng)用[2]。在隧道工程領(lǐng)域，我國(guó)巖溶發(fā)育的隧址地區(qū)，特別是高山險(xiǎn)峻的西部，突涌水災(zāi)害由于高水壓、高地應(yīng)力因素使破壞力激增，災(zāi)變機(jī)理更加復(fù)雜。利用微震監(jiān)測(cè)技術(shù)對(duì)防突結(jié)構(gòu)的微震活動(dòng)進(jìn)行監(jiān)測(cè)和分析，實(shí)現(xiàn)防突結(jié)構(gòu)整體穩(wěn)定性評(píng)價(jià)，是突涌水災(zāi)害監(jiān)測(cè)的有效方法。因此，有必要開展一系列微震監(jiān)測(cè)技術(shù)在突涌水災(zāi)害監(jiān)測(cè)方面的應(yīng)用研究，建立防突結(jié)構(gòu)穩(wěn)定性評(píng)價(jià)與微震結(jié)果之間的聯(lián)系，為災(zāi)害監(jiān)測(cè)預(yù)警奠定基礎(chǔ)。

在微震監(jiān)測(cè)技術(shù)應(yīng)用過程中，震源定位和震源機(jī)制分析的關(guān)鍵是P 波初至?xí)r間拾取。突涌水過程存在水力壓裂、通道流體撞擊等微小振動(dòng)信號(hào)，需要從中提取清晰的表征微震事件并實(shí)現(xiàn)精準(zhǔn)定位。理論方法在理想信號(hào)或高信噪比信號(hào)中效果較好，但在低信噪比信號(hào)中具有一定難度，誤差較大；算法選擇、拾取閾值等人為因素也會(huì)導(dǎo)致誤差加大。

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

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

1 單自由度體系振動(dòng)規(guī)律分析

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

在單自由度體系中，阻尼主要由2 部分因素構(gòu)成：一是外部介質(zhì)的摩擦阻尼力；二是結(jié)構(gòu)內(nèi)部變形的內(nèi)耗。阻尼力的性質(zhì)比較復(fù)雜，為了進(jìn)行阻尼結(jié)構(gòu)的振動(dòng)計(jì)算，采用線性粘滯阻尼理論求解。最簡(jiǎn)單的單自由度阻尼結(jié)構(gòu)模型見圖1。

圖1 單自由度阻尼結(jié)構(gòu)模型

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

式中：D(t)為阻尼力；y為質(zhì)量塊的運(yùn)行速度；c為結(jié)構(gòu)系統(tǒng)的阻尼系數(shù)。結(jié)構(gòu)的阻尼力由多種因素引起，機(jī)理復(fù)雜，多數(shù)情況將阻尼系數(shù)理解為實(shí)際結(jié)構(gòu)中各種因素的綜合阻尼系數(shù)，不影響最終計(jì)算結(jié)論。

設(shè)ξ為阻尼比，定義為：

式中：ω為結(jié)構(gòu)自振角頻率。

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

式（4）方程特征根為：

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

圖2 典型含阻尼自由振動(dòng)位移-時(shí)間曲線

為了理解微震信號(hào)能量的輸入過程，并理解體系內(nèi)能量轉(zhuǎn)化，對(duì)以下關(guān)鍵參量進(jìn)行說明：

（1）振幅衰減。

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

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

（2）頻率減小。

為進(jìn)一步表明結(jié)構(gòu)體固有阻尼參數(shù)和自振周期對(duì)輸入能衰減過程不同的反應(yīng)，基于Matlab 完成相同初始位移、初始速度及時(shí)間步長(zhǎng)條件下不同阻尼比的單自由度體系位移-時(shí)間曲線計(jì)算（見圖3）。

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

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

2 單自由度體系振動(dòng)能量反應(yīng)

研究單自由度體系的微震能量輸入，通過結(jié)構(gòu)彈塑性地震反應(yīng)的時(shí)程分析法，計(jì)算結(jié)構(gòu)在微震作用下的振動(dòng)輸入能、阻尼耗能和滯回耗能時(shí)程曲線[21]，為微震事件信號(hào)中P 波初至?xí)r間拾取提供依據(jù)。單自由度體系在微震作用下的動(dòng)力方程為：

將式（7）進(jìn)行從y(0)到y(tǒng)(t0)的積分：

上式簡(jiǎn)化為：

式中：EK為單自由度體系的相對(duì)動(dòng)能；ED為單自由度體系的阻尼耗能；EA為單自由度體系的變形能，由2 部分組成，一是可恢復(fù)的彈性應(yīng)變能，另一部分是不可恢復(fù)的塑性累計(jì)滯回耗能；EI為微震事件輸入能。事實(shí)上，如將包含微震事件信號(hào)的加速度振動(dòng)譜作為輸入能，其因阻尼逐漸發(fā)生能量耗散，直至結(jié)構(gòu)停止振動(dòng)，即體系EK=0。

3 微震信號(hào)能量辨識(shí)與初至拾取

計(jì)算截取2 種典型微震監(jiān)測(cè)的現(xiàn)場(chǎng)數(shù)據(jù)作為輸入能，基于Matlab 計(jì)算結(jié)構(gòu)體系相對(duì)于地面的振動(dòng)輸入能、動(dòng)能、阻尼耗能和變形能等各種能量時(shí)程曲線，解析微震信號(hào)時(shí)間序列上的能量變化過程。

將微震彈性波能量作為單自由度體系輸入的總能量，受體系內(nèi)阻尼力和變形作用，總能量以阻尼能和變形能的形式耗散，轉(zhuǎn)化過程可由Matlab 計(jì)算實(shí)現(xiàn)。單自由度體系能量響應(yīng)計(jì)算流程見圖4。

圖4 單自由度體系能量響應(yīng)計(jì)算流程

采集某隧道微震監(jiān)測(cè)現(xiàn)場(chǎng)數(shù)據(jù)，微震信號(hào)加速度時(shí)間序列見圖5，輸入能時(shí)程曲線見圖6。

圖5 微震信號(hào)加速度時(shí)間序列

圖6 輸入能時(shí)程曲線

由圖5 可見，高信噪比微震信號(hào)P 波振幅突跳明顯，波形輪廓清晰；低信噪比微震信號(hào)P 波突跳振幅不明顯，P 波波至前存在噪聲干擾。針對(duì)圖5 的2 種現(xiàn)場(chǎng)典型微震信號(hào)進(jìn)行分析，可得微震信號(hào)計(jì)算處理結(jié)果，輸入能在體系能轉(zhuǎn)換后的動(dòng)能響應(yīng)時(shí)程曲線見圖7，阻尼耗能時(shí)程曲線見圖8。

圖7 動(dòng)能響應(yīng)時(shí)程曲線

圖8 阻尼耗能時(shí)程曲線

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

低信噪比的微震信號(hào)處理是研究難點(diǎn)，微震事件的初至受噪聲影響容易被掩蓋。如圖7 所示，微震信號(hào)來自同一次微震事件的加速度波形，但圖7（b）相對(duì)于圖7（a），有效信號(hào)在抵達(dá)之前存在一個(gè)噪聲段，幅值相當(dāng)于微震事件信號(hào)段峰值的1/2，容易造成算法識(shí)別錯(cuò)誤。

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

圖9 阻尼耗能激增點(diǎn)位置

目前，常用的STA/LTA、AIC 等算法對(duì)微震初至的拾取，要求設(shè)定檢測(cè)間隔或觸發(fā)門檻，一定程度上增加了研究人員的主觀因素[22-23]。將微震引起的質(zhì)點(diǎn)運(yùn)動(dòng)加速度記錄轉(zhuǎn)化為體系能量域，使得在時(shí)間序列中的波至特征在體系輸入能、阻尼耗能及彈性動(dòng)能序列上得到凸顯，以能量出現(xiàn)特征反應(yīng)的時(shí)刻作為震波初至?xí)r間，進(jìn)行定位計(jì)算。

通過上述計(jì)算結(jié)果分析可知，阻尼耗能在P 波相位到達(dá)之前為0 或接近0，并在P 波相位之后迅速積累。在阻尼作用下，隨時(shí)間推移產(chǎn)生光滑的能量耗散累積曲線。因?yàn)樽枘崮芎瘮?shù)在低信噪比信號(hào)中反應(yīng)仍然敏感，所以阻尼能變化過程可作為追蹤和檢測(cè)P 波初至?xí)r間的指標(biāo)。根據(jù)阻尼耗能變化曲線，取耗能激增點(diǎn)位置0.16 s 為微震初至?xí)r刻，比采用長(zhǎng)短時(shí)均值比檢測(cè)延后了0.03 s（見圖10（a））。針對(duì)低信噪比的微震信號(hào)，微震信號(hào)波至在阻尼耗能曲線上會(huì)出現(xiàn)第2 次增幅更大的跳躍，由此將微震初至定為1.85 s，與長(zhǎng)短時(shí)均值比檢測(cè)結(jié)果相差較大（見圖10（b））。這是因?yàn)殚L(zhǎng)短時(shí)均值比的檢測(cè)方式一旦檢測(cè)閾值確定，對(duì)于低信噪比的微震信號(hào)來說，其容易受噪聲干擾，這就導(dǎo)致誤將噪聲識(shí)別為有效波至。

4 拾取誤差分析

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

圖10 拾取結(jié)果對(duì)比

圖11 波至到時(shí)差曲線模型

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

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

圖12 拾取結(jié)果誤差分析

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

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

5 結(jié)束語(yǔ)

基于單自由度振動(dòng)體系的能量分析過程實(shí)現(xiàn)有效波動(dòng)初至的拾取，并采用雙曲線定位模型進(jìn)行驗(yàn)證，研究結(jié)果表明：

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

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

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