李利平,陳彥好,靳昊,成帥,孫子正,胡慧江,陳迪楊
(山東大學(xué) 齊魯交通學(xué)院,山東 濟南 250002)
隧道突涌水災(zāi)害是嚴重威脅我國交通基礎(chǔ)設(shè)施建設(shè)安全的地質(zhì)災(zāi)害之一,具有強破壞性、突發(fā)性等特點[1],利用先進技術(shù)實現(xiàn)對地質(zhì)災(zāi)害的時空監(jiān)測預(yù)警是當(dāng)前研究的熱門方向。微震監(jiān)測技術(shù)因其遠程無損監(jiān)測的優(yōu)勢得以在隧道工程、油田服務(wù)、礦山安全監(jiān)測等領(lǐng)域廣泛應(yīng)用[2]。在隧道工程領(lǐng)域,我國巖溶發(fā)育的隧址地區(qū),特別是高山險峻的西部,突涌水災(zāi)害由于高水壓、高地應(yīng)力因素使破壞力激增,災(zāi)變機理更加復(fù)雜。利用微震監(jiān)測技術(shù)對防突結(jié)構(gòu)的微震活動進行監(jiān)測和分析,實現(xiàn)防突結(jié)構(gòu)整體穩(wěn)定性評價,是突涌水災(zāi)害監(jiān)測的有效方法。因此,有必要開展一系列微震監(jiān)測技術(shù)在突涌水災(zāi)害監(jiān)測方面的應(yīng)用研究,建立防突結(jié)構(gòu)穩(wěn)定性評價與微震結(jié)果之間的聯(lián)系,為災(zāi)害監(jiān)測預(yù)警奠定基礎(chǔ)。
在微震監(jiān)測技術(shù)應(yīng)用過程中,震源定位和震源機制分析的關(guān)鍵是P 波初至?xí)r間拾取。突涌水過程存在水力壓裂、通道流體撞擊等微小振動信號,需要從中提取清晰的表征微震事件并實現(xiàn)精準定位。理論方法在理想信號或高信噪比信號中效果較好,但在低信噪比信號中具有一定難度,誤差較大;算法選擇、拾取閾值等人為因素也會導(dǎo)致誤差加大。
對此,各國學(xué)者開展了大量研究。針對瞬態(tài)非平臺的復(fù)雜微震信號,常規(guī)識別檢測方法有時域法、頻域法、時頻域法及綜合法[3]?;陬l域、時頻域的算法也被相繼提出[4],如時間-頻率域能量比[5]、時頻分析[6]等。此外,運用分形理論[7-8]、互相關(guān)[9]、神經(jīng)網(wǎng)絡(luò)[10]、數(shù)字圖像[11]等技術(shù)亦可進行地震波初至?xí)r間拾取。通過對波形數(shù)據(jù)的統(tǒng)計分析,獲取信號的偏斜度和峰度,進而拾取初至?xí)r間的PAI-S/K 法,可獲得較準確的初至?xí)r間。在時域法——數(shù)值P 波初至拾取中,長短時均值比拾取技術(shù)是應(yīng)用最廣泛的技術(shù)之一[12],該技術(shù)利用信號特征函數(shù)STA 和LTA 之比識別P 波初至相位,當(dāng)STA/LTA 超過預(yù)定或動態(tài)指定的閾值時觸發(fā)信號識別。隨著現(xiàn)代計算機技術(shù)的迅猛發(fā)展,人工神經(jīng)網(wǎng)絡(luò)[13]、小波變換[14-15]、波偏振分析[16]等也逐漸應(yīng)用到P 波初至?xí)r間拾取中,一定程度上給P波初至檢測帶來新的突破口,但這些方法需要1 個或多個預(yù)設(shè)標(biāo)準,如檢測間隔、閾值設(shè)置,不同的研究人員對此難以形成共識標(biāo)準,因此分析結(jié)果存在一定差異。
在單自由度振動體系結(jié)構(gòu)抗震研究領(lǐng)域中,已有研究以結(jié)構(gòu)的總地震輸入能及滯回耗能占總耗能的比例評估結(jié)構(gòu)滯回耗能水平,相關(guān)理論研究在該領(lǐng)域已趨于成熟[17],并且用于抗震設(shè)計的彈塑性單自由度體系的總能量譜研究也有較大進展[18],許多學(xué)者采用以“單”代“多”的方法,研究單自由度體系和多自由度體系輸入能的相關(guān)性,通過單自由度體系等效疊加得到多自由度體系的相關(guān)結(jié)論[19]。以上研究成果在微震信號分析方面尚無成熟應(yīng)用,采用單自由度體系這一相對穩(wěn)定的結(jié)構(gòu)處理微震信號,利用體系能量譜的形式估計振動輸入能,嘗試實現(xiàn)對微震信號P 波初至特征的捕捉,提供了一種能量分析拾取算法的新角度。
單自由度體系分為無阻尼和有阻尼2 種,前者振動總是以動能和勢能交換為特征,并未考慮體系能量的耗散,即結(jié)構(gòu)體系在振動過程中總能量保持不變,因而與能量大小密切相關(guān)的振幅始終不變[20]。實際上,單個微震事件信號振幅逐漸衰減,最終消失于背景噪聲中,在單自由度體系中表現(xiàn)為質(zhì)量塊m 逐漸靜止在靜力平衡位置,即質(zhì)量塊m 從微震波動能量輸入后進行的是有阻尼振動。
在單自由度體系中,阻尼主要由2 部分因素構(gòu)成:一是外部介質(zhì)的摩擦阻尼力;二是結(jié)構(gòu)內(nèi)部變形的內(nèi)耗。阻尼力的性質(zhì)比較復(fù)雜,為了進行阻尼結(jié)構(gòu)的振動計算,采用線性粘滯阻尼理論求解。最簡單的單自由度阻尼結(jié)構(gòu)模型見圖1。
圖1 單自由度阻尼結(jié)構(gòu)模型
假定阻尼力與變形速度成正比,但是方向與速度方向相反,公式為:
式中:D(t)為阻尼力;y為質(zhì)量塊的運行速度;c為結(jié)構(gòu)系統(tǒng)的阻尼系數(shù)。結(jié)構(gòu)的阻尼力由多種因素引起,機理復(fù)雜,多數(shù)情況將阻尼系數(shù)理解為實際結(jié)構(gòu)中各種因素的綜合阻尼系數(shù),不影響最終計算結(jié)論。
設(shè)ξ為阻尼比,定義為:
式中:ω為結(jié)構(gòu)自振角頻率。
將式(3)代入式(2),設(shè)方程解為y=ert,則:
式(4)方程特征根為:
式(6)表達了單自由度含阻尼體系質(zhì)量塊在外界振動能量輸入條件下,位移y隨時間t的變化規(guī)律。當(dāng)阻尼比設(shè)置為0.05,自振頻率為1 時,含阻尼自由振動位移-時間曲線見圖2。
圖2 典型含阻尼自由振動位移-時間曲線
為了理解微震信號能量的輸入過程,并理解體系內(nèi)能量轉(zhuǎn)化,對以下關(guān)鍵參量進行說明:
(1)振幅衰減。
由式(6)可知,含阻尼自由振動的振幅為Ae-ξωt,且振幅值按照指數(shù)e-ξωt的規(guī)律迅速衰減。如圖2 所示,對于每1 個周期振幅的衰減情況,采用衰減率η表示振幅的衰減:
由此可見,阻尼造成振幅不斷衰減,從初始時刻開始外界振動輸入能逐漸耗散。
(2)頻率減小。
為進一步表明結(jié)構(gòu)體固有阻尼參數(shù)和自振周期對輸入能衰減過程不同的反應(yīng),基于Matlab 完成相同初始位移、初始速度及時間步長條件下不同阻尼比的單自由度體系位移-時間曲線計算(見圖3)。
圖3 不同阻尼比的單自由度體系位移-時間曲線
如圖3 所示,在輸入體系能量一致的情況下,含有不同阻尼比的體系對輸入能的振動反應(yīng)有所不同。阻尼比越大,原始振動振幅和頻率衰減越明顯,最終趨于靜止,阻尼比對輸入能的耗散作用明顯,同時也反映了單自由度模型對復(fù)雜振動過程的表征能力。本能量分析拾取法即利用單自由度體系的固有阻尼對微震信號輸入能的耗散作用來表征有效微震信號的能量變化過程,從而拾取P 波初至?xí)r間。
研究單自由度體系的微震能量輸入,通過結(jié)構(gòu)彈塑性地震反應(yīng)的時程分析法,計算結(jié)構(gòu)在微震作用下的振動輸入能、阻尼耗能和滯回耗能時程曲線[21],為微震事件信號中P 波初至?xí)r間拾取提供依據(jù)。單自由度體系在微震作用下的動力方程為:
將式(7)進行從y(0)到y(tǒng)(t0)的積分:
上式簡化為:
式中:EK為單自由度體系的相對動能;ED為單自由度體系的阻尼耗能;EA為單自由度體系的變形能,由2 部分組成,一是可恢復(fù)的彈性應(yīng)變能,另一部分是不可恢復(fù)的塑性累計滯回耗能;EI為微震事件輸入能。事實上,如將包含微震事件信號的加速度振動譜作為輸入能,其因阻尼逐漸發(fā)生能量耗散,直至結(jié)構(gòu)停止振動,即體系EK=0。
計算截取2 種典型微震監(jiān)測的現(xiàn)場數(shù)據(jù)作為輸入能,基于Matlab 計算結(jié)構(gòu)體系相對于地面的振動輸入能、動能、阻尼耗能和變形能等各種能量時程曲線,解析微震信號時間序列上的能量變化過程。
將微震彈性波能量作為單自由度體系輸入的總能量,受體系內(nèi)阻尼力和變形作用,總能量以阻尼能和變形能的形式耗散,轉(zhuǎn)化過程可由Matlab 計算實現(xiàn)。單自由度體系能量響應(yīng)計算流程見圖4。
圖4 單自由度體系能量響應(yīng)計算流程
采集某隧道微震監(jiān)測現(xiàn)場數(shù)據(jù),微震信號加速度時間序列見圖5,輸入能時程曲線見圖6。
圖5 微震信號加速度時間序列
圖6 輸入能時程曲線
由圖5 可見,高信噪比微震信號P 波振幅突跳明顯,波形輪廓清晰;低信噪比微震信號P 波突跳振幅不明顯,P 波波至前存在噪聲干擾。針對圖5 的2 種現(xiàn)場典型微震信號進行分析,可得微震信號計算處理結(jié)果,輸入能在體系能轉(zhuǎn)換后的動能響應(yīng)時程曲線見圖7,阻尼耗能時程曲線見圖8。
圖7 動能響應(yīng)時程曲線
圖8 阻尼耗能時程曲線
在圖6 所示各項能量中,輸入能和阻尼耗能的表征曲線整體呈遞增狀。如圖5(a)所示,動能曲線在200 ms 前遞增明顯,但維持時間較短,到達頂峰值后迅速衰減,在400 ms 后微震信號持續(xù)時間內(nèi)在零線之上微小波動。從微震信號原始加速度波形分析,體系質(zhì)量塊受微震輸入能影響,質(zhì)點從靜止開始振動,阻尼能量耗散作用開始,相應(yīng)動能迅速增加,與微震信號轉(zhuǎn)化為單自由度能量響應(yīng)域的理論預(yù)期相符。如圖6(a)所示,較高信噪比微震信號可清晰識別P 波初至振幅突跳位置,在輸入能時程曲線相應(yīng)時刻出現(xiàn)陡增。阻尼耗能隨體系輸入能一開始即發(fā)生能量耗散作用,所以阻尼耗能與輸入能時程曲線變化趨勢相對一致。動能是體系質(zhì)量塊運動狀態(tài)的直接反映,體系動能曲線相對阻尼耗能變化復(fù)雜。在變化趨勢上,動能時程曲線與微震信號加速度譜幾乎同步達到峰值,然后迅速衰減。
低信噪比的微震信號處理是研究難點,微震事件的初至受噪聲影響容易被掩蓋。如圖7 所示,微震信號來自同一次微震事件的加速度波形,但圖7(b)相對于圖7(a),有效信號在抵達之前存在一個噪聲段,幅值相當(dāng)于微震事件信號段峰值的1/2,容易造成算法識別錯誤。
如計算結(jié)果所見,體系輸入能和阻尼耗能隨時間累積,整體保持遞增趨勢。阻尼耗能激增點位置見圖9。動能時程曲線隨微震信號時間序列出現(xiàn)噪聲信號段平穩(wěn)、微震信號段內(nèi)激增后迅速衰減的規(guī)律。如圖7(b)所示,廣泛應(yīng)用的STA/LTA 算法拾取結(jié)果受噪聲干擾,初至識別位置有所提前。這是因為算法觸發(fā)機制設(shè)定在長短時均值比值超過預(yù)設(shè)閾值時發(fā)生拾取,而一般噪聲在加速度時間序列內(nèi)有著和微震事件相似的地震屬性,存在誤判的可能。當(dāng)微震加速度譜轉(zhuǎn)入單自由度體系作為能量變化考慮時,由于噪聲信號具有隨機特性,頻散嚴重,計算的能量累積處于穩(wěn)定水平,如圖7(b)噪聲段的動能時程曲線部分,但微震有效信號具有特定主頻段和震源極性,所以噪聲和微震事件信號在體系各部分能量時程曲線(見圖8(b)、圖9)變化上有清晰的辨識特征。
圖9 阻尼耗能激增點位置
目前,常用的STA/LTA、AIC 等算法對微震初至的拾取,要求設(shè)定檢測間隔或觸發(fā)門檻,一定程度上增加了研究人員的主觀因素[22-23]。將微震引起的質(zhì)點運動加速度記錄轉(zhuǎn)化為體系能量域,使得在時間序列中的波至特征在體系輸入能、阻尼耗能及彈性動能序列上得到凸顯,以能量出現(xiàn)特征反應(yīng)的時刻作為震波初至?xí)r間,進行定位計算。
通過上述計算結(jié)果分析可知,阻尼耗能在P 波相位到達之前為0 或接近0,并在P 波相位之后迅速積累。在阻尼作用下,隨時間推移產(chǎn)生光滑的能量耗散累積曲線。因為阻尼能函數(shù)在低信噪比信號中反應(yīng)仍然敏感,所以阻尼能變化過程可作為追蹤和檢測P 波初至?xí)r間的指標(biāo)。根據(jù)阻尼耗能變化曲線,取耗能激增點位置0.16 s 為微震初至?xí)r刻,比采用長短時均值比檢測延后了0.03 s(見圖10(a))。針對低信噪比的微震信號,微震信號波至在阻尼耗能曲線上會出現(xiàn)第2 次增幅更大的跳躍,由此將微震初至定為1.85 s,與長短時均值比檢測結(jié)果相差較大(見圖10(b))。這是因為長短時均值比的檢測方式一旦檢測閾值確定,對于低信噪比的微震信號來說,其容易受噪聲干擾,這就導(dǎo)致誤將噪聲識別為有效波至。
為了驗證微震初至拾取方法的精度,通過建立震源定位雙曲線控制的定位誤差計算模型[24-25],基于核密度估計函數(shù)得到能量分析拾取法與常用STA/LTA 在同一微震定位模型下的定位誤差分布密度函數(shù),從而判斷2 種方法的拾取精度。以震波從震源到2 個傳感器i和j為例,波至到時差曲線模型見圖11。
圖10 拾取結(jié)果對比
圖11 波至到時差曲線模型
微震彈性波從震源到達傳感器i和j之間的到時差(二維平面)表示為:
式中:vP表示彈性波P 波波速。式(10)是震源定位雙曲線(平面范圍)控制方程,形成的曲線為ti-tj到時差軌跡線,線上每點到2 個傳感器的距離為vP(ti-tj)。在概率論方法中,概率密度函數(shù)是一個描述此變量的輸出值,在某個確定取值點附近的可能性的函數(shù)。對模型臺網(wǎng)內(nèi)所有臺站添加了滿足相同的正態(tài)分布ξ~N(0,σ2I)的旅行時隨機誤差,I為單位矩陣,σ為隨機誤差的方差。利用能量分析拾取法與STA/LTA 計算同一觀測系統(tǒng)下的定位雙曲線,檢驗2 種算法在含隨機誤差情況下的拾取結(jié)果,計算定位誤差概率分布。通過獲得的誤差概率密度分布函數(shù)判斷不同拾取結(jié)果下的定位精度,從而實現(xiàn)對2 種方法P 波初至拾取精度的比較。經(jīng)Matlab 計算,進行2 種方法的反演定位結(jié)果與精確值對比,可得P 波初至拾取結(jié)果誤差分析(見圖12)及各自定位誤差概率密度分布(見圖13)。
圖12 拾取結(jié)果誤差分析
圖13 水平向定位誤差概率密度分布曲線
如圖13 所示,能量分析拾取法的結(jié)果在水平向上誤差控制更好,誤差分布帶寬更窄,定位誤差在7.8 m范圍內(nèi)波動,STA/LTA 的定位誤差則在13.3 m 范圍內(nèi)波動。較小的誤差分布帶寬,表示在零誤差范圍內(nèi)具有較高集中度,說明基于能量分析的拾取方法在反演計算微震事件位置時,對于定位誤差具有更好的控制力。
基于單自由度振動體系的能量分析過程實現(xiàn)有效波動初至的拾取,并采用雙曲線定位模型進行驗證,研究結(jié)果表明:
(1)不同信噪比的微震信號在單自由度振動體系能量分析中可以捕獲明顯特征點。高信噪比信號段在阻尼耗能、動能曲線上激增點凸顯,輸入能曲線激增幅度大。低信噪比信號在阻尼耗能曲線上出現(xiàn)階段性遞增,為微震事件波至拾取提供參考依據(jù)。
(2)由于能量特征點的檢測不需要預(yù)設(shè)閾值,避免了其他預(yù)設(shè)閾值方法中的人為誤差,具有更好的誤差控制,通過數(shù)值計算驗證,能量分析拾取法可獲得微震反演更窄的定位誤差分布。
隧道工程監(jiān)測領(lǐng)域所研究的微震監(jiān)測對事件定位精度比天然地震監(jiān)測的精度要求更高,眾多研究旨在解決這項技術(shù)在尺度變化下的適用性問題。本研究的微震事件初至拾取作為復(fù)雜系統(tǒng)的一環(huán),力求提供一種新方法和思路。要做好定位精度的控制仍需考慮前后環(huán)節(jié)的銜接處理和相互影響,如數(shù)字濾波技術(shù)和定位計算,未來很多相關(guān)研究工作必定繼續(xù)開展。