井工煤礦全波形反演的初始模型構建探究

梁 旭，徐 裴，白喜成，李志勇，孫瑞雪，聶天文

(1.陜西小保當礦業有限公司，陜西省榆林市，719399；2.北京中礦大地地球探測工程技術有限公司，北京市海淀區，100083)

0 引言

隨著煤炭開采深度的加大，地下構造變得越來越復雜，煤礦開采對煤層頂底板成像精度的要求越來越高[1]。而全波形反演技術是基于廣義反演理論下對地震數據的全部波場信息的反演擬合，是一種高精度高分辨率的物探方法，對于透明工作面構建有指導作用。井工煤礦工作面的全波形反演中，良好的初始模型是關鍵，準確的速度場又會降低波形反演的周跳現象，從而提升全波形反演的精度。一般可以通過鉆井數據建立簡單層狀速度模型來進行反演計算，但是在對于三維方向速度精度要求較高的運算中，就必須要對模型校準調整，從而全波形反演準確。針對此問題，前人做了大量研究，包含走時層析成像[2-3]、偏移速度分析[4-5]和不同空間域反演[6-7]等領域。石油工業上主要以走時層析成像(CT)作為全波形反演的主要初始模型來源，而在煤礦工作面勘探中，還未有相關研究的應用探究，只限于地震CT的單一應用。JACKSON M J等[8]將三維地震CT技術引入到煤層的構造探測中；SALNIKOV A S等[9]提出透射波地震數據成像方法，應用在煤礦透射勘探中；彭蘇萍等[10]首次采用地震CT技術探測放頂煤工作面構造，表明成像結果可以凸顯綜放面內的隱伏斷層；胡澤安等[11]使用2.5維地震速度層析成像方法，通過加入高程參數，消除了雙巷道高程的影響，使得射線追蹤的旅行路徑更加準確，提升了速度反演結果。

本文主要采用三維地震CT技術的地震透射數據反演，提升“圍巖-煤層-圍巖”的速度模型，從而探究幫助改進全波形反演計算的方法。首先基于程函方程求解時間場[12-13]、梯度射線追蹤法獲得走時[14]，與實際初至走時做殘差得到走時數據不適定變量，構建走時對于慢度的敏感矩陣，采用帶阻尼的最小二乘反演算法[15-16]，校正三維速度場(慢度)，從而確定作為波形反演的初始大概速度輪廓，為改善全波形反演效果做輔助。

1 走時正、反演技術

1.1 正演計算

1.1.1 有限差分格式

地震波的瞬間走時可以通過程函方程計算，如式(1)所示，描述了走時梯度和地下介質慢度的關系，然后根據VIDALE J E[4]推導出的有限差分公式遞推出地下網格上的每一個點的走時。

(1)

式中：t——地震走時；

x，y，z——卡尼吉亞坐標軸不同方向；

s——慢度(速度的倒數)。

為了求解式(1)，將網格上的點分成3種類型，如圖1所示[4]。在圖1(a)中，若已知0～6號共7個角點的走時，可以利用式(1)的有限差分形式獲得式(2)，從而求解得到角點7的走時。

圖1 計算三維均一介質的三種基本點模板

(2)

式中：h——網格間距；

ti——第i個點的走時。

圖1(b)是第2種網格點類型，被用于式(2)沒有計算出來的點。5號點的走時可以通過式(3)由點0～4推導得到。其中式中的慢度s是2個立方體的中心慢度值，也就是2個立方體平均慢度的平均值，而每個立方體的平均慢度是該立方體的8個角點慢度的平均值。圖1(c)是第三類點，被用于計算剩余未被計算出來的網格點走時。圖中5號點的走時可以通過式(4)由0～4號點計算得到。s是圖中所有網格點上的慢度的平均值。

由以上式子先計算出圍繞震源點的26個網格點的走時，然后按照快速行進法推進，沿波前面計算其余網格點的走時。

1.1.2 迎風格式與快速行進方法

SETHIAN J A[13]提出程函方程迎風解的快速行進法加快了計算效率。首先需要推導程函方程公式(1)的迎風格式近似于其多維度的走時場梯度值，如式(5)所示。

(5)

式中：D-、D+——正、反向算子，代表在x、y、z各方向上的導數；

sijk——網格點(i，j，k)的慢度。

簡化后，式(5)可以寫成以下形式：

(6)

快速行進法的核心思想是構建程函方程的迎風格式來解方程求走時場。本質上，式(6)意味著地震波信息從較小的走時傳播到更大的走時。因此，快速行進法就是建立一個解方程的方案，從最小的走時向外擴散。該算法通過構建一個波前窄帶讓計算加快，也就是先考慮當下的波前窄帶里的點，利用迎風格式去掃描，試圖推進這個窄帶向前，凍結當下的點，然后將新的點帶入窄帶結構[5]，如圖2所示。

圖2 迎風格式示意

這里的關鍵在于更新窄帶結構的網格點的選擇方式。具體的算法步驟如下所述。

(1)標記初始條件下的點為“可接受的值”，然后標記波前的窄帶點為“近點”，最后，標記其他點為“遠點”。

(2)開始循環，讓在“近點”中的點取時間最小值的點作為試驗點。

(3)將試驗點填入到“可接受的值”群里，將之從“近點”群里移除。

(4)將所有試驗點的鄰點標記為“不可接受點”，如果有鄰點在“遠點”群中，將這個鄰點從“遠點”中移除，并將之添加到“近點”中。

(5)按照迎風方程(6)重新計算所有窄帶鄰點的時間。

(6)返回到步驟(2)。

這種算法的合理性在于重新計算迎風區的相鄰點的走時無法產生比接受的點更小的值了。所以才可以向外推進走時解，總是選擇最小的走時作為窄帶中的網格點值，并且重新調整鄰近點。

通過建立一個三維簡單層狀模型來測試本文正演算法，得到走時等值線，如圖3所示，模型大小x×y×z=100 m×100 m×80 m，震源設置在x=50 m，y=2 m，z=2 m的位置，深度z為0～30 m地層中速度為1 800 m/s，30～80 m地層中速度為5 300 m/s。

圖3 簡單層狀結構速度模型三方向剖面及旅行時等值線

在模擬走時場時，震源作為走時起始位置，本程序中只做了簡單設置為0的處理，也只做了二階的偏微分方程求解，因此在走時傳播時，存在誤差，導致走時等值線有一定的抖動現象，在后續研究中，會改進算法，提升精度，本文重點對初級版本的程序做了應用探究。

1.1.3 射線追蹤

在走時場t(x，y，z)中，走時場的方向余弦為l=(cosα，cosβ，cosγ)，則走時場在方向l上的方向導數式為：

(7)

(8)

由式(8)可知，當走時場的方向向量和梯度向量的方向一致時，方向導數取最大值，與梯度的模相等，這說明走時變化率最大方向就是梯度的方向，由費馬最短路徑原理知，地震波傳播的射線方向應該沿走時場的梯度方向變化，那么每條射線的坐標更新式如式(9)～式(11)所示。

式中：(xcurrent，ycurrent，zcurrent)——射線起始點；

(xnext，ynext，znext)——射線下一個結點；

dx，dy，dz——射線追蹤的不同方向上的步長。

1.2 反演計算

使用線性最小二乘法[6]解決走時反演問題。走時和慢度之間的關系如式(12)所示，它是非線性關系，而做反演時，需要把此關系式轉換為泰勒級數展開，然后通過迭代的方式求解方程。考慮一個離散的走時反演問題，此處用絕對走時。慢度模型包含離散的慢度單元，而數據包含有限個走時觀測值。走時t和慢度s的關系可以表示為以下的矩陣方程：

t=Ls

(12)

式中：L——一個穿過每一個慢度單元的射線路徑的長度矩陣。

由于射線路徑定義矩陣L，取決于慢度s的非線性關系，式(12)的解可以通過一個初始模型開始，迭代擬合到一個正確的慢度模型和射線路徑。這往往被當作一個非線性的走時反演問題。為了提升反演效率，在有限的計算空間下，采取將非線性問題線性化的方式去解決這個問題。走時和慢度的方程描述如下：

t=Fs

(13)

式中：F——非線性方程，代表走時t與慢度s的計算關系。

可以將F從初始慢度模型s0擴展開來，展開式如下：

(14)

式中：D——走時對慢度的偏導數；

δs——慢度的變化值矩陣。

忽略慢度變化值的高階項，重組方程，可以得到：

t-Fs0≈Dδs

(15)

式(15)左邊是走時剩余值，右邊是一個簡單的矩陣乘積項。矩陣D可以通過射線追蹤后得到的有效射線長度的累計和求取獲得。式(15)可以通過標準的最小二乘方法求解。對式(15)做迭代，得到慢度校正量δs，更新當前慢度模型。對模型光滑度的限制也被引入到模型校正計算過程中。一個直接的方法就是應用限制到模型光滑度上。在式(15)兩側加上Ds0，得到式(16)：

t-Fs0+Ds0≈Ds

(16)

式(16)和式(15)相同，卻更容易引入反演慢度的光滑度。使用式(15)做迭代收斂到一個慢度值，使得走時剩余值r最小。將式(16)修改，包含模型的拉普拉斯差分。光滑度限制方程與走時方程相關，求解一個與走時觀測值一致的光滑模型。推導之后的迭代方程如式(17)所示：

(17)

式中：si——第i次迭代的模型解；

si-1——前一次迭代的解，雅可比矩陣D可以使用si-1。

加權因子λ控制走時擬合和模型光滑的平衡，Δ是構建的慢度圖像的拉普拉斯矩陣，添加一個方程系統到原始的走時反演問題中。

式(18)為慢度圖像的拉普拉斯矩陣構建，相當于Sx，y，z這個中心點慢度分別離散成周邊6個點的慢度值來表示，從而起到光滑圖像的效果。式(19)表示下一個點距離上一個點的距離累加，從而計算隨著網格行進傳遞的行走距離，用來表示不同射線對應的雅可比矩陣。

基于以上的公式可以構建靈敏度矩陣A：

(20)

選用最小二乘法LSQR解上述反演方程：

已知靈敏度矩陣和走時殘差，得到模型校正量ds。

si=si-1+ds

(23)

式中：si——更新后的慢度模型，從而不斷迭代更新速度模型。

2 合成數據測試

建立理論速度模型，如圖4(a)所示。以圖5作為觀測系統，包括23個地面檢波器，并設置3個震源，觀測系統中央位置深度分別為4 500、4 100、3 800 m，它們的北方向為3 000 m，東方向為3 000 m處。圖4(a)為速度模型的三維示意圖，北-東方向為北-深度剖面切片向東方向的等值延展，不同顏色代表不同巖性，由不同速度值構成，其中，紅色部分為4 800 m/s，深紅色部分為5 500 m/s，藍色部分為3 000、2 500、2 000 m/s(由深部到淺部依次)。用上述正演方法計算相應的初至走時。將計算的數據當作真實值，再平滑，作為初始模型，利用以上反演方法初始擬合真實數據，不斷迭代，縮小擬合誤差。本次實驗經過20次迭代收斂到合理的誤差，如圖4(f)所示，反演的速度模型結果如圖4(d)所示，該北視圖剖面顯示傾斜斷層，與真實模型近似，而真實走時和計算走時擬合效果良好，如圖4(e)所示。結果說明該方法可行，只是在震源處存在假異常現象，是由于觀測系統局限性導致，本研究在應用時將震源異常做扣除和梯度均衡處理。

圖4 復雜三維構造理論數據正反演結果

圖5 理論模型觀測系統

3 小保當煤礦應用案例

基于上述對透射地震數據的三維走時正反演框架，嘗試引入到煤礦工作面的雙巷道透射勘探中，通過拾取透射炮集的縱橫波折射初至，反演擬合初至走時，從而勾勒出初始的三維速度輪廓。最終使用不同初始模型對比測試全波形反演的梯度形態來說明應用效果。

3.1 工作面概況

以小保當二號煤礦工作面132204為例，該工作面寬450 m，長1 000 m，煤層頂底板情況見表1，頂底板及圍巖以砂巖和泥巖為主。地質資料顯示，煤層平均厚度為2～3 m，工作面起伏較平緩，沒有明顯大構造揭露。

表1 煤層頂底板情況

工作面透射觀測系統震源和檢波器布置如圖6所示，震源和檢波器均大體在煤層中間層位，分別沿輔運巷和運輸巷布設，檢波器間距為5 m，輔運巷布設201個檢波器，運輸巷布設201個檢波器，炮間距為10 m，輔運巷100炮，運輸巷100炮。檢波器為三分量采集形式，震源為炸藥。

圖6 相對坐標系統中的震源和檢波器位置

3.2 實際地震數據的分析計算

以132204工作面的運輸巷放炮輔運巷接收的第15炮三分量地震記錄為例，依據其透射最短走時，估計圍巖-煤層速度見表2，得到彩色走時曲線如圖7所示。煤層直達波在本勘探數據中不明顯，初至走時是煤層到圍巖的折射波走時，所以本文主要利用折射走時反演計算，而得到圍巖的橫縱波速度響應，從而刻畫煤層-圍巖界面，提高初始模型準確度。圍巖橫縱波速度可以根據走時曲線擬合地震記錄相位大體位置，初步估計為1 700 m/s和3 800 m/s。該工作面地震數據第15炮運輸巷到輔運巷透射槽波第121道的頻散譜如圖8所示。由圖8可知，通過比對頻譜拾取點與地震記錄對應相位，可以判定頻散曲線的埃里相速度位置大體在923 m/s，地質資料顯示本工作面屬于延安2-2型煤層，與常規槽波速度范圍相符[16]，其常規橫波速度在1 000～1 700 m/s范圍，通過三維彈性波場正演模擬計算，可以將煤層橫波速度大體定為1 000 m/s，而根據所產生的槽波頻譜可初步判定該組數據槽波只有勒夫型槽波一種類型，所以可以簡單判定煤層的縱波速度大于圍巖的橫波速度，以常規縱橫波比1.7為參考值，估計煤層縱波速度大體為2 200 m/s。最終設定初始速度模型參數初步見表2。

表2 初始速度模型參數

注：圖中的走時曲線中，粉色表示縱波圍巖折射初至，綠色表示煤層直達縱波，藍色表示橫波圍巖折射到時，黃色表示煤層直達橫波。

圖8 實際透射數據的某一道頻散譜

受透射法探測距離所限，工作面實際觀測系統沿x方向分為兩段0～350 m和320～1 000 m進行布設，分別建立x=350 m，y=451 m，z=28 m和x=680 m，y=451 m，z=28 m的3層模型，在深度(z)10～13 m處設置煤層，大體檢波器布置到煤層中間層位，道間距為5 m，炮間距為10 m，有效炮共有192炮。按照表2為走時反演賦初值建立三層模型，利用高斯平滑將三層模型光滑，提升射線追蹤效果。拾取每一炮的走時信息，如圖7中粉色和藍色走時曲線標定處，則所有走時拾取點如圖9所示。將拾取橫縱波的走時分別求解程函方程的反演計算系統中，縱波初至經過迭代100次，橫波初至經過迭代50次，走時擬合剩余殘差均方差曲線如圖10所示，第15炮的合成走時和實際走時擬合情況如圖11所示，得到最佳速度模型，如圖12所示。

圖9 縱波和橫波初至拾取

圖10 走時擬合剩余殘差曲線

圖11 第15炮單炮走時曲線(實際與合成)

圖12 縱波和橫波速度走時反演速度結果

將簡單三層速度模型(表2)和反演得到速度模型(圖12)分別輸入到三維全波形反演系統中，由實際數據提取的地震子波做震源激發正演波場記錄，與實際波場記錄做伴隨互相關得到相應模型梯度，如圖13所示(在z=14 m、y=229 m、x=692 m處取切片)，比較x-y-z三方向切片，走時反演改進后的速度模型使得合成數據和實際數據的相關性更好，梯度結果趨于穩定，體現出煤層-圍巖界面的連續變化，證明本方法可有效提升煤礦工作面中全波形反演效果。

圖13 不同初始速度對應的全波形反演梯度

3 結論

(1)探究了對于井工煤礦工作面中三維全波形反演初始模型校正的方法，旨在緩解鉆井資料建立速度模型在波形反演中的不適定性，經過三維透射走時成像大體刻畫出初始速度場的輪廓，從而輔助波動反演成像的計算效果，在煤礦領域初次實踐了走時CT結合波動反演的技術策略。

(2)通過結合程函方程求解時間場和反演解最小二乘目標函數的思路，初步提升三維速度場，對復雜模型的正反演理論測試驗證該算法可行。

(3)實際數據結果顯示，初至走時擬合而反推的速度獲得的波場與簡單一維層狀速度模型得到的波場，分別與實際波場伴隨互相關，對比其對應的梯度，前者梯度變換更連續穩定，更符合煤礦工作面的地質意義。因此，可知本研究結果可提升煤礦井下工作面中全波形反演效果的價值，可以促進煤巖分界面的刻畫。