基于機(jī)器學(xué)習(xí)的時間反轉(zhuǎn)散射體檢測方法

韓令賀， 狄?guī)妥專?胡自多， 劉威， 王國慶， 徐中華

1 中國石油大學(xué)(北京)油氣資源與探測國家重點實驗室， 北京 102249 2 中國石油勘探開發(fā)研究院西北分院， 蘭州 730020 3 中國石油天然氣集團(tuán)有限公司油藏描述重點實驗室， 蘭州 730020

0 引言

傳統(tǒng)地震勘探通常使用反射波信息估算地下介質(zhì)的速度，并通過不同的偏移方法對地下連續(xù)反射層和地質(zhì)構(gòu)造進(jìn)行成像(李振春，2014；陳生昌和周華敏，2016).雖然近年來速度建模和偏移方法的進(jìn)步大大提高了地下地質(zhì)體的成像精度(Wang and Tsvankin，2013; Fletcher et al.，2019)，但反射波勘探由于其本身的局限性，仍無法準(zhǔn)確識別出地下的不連續(xù)地質(zhì)體，如斷層、尖滅點、溶洞、裂縫等小尺度散射體.而這些與地震主波長尺度相當(dāng)或更小的局部構(gòu)造和巖性信息對于提高地震勘探的分辨率具有重要意義，也有利于提高小尺度地質(zhì)體的預(yù)測精度和鉆探目標(biāo)靶點選擇的可靠性.

地下小尺度散射體的地震響應(yīng)常表現(xiàn)為繞射波或散射波(繞射波即狹義上的散射波)，因此，繞射波成像成為檢測地下散射體最常用的方法.在地震原始數(shù)據(jù)中，繞射波能量一般比反射波能量弱一到兩個量級(Landa et al.,2008)，往往淹沒在較強(qiáng)的反射波信號中，準(zhǔn)確識別難度較大，而在常規(guī)資料處理(動校正、疊加)中，繞射波能量也會被壓制.因此，繞射波成像首先要對繞射波和反射波進(jìn)行較好的分離，再單獨對繞射波進(jìn)行偏移成像.國內(nèi)外學(xué)者根據(jù)繞射波與反射波在運動學(xué)和動力學(xué)特征上的差異，在繞射波分離與成像方面提出了很多不同的方法：Landa等(1987)提出在共繞射點剖面上利用相關(guān)來增強(qiáng)繞射點位置的地震信號振幅，以檢測地下局部非均勻地質(zhì)體；Khaidukov等(2004)根據(jù)反射波和繞射波的聚焦特性不同，采用聚焦-反聚焦方法提取出主要包含繞射波信息的炮集并進(jìn)行成像；Koren和Ravve(2010)則利用局部角度域成像構(gòu)建全方位方向角道集，并通過鏡像加權(quán)疊加和繞射加權(quán)疊加分別得到反映地下連續(xù)地質(zhì)體和小尺度地質(zhì)體的成像結(jié)果；朱生旺等(2013)利用局部傾角濾波和預(yù)測反演聯(lián)合分離繞射波，可以提取出完整的繞射波信息.除此之外，還可以利用本征矢量法(Bansal and Imhof, 2005)、Radon方法(Nowak and Imhof,2004)、平面波分解(PWD)方法(黃建平等，2012)和多聚焦方法(Berkovitch et al., 2009)進(jìn)行繞射波分離與成像.

地下小尺度散射體的物性參數(shù)通常與圍巖有較大差異，在速度模型中表現(xiàn)為速度突變或速度異常區(qū)，對這些速度異常區(qū)的檢測和定位是無損檢測技術(shù)主要解決的問題.同時地下小尺度散射體也可以看作是繞射源，對小尺度散射體的檢測和地球物理中廣泛存在的震源定位問題有一定的相似性，如天然地震震源定位、微地震震源檢測、火山震顫等.

近年來，時間反轉(zhuǎn)原理逐漸被應(yīng)用于震源定位問題，該原理基于無耗散介質(zhì)中波動方程具有時間不變性和互易性的性質(zhì)，通過將記錄的波場數(shù)據(jù)時間反轉(zhuǎn)、反向傳播并聚焦到初始的震源位置.時間反轉(zhuǎn)是一種基本的物理規(guī)律，具有物理可實現(xiàn)性.法國巴黎大學(xué)的Fink(1992)在超聲學(xué)領(lǐng)域提出時間反轉(zhuǎn)鏡(TRM，Time Reverse Mirror)技術(shù)，并將其成功應(yīng)用于粉碎腎結(jié)石的實驗.隨后，時間反轉(zhuǎn)在水下聲學(xué)、通信、無損檢測、碎石術(shù)和癌癥治療等很多領(lǐng)域得到推廣應(yīng)用(Kuperman et al.,1998; Chakroun et al., 1995, Thomas and Fink, 1996).如果給定準(zhǔn)確的地下介質(zhì)模型，時間反轉(zhuǎn)可以通過波動方程數(shù)值實現(xiàn).在勘探地震學(xué)領(lǐng)域，時間反轉(zhuǎn)原理主要用于被動源的震源定位.McMechan等(1985)提出可以利用波場反向外推進(jìn)行震源分析，并利用該技術(shù)對1983年加州一次天然地震的震源進(jìn)行了分析.Lu等(2008)利用時間反轉(zhuǎn)原理定位微地震震源，以在開發(fā)過程中監(jiān)控儲層變化.Zhu(2014)對黏聲介質(zhì)的聲波方程進(jìn)行修改，利用時間反轉(zhuǎn)原理實現(xiàn)了井間地震震源的定位.Lellouch 和Landa(2018)利用時間反轉(zhuǎn)聚焦對井間地震數(shù)據(jù)進(jìn)行速度建模.

對于被動源的震源定位，利用時間反轉(zhuǎn)原理搜索模型中所有時刻反傳波場的最大振幅值即可定位出震源位置及震源激發(fā)時刻.而對于主動源激發(fā)并由地下散射體產(chǎn)生的繞射波，由于繞射波能量相對于反射波能量較弱，同時由于散射體的多樣性導(dǎo)致繞射波形態(tài)也較為復(fù)雜，因此，利用時間反轉(zhuǎn)原理檢測地下散射體的位置面臨很多困難.本文通過分析主動源模型中反傳波場在繞射點的聚焦?fàn)顟B(tài)，提取不同的特征反映波場聚焦的程度，并采用機(jī)器學(xué)習(xí)中的樸素貝葉斯分類算法計算出模型中每個點成為繞射點的概率，最終檢測出地下散射體最有可能存在的位置.最后通過散射體模型和Sigsbee2a模型數(shù)據(jù)的試算，驗證了本文方法的有效性和正確性.

1 時間反轉(zhuǎn)原理

無耗散介質(zhì)的聲波方程可以寫為

(1)

二維情況下，式(1)中P(x,z,t)為壓力場，v(x,z,t)和ρ(x,z,t)分別為速度和密度.該方程中只包含了波場P對時間的二階偏導(dǎo)數(shù)算子，具有時間不變性.如果P(x,z,t)是方程的一個解，那么P(x,z,-t)也是方程的一個解.考慮到時間的因果性關(guān)系，用P(x,z,T-t)表示方程的解，其中T表示記錄的總時間，并假設(shè)該時間T足夠長且P(t>T)=0.

時間反轉(zhuǎn)原理基于無耗散介質(zhì)中波動方程具有時間不變性和炮檢互易性的性質(zhì)，并具有物理可實現(xiàn)性.理想情況下的時間反轉(zhuǎn)腔體實驗由Fink在超聲實驗室完成，如圖1所示.假設(shè)包圍震源的封閉表面上分布著足夠密集的傳感器，可以接收壓力場及其法向?qū)?shù)，根據(jù)惠更斯原理即可以計算出封閉表面內(nèi)任一點的波場值.時間反轉(zhuǎn)腔體實驗包括兩步，第一步：一個點狀震源在封閉表面內(nèi)激發(fā)出球面波，波前在經(jīng)過非均勻介質(zhì)時發(fā)生扭曲，并被封閉表面上的傳感器接收，如圖1左側(cè)所示.第二步：將傳感器接收到的信號沿時間方向反轉(zhuǎn)，并從每個傳感器將時間反轉(zhuǎn)后的信號發(fā)出，波場將沿原路徑返回，并最終在震源點位置聚焦，如圖1右側(cè)所示.

圖1 時間反轉(zhuǎn)腔體實驗(Fink，1999)Fig.1 Time-reversal cavity experiment (Fink，1999)

時間反轉(zhuǎn)腔體只是一種理想情況下的實驗，實際地震勘探中只能在地表放置檢波器，也無法觀測到波場的法向?qū)?shù).因此，實際應(yīng)用中只能采用有限孔徑的時間反轉(zhuǎn)鏡(TRM)來實現(xiàn)波場在震源點位置的聚焦，TRM可以是任意的形態(tài)和維度.醫(yī)學(xué)和無損檢測等領(lǐng)域的研究已證實，在非均勻介質(zhì)中由于存在多次聚焦，TRM對震源點仍具有較好的聚焦能力(Fink,1999).

時間反轉(zhuǎn)原理的數(shù)值實現(xiàn)在地震勘探中已得到應(yīng)用，通常用于被動源的震源定位，其數(shù)值實現(xiàn)包括兩步：(1)將檢波點記錄的地震數(shù)據(jù)沿時間方向反轉(zhuǎn)；(2)將時間反轉(zhuǎn)后的地震數(shù)據(jù)進(jìn)行反向傳播，反向傳播通常采用波動方程有限差分法實現(xiàn).時間反轉(zhuǎn)也是逆時偏移的重要組成部分，逆時偏移在時間反轉(zhuǎn)的基礎(chǔ)上增加了兩步：炮點波場的正向傳播及成像條件的應(yīng)用.成像條件通常選擇炮點正傳波場和檢波點反傳波場的零延遲互相關(guān)，但這兩個波場在波阻抗反射界面上會產(chǎn)生較強(qiáng)的逆散射能量，經(jīng)過互相關(guān)后會產(chǎn)生嚴(yán)重的低頻噪聲.同時，錯誤的子波估計也會導(dǎo)致逆時偏移的成像結(jié)果產(chǎn)生一些假象.

圖2和圖3是對被動源模型進(jìn)行的時間反轉(zhuǎn)試驗，在圖2a的兩層介質(zhì)模型中，被動源位于(1000 m,800 m)處，被動源可以是天然地震或微地震震源，地表接收到的單炮記錄如圖2b所示，記錄長度為1.2 s，波場包含直達(dá)波和層界面產(chǎn)生的反射波，對此單炮記錄進(jìn)行時間反轉(zhuǎn)，并采用聲波方程高階有限差分法進(jìn)行反向傳播，圖3a—d為不同時刻的反傳波場.可以看到，隨著反傳時間從1.2 s減小到0 s，反傳波場在震源點位置逐漸聚焦，在t=0 s即震源激發(fā)時刻震源點位置的能量最強(qiáng).因此，對于被動源的震源定位，利用時間反轉(zhuǎn)原理搜索模型中所有時刻反傳波場的最大振幅值即可定位出震源位置及震源激發(fā)時刻.

圖2 被動源模型時間反轉(zhuǎn)試驗 (a) 速度模型； (b) 單炮記錄.Fig.2 Time-reversal experiment of passive source (a) Velocity model； (b) One shot record.

圖3 不同時刻反傳波場快照 (a) 0.6 s； (b) 0.4 s； (c) 0.2 s； (d) 0 s.Fig.3 Back propagation wavefield snapshot at different times

圖4和圖5是對主動源模型進(jìn)行的時間反轉(zhuǎn)試驗，實際地震勘探都采用人工震源，屬于主動源模型.在圖4a的均勻介質(zhì)模型中有一個散(繞)射點，均勻介質(zhì)的速度為2000 m·s-1,散射點的速度為4000 m·s-1.震源位于(1000 m,0 m)處，散射點位于(1000 m,800 m)處，檢波點位于地表，記錄長度為1.2 s，波場包含直達(dá)波和散射點產(chǎn)生的繞射波.由于繞射波能量一般相對于直達(dá)波小2個量級，相對于反射波小1個量級，因此反傳時去掉直達(dá)波，去除直達(dá)波后的單炮記錄如圖4b所示.圖5a—d為不同時刻的反傳波場.可以看到，在不同時刻的反傳波場呈現(xiàn)出不同的聚焦?fàn)顟B(tài)，在t=0.4 s即震源激發(fā)的波場到達(dá)散射點的時刻，反傳波場恰好在散射點位置聚焦(此時散射點相當(dāng)于二次震源)，在其他時刻反傳波場呈現(xiàn)出非聚焦?fàn)顟B(tài).因此，對于主動源模型，反傳波場理論上應(yīng)該在散射點的激發(fā)時間聚焦，但對于實際資料，速度不可能完全準(zhǔn)確，無法準(zhǔn)確計算出散射點的激發(fā)時刻，只能以某種方法判斷反傳波場的聚焦?fàn)顟B(tài)，達(dá)到檢測地下散射體位置的目的.

圖4 主動源模型時間反轉(zhuǎn)試驗 (a) 速度模型； (b) 去除直達(dá)波后的單炮記錄.Fig.4 Time-reversal experiment of active source (a) Velocity model； (b) One shot record after removing direct wave.

圖5 不同時刻反傳波場快照 (a) 0.8 s； (b) 0.6 s； (c) 0.4 s； (d) 0.2 s.Fig.5 Back propagation wavefield snapshot at different times

2 時間反轉(zhuǎn)散射體檢測方法

實際地震勘探中人們關(guān)注的地下小尺度散射體可以看作是主動源模型中的散射點，根據(jù)時間反轉(zhuǎn)原理，通過分析每個時刻反傳波場的聚焦?fàn)顟B(tài)，可以判斷模型中的每一個點是散射點還是非散射點，理論上可以檢測出地下散射體的位置.我們需要計算模型中每個點的聚焦質(zhì)量來評估其聚焦?fàn)顟B(tài)，對二維空間中任一點(x,z)和任一反傳時刻t，需要計算在一個時空窗內(nèi)的聚焦函數(shù)F(x,z,t)，聚焦函數(shù)的大小正比于該點成為散射點的概率.由于存在多次聚焦現(xiàn)象，波場會在不同的時刻在散射點聚焦多次，將所有時刻的聚焦函數(shù)累加得到每個點最終的聚焦函數(shù)，累加值越大，該點為散射點的概率越大，結(jié)果也越可靠.關(guān)于聚焦函數(shù)的選取可以借鑒圖形圖像處理方面的算法(Pertuz et al.,2013)，但對于地震反傳波場，由于波場形態(tài)的復(fù)雜性和散射體的多樣性，同時缺少準(zhǔn)確的震源子波和模型速度信息，很難給出一個準(zhǔn)確的聚焦函數(shù)分析波場的聚焦?fàn)顟B(tài).

2.1 樸素貝葉斯分類算法原理

考慮到時間反轉(zhuǎn)散射體檢測需要對模型中任一點分析其反傳波場的聚焦?fàn)顟B(tài)，從而判斷該點是散射點還是非散射點，因此，這是一個典型的分類問題，更適合采用基于機(jī)器學(xué)習(xí)的分類算法.目前常用的分類算法包括樸素貝葉斯分類、決策樹、支持向量機(jī)、邏輯回歸、隨機(jī)森林、K最近鄰、人工神經(jīng)網(wǎng)絡(luò)和深度學(xué)習(xí)等算法，其中樸素貝葉斯分類算法具有算法簡單、分類準(zhǔn)確率高、速度快的優(yōu)點，在工業(yè)界得到了廣泛應(yīng)用.

樸素貝葉斯分類是一種基于概率統(tǒng)計的分類方法，以貝葉斯定理為理論基礎(chǔ)，并假設(shè)特征條件之間相互獨立.該方法是一種監(jiān)督學(xué)習(xí)方法，先通過已給定的訓(xùn)練數(shù)據(jù)，以特征條件之間獨立作為前提假設(shè)，學(xué)習(xí)從輸入到輸出的概率密度函數(shù)，再將學(xué)習(xí)到的模型應(yīng)用到新數(shù)據(jù)，輸出新數(shù)據(jù)的分類結(jié)果.

假設(shè)有變量x和y，x表示特征變量，且x有n種取值，y表示分類變量，且y有m種取值，那么對于貝葉斯分類問題，根據(jù)貝葉斯公式和全概率公式，可通過式(2)計算出類別y的后驗概率p(y|x)，即特征x屬于y類的概率：

(2)

式中，p(x,yi)表示聯(lián)合概率，p(yi)表示類別yi的先驗概率，p(x|yi)表示似然函數(shù)，即特征x在類別yi下的條件概率，p(x)表示證據(jù)因子，對于給定的特征x，p(x)為固定值，且與類別y無關(guān).

樸素貝葉斯分類基于貝葉斯定理并假設(shè)特征條件之間相互獨立，則式(2)可改寫為：

(3)

式(3)即為樸素貝葉斯分類器的表達(dá)式，對于訓(xùn)練樣本來說，需要提取樣本的特征屬性，并計算先驗概率p(yi)，即每個類別yi在訓(xùn)練樣本中的出現(xiàn)頻率，及每個特征屬性對每個類別的條件概率p(xj|yi).

2.2 樸素貝葉斯分類實現(xiàn)過程

本文將樸素貝葉斯分類算法應(yīng)用于時間反轉(zhuǎn)散射體檢測中，對模型中每個時刻每個點的反傳波場進(jìn)行分類，判斷該點是散射點還是非散射點.根據(jù)對波場分類的任務(wù)特點，樸素貝葉斯分類算法的實現(xiàn)過程包括三大步：

2.2.1 數(shù)據(jù)準(zhǔn)備階段

根據(jù)分類任務(wù)的具體情況確定特征屬性，對每個特征屬性進(jìn)行適當(dāng)劃分并分類，形成訓(xùn)練樣本集合.這一階段是整個樸素貝葉斯分類中唯一需要人工完成的階段，其質(zhì)量對整個過程將有重要影響，分類器的質(zhì)量很大程度上由特征屬性、特征屬性劃分及訓(xùn)練樣本質(zhì)量決定.

對于模型中每個時刻每個點的反傳波場來說，特征屬性應(yīng)該能反映反傳波場的聚焦特征，具有物理意義的特征屬性能提高訓(xùn)練效果，本文對一個時空窗內(nèi)的反傳波場計算如下兩種特征屬性：

(1)聚焦點到時空窗中心的加權(quán)距離：

(4)

其中，f1(x,z,t)表示特征屬性，w為模型中任一點的時空窗，時空窗的大小取為時間和空間域的平均波長，di=‖li-lcenter‖表示時空窗內(nèi)任一點到中心點的距離，li表示將時空窗按像素點排列后第i個像素點，Ai為第i個像素點的振幅.如果聚焦點恰好位于時空窗的中心點，f1的值為最小.

(2)反傳波場和高斯函數(shù)的相關(guān)系數(shù)：

(5)

其中，X、Y分別為時空窗的反傳波場和高斯函數(shù)，cov(X,Y)表示兩者的協(xié)方差，σX、σY分別為X、Y的方差.相關(guān)系數(shù)越大，表明時空窗內(nèi)反傳波場的形態(tài)和高斯函數(shù)的形態(tài)越接近，反傳波場聚焦特征越好.

2.2.2 訓(xùn)練階段

這一階段主要工作是計算每個類別在訓(xùn)練樣本中的出現(xiàn)頻率及每個特征屬性對每個類別的條件概率，輸入是特征屬性和訓(xùn)練樣本，輸出是分類器.本文采用高斯樸素貝葉斯分類方法，即假設(shè)每一類的數(shù)據(jù)都服從高斯分布，那么每一類的概率密度函數(shù)為：

p(xj|yi)=

(6)

其中，d為特征屬性的維度，本文計算兩種特征屬性，因此d=2；μi為第i個類別所有特征的均值，Σ為協(xié)方差矩陣.

訓(xùn)練階段通過式(6)計算每一類的高斯概率密度函數(shù)，即可得到每個特征屬性對每個類別的條件概率，同時統(tǒng)計每個類別在訓(xùn)練樣本中的出現(xiàn)頻率，利用式(3)即可得到樸素貝葉斯分類器.

2.2.3 分類階段

分類階段輸入新數(shù)據(jù)，并計算新數(shù)據(jù)的特征屬性，利用訓(xùn)練階段得到的樸素貝葉斯分類器對其進(jìn)行分類.輸入待分類的每炮數(shù)據(jù)的反傳波場，輸出每個反傳時刻模型中每個點成為散射點的概率，因存在多次散射現(xiàn)象，將所有時刻所有炮的計算結(jié)果累加并進(jìn)行歸一化，最后得到模型中每個點成為散射點的概率.

圖6為應(yīng)用樸素貝葉斯分類進(jìn)行時間反轉(zhuǎn)散射體檢測的算法框架圖.在訓(xùn)練階段需分別對無繞射波的反傳波場和有繞射波的反傳波場數(shù)據(jù)提取兩種特征屬性進(jìn)行訓(xùn)練，生成樸素貝葉斯分類器，無繞射波的反傳波場可由平滑速度模型和密度模型進(jìn)行正演模擬得到，有繞射波的反傳波場可由繞射源模型進(jìn)行正演模擬得到；在分類階段，輸入實際數(shù)據(jù)、速度模型和密度模型即可利用分類器進(jìn)行分類，計算模型中每個點成為散射點的概率.

圖6 時間反轉(zhuǎn)散射體檢測算法框架Fig.6 Time-reversal scatterer detection algorithm framework

3 模型算例

3.1 簡單散射點模型

設(shè)計如圖7所示的模型，模型大小為100 m×100 m，網(wǎng)格間距為1 m，速度為常數(shù)700 m·s-1，密度模型為在常密度背景模型上設(shè)計了2個反射層和3個散射點，3個散射點分別位于(30 m，45 m)、(50 m，45 m)、(70 m，45 m).采用二維聲波方程數(shù)值模擬方法計算得到20炮正演數(shù)據(jù)，采用的觀測系統(tǒng)如下：檢波點位于地表，從x=1 m到x=99 m每隔1 m放置一個檢波點，每炮99道接收；炮點位于地表，從x=1 m到x=96 m每隔5 m放置一個炮點，共20炮，主頻80 Hz，采樣間隔1 ms，采樣長度500 ms.正演的第11炮和第12炮數(shù)據(jù)如圖8a所示，為減小繞射波和反射波的能量差異，對正演數(shù)據(jù)加增益(AGC, Automatic Gain Control)并添加高斯白噪后單炮如圖8b所示.從單炮記錄上可以識別出2個反射層的1次反射波及3個散射點的1次散射波，同時由于反射層界面上下較強(qiáng)的阻抗差異，導(dǎo)致單炮發(fā)育幾套短周期多次反射波及多次散射波.本文的時間反轉(zhuǎn)散射體檢測方法不需對反射波和繞射波進(jìn)行分離，只需對單炮數(shù)據(jù)加增益，使兩者的能量處于同一量級即可.

圖7 散射點密度模型Fig.7 Density of scatterer model

圖8 正演炮記錄 (a) 第11和12炮記錄； (b) 加AGC和高斯白噪后炮記錄.Fig.8 Forward modeling shot records (a) The 11th and 12th shot records； (b) Shot records after adding AGC and Gaussian white noise.

圖9為采用圖6所示的時間反轉(zhuǎn)散射體檢測方法對有繞射波和無繞射波的兩類反傳波場進(jìn)行訓(xùn)練的結(jié)果，縱、橫坐標(biāo)分別為利用式(4)、(5)計算的兩種特征屬性，每類訓(xùn)練樣本為100000個，其中紅色圓圈代表有繞射波數(shù)據(jù)的訓(xùn)練結(jié)果，藍(lán)色圓圈代表無繞射波數(shù)據(jù)的訓(xùn)練結(jié)果.從圖中可以看到，兩類數(shù)據(jù)的特征屬性具有明顯差異，有繞射波數(shù)據(jù)和高斯函數(shù)的相關(guān)系數(shù)較高，中心點加權(quán)距離較小；而無繞射波數(shù)據(jù)和高斯函數(shù)的相關(guān)系數(shù)較低，中心點加權(quán)距離較大.利用這兩種特征屬性可以對有繞射波數(shù)據(jù)和無繞射波數(shù)據(jù)進(jìn)行較好的區(qū)分，經(jīng)訓(xùn)練后產(chǎn)生的樸素貝葉斯分類器可用于對模型中的散射體和非散射體進(jìn)行分類.

圖9 兩類訓(xùn)練數(shù)據(jù)特征屬性的交會圖 紅色、藍(lán)色圓圈分別代表有繞射波數(shù)據(jù)和無繞射波數(shù)據(jù)的訓(xùn)練結(jié)果.Fig.9 Crossplot of two kinds of training data feature attributes The red and blue circles represent the training results with and without diffraction data respectively.

應(yīng)用本文方法對20炮正演數(shù)據(jù)進(jìn)行試算，將所有時刻的結(jié)果累加可以得到每炮的散射體檢測結(jié)果.對第11炮正演數(shù)據(jù)檢測結(jié)果進(jìn)行分析，圖10為3個散射點處的檢測概率隨反傳時間的分布圖，可以看出，由于多次散射的存在，在幾個不同時刻都會檢測到散射體聚焦的現(xiàn)象，并且由于炮點位置的不同，不同散射體的聚焦次數(shù)也不同，對第11炮數(shù)據(jù)而言，左側(cè)散射點檢測到2次聚焦，中間散射點檢測到1次聚焦，右側(cè)散射點檢測到4次聚焦，但t=0.25 s處的檢測概率很小，在對所有時刻的結(jié)果累加時會選擇一個門檻值，將較小的概率檢測結(jié)果剔除，并用聚焦次數(shù)對累加結(jié)果進(jìn)行歸一化，得到每炮最終的檢測結(jié)果.圖11a為第11炮數(shù)據(jù)的散射體檢測結(jié)果，對20炮數(shù)據(jù)的散射體檢測結(jié)果進(jìn)行累加，并以總炮數(shù)進(jìn)行歸一化，得到最終的散射體檢測結(jié)果，如圖11b所示.從最終的散射體檢測結(jié)果中可以看出，在3個散射點位置都能檢測到較大的概率，非散射點位置的檢測概率相對很小，檢測結(jié)果能較好的反映散射點的位置.

圖10 散射點位置檢測概率隨反傳時間的分布 (a) 左側(cè)散射點(30 m，45 m)； (b) 中間散射點(50 m,45 m)； (c) 右側(cè)散射點(70 m,45 m).Fig.10 Distribution of detection probability of scattering point location with back propagation time (a) Left scatterer (30 m，45 m)； (b) Middle scatterer (50 m，45 m)； (c) Right scatterer (70 m，45 m).

圖11 散射體檢測結(jié)果 (a) 第11炮數(shù)據(jù)； (b) 20炮數(shù)據(jù).Fig.11 Scatterer detection results (a) The 11th shot； (b) 20 shots.

3.2 多個散射體模型

為檢驗本文算法對地下散射體的縱橫向分辨能力，設(shè)計如圖12所示的多個散射體模型，模型大小為2000 m×1400 m，縱橫向網(wǎng)格間距均為10 m.模型從上至下三層介質(zhì)的速度分別為2500 m·s-1、3000 m·s-1、4000 m·s-1,在第2層介質(zhì)中分布有12個縱橫向不同間距的高速散射體，散射體的大小均為10 m×10 m、速度為4500 m·s-1.第1和第2個散射體與第1層反射界面重疊，第3和第4個散射體的縱向距離為100 m，第5和第6個散射體的橫向距離為100 m，第7和第8個散射體的縱向距離為50 m，第9和第10個散射體的縱向距離為20 m，第11和第12個散射體的橫向距離為50 m.設(shè)計觀測系統(tǒng)如下：檢波點位于地表，每個網(wǎng)格點均放置一個檢波點，每炮201道接收；炮點位于地表，從第2個網(wǎng)格點，每隔6個網(wǎng)格放置一個炮點，共34炮，主頻30 Hz，采樣間隔1 ms，采樣長度1400 ms.采用聲波方程數(shù)值模擬方法進(jìn)行正演，將正演數(shù)據(jù)的直達(dá)波切除并加增益作為本文算法的輸入，其中第16炮數(shù)據(jù)如圖13所示，從正演單炮中可以看到兩層的反射波及散射體產(chǎn)生的大量繞射波.應(yīng)用本文方法對34炮正演數(shù)據(jù)進(jìn)行試算，將所有時刻的結(jié)果累加可以得到每炮的散射體檢測結(jié)果，將每炮檢測結(jié)果累加并進(jìn)行歸一化后得到最終的散射體檢測結(jié)果，如圖14所示.從檢測結(jié)果中可以發(fā)現(xiàn)，對于第1和第2個散射體，雖然兩者和第1層反射界面重疊在一起，輸入數(shù)據(jù)也沒有將反射波和繞射波分離，但本文算法不受反射層的影響，仍能清楚地檢測到散射體的位置，這對于塔里木盆地檢測某些和反射層位置重疊的溶洞等散射體具有重要的參考價值.

圖12 多個散射體模型Fig.12 Multiple scatterer model

圖13 多個散射體模型正演單炮Fig.13 One forward modeling shot of multiple scatterer model

圖14 34炮數(shù)據(jù)的散射體檢測結(jié)果Fig.14 Scatterer detection result of 34 shots

該模型也可以驗證本文方法在檢測散射體時的縱橫向分辨率.第2層介質(zhì)的速度為3000 m·s-1，震源主頻為30 Hz，因此第2層內(nèi)震源子波的波長為100 m.第3和第4個散射體的縱向距離為1個波長，第7和第8個散射體的縱向距離為半個波長，從檢測結(jié)果上可以看到，第3和第4個散射體、第7和第8個散射體都能區(qū)分開，只是由于下方散射體產(chǎn)生的繞射波在向上傳播時會受到上方散射體的影響，所以下方散射體的檢測概率要小于上方散射體；第9和第10個散射體的縱向距離為20 m，小于1/4個波長，兩個散射體的檢測結(jié)果相互重疊為一個檢測概率較大的散射體，因此兩個散射體不能被有效區(qū)分開.

從檢測結(jié)果中可以看到，第11和第12個散射體的橫向距離為50 m，兩個散射體無法被檢測到，究其原因，筆者認(rèn)為這與算法選取的時空窗大小有關(guān)，由于時空窗的大小至少要包括整個地震子波的長度，因此時空窗的大小通常選取為1個波長.當(dāng)兩個散射體的橫向距離小于或等于半個波長時，兩個散射體的繞射波相互疊合導(dǎo)致在1個時空窗內(nèi)的聚焦特征較差，被檢測到的概率較小.當(dāng)兩個散射體的橫向距離大于半個波長時，才能被檢測到并有效區(qū)分開，如第5和第6個散射體的橫向距離為100 m，兩者均可以被清晰地檢測到.

3.3 Sigsbee2a模型

為驗證本文算法對復(fù)雜構(gòu)造的適用性，選取部分Sigsbee2a模型進(jìn)行測試，速度模型如圖15所示，模型中下部含有兩個高速散射體.模型大小為276×226個網(wǎng)格點，縱橫向網(wǎng)格間距均為15.24 m.設(shè)計觀測系統(tǒng)如下：檢波點位于地表，每個網(wǎng)格點均放置一個檢波點，每炮276道接收；炮點位于地表，從第15個網(wǎng)格點，每隔20個網(wǎng)格放置一個炮點，共13炮，主頻20 Hz，采樣間隔2 ms，采樣長度4000 ms.采用聲波方程數(shù)值模擬方法進(jìn)行正演，將正演數(shù)據(jù)的直達(dá)波切除并加增益作為本文算法的輸入，其中第7炮數(shù)據(jù)如圖16所示.應(yīng)用本文算法對13炮正演數(shù)據(jù)進(jìn)行試算，將每炮檢測結(jié)果累加并進(jìn)行歸一化后得到最終的散射體檢測結(jié)果，如圖17所示.可以看到，檢測結(jié)果在模型下部的兩個高速散射體位置能得到較大的概率，同時檢測結(jié)果對于模型3條斷裂中速度突變的斷點位置也能得到較好的反映，將檢測到的斷點位置用紅色虛線連接起來得到斷裂檢測結(jié)果如圖18所示，和模型中的3條斷裂具有較高的吻合度.同時由于斷點位置的速度突變?nèi)跤谏⑸潴w位置，斷點產(chǎn)生的繞射波能量也弱于散射體產(chǎn)生的繞射波能量，因此，斷點位置的檢測概率也小于散射體的檢測概率.

圖15 部分Sigsbee2a速度模型Fig.15 Partial Sigsbee2a velocity model

圖16 Sigsbee2a模型正演單炮Fig.16 One forward modeling shot of Sigsbee2a model

圖17 13炮數(shù)據(jù)的散射體檢測結(jié)果Fig.17 Scatterer detection result of 13 shots

圖18 斷裂檢測結(jié)果Fig.18 Fracture detection result

4 結(jié)論與討論

時間反轉(zhuǎn)原理基于無耗散介質(zhì)中波動方程具有時間不變性和互易性的性質(zhì)，已被廣泛應(yīng)用于被動源震源定位問題，同時也是逆時偏移的重要組成部分.本文將時間反轉(zhuǎn)原理引入到地下散射體檢測和定位中，并針對地下反傳波場形態(tài)的復(fù)雜性和散射體的多樣性，引入機(jī)器學(xué)習(xí)中的樸素貝葉斯分類算法，計算兩種反映波場聚焦?fàn)顟B(tài)的特征屬性，給出適用于時間反轉(zhuǎn)散射體檢測的分類算法框架，計算得到模型中每個點成為散射體的概率，最終檢測出地下散射體最有可能存在的位置.散射體模型和Sigsbee2a模型的試算結(jié)果表明，本文方法在不需對反射波和繞射波分離的情況下，即可實現(xiàn)對地下散射體的檢測和定位，即使散射體和反射層重疊在一起也不受反射層的影響，同時由于考慮了多次散射的影響，檢測結(jié)果能準(zhǔn)確反映地下散射體的位置.利用多個散射體模型對本文方法的縱橫向分辨能力進(jìn)行了探討，結(jié)果表明受選取時空窗大小的影響，本文方法檢測散射體的縱向分辨率要大于橫向分辨率.

訓(xùn)練樣本的質(zhì)量對樸素貝葉斯分類算法的效果具有重要影響，本文給出的時間反轉(zhuǎn)散射體檢測方法流程中，對于有繞射波的訓(xùn)練樣本假定散射體是點源形態(tài)，實際地下散射體的形態(tài)比較多樣，在后續(xù)研究中利用多種散射體產(chǎn)生有繞射波的訓(xùn)練樣本，可能會提高樸素貝葉斯分類算法的分類效果.同時，對于三維情況下的散射體檢測，如何提取能反映三維反傳波場聚焦?fàn)顟B(tài)的特征屬性是一個較大的挑戰(zhàn)，也是作者后續(xù)的研究重點.