李德偉楊瑞召孟令彬

中國礦業(yè)大學(xué)(北京)地球科學(xué)與測(cè)繪工程學(xué)院,北京 100083

微地震事件定位面臨兩大難題,一是發(fā)震時(shí)刻難以確定,二是微地震事件的信號(hào)難以識(shí)別[1]。相對(duì)于井中采集的信號(hào),地面檢波器陣列通常信噪比較低,原因主要是震源到地面檢波器距離較遠(yuǎn),導(dǎo)致子波衰減嚴(yán)重,其次是地面觀測(cè)陣列更容易受噪聲(施工噪聲、環(huán)境噪聲等)影響[2]。

為了克服這兩大難題,微地震震源定位發(fā)展了兩類方法,一類是震源成像的方法,另一類是基于微地震事件振幅疊加方法[3-4]。振幅疊加或網(wǎng)格搜索方法,是將儲(chǔ)層一定區(qū)域內(nèi)劃分為大小相同的網(wǎng)格,根據(jù)速度模型計(jì)算每個(gè)網(wǎng)格到所有檢波器的走時(shí),得到所有網(wǎng)格點(diǎn)到全部檢波器的走時(shí)表[5]。根據(jù)走時(shí)表的時(shí)差關(guān)系或利用波形的互相關(guān),提取走時(shí)差進(jìn)行偏移疊加數(shù)據(jù)道波形,從而反演儲(chǔ)層區(qū)域的成像結(jié)果,成像結(jié)果的最大值位置即為震源位置。該方法的優(yōu)點(diǎn)是,無須識(shí)別、拾取初至即可進(jìn)行對(duì)地下目標(biāo)范圍和地震記錄進(jìn)行掃描,完成震源的成像結(jié)果[6]。但通常疊加的數(shù)據(jù)不會(huì)進(jìn)行篩選,導(dǎo)致噪聲干擾數(shù)據(jù)參與運(yùn)算,降低了成像效果。同時(shí),不同的震源機(jī)制(ISO 源、DC 源、CLVD 源)會(huì)輻射不同的波場(chǎng)特征,地面不同位置檢波器觀測(cè)到的P 波和S 波特征也不相同[7]。在疊加成像過程中,低信噪比和信號(hào)特征差異都會(huì)對(duì)成像結(jié)果造成影響。

本文應(yīng)用了一種聚類分析方法消除這種影響。聚類分析又稱為群分析,它是研究(樣品和變量)分類問題的一種多元統(tǒng)計(jì)方法,并對(duì)相似屬性特征的對(duì)象進(jìn)行分組的數(shù)據(jù)挖掘方法[8]。正如統(tǒng)計(jì)學(xué)中樣本之間有不同的分類方法,信號(hào)也可以看作是等間隔、連續(xù)采樣的樣本。通過聚類分析可以定義信號(hào)與信號(hào)之間的距離,距離越小,信號(hào)差異越??；反之,差異越大。根據(jù)不同信號(hào)的特征對(duì)信號(hào)進(jìn)行聚類計(jì)算,篩選出信號(hào)特征接近的數(shù)據(jù)集進(jìn)行成像,改善成像的質(zhì)量。相比于傳統(tǒng)波形互相關(guān)計(jì)算,聚類分析對(duì)地震道噪聲敏感性更低,不易受噪聲影響。波形互相關(guān)計(jì)算需精確提取事件的完整波形才能計(jì)算出較為精確的相似性結(jié)果,而聚類分析只需包含事件本身就可對(duì)其進(jìn)行分類。

1 波形聚類方法

波形聚類方法主要包括數(shù)據(jù)道之間的距離計(jì)算方法和聚類算法。

學(xué)者提出多種距離計(jì)算方法,如歐氏(Euclidean)距離、切比雪夫(Chebychev)距離、蘭氏(Canberra)距離、馬氏(Mahalanobis)距離、夾角余弦(Cosine)距離以及漢明(Hamming)距離等[9]。每種距離計(jì)算適用不同的樣本特點(diǎn),嘗試以上方法對(duì)微地震數(shù)據(jù)進(jìn)行了測(cè)試,發(fā)現(xiàn)歐式距離能獲得不錯(cuò)的聚類結(jié)果。歐氏距離計(jì)算方法如下:

式中,dij為數(shù)據(jù)道Xi與數(shù)據(jù)道Xj的距離；P為數(shù)據(jù)道中采樣點(diǎn)總數(shù)；xia、xja分別為數(shù)據(jù)道Xi與Xj中第a個(gè)采樣點(diǎn)位置的振幅值。

定義dij有如下性質(zhì):dij≥0,dij=dji且dij≤dik+dkj。歐式距離計(jì)算結(jié)果與波形數(shù)據(jù)的量綱有關(guān),因此在計(jì)算前需要對(duì)各道數(shù)據(jù)進(jìn)行歸一化處理。通過式(1)可得距離矩陣D(0):

式(2)中,d11=d22=…=dnn=0,D是一個(gè)對(duì)稱矩陣,在聚類時(shí)只需計(jì)算上三角或下三角形部分即可。根據(jù)D(0)可對(duì)n道數(shù)據(jù)進(jìn)行分類,距離近的數(shù)據(jù)歸為一類,距離遠(yuǎn)的數(shù)據(jù)歸為另一同類。

聚類也有多種方法,如最短距離法、最長(zhǎng)距離法以及重心法等。本文采用最短距離法來說明聚類的過程:

(1) 確定距離的計(jì)算方法。計(jì)算波形數(shù)據(jù)之間的距離,得到距離矩陣D(0),定義類與類之間的距離為Dij,開始時(shí)每道數(shù)據(jù)自成一類,顯然此時(shí)Dij=dij。

(2) 找出D(0)的非對(duì)角線最小元素,設(shè)為dpq,則將波形數(shù)據(jù)道X(p)和X(q)合并成一個(gè)新的類,記為Gr={X(p),X(q)}。

(3) 計(jì)算所有剩余類Gk與新類Gr的距離:

(4) 將D(0)中p、q行及p、q列刪除,然后在D(0)的第一行第一列位置插入Gr,得到一個(gè)n-1維距離矩陣,記為D(1)。

(5) 對(duì)D(1)重復(fù)步驟(2)(3)(4)得到D(2),直到所有的元素歸為一類(距離矩陣中剩下一行一列)時(shí),聚類結(jié)束。

2 聚類分析的應(yīng)用

2.1 正演數(shù)據(jù)聚類

現(xiàn)以合成波形數(shù)據(jù)為例,說明波形聚類的過程。圖1(a)(c)(e)(g)分別為Ricker 子波、Morlet子波、偏移的Morlet 子波以及以上3 種子波疊加的波形。微地震事件在發(fā)震時(shí)頻率很高,傳播過程中衰減嚴(yán)重[10-12],因此分別對(duì)4 種子波疊加高斯白噪聲模擬真實(shí)的信號(hào),如圖1(b)(d)(f)(h)所示。可以看出,疊加噪聲后的合成信號(hào)信噪比很低,波形特征較為接近。

圖1 4 種類型子波與對(duì)應(yīng)加噪后合成波形Fig.1 Four types of wavelets and the corresponding synthesized waveform with noise

對(duì)4 種不同的子波分別疊加3 次高斯白噪聲,如圖2 所示。每次疊加噪聲均為隨機(jī)疊加,因此疊加噪聲后波形各不相同,將12 道合成數(shù)據(jù)隨機(jī)打亂,其中1、5、9 道為Ricker 子波合成波形,2、3、8道為Morlet 子波合成波形,4、6、10 道為偏移Morlet子波合成波形,7、11、12 道為合成子波合成波形。

圖2 隨機(jī)打亂后的合成數(shù)據(jù)剖面Fig.2 Synthetic data traces after random order

聚類結(jié)果如圖3 所示。根據(jù)式(1)計(jì)算各道之間歐式距離,得到距離矩陣,在初始狀態(tài)下每一道各為一類。第1 步,查找各類之間的最小距離,將1、9 合并為一個(gè)類；第2 步,重復(fù)計(jì)算新類和剩余類的距離,將最小距離的5 和1、9 合并為新的類1、9、5。重復(fù)上述2 個(gè)步驟后,2、8、3 與4、10、6 各分為同一類。分類結(jié)果中只有7 沒有與11、12 直接歸為同一類,這是由于7 前面組合的新類較為相似。分類結(jié)果表明,聚類方法基本能將12 道合成波形數(shù)據(jù)正確分類。

圖3 合成數(shù)據(jù)聚類結(jié)果Fig.3 Clustering results of synthetic data

2.2 微地震信號(hào)聚類

將聚類分析方法在事件微地震監(jiān)測(cè)數(shù)據(jù)中進(jìn)行了應(yīng)用。M 井組位于山西省臨汾市沁水縣馬壁鄉(xiāng)(沁水盆地南部),該井組為水平井,水平井段垂直深度在920 ～1 030 m 之間。該地區(qū)主要煤層為3 號(hào)煤,煤層深度由東南向西北方向逐漸變淺,煤層厚度2.4 ～7.5 m,水平井段附近無顯著大斷層發(fā)育[13-14]。根據(jù)聲波測(cè)井曲線計(jì)算,儲(chǔ)層P 波速度為3 300 ～3 500 m/s,地表速度較高,P 波速度為3 900 ～4 200 m/s。煤樣測(cè)試表明,3 號(hào)煤層滲透率在0.025 ～0.029 mD 范圍,整體滲透性較差。

M 井組進(jìn)行了水力壓裂改造,以擴(kuò)大煤儲(chǔ)層的煤層氣產(chǎn)量,同時(shí)地面進(jìn)行了微地震監(jiān)測(cè)。微地震監(jiān)測(cè)觀測(cè)系統(tǒng)地面投影如圖4 所示,其中黃色圖釘表示檢波器位置,紅色圖釘為水平井段起始位置,綠色圖釘為井底,紅色實(shí)線為井軌跡。

圖4 M 井組地面微地震觀測(cè)陣列與井軌跡地面投影Fig.4 Geophone array of well M group and ground projection of well trajectory

M-1 井第5 段水力壓裂地面微地震監(jiān)測(cè)事件A 波形如圖5 所示,數(shù)據(jù)記錄共43 道,時(shí)窗內(nèi)共800 采樣點(diǎn),采樣頻率2 ms,波形顯示時(shí)窗的長(zhǎng)度為1.6 s。從波形圖中可以看出,2、14、25 為明顯的噪聲干擾道。此外,由于地面檢波器分布位置不同,波形特征存在明顯的差異。

圖5 地面監(jiān)測(cè)微地震事件A 波形Fig.5 Waveform of surface microseismic event A

對(duì)圖5 的微地震進(jìn)行聚類分析,結(jié)果如圖6所示。橫軸坐標(biāo)對(duì)應(yīng)圖5 中的數(shù)據(jù)序號(hào),縱坐標(biāo)表示分類距離。由于數(shù)據(jù)較多,為了能清晰地看出聚類結(jié)果,分類距離從左至右依次增大。根據(jù)分類結(jié)果中距離的變化趨勢(shì)可以看出,矩形內(nèi)數(shù)據(jù)道之間的距離差值較小,變化梯度低；矩形框外的數(shù)據(jù)道相對(duì)距離變化較大,因此矩形框內(nèi)的數(shù)據(jù)道可以看成一類事件。選取圖6 中藍(lán)色虛線矩形框內(nèi)的數(shù)據(jù)進(jìn)行疊加定位。保留的數(shù)據(jù)波形都具有相似的波形特征,如較強(qiáng)的P 波或S波。去除數(shù)據(jù)道波形為噪聲道或信號(hào)差異較大的數(shù)據(jù)波形。

圖6 事件A 聚類結(jié)果Fig.6 Clustering result of event A

從圖5 的波形數(shù)據(jù)中提取3 類典型特征的微地震信號(hào),如圖7 所示。

圖7 事件A 中3 種不同信號(hào)特征Fig.7 Three different signal characteristics in event A

Ⅰ類信號(hào)特征具有明顯P 波和S 波,且P 波波動(dòng)明顯；Ⅱ類信號(hào)無明顯P 波,但S 波振幅較強(qiáng)；Ⅲ類信號(hào)具有較強(qiáng)的P 波,而無明顯的S 波。3 類信號(hào)主要由震源機(jī)制與地面采集位置確定。由于基于振幅疊加的微地震事件成像方法與波形特征密切相關(guān),噪聲干擾道以及不同信號(hào)特征會(huì)影響最終的成像結(jié)果,因此在波形疊加過程中,一般采用絕對(duì)值振幅或平方振幅值[15],Ⅰ類和Ⅲ類可避免極性相反造成的振幅相互抵消,但Ⅱ類在P 波振幅疊加中不可避免地帶來振幅削弱的影響。利用波形之間的相關(guān)性,去除了Ⅱ類數(shù)據(jù),最終保留的數(shù)據(jù)道如圖8 所示。

圖8 事件A 聚類后保留數(shù)據(jù)與刪除數(shù)據(jù)波形Fig.8 Retained traces and deleted traces after event A clustering

2.3 信號(hào)聚類對(duì)定位成像的優(yōu)化

分別對(duì)聚類篩選前后的數(shù)據(jù)道進(jìn)行繞射疊加成像。繞射疊加的原理是,通過射線追蹤算法計(jì)算同一個(gè)震源點(diǎn)到每一個(gè)檢波器的時(shí)間偏移量,當(dāng)不同檢波器校準(zhǔn)偏移量為同一發(fā)震時(shí)刻時(shí),振幅疊加必然最大。其定位基本流程是,將目標(biāo)區(qū)域按定位精度劃分成多個(gè)小體元(假設(shè)震源點(diǎn)),計(jì)算體元到每個(gè)檢波器的偏移量,然后沿著偏移量進(jìn)行繞射疊加,直至遍歷所有體元。其成像公式[16]為

式中,N為檢波器個(gè)數(shù)；ui為第i個(gè)檢波器記錄的地震記錄,即原始波形信息,通常為波形振幅的絕對(duì)值或振幅平方；ti(x)為從體元到第i個(gè)檢波器的理論到時(shí)；τ為發(fā)震時(shí)刻。

整個(gè)成像過程是將所有檢波器波形記錄上的振幅按照計(jì)算的走時(shí)進(jìn)行疊加,在震源位置的網(wǎng)格點(diǎn)處疊加能量最強(qiáng)。

分別對(duì)聚類前后的數(shù)據(jù)進(jìn)行繞射疊加成像。由于已知微地震事件所在位置,疊加時(shí)選取600 ms 時(shí)窗長(zhǎng)度進(jìn)行繞射疊加。全部波形數(shù)據(jù)成像結(jié)果如圖9(a)所示,聚類后保留數(shù)據(jù)成像結(jié)果如圖9(b)所示。背景顏色為相似性疊加值,紅色表示高相似性值,黃色、綠色依次降低,藍(lán)色最低。根據(jù)疊加成像的原理,一個(gè)震源點(diǎn)釋放的地震波信號(hào)由地面各站點(diǎn)接收后雖然有時(shí)差、振幅差或畸變[17],但是具有內(nèi)在相關(guān)性。繞射疊加成像結(jié)果中最大能量值位置即是震源發(fā)生的位置[18-19]。聚類分析的本質(zhì)是增加事件的成像值,改善定位結(jié)果。

圖9 事件A 成像結(jié)果對(duì)比Fig.9 Comparison of imaging results of event clustering

對(duì)比聚類前后可以看出,圖9(b)成像結(jié)果相比圖9(a)更加聚焦,紅色區(qū)域過渡更集中,并且聚類保留數(shù)據(jù)的相似性值域明顯提高(全部數(shù)據(jù)疊加相似性值域?yàn)?.33 ～0.47,篩選后數(shù)據(jù)疊加相似性值域?yàn)?.40 ～0.54),從而提高了定位結(jié)果的可信度。在成像結(jié)果中,只需提取對(duì)應(yīng)的最大相似性值位置或拓?fù)浣Y(jié)構(gòu)的中心,便可確定震源點(diǎn)[20]。因此,聚類篩選后的數(shù)據(jù)在成像效果和相似性值域均得到了明顯改善,有利于微地震事件的進(jìn)一步分析。

事件B 與事件C 的微地震事件成像結(jié)果如圖10、圖11 所示?？梢钥闯?聚類分析方法改善了成像的質(zhì)量,提高了相似性值域范圍。

圖10 事件B 成像結(jié)果對(duì)比Fig.10 Comparison of imaging results of event B

圖11 事件C 成像結(jié)果對(duì)比Fig.11 Comparison of imaging results of event C

3 結(jié) 論

本文將統(tǒng)計(jì)學(xué)中的聚類分析方法在微地震成像中進(jìn)行了應(yīng)用,將每道信號(hào)看作一個(gè)樣本,根據(jù)信號(hào)的特征進(jìn)行分類,再通過各道相關(guān)性去除P波振幅較弱的數(shù)據(jù)道。

聚類算法在實(shí)際應(yīng)用中計(jì)算效率高,識(shí)別較為準(zhǔn)確,能夠快速識(shí)別和提取距離相近的數(shù)據(jù)道,無需人工初至拾取進(jìn)行微地震事件成像,適合于微地震監(jiān)測(cè)實(shí)時(shí)定位處理。

篩選后的地震數(shù)據(jù)成像結(jié)果明顯改善,且疊加的相似性值域明顯增大,提高了定位的可靠度。

致 謝

感謝深聽（北京） 科技有限責(zé)任公司DeepListen 解釋軟件，為本文提供了數(shù)據(jù)處理的基礎(chǔ)。