趙驚濤彭蘇萍陳宗南柳倩男

1. 煤炭資源與安全開采國家重點實驗室,北京 100083；2. 中國礦業(yè)大學(xué)(北京)地球科學(xué)與測繪工程學(xué)院,北京 100083；3. 空軍研究院工程設(shè)計研究所,北京 100068

在煤礦開采過程中,隱蔽小尺度斷層和陷落柱等不連續(xù)地質(zhì)體破壞了煤層的連續(xù)性,易誘發(fā)透水與瓦斯突出等事故,嚴(yán)重威脅煤礦安全生產(chǎn)[1]。斷層易使頂板巖層破碎,支架失穩(wěn),陷落柱迫使采煤機械停產(chǎn)搬家,由此造成的儲量損失高達15% ～30%。

長期以來,我國地質(zhì)、采礦界科技人員對此類問題十分重視,進行了大量的研究工作,但由于小尺度斷層和陷落柱隱伏性強、分布規(guī)律性差,僅靠地質(zhì)和鉆井資料無法滿足煤礦生產(chǎn)需求,礦井設(shè)計存在盲目性,制約礦井安全高效生產(chǎn)。三維地震勘探方法廣泛應(yīng)用在斷層和陷落柱等不連續(xù)地質(zhì)體空間分布的研究上。但傳統(tǒng)三維地震勘探以反射波為主。在此基礎(chǔ)上開展的不連續(xù)地質(zhì)體檢測方法,如相干技術(shù)[2]、地震屬性分析方法[3]和螞蟻體追蹤算法[4]等,主要適用于大尺度斷層和陷落柱,分辨率難以滿足煤礦安全高效生產(chǎn)需求。一般地,煤礦隱蔽性小尺度地質(zhì)異常體主要指長軸直徑小于20 m 的陷落柱和落差小于5 m 的斷層。據(jù)統(tǒng)計,傳統(tǒng)三維地震勘探技術(shù),對于長軸直徑20 m 以上陷落柱解釋的準(zhǔn)確率約為40% ～50% ,在探明落差小于5 m 斷層上仍存在極大挑戰(zhàn)[5]。因此,如何提高隱蔽小尺度斷層和陷落柱探測精度,一直都是地震勘探的難題。

從地震波場看,繞射波是地質(zhì)體尺度小于一個菲涅爾帶的地震響應(yīng),攜帶了關(guān)于斷層、尖滅點、陷落柱和巖溶洞穴等小尺度地質(zhì)體信息,因此有效利用地震數(shù)據(jù)中的繞射波,對提高三維地震成像分辨率至關(guān)重要。但由于繞射波能量較反射波弱1 ～2個數(shù)量級,因此繞射波利用難度較大,需單獨分離繞射波。近年來,利用繞射波識別小尺度不連續(xù)地質(zhì)體(斷層、巖溶洞穴和陷落柱等)的研究日益增多,研究熱點集中在三維疊前繞射波分離與成像方法。繞射波為弱信號,通常淹沒在強反射波背景中,因此分離繞射波是研究難點,在方法上主要分為數(shù)據(jù)域和成像域兩類。

數(shù)據(jù)域分離方法研究較早,是指在地震數(shù)據(jù)偏移之前通過信號變換等方法分離繞射波,主要利用數(shù)據(jù)域反射波和繞射波幾何形態(tài)差異,通過壓制平滑連續(xù)的反射波,分離曲率大、動力學(xué)特征連續(xù)性差的雙曲繞射波。Asgedom 等[6]研究基于相似度和多重信號分類準(zhǔn)則的繞射波分離方法,用于增強繞射波能量。Figueiredo 等[7]利用模式識別技術(shù)研究繞射波自動成像方法,通過k最近鄰值技術(shù)區(qū)分反射波與繞射波。為克服地震采集數(shù)據(jù)有限孔徑和帶寬限制,Gelius 等[8]利用共反射面元疊加和多重信號分類方法,探討了窗口導(dǎo)向的高分辨率繞射波成像方法。朱生旺等[9]利用局部傾角濾波和預(yù)測反演方法實現(xiàn)了炮集數(shù)據(jù)繞射波分離。蔣波等[10]利用反射層拉平方法研究了繞射波分離問題。劉太臣等[11]研究了基于奇異值分解的共偏移距繞射波分離方法。Zhao 等[12-13]利用一致漸近繞射理論和雙指數(shù)擬合最優(yōu)化方法,研究了炮集數(shù)據(jù)域繞射波分離與成像方法,以及基于馬氏距離和繞射波振幅衰減的繞射波自動分離與成像方法。Yu等[14]通過正則化平面波破壞方法,研究了共偏移距數(shù)據(jù)繞射波分離方法,解決了在噪聲背景下反射波傾角估計穩(wěn)定性問題。

成像域分離方法需要提供偏移速度模型,利用反射波和繞射波在成像域的運動學(xué)和動力學(xué)特征差異來分離繞射波。孫贊東和白英哲[15]提出了一種基于傾角域反射波廣義拉東譜的繞射波分離方法。Zhang 等[16]通過在角度域壓制菲涅爾帶,實現(xiàn)了疊前時間域繞射波成像。劉斌等[17]在傾角域研究了基于局部傾角估計的繞射波分離與成像方法。徐輝[18]利用張角和地層傾角成像道集中反射波和繞射波特征差異,通過傾角信息分離了繞射波。Silvestrov 等[19]利用逆時偏移方法產(chǎn)生共成像點道集,并采用Radon 變換分離繞射波,但由于逆時偏移對速度模型依賴性極強,該方法應(yīng)用難度較大。

綜上所述,數(shù)據(jù)域分離方法難以有效處理復(fù)雜波場情況下繞射波分離問題,當(dāng)?shù)卣鸩ㄏ嘟缓拖嗲袝r,很難保證繞射波極性和振幅特征完整性,致使在后續(xù)繞射波成像過程中小尺度地質(zhì)信息極易缺失；成像域分離方法雖能夠利用偏移算子處理地震波場相交情況下繞射波分離問題,但對速度模型依賴性強,復(fù)雜地質(zhì)條件適用性不強。

本文基于曲波稀疏變換和平面波分解,研究三維疊前炮域多參數(shù)稀疏優(yōu)化繞射波分離方法,解決在地震波場相交和相切情況下繞射波分離的難題,使得分離的繞射波具有波形一致性和完整性,進而實現(xiàn)繞射波單獨成像,提高不連續(xù)地質(zhì)體成像分辨率。該研究一定程度上解決了煤礦小尺度斷層和陷落柱等隱蔽致災(zāi)地質(zhì)體探測難題,為我國煤礦安全高效開采提供技術(shù)支撐。

1 繞射波弱信號提取方法

傳統(tǒng)繞射波分離方法一般利用平面波破壞濾波器,在共偏移距數(shù)據(jù)分離繞射波。該方法要求觀測系統(tǒng)規(guī)則化,但一般近偏移距和遠偏移距中地震數(shù)據(jù)分布不規(guī)則,繞射波難以有效分離,且在地震波場交叉或相切情況下,該方法也不能保證繞射波分離的完整性。

本研究以曲波變換表示繞射波,結(jié)合平面波分解,構(gòu)建三維疊前炮域多參數(shù)稀疏優(yōu)化繞射波分離方法,解決在地震波場相交和相切情況下繞射波分離難題。首先,對三維疊前炮集數(shù)據(jù)動校正,消除偏移距造成的時差,使得反射波表現(xiàn)為線性特性,而繞射波為雙曲線特征；其次,利用曲波變換和平面波分解,構(gòu)建多參數(shù)稀疏最小二乘目標(biāo)函數(shù),其中曲波變換系數(shù)為稀疏約束,對應(yīng)于繞射波,而平面波分解算子對應(yīng)的傾角場為平滑約束,對應(yīng)于反射波；最后,通過反動校正算子,完成三維疊前炮域繞射波分離。

1.1 多尺度多方向曲波變換方法

基于Wrapping 算法的曲波變換,在每個象限每個角度的網(wǎng)格是一樣的,每個曲波變換都有正確的方向,定義笛卡爾坐標(biāo)下的曲波系數(shù)為

式中,?(ω)為信號f的傅里葉變換；為剪切矩陣Sθl的逆矩陣；j、l、k分別為尺度、角度和角度平移量；ω為角頻率。

定義原點周圍L1,j×L2,j的區(qū)域:

基于Wrapping 算法的curvelet 變換方法,改變了數(shù)據(jù)的符號,數(shù)據(jù)d[n1,n2]的Wrapping 可以簡化為Wd[n1modL1,j,n2modL2,j]=d[n1,n2],(n1,n2)∈Pj,l。

這種表示將原來的點(n1,n2)映射到原點附近新的位置。對于角度θ∈(π/4,3π/4),Wrapping算法是類似的,需要交換兩個坐標(biāo)的位置。

基于Wrapping 算法的快速曲波變換如下:

第一步:對f[t1,t2]∈L2(R)采用二維快速傅里葉變換得到?[n1,n2],其中t1,t2為空間變量；n1,n2為頻率變量,-n/2＜n1,n2＜n/2。

第二步:對每個尺度j和方向參數(shù)l,構(gòu)造乘積,從而實現(xiàn)局部化。

這個算法的計算復(fù)雜度為O(n2logn)。在實際計算中,計算量不會超出二維快速傅里葉變換的6 ～10 倍。該變換的逆變換為其轉(zhuǎn)置,屬于一種連續(xù)變換。

1.2 平面波分解方法

平面波分解利用局部平面波方程表征地震數(shù)據(jù)構(gòu)造特征,是一種基于地震局部傾角估計的預(yù)測方法。該算法依次從相鄰地震道預(yù)測當(dāng)前地震道,預(yù)測算子按照局部傾角模式對原始輸入地震數(shù)據(jù)進行平滑估計,因此可以預(yù)測地震反射波。平面波分解方法可定義為線性算子,表達式如下:

式中,r為預(yù)測的地震反射波數(shù)據(jù)；s為動校正后的炮集地震道；L2(σ)為線性平面波算子；σ為地震傾角斜率。

式中,I為單位矢量；σ為局部傾角斜率矢量；Pi,j(σi)表示由第i道地震數(shù)據(jù)預(yù)測第j道數(shù)據(jù)。

1.3 多參數(shù)稀疏優(yōu)化繞射波分離方法

基于曲波稀疏變換和平面波分解構(gòu)建的多參數(shù)稀疏優(yōu)化繞射波分離模型如下:

式中,‖·‖1為l1稀疏約束范數(shù)；d為三維疊前炮集動校正數(shù)據(jù)；L1為曲波變換算子；m1為其變換系數(shù)；m2為其傾角場；α為正則化參數(shù)。

雖然基于l1范數(shù)誤差函數(shù)具有非平滑性質(zhì),但數(shù)值極小值難以計算,并且l1范數(shù)誤差函數(shù)在小誤差處也會按照大誤差處進行“嚴(yán)重對待”,因此適用性有限。Huber 范數(shù)收斂更快,且由于與較少參數(shù)有關(guān)而更容易操作,可用于近似表示l1范數(shù):

式中,ξ為可調(diào)節(jié)的正則化參數(shù)；hξ為分段函數(shù),稀疏約束項。

引入拉格朗日算子后,多參數(shù)稀疏優(yōu)化繞射波分離模型等價于求解如下無約束問題:

式中,β為正則化因子。

該模型可通過交替方向迭代法求解,步驟如下:

由于三維疊前炮集地震數(shù)據(jù)規(guī)模較大,因此采用擬牛頓法求解子問題,擬牛頓法只需計算一階偏導(dǎo)數(shù),計算簡單且運算量小,比解析求得Hessian矩陣更加高效穩(wěn)定。

基于curvelet 稀疏變換的三維疊前炮域多參數(shù)稀疏優(yōu)化繞射波分離方法的流程如圖1 所示。

圖1 三維疊前炮域多參數(shù)稀疏優(yōu)化繞射波分離流程Fig.1 Flow chart of three-dimensional pre-stack diffraction separation with multi-parameter sparse optimization

2 數(shù)值模型實驗

三維地質(zhì)模型由水平反射界面、斷層和散射點組成,如圖2 所示,圖中A 和B 標(biāo)記為散射點,通常描述地質(zhì)上的巖溶洞穴或小尺度陷落柱。水平反射界面將產(chǎn)生較強能量的反射波,而斷層邊界和散射點產(chǎn)生的地震響應(yīng)分別為邊緣繞射波與散射波。采用三維Kirchhoff 正演算法合成地震記錄,射線追蹤路徑如圖3 所示,用于模擬地震波傳播規(guī)律。激發(fā)震源采用主頻為30 Hz 的Ricker 子波,其中全波場(反射波與繞射波)單炮記錄的一條接收線如圖4(a)所示,在該圖中強反射波對應(yīng)于水平反射界面,1 和2 為斷層上下邊界產(chǎn)生的邊緣繞射波。通常,斷層產(chǎn)生的繞射波能量較弱于反射波,但強于散射點,因此在圖4(a)中,更為微弱的散射波幾乎不可見。利用本文提出的三維疊前炮域多參數(shù)稀疏優(yōu)化繞射波分離方法去除強反射后,得到的繞射波炮集記錄如圖4(b)所示,該圖中強反射波去除干凈,1 和2 對應(yīng)的繞射波特征清晰,A 和B對應(yīng)的散射波弱信號突顯。在反射波與繞射波相切或交叉情況下,傳統(tǒng)平面波濾波器的繞射波分離方法[9]由于波場估計不準(zhǔn),使得分離的繞射波不完整,出現(xiàn)繞射波頂部殘缺或者交叉位置波場不連續(xù)等問題。本文提出的方法在繞射波分離結(jié)果圖4(b)中較好地解決了此類問題。利用三維疊前時間偏移算法分別處理全波場和繞射波地震數(shù)據(jù),獲得成像結(jié)果如圖5 所示。圖5(a)能夠清晰刻畫水平反射界面和斷層邊界,但散射點無法表征；圖5(b)強反射去除干凈,散射點和斷層斷點成像清楚。由于圖5(b)中的斷層斷點成像存在左右極性反轉(zhuǎn),而散射點無此特征,因此該現(xiàn)象可作為區(qū)分?jǐn)鄬雍拖萋渲鶅深愋〕叨入[蔽地質(zhì)體的關(guān)鍵因素。

圖2 三維地質(zhì)模型Fig.2 3D geological model

圖3 三維模型單炮射線追蹤路徑Fig.3 3D ray tracing path exciated by one shot

圖4 模擬數(shù)據(jù)炮集Fig.4 Synthetic shot gather

圖5 模擬數(shù)據(jù)偏移剖面Fig.5 Profile of synthetic migation result

數(shù)值模型測試表明,本文提出的基于曲波和平面波分解的三維疊前炮域多參數(shù)稀疏優(yōu)化繞射波分離方法,能夠有效去除反射波,突顯弱繞射信息,有助于實現(xiàn)小尺度不連續(xù)地質(zhì)體的高精度探測。

3 實際資料應(yīng)用

3.1 研究區(qū)基本地震地質(zhì)條件

在煤炭生產(chǎn)中,斷層、陷落柱及采空區(qū)等不良地質(zhì)體,是導(dǎo)致煤礦事故發(fā)生的主要災(zāi)害源。因此,在礦井采掘活動前,需要查明開采區(qū)域地質(zhì)構(gòu)造情況,以確保礦井安全生產(chǎn)。研究區(qū)為山西陽泉煤業(yè)二礦九采擴區(qū),含煤地層從老到新依次為石炭系上統(tǒng)太原組、二疊系下統(tǒng)山西組；下伏地層為石炭系中統(tǒng)本溪組,上覆地層為二疊系下統(tǒng)下石盒子組。構(gòu)造基本形態(tài)為一走向北西、向南西傾伏的單斜構(gòu)造,以褶皺構(gòu)造為主,斷裂構(gòu)造次之。陷落柱剖面形狀多為上小下大的不規(guī)則柱體,局部呈葫蘆狀,柱體內(nèi)巖性雜亂。在陷落柱邊緣,煤層及頂板裂隙密集,有時伴有小斷層或煤層底板牽引下彎現(xiàn)象。本次主要勘探目的層為8 號和15 號煤層,煤層發(fā)育較穩(wěn)定、結(jié)構(gòu)較簡單,與圍巖物性差異大,煤巖之間能構(gòu)成較好的波阻抗界面,可形成能量較強、連續(xù)性較好的反射波,具有較好的地震地質(zhì)條件。

3.2 實際資料繞射成像結(jié)果

野外三維地震采集為中點激發(fā)12 線16 炮線束狀觀測系統(tǒng),接收道數(shù)840,道間距10 m,接收線距40 m,最小炮檢距5 m,最大炮檢距458 m。經(jīng)過觀測系統(tǒng)加載、靜校正、噪聲去除和反褶積等一系列常規(guī)處理后,得到的全波場炮集記錄的部分?jǐn)?shù)據(jù)如圖6(a)所示,其中強振幅連續(xù)好的雙曲線同相軸對應(yīng)于較為平緩的強波阻抗反射界面,如煤層,而斷層和陷落柱則表現(xiàn)為弱繞射波信號,通常在炮集中很難被直接觀察到。利用本文提出的三維疊前炮域多參數(shù)稀疏優(yōu)化繞射波分離方法,消除強反射信號屏蔽作用后,得到的繞射波炮集記錄部分如圖6(b)所示,該記錄去除了雙曲形態(tài)的強反射信號。為進一步對分離的繞射波信號進行質(zhì)量控制,將繞射波疊加剖面圖與常規(guī)全波場疊加剖面對比(圖7),繞射波疊加剖面中呈現(xiàn)為線性的反射波同相軸去除較好,而雙曲形態(tài)的繞射波弱信號得到了突顯。分別對全波場和繞射波炮集記錄進行三維疊前偏移處理,成像結(jié)果如圖8 所示。

圖6 繞射和全波場炮集Fig.6 Diffraction and full-wave shot gather

圖7 繞射和全波場疊加Fig.7 Diffraction and full-wave stacked section

圖8 繞射和全波場疊前偏移Fig.8 Diffraction and full-wave migration

陷落柱塌陷引起地層不連續(xù)性,在三維地震數(shù)據(jù)體上通常表現(xiàn)為同相軸扭曲和分叉、能量變?nèi)?、連續(xù)性變差。常規(guī)三維疊前偏移剖面解釋出的DX6、DX7 和DX8 三個陷落柱,發(fā)育在8 號和15 號煤層中,形狀近橢圓狀,但陷落柱邊界形態(tài)不清楚,尤其是DX7 陷落柱。在繞射波三維疊前偏移剖面中,3 個陷落柱邊界形態(tài)斷點清晰,且陷落柱內(nèi)部特征豐富,陷落柱邊緣伴有小斷裂,說明隱蔽地質(zhì)體探測精度得以提高。三維疊前繞射波偏移剖面顯示DX6 陷落柱在8 號煤相交位置處存在短弱同相軸錯斷,在8 號煤層以下還存在短反射同相軸；在15 號煤層以下發(fā)現(xiàn)了斷裂。常規(guī)三維疊前偏移剖面解釋出的DX7 陷落柱頂點位置在繞射波成像剖面中存在一段弱同相軸。繞射波成像剖面中給出的DX8 陷落柱與8 號和15 號煤的相交位置更為清晰,內(nèi)部的同相軸下彎變形嚴(yán)重,一定程度上反映了內(nèi)部的塌陷情況。

3.3 繞射波沿層屬性分析

地震數(shù)據(jù)體中包含大量的地質(zhì)信息,不同的地震屬性可能與地質(zhì)情況存在關(guān)聯(lián),能夠增加對斷層、陷落柱和撓曲等地質(zhì)構(gòu)造的識別精度。沿著15 號煤層,對常規(guī)三維疊前偏移數(shù)據(jù)進行方差屬性分析,計算參數(shù)選取為3 道×3 道乘法模式,時窗長度為15 ms,方差屬性平面如圖9 所示。圖9 中斷層和陷落柱等構(gòu)造引起的地層出現(xiàn)不連續(xù),大斷裂形態(tài)清晰,但隱蔽陷落柱和斷裂刻畫精度不夠。三維疊前繞射波偏移數(shù)據(jù)方差屬性平面如圖10 所示。圖中隱蔽斷裂清晰,小尺度陷落柱刻畫精度提高；主測線201 和聯(lián)絡(luò)線251 交叉位置顯示出若干環(huán)形區(qū)帶,該特征符合陷落柱平面展布規(guī)律,而附近存在的小尺度線性斷裂,在常規(guī)偏移數(shù)據(jù)方差平面圖上并無顯示。同時,還出現(xiàn)了一系列與主斷層具有相同走向的線性小尺度斷裂,呈現(xiàn)北東向間斷分布特征。這種線性斷裂進一步揭示了巖溶發(fā)育與斷層活動的伴生關(guān)系。煤礦應(yīng)用效果表明,本文提出的三維疊前炮域多參數(shù)稀疏優(yōu)化繞射波分離方法,能夠有效去除煤層強反射信息屏蔽作用,分離出完整的繞射弱信號,疊加剖面中繞射波的雙曲形態(tài)進一步表明分離出繞射波的可靠性。在繞射波地震解釋上,通常需要全波場成像結(jié)果提供宏觀的地質(zhì)背景,在此基礎(chǔ)上進一步利用繞射波成像結(jié)果揭示隱蔽地質(zhì)體,如陷落柱和斷層等。從成像剖面上看,繞射波成像結(jié)果對陷落柱的邊界形態(tài)刻畫更為清晰,內(nèi)部細節(jié)描述更為豐富。從沿層屬性上看,繞射波成像數(shù)據(jù)對隱蔽斷裂的揭示更為細膩,對陷落柱的探測能力有很大提升。

圖9 15 號煤常規(guī)偏移數(shù)據(jù)方差平面Fig .9 Variance map of conventional migration data of No.15 coal

圖10 15 號煤三維疊前繞射波偏移數(shù)據(jù)方差平面Fig.10 Variance map of 3d pre-stack diffraction migration data of Coal No.15

4 結(jié) 論

(1) 常規(guī)疊前地震繞射波分離方法利用平面波破壞濾波器在共偏移距進行,該方法要求觀測系統(tǒng)規(guī)則化,并且在地震波場交叉或相切情況下很難保證繞射波分離的完整性?；谇ㄗ儞Q和平面波分解,本文構(gòu)建了三維疊前炮域多參數(shù)稀疏優(yōu)化繞射波分離方法,解決了地震波場相交和相切情況下繞射波分離問題。

(2) 在地震高分辨率探測上,本文提出的方法能夠有效消除強反射波,分離出完整的繞射波/散射波弱信號,提高斷層斷點和散射點成像精度,有助于實現(xiàn)小尺度不連續(xù)地質(zhì)體的高精度探測。

(3) 在繞射波成像結(jié)果解釋上,需要結(jié)合全波場成像結(jié)果的宏觀地質(zhì)背景,進一步揭示隱蔽致災(zāi)地質(zhì)體,其中極性反轉(zhuǎn)特征可用于區(qū)分?jǐn)鄬訑帱c和散射點。

(4) 繞射波成像雖然具有高分辨率探測能力,但在地質(zhì)認(rèn)識上還存在很多研究難點,如繞射波高分辨率成像結(jié)果的地質(zhì)可靠性分析、繞射波成像智能解釋等。利用繞射波多方位性波場特征,結(jié)合人工智能技術(shù)建立多屬性模型是未來的研究方向,將有助于進一步提升煤礦小尺度不連續(xù)隱蔽致災(zāi)體探測能力。