劉 超，劉 鵬，張慧峰，白 坤，李志勇，潘文勇

(1.陜西小保當?shù)V業(yè)有限公司，陜西 神木 719399；2.北京中礦大地地球探測工程技術有限公司，北京 100038；3.中國科學院地質與地球物理研究所 油氣資源研究院重點實驗室，北京 100029)

為實現(xiàn)煤炭資源高效、安全、智能化開采，提高資源利用率，獲取精準的三維地下地質模型尤為重要，是采煤工作面信息透明化管理的重要基礎[1-4]。然而，常規(guī)地震探測技術無法滿足當前煤礦開采的成像精度和分辨率需求，亟需全新的地震反演與成像技術。全波形反演(Full-waveform inversion，即FWI)是近年來興起的“新一代”高分辨率地震成像技術，通過對彈性波動方程的數(shù)值求解，充分利用全波形信息進行反演成像，可構建地下復雜構造的彈性參數(shù)模型，實現(xiàn)三維地質模型的精細刻畫[5-8]。地震全波形反演技術在油氣勘探領域已獲得大量成功應用，但在煤炭行業(yè)應用較少，國外僅基于2D地質模型做理論數(shù)據(jù)的FWI，幫助煤礦勘探[9]，國內研究機構只針對二維工作面進行合成數(shù)據(jù)的FWI擬合[10]，受限于精度要求高，計算成本過大的實際操作問題，至今國內外還沒有三維FWI在煤礦上的實際應用研究。而為了精細刻畫井下煤礦模型，對陷落柱、斷層以及煤巖分界面清晰成像，使用三維FWI算法可同時降低計算成本，并確保反演結果精確度。

本研究，首先需要通過求解三維彈性波動方程進行正演模擬獲得模擬數(shù)據(jù)，在計算梯度時還需要求解伴隨波場。常規(guī)的地震波正演模擬方法包括有限差分法、有限元法、偽譜法和邊界元法等[11-13]，與這些方法相比，譜元法(Spectral-Element Method，即SEM法)，具有偽譜法高精度優(yōu)點，同時可靈活處理不規(guī)則邊界，非常適合GPU并行計算。本研究采用SEM法模擬地震波在復雜三維彈性介質中的傳播[14]，使用CPU/GPU異構并行技術提高計算效率。

以陜西小保當煤礦為例，采用三維彈性波全波形反演技術，對礦區(qū)地下三維速度結構進行高分辨反演成像，預測煤層頂?shù)捉缑?，實現(xiàn)采煤工作面的透明化，達到整個礦區(qū)的信息透明化管理和安全高效采煤。根據(jù)煤礦巷道勘探的特點，發(fā)展了多尺度、多目標函數(shù)和多波反演方法與技術。由于地震數(shù)據(jù)缺失低頻信息，F(xiàn)WI具有非常強的非線性，容易陷入局部最小值，即周期跳躍問題。采用多尺度反演策略[15]，首先使用低頻信息反演速度的長波長部分，然后使用高頻信息反演速度結構短波長部分[15]；傳統(tǒng)的波形差目標函數(shù)非常容易受周期跳躍問題的影響，發(fā)展多目標函數(shù)反演策略，依次使用互相關旅行時[16]、包絡[17]和波形差目標函數(shù)進行反演，可有效克服周期跳躍問題；在煤礦巷道采集的地震數(shù)據(jù)中，包含P波、S波和槽波，可使用P波和S波反演背景速度結構，槽波在煤層的頂?shù)捉缑娼涍^多次全反射，在煤層中具有較強的穿透性，對煤層頂?shù)捉缑娉上裼兄匾饔谩２捎枚嗖ǚ囱莶呗?，可更好地構建煤層的三維空間結構。

為實現(xiàn)三維彈性波FWI在實際地震數(shù)據(jù)處理中的成功應用，地震數(shù)據(jù)預處理至關重要，對小保當?shù)V區(qū)采集的實際地震數(shù)據(jù)進行了一些機預處理，具體包括：重新采樣、去除壞道、帶通濾波和振幅歸一化等。使用已知的測井數(shù)據(jù)(包括煤層頂?shù)捉缑嫖恢煤兔簩雍穸?作為FWI的先驗信息進行正則化約束，可減弱反演的多解性。從每炮地震數(shù)據(jù)的直達波中提取子波函數(shù)，作為正演模擬的震源函數(shù)。在反演中，對梯度進行異常值去除和三維高斯平滑，降低反演的病態(tài)性。最終，獲得高分辨率的三維彈性速度結構模型，預測煤層的三維空間展布和頂?shù)捉缑?，使用測井數(shù)據(jù)驗證反演結果的準確性。

1 三維彈性波全波形反演基本理論

在FWI中，傳統(tǒng)波形差目標函數(shù)測量了觀測數(shù)據(jù)di與模擬數(shù)據(jù)ui之間的直接差值，通過求取目標函數(shù)的最小值，使用最優(yōu)化算法對模型參數(shù)進行迭代反演，最終使得模擬數(shù)據(jù)與觀測數(shù)據(jù)達到最佳匹配，反演的模型與地下真實模型最大程度接近。波形差目標函數(shù)可以寫為[18]：

式中，i，j，k，l∈[1，2，3]或[x，y，z]；ui為位移場；ρ為密度式中，fi為震源式中，cijkl和ekl分別為彈性剛度張量和應變張量。根據(jù)Voigt標記法，cijkl可簡化為CIJ(I，J∈[1，2，3，4，5，6])。在三維各向同性彈性介質中，獨立的彈性系數(shù)簡化為C33和C44，與速度的關系可表示為：

式中，VP和VS分別為P波和S波速度?；谧V元法模擬地震波在復雜三維彈性介質中的傳播[14]，與傳統(tǒng)有限差分法相比，精度更高[18，19]。使用國家“天河”超算平臺的GPU節(jié)點(顯卡型號為NVIDIA A100，40 GB顯存，6192 CUDA核)加速計算，與常規(guī)多線程CPU節(jié)點(18核)相比，計算效率提高50～60倍。

FWI使用局部最優(yōu)化算法(如最速下降法和牛頓法等)求取目標函數(shù)的最小值，需要計算目標函數(shù)對模型參數(shù)的一階偏導，即梯度(或敏感核)。為避免大型雅克比矩陣的直接求取，基于伴隨狀態(tài)法，彈性速度的敏感核可以通過正演模擬波場和伴隨波場的互相關計算[20-22]：

mn+1=mn+αnΔmn(6)

式中，αn為第n次迭代的步長，可使用線性搜索方法計算獲得；Δmn為搜索方向，使用擬牛頓L-BFGS優(yōu)化算法計算[23，24]，并實現(xiàn)迭代反演。本研究使用mask函數(shù)去除梯度異常值，通過三維高斯函數(shù)對梯度進行平滑，提高反演的穩(wěn)定性。

由于地震波形與模型參數(shù)之間具有非常強的非線性[25-27]，使用傳統(tǒng)波形差目標函數(shù)容易陷入局部最小值，得到不準確的反演結果。針對該問題，發(fā)展了多尺度、多目標函數(shù)和多波反演以及正則化約束等方法：

1)在多尺度反演方法中，對數(shù)據(jù)進行帶通濾波，使用低頻段數(shù)據(jù)構建模型的長波長部分，拓展頻率范圍，基于低頻反演結果，使用高頻信息構建短波長部分，獲得最終反演結果。

2)采用多目標函數(shù)反演策略，依次使用互相關旅行時、包絡和波形差目標函數(shù)進行反演。旅行時目標函數(shù)為：

式中，Δti為旅行時差。包絡目標函數(shù)為：

3)發(fā)展了煤礦巷道多波反演方法，首先使用P波和S波分別反演目標區(qū)域的VP和VS結構。由于槽波在煤層頂?shù)捉缑娼涍^多次全反射，引入槽波使用整個炮集記錄進行反演，更加有助于對煤層頂?shù)捉缑孢M行成像。

4)基于已有測井數(shù)據(jù)(包括煤層頂?shù)捉缑嫖恢煤兔簩雍穸?構建先驗模型mprior，在FWI目標函數(shù)中加入懲罰項，建立正則化目標函數(shù)：

式中，β為權重算子；W為正則化參數(shù)；χ(m)為常規(guī)目標函數(shù)。使用該目標函數(shù)可壓制不符合先驗信息的解，降低反演的多解性。

2 小保當?shù)V區(qū)概況及地質背景

小保當?shù)V區(qū)位于陜西省神木市西南部，侏羅紀煤田榆神礦區(qū)西部鄂爾多斯臺地向斜東翼，面積約220 km2。小保當一號煤礦開采2-2煤層，112201和112202為傾角較小的單斜構造，三維地震未發(fā)現(xiàn)較大斷層和明顯的褶皺構造。煤層存在寬緩的波狀起伏，其軸向大致平行于勘探線方向，主要表現(xiàn)為西部兩翼波長較短，東部波長較長，由西向東呈過渡狀態(tài)，在平面上形成中部為喇叭狀凹型地勢。此構造形態(tài)在各個煤層底板等高線中都有表現(xiàn)，說明該構造形態(tài)與含煤地層沉積基底有關。

FWI的目標區(qū)域為112204工作面，位于井田西南側，長約650 m，寬約350 m，如圖1所示，圖1中紅色和藍色圓球分別為震源和接收器。煤層賦存于延安組第四段頂部，埋深為305～387 m，煤層底板標高924～990 m，煤厚6.30～7.24 m，平均煤厚6.50 m。煤層由北向南逐漸變厚，地質構造簡單，存在寬緩的波狀起伏。為實現(xiàn)采煤工作面的地質信息透明化，達到整個礦區(qū)的透明化管理和安全高效采煤，需通過三維彈性波全波形反演，構建目標地層的精細三維地質結構模型。

圖1 FWI工區(qū)平面

3 實際地震數(shù)據(jù)應用

3.1 實際地震數(shù)據(jù)分析與預處理

首先對采集的實際地震數(shù)據(jù)進行了分析和預處理。巷道地震數(shù)據(jù)包括三分量透射接收數(shù)據(jù)和反射接收數(shù)據(jù)，震源類型為炸藥，共248炮，采樣間隔為0.25 ms，最大記錄時間為2 s，采樣點數(shù)為8000。在透射接收的原始地震數(shù)據(jù)中(如圖2所示)，地震波形(包括P波、S波和槽波)連續(xù)且清晰，信噪比較高，但也存在一定的隨機噪音和一定數(shù)量的壞道。圖3為原始地震數(shù)據(jù)的振幅譜，可看到有效信號的頻率范圍為40 Hz至200 Hz。由于反射接收的原始地震數(shù)據(jù)存在非常強的高頻和低頻噪音，有效反射信號非常弱，幾乎無法識別。因此，只使用透射接收的地震數(shù)據(jù)進行反演和成像。對透射接收的地震數(shù)據(jù)進行了預處理，包括去壞道、濾波、去噪和能量均衡，如圖4所示?？梢钥吹?，x和z分量質量較高，y分量數(shù)據(jù)質量較差?？傮w而言，處理后的數(shù)據(jù)集可作為FWI的輸入，并達到預期反演效果。

圖2 工作面運輸巷透射接收的原始地震數(shù)據(jù)

圖3 工作面運輸巷透射接收數(shù)據(jù)頻譜

圖4 處理后的工作面運輸巷透射接收數(shù)據(jù)

為實現(xiàn)有效反演，設計了專門針對巷道全波形反演的地震數(shù)據(jù)預處理流程，首先根據(jù)正演模擬得到相應的采樣率，地震數(shù)據(jù)再按此采樣率進行重新采樣，選取合適的時間窗口分別提取P波、S波和槽波，去除數(shù)據(jù)中的壞道，通過帶通濾波控制反演的頻率范圍，并使用能量均衡將振幅歸一化。具體數(shù)據(jù)處理流程為：輸入數(shù)據(jù)—重新采樣——時間窗口選取—去除壞道—帶通濾波—能量均衡—輸出數(shù)據(jù)。經數(shù)據(jù)處理后獲得的直達P波數(shù)據(jù)如圖5所示。在迭代反演中，需要對地震進行多次處理，為提高效率，開發(fā)了地震數(shù)據(jù)批量處理腳本程序，并與正演模擬和反演程序融合，形成了一套高效的反演算法流程。

圖5 處理后的地震直達P波數(shù)據(jù)

3.2 三維初始模型和網格剖分

根據(jù)勘探地區(qū)的實際情況以及震源和接收器的三維空間絕對位置，經過坐標旋轉和相對坐標轉換，建立了正演模擬和反演所需的相對坐標系統(tǒng)。在工作面運輸巷和工作面回風巷中震源和接收器的分布如圖6所示。該模型X方向長度為650 m，Y方向長度為350 m，Z方向的深度為-30 m(對應海拔為915 m至945 m)。

圖6 旋轉后的相對坐標系統(tǒng)以及觀測系統(tǒng)分布

FWI需要初始彈性參數(shù)模型作為輸入，然后進行迭代更新。根據(jù)地震數(shù)據(jù)中P波和S波的旅行時，獲得了初始模型的P波速度為4000 m/s，S波速度為1750 m/s，密度為2465 kg/m3，作為FWI的初始輸入模型。在X、Y和Z方向的網格數(shù)量分別為120，72和60，其中垂直Z向的網格大小為0.5 m。

3.3 三維彈性速度反演

使用全波形反演算法與流程以及反演參數(shù)對小保當目標區(qū)域的三維速度模型進行構建。首先，在頻率段[30 Hz，60 Hz]內，使用旅行時目標函數(shù)，對三維P波和S波速度結構進行反演。由于旅行時與速度參數(shù)的非線性關系較低，可避免周期跳躍現(xiàn)象。將P波和S波分開反演，可減弱由多參數(shù)串擾導致的噪音。隨后，在該頻率段內分別使用包絡和波形差目標函數(shù)，進行反演，可以進一步提高反演的分辨率。最后使用P波、S波和槽波對速度結構進行反演。加入槽波信息對煤層的頂?shù)捉缑娴姆囱莞佑行А８鶕?jù)多尺度反演策略，繼續(xù)在頻率段[30 Hz，90 Hz]、[30 Hz，120 Hz]和[30 Hz，150 Hz]，重復使用同樣的反演流程與算法，從低頻到高頻逐步構建小保當?shù)V區(qū)地下的速度結構。低頻信息可反演速度的長波長背景結構，高頻信息可構建地層的高頻短波長精細結構。

最終反演得到三維P波速度結構模型如圖7所示，P波速度結構切片如圖8所示。最終反演得到的三維S波速度模型如圖9所示，S波速度結構切片如圖10所示。從速度反演結果可以看出，煤層的三維空間結構和頂?shù)捉缑婵梢郧逦上?，地下煤層分布可清晰識別，總體較平緩，無復雜構造(如斷層和褶皺等)。而基于傳統(tǒng)槽波的探測方法，只能給出工作面內的二維能量成像圖或煤層厚度的等值線圖[28，29]，無法對煤層的三維空間變化進行成像和預測。由此說明，全波形反演對地下三維空間結構的成像結果可達到煤層透明工作面的實際需求。

圖7 反演得到的三維VP速度模型

圖8 VP速度模型切片

圖9 反演得到的三維VS速度模型

圖10 VS速度模型切片

3.4 煤層頂?shù)捉缑婵臻g分布預測

基于三維全波形反演得到的三維彈性速度結構反演結果，對煤層的頂?shù)捉缑娴目臻g分布進行預測，結果如圖11所示。從圖11中可看出，煤層的總體厚度約為6～7 m，空間變化較為平緩，無明顯的特殊構造。為驗證預測結果的準確性，將在測井數(shù)據(jù)中獲得煤層頂?shù)装褰缑嫖恢门c全波形反演預測結果進行了對比，見表1。測井位于工作面運輸巷內，隨機選取了5口測井，在X方向的位置分別為100，230，450，575和625 m?？梢?，全波形反演結果獲得煤層頂板誤差較小，而地板誤差較大，預測精度達到了煤層透明工作面的實際需求，為高效和安全的煤礦開采和生產提供技術支撐和重要信息。

表1 測井數(shù)據(jù)煤層頂?shù)捉缑媾c全波形反演結果對比 m

圖11 煤層頂面(藍色)和底面(紅色)分布預測

4 結 語

針對我國煤礦開采的實際需求，研發(fā)了煤礦巷道三維彈性波全波形反演算法與技術流程，并應用于陜西小保當一號煤礦112204工作面的實際三維地震數(shù)據(jù)處理。在反演中，采用了多尺度、多目標函數(shù)和多波的反演策略，并使用測井先驗信息進行正則化約束，與標準的全波形反演方法相比，可減弱反演的非線性和多解性，獲得更加準確的結果。最終構建了112204采煤工作面精細三維彈性速度結構模型，對煤層的頂?shù)装暹M行了高精度成像，與傳統(tǒng)槽波探測方法相比，可對煤層三維空間分布進行精細刻畫，滿足了采煤工作面的透明化的實際需求。