基于時間域的礦山地質(zhì)結(jié)構(gòu)變化評估研究

車正聲

（浙江省有色金屬地質(zhì)勘查局，浙江 紹興 312000）

1 研究動機(jī)與目的

地質(zhì)變化的發(fā)展過程是一個復(fù)雜的行為，諸多地質(zhì)結(jié)構(gòu)源參數(shù)可以幫助了解地質(zhì)結(jié)構(gòu)的破裂行為，而其中地質(zhì)結(jié)構(gòu)發(fā)展時間函數(shù)是了解結(jié)構(gòu)變化機(jī)制的重要參數(shù)之一。地質(zhì)結(jié)構(gòu)變化時間函數(shù)除反應(yīng)了地震從發(fā)生到結(jié)束的時間，也能描繪整個破裂過程的復(fù)雜變化。

傳統(tǒng)上，利用水平(water-level)頻率域解回旋求得結(jié)構(gòu)變化時間函數(shù)，但經(jīng)常會因頻譜相除而造成不穩(wěn)定的現(xiàn)象，近年來，非負(fù)解(non-negative solution)時間域解回旋(time-domain deconvolution)被利用去調(diào)查地質(zhì)結(jié)構(gòu)變化及災(zāi)害時間函數(shù)(例如Ammon et al.,2006)[1]，能夠解出大地震的破裂特征。

本研究利用Hwang(2016)[2]所發(fā)展的非負(fù)解阻尼平滑化時間域解回旋技術(shù)調(diào)查是否能穩(wěn)定求得地質(zhì)結(jié)構(gòu)變化的歷時，希冀未來重新建立地質(zhì)結(jié)構(gòu)變化歷時與結(jié)構(gòu)矩陣的關(guān)系 (例如 Furumoto and Nakanishi,1983)[3]。

本研究先以四個地質(zhì)結(jié)構(gòu)變化為例(表1)，將其結(jié)果與先前的研究比較，探討非負(fù)解平滑化時間域解回旋方法的可用性。

表1 研究所分析地質(zhì)結(jié)構(gòu)變化歷時

2 遠(yuǎn)場P波資料

本研究使用的遠(yuǎn)場P波數(shù)據(jù)來自于美國IRIS地質(zhì)變化數(shù)據(jù)中心。

首先，挑選地質(zhì)變化中心距離30至90的測站數(shù)據(jù)，避免PcP波及PP波的干擾；接著，依據(jù)各個地震的余震分布，判斷此地震的破裂方向；最后，挑選垂直破裂方向且具清晰的初達(dá)P波進(jìn)行分析。而挑選垂直破裂方向的數(shù)據(jù)，主要在避免破裂方向性的影響。最后每筆數(shù)據(jù)在經(jīng)儀器修正后，再經(jīng)0.01至0.5 Hz濾波，而截切P波的時間長度先依各地質(zhì)結(jié)構(gòu)變化的矩陣配合Furumoto and Nakanishi(1983)的經(jīng)驗公式所得震源歷時的兩倍計算，最后再繪圖判讀合適的截切時間長度。

3 方法與原理

本研究所探討的地質(zhì)結(jié)構(gòu)變化災(zāi)害皆為單向地質(zhì)破裂。對于單向地質(zhì)破裂的情況而言，在某一方向觀測到的結(jié)構(gòu)變化歷時可表示為：

TSPT為地質(zhì)結(jié)構(gòu)破裂歷時，L是破裂長度，Vr為破裂速度，C為在地質(zhì)結(jié)構(gòu)災(zāi)害源區(qū)的表面波相速，θ為站方位與災(zāi)害源破裂方向的夾角。

為了移除破裂方向?qū)Y(jié)構(gòu)變化歷時的影響，選取與破裂方向近乎垂直的測站進(jìn)行研究。

首先對這些測站觀測到的P波進(jìn)行濾波得到測站的觀測波，然后依照Ammon et.al.(2006)的方式產(chǎn)生合成P波，合成P波up（t）可由下式求得:

其中M0為結(jié)構(gòu)矩陣αh、ρh、βh分別為P波的速度、密度以及S波的速度，g（△）為幾何擴(kuò)散因子（geometrical spreading factor），a為 地 球 半 徑，Rp、RpP、RsP分 別P波、pP波及sP波的地質(zhì)結(jié)構(gòu)變化范圍模式（radiation pattern），VP、VpP、VsP分 別 為P波、pP波 及sP波 與 地表反射的反射系數(shù)（coefficients of reflection），ih、jh分別為P波、S波，出射角（take-off angle）CP（ih）為接收站的地表反應(yīng)，f（t）為災(zāi)害源時間函數(shù)，tP、tpP、tsP分別為P波、pP波、sP波的走時（travel time），Q（t）為地質(zhì)結(jié)構(gòu)災(zāi)害衰減項，I（t）為儀器反應(yīng)相。其中關(guān)于與擴(kuò)散幾何因子、反射系數(shù)、輻射模式、接收站的地表反應(yīng)等等的計算，可參考Kanamori and Stewart(1976)[4]、Aki and Richard（1980）[5]等等，災(zāi)害波動衰減項目可參考Yoshid（1988）[6]。速度模型采用IASP91速度模型（Kennett and Engdahl,1991）。最后進(jìn)行利用觀測波的數(shù)值與合成波的數(shù)值進(jìn)行解回旋以得到地質(zhì)結(jié)構(gòu)變化歷時[7，8]（如圖1所示）。

圖1 由上往下依次是測站BKE的觀測波圖、合成波圖以及觀測波與合成波解回旋后的圖

若g（t）是格林函數(shù)，是s（t）是災(zāi)害源時間函數(shù)，b（t）是觀測災(zāi)害波動范圍，則g＊s=b，＊表示回旋，其積分型態(tài)為：

上述公式得出，L為基本形態(tài)。

在離散（disrtet）且有限的長度下，其離散形式為：

△t是取樣率，m是災(zāi)害源時間函數(shù)s（t）的總點數(shù)。（2）式可改寫成矩陣形式如下：

令G、S、b分別等于上式三項中的一項，表示為：

從（4）式中求出S，一般是利用最小平方法逆推（如Menke，1989），其觀念是使觀測值與預(yù)測值得誤差平方總和（，eT是e的轉(zhuǎn)置矩陣）最小，即使

圖2 觀測站解回旋后的波形圖

圖3 地質(zhì)結(jié)構(gòu)變化災(zāi)害分布四種分布位置

最后得到的解（Sest）可表示為：

GT為G的轉(zhuǎn)置矩陣。然而（5）式所得結(jié)果并不穩(wěn)定所以在這一過程中加入平滑因子（W）與阻尼值（λ2）最終的到：

Sest即為索要求的震源時間函數(shù)。

4 結(jié)果與討論

表2 震源歷時與WCMT、GCMT的結(jié)果比較[9]

5 結(jié)論

本研究利用非負(fù)解阻尼平滑化時間域解回旋法調(diào)查地質(zhì)結(jié)構(gòu)災(zāi)害變化時間函數(shù)，并在主要變化位置產(chǎn)生合成波用以代替地質(zhì)變化進(jìn)行時間域解回旋，解決了因在地質(zhì)結(jié)構(gòu)變化附近沒有機(jī)構(gòu)機(jī)制相似的災(zāi)害而無法求出地質(zhì)結(jié)構(gòu)變化源時間函數(shù)的問題。通過加入非負(fù)解與阻尼值來使所得的結(jié)果變得穩(wěn)定，嘗試多個阻尼值最終確定當(dāng)阻尼值為100時的波形圖是較為理想的，說明此法是相對較穩(wěn)定的。