煤礦井下全空間地震波時頻域極化分析及其試驗研究

王 勃，曾林峰，張 衍，劉盛東，章 俊，陳泓云

(1.中國礦業大學 深部巖土力學與地下工程國家重點實驗室，江蘇 徐州 221116；2.中國礦業大學 資源與地球科學學院，江蘇 徐州 221116；3.中國礦業大學 力學與土木工程學院，江蘇 徐州 221116；4.深地科學與工程云龍湖實驗室，江蘇 徐州 221116)

“雙碳”目標下我國能源資源稟賦特征及當前復雜國際能源形勢，決定了煤炭在能源中的基礎和兜底保障作用。煤炭綠色開采、智能精準開采等對煤礦安全高效開采地質保障系統提出了更高的要求，礦井地質透明化是當前礦井地質保障系統發展的努力方向。礦井地震勘探具有不受復雜地面條件影響、距離目標體近、地震波能量和高頻成分衰減少、分辨率高、探采對比易驗證等優點，進而被廣泛運用于地質精細探測及透明地質建模。

不同于地面半空間地震勘探，礦井地震屬于全空間條件下勘探。震源在井下激發后，地震波向四周傳播，檢波點會接收到全方位的地震信號。單分量信號難以確定地震波的方向，雙分量地震信號僅能算出二維平面上地震波傳播方向，三維空間條件下準確判定地震波傳播方向則依賴三分量信號的矢量特性，并通過極化分析實現。極化分析在地震信號處理中主要用于波場分離，根據各類地震波的極化屬性差異設置濾波器來分離提取有效波。極化濾波方法在天然地震、石油地震勘探領域研究較多，DIALLO 等對多分量地震數據進行小波變換，將自適應瞬時極化濾波引入時頻域，進行面波壓制及轉換波分離。KULESH 等提出了基于自適應協方差矩陣的小波域時頻極化分析方法，在時頻域實現了波場分離。PINNEGAR采用S變換在時頻域計算三分量地震信號極化參數，設置極化率濾波器壓制極化橢圓性地震波。程冰潔等在小波域進行能量分類約束從而實現極化濾波。

在煤炭領域，張平松等在時間域構建協方差矩陣，采用極化濾波方法實現了巷道超前探測多波有效分離。王勃利用三分量信號的矢量特征，提出了一種全空間條件下集波場分離、偏移成像于一體的極化偏移方法。胡澤安等在時間域對礦井槽波地震數據進行了極化分析，實現了波場分離與噪聲壓制。金丹等利用槽波信號與干擾波在偏振度上的差異，改進頻率域的極化濾波權函數，將頻率域極化濾波用于槽波記錄的噪聲壓制。馮磊等對二分量槽波數據進行S變換，采用時頻域自適應協方差矩陣極化濾波方法實現了兩類槽波分離。劉盛東等利用三分量地震記錄，通過時窗自適應的極化分析方法獲取極化率，提取了縱橫波及勒夫型槽波等線性極化波。

上述方法主要用于地面半空間或煤礦井下全空間的極化濾波，針對礦井全空間條件下三分量地震波傳播方向研究較少，特別在礦井地震近場勘探多類型波場混疊條件下。筆者提出了一種礦井全空間三分量地震波時頻域極化分析方法，首先通過時間域混疊合成信號驗證地震波優勢方位角、傾角準確度，然后開展全空間條件下中心激發-全方位接收的三維三分量數值模擬實驗對比縱波、橫波定向精度，最后通過礦井地震勘探常用的反射波現場試驗驗證極化參數求解的可靠性。

1 全空間地震波時頻域極化分析方法

礦井地震屬于全空間條件下近場勘探，多類地震波混疊，需從時間、頻率2個維度聯合計算分析。為準確計算地震波傳播方向，筆者首先構建基于基本小波函數形態不固定的廣義S變換的時頻域復協方差矩陣，然后求解三分量地震信號的優勢能量方向。同時，為進一步驗證方向計算的準確性，開展已知震源點聚焦定位對比研究，具體理論方法如下。

1.1 基于廣義S變換的時頻域復協方差矩陣構建

設時間域信號()的傅里葉變換為

(1)

式中，j為虛數單位；為頻率；為時間。

對時間序列()乘以一個窗函數()得

(2)

設()為歸一化的高斯窗，且利用參數控制其時窗寬度，利用參數控制其時窗位置，則有

(3)

將式(3)代入式(2)，可得到對時間序列()在時刻加高斯窗的譜為

(4)

將窗寬控制參數設置成與頻率成反比的關系，以此讓高斯窗的寬度自適應于頻率，則有

(5)

和為控制時窗寬度變化的2個參數，進一步獲得廣義S變換的表達式為

(6)

通過廣義S變換計算時頻譜，再利用Hilbert變換構建時頻譜的解析信號。在時頻域內時刻，頻率處的復協方差矩陣(,)可描述為

(7)

矩陣中各元素定義為

(8)

式中，|SC(,)|和|SC(,)|為2組解析信號的瞬時振幅；和為瞬時頻率；arg為瞬時相位；(,)為均值，其定義為

(,=,,)

(9)

式中，R(·)為復數的實部，函數sin()為辛格函數，其定義為

(10)

(,)為自適應時窗長度，其定義為

(11)

式中，為整數，是一個經驗參數，用于刻畫不同極化屬性，取較大值時可刻畫三維復雜極化屬性,一般取1或2。

當3個分量時頻譜中對應時頻點的瞬時頻率相等時，(,)=(,)=(,)，(,)可簡化為2π(,)，復協方差矩陣可以化簡為

(,)=|SC(,)||SC(,)|×

cos[arg SC(,)-arg SC(,)]=

(12)

式中，為復共軛。

1.2 極化參數計算

地震波矢量特征可以用極化參數描述，這些參數可以通過在時頻域中求取復協方差矩陣的最大特征值及其對應的歸一化特征向量(,,)獲得。

優勢極化方位角計算公式為

()=arctan[R()R()]

(13)

其中-90°≤()≤90°，為極化主軸在面的投影與軸的夾角，當極化主軸偏向軸正方向時，()>0°；當極化主軸偏向軸負方向時，()<0°。

優勢極化傾角計算公式為

(14)

其中，-90°≤()≤90°，為極化主軸與面的夾角，當極化主軸偏向軸正方向時，()>0°；當極化主軸偏向軸負方向時，()<0°。

1.3 聚焦定位方法

與微震定位不同，地震勘探震源點位置已知，為了驗證上述極化參數準確性，利用檢波點位置及其極化傾角、方位角反向求解震源點位置，其計算方法如下：設檢波點的空間坐標和主極化軸歸一化矢量分別為=(,,)和=(,,)，設空間任意一點=(,,)。向量與向量的向量積代表向量和向量共起點的情況下所構成平行四邊形的面積。對該向量積除以空間矢量的模可獲得點到空間矢量的最短距離：

(15)

(16)

求取的最小值所對應的(,,)為震源點在三維空間內的位置。

2 三維數值模擬及算法驗證

2.1 合成信號及極化特征分析

為了驗證極化分析對三分量信號在時頻域上區分多類型波和計算極化參數的準確性，開展正弦波合成信號分析實驗。合成信號由A,B,C,D,E,F共6個信號組成，合成信號具體參數見表1。其中B，C信號在時間域混疊，E，F信號在時間域混疊，在時間域可分為4段信號，如圖1所示。其中每組信號均為501個采樣點，采樣頻率為2.5 kHz，對合成信號進行時頻域極化分析，結果如圖2所示。時頻譜上可清晰區分A，B，C，D，E，F信號，圖2(a)顯示的6個信號的方位角與理論方位角一致，正、負方位角均無誤差；圖2(b)顯示的6個信號的傾角與理論傾角一致，正、負傾角均無誤差。

表1 合成信號基本參數

圖1 合成三分量信號

圖2 合成三分量信號方位角及傾角分布

2.2 正演模擬信號極化特征分析及聚焦定位驗證

為研究全空間條件下地震波矢量特征，建立三維數值模型。模型中心為震源點，在震源點的左側、右側、前方、后方、頂部及底部布置檢波器，形成中心激發全方位接收的三維三分量地震觀測系統，如圖3所示。三維模型在,,方向的大小分別為400 m×400 m×400 m，震源點位于(0，0，0)原點處，三分量檢波器在,,面上以原點為中心直徑為200 m的圓形測線上以15°間隔布置，累計66個檢波器。模型在,,方向上進行網格化，網格間距均為0.3 m。模擬采用主頻為125 Hz的零相位雷克子波，采樣頻率為2.5 kHz。模型添加PML吸收邊界，采用三維時空域高階有限差分法進行數值模擬。

圖3 三維觀測系統

在數值模擬地震數據中選擇一道(第64道)數據進行時頻域極化分析，通過對比理論、計算的極化參數驗證上述方法效果。利用直達波的理論到達時間、震源主頻確定直達波的時頻范圍，在圖4中可分析直達縱波、橫波的極化參數，從圖5可見直達縱波、橫波質點振動軌跡。圖4為第64道數據求得的極化方位角、極化傾角。由圖4(a)可看出，在時間0.03 s 和0.06 s附近各有一團能量，分別為縱波、橫波能量。直達縱波、橫波的方位角均為3.3°，與理論值0°存在部分偏差。由圖4(b)可看出，時間0.03 s附近的直達縱波傾角為45°，與理論值45°無偏差。時間0.06 s附近處的直達橫波傾角為-45°，與理論值-45°無偏差。計算其他65道地震波時頻域極化參數，并利用直達縱波的極化信息對震源點進行反向聚焦定位，結果如圖6所示，在200 m直徑范圍內，求解震源點與已知震源點的直線距離誤差僅為 0.548 6 m。

圖4 計算所得極化方位角與極化傾角

圖5 第64道三分量信號質點振動軌跡

圖6 檢波點反向聚焦定位震源點(右上角為原點位置局部放大)

3 現場試驗

實際礦井地震勘探過程中，除利用透射波之外，還常采用反射波，特別是煤層條件下反射槽波勘探是井下通用方法。為此，開展了現場實測并針對反射槽波信號進一步驗證方法有效性。

安徽某礦1034工作面位于三采區深部，煤層頂板標高-603.8～-466.0 m，煤厚3.0～5.2 m，平均3.9 m，煤層傾角5°～18°，平均10°，頂板為灰白色中粒砂巖，底板為粉砂巖。1034回風巷外幫發育有F13正斷層，與回風巷相距57～108 m，走向NE，傾向NW，傾角60°～70°，落差100～200 m，其在礦區內延展長度5.20 km，利用地面、井下鉆探進行斷層探查，斷層控制程度可靠。以R1檢波點為原點，每隔10 m布置1個檢波點，震源點設置在132.5 m處，三分量地震勘探觀測系統如圖7所示，S1炮激發R1～R27道接收的三分量地震記錄如圖8所示。

圖7 某礦1034工作面地震觀測系統

圖8 三分量地震記錄

觀察三分量地震記錄，圖中紅圈處(R1～R8道記錄)存在明顯反射槽波信號。選取R4道進行分析，S1炮、F13斷層與R4檢波點的反射路徑距離約213.28 m，R4檢波點接收的反射槽波理論方位角約為27°，傾角則等同于煤層傾角，約10°。S1炮激發、R4三分量檢波器接收的地震信號進行時頻域極化分析，時頻域方位角及時頻域傾角如圖9所示。根據已知斷層位置計算反射槽波理論傳播路徑，時間0.18～0.25 s、頻率120～180 Hz區域為反射槽波。圖9(a)中該區域的方位角計算結果可見，反射槽波的方位角較穩定，約為28°，分析其主要為Rayleigh型反射槽波，與理論方位角偏差1°。圖9(b)顯示反射槽波的傾角為10°～12°，與實際煤層10°傾角吻合。

圖9 R4道數據時頻域極化方位角傾角分布

4 結 論

(1)對三分量合成信號進行時頻域極化分析，得到信號的時頻位置與理論時頻位置相對應，求解的方位角、傾角與理論值吻合；在時間和頻率2個維度上，時頻域極化分析可以精準確定混疊情況下多類型信號極化方向。

(2)針對中心激發-全方位接收的三維時空域高階三分量模擬信號，直達波的方位角、傾角誤差分別為3.3°，0°；在200 m直徑范圍內，根據檢波點位置及地震波方位角、傾角聚焦定位的震源點與已知震源點的距離誤差僅為0.548 6 m，驗證了全空間條件下三分量地震波的時頻域極化參數計算的精度。

(3)針對礦井實測數據進行廣義S變換的時頻域極化分析處理，反射槽波的時頻域方位角、傾角參數與理論值一致，證實了礦井全空間三分量地震波時頻域極化分析方法的有效性。

(4)全空間條件下時頻域極化分析方法確定的地震波傳播方向可為礦井地震勘探精細成像提供基礎性支撐，并為井下地質構造的透射波層析成像、反射成像等提供基礎約束條件。