共聚焦點層析成像方法

沈天晶 胡葉正 黃旭日 曹衛(wèi)平 徐云貴

(西南石油大學地球科學與技術(shù)學院，四川成都 610000)

0 引言

地震層析成像是了解地下地質(zhì)結(jié)構(gòu)和巖層物理性質(zhì)的主要方法，根據(jù)理論不同可分為射線類和波形類。射線類層析成像由于理論簡單和計算效率高，已廣泛應(yīng)用于實際的資料處理[1-3]；而波形類反演需依賴于初始的速度模型和觀測數(shù)據(jù)的低頻信息[4]。

1992年，Stork[5]利用最小化時間殘差實現(xiàn)層析成像，奠定了層析反演的理論基礎(chǔ)。一般來說，偏移與層析成像交替進行，分別恢復(fù)速度場的高波數(shù)信息(速度界面)以及低波數(shù)信息(速度)[6]。在偏移速度分析中主要有基于共成像點道集(CIG)的曲率分析[7]和深度聚焦分析[8]兩種方法。由于聚焦深度很難拾取，所以深度聚焦分析方法使用較少。當CIG包含多路徑波場傳播信息，曲率分析易造成假象，甚至速度準確時CIG也不能完全拉平。同時，在求解過程中需要拾取剩余曲率以獲得旅行時殘差。總的來說，層析成像是一個迭代求解模型參數(shù)的過程，在每次迭代中都需要進行偏移，計算量較大。同時在射線追蹤的過程中，需要模擬復(fù)雜地質(zhì)結(jié)構(gòu)的反射，增加了正演的難度[9]。

為了解決這些問題，Kabir 等[10]提出了一種基于共聚焦點(CFP)技術(shù)[11]的速度反演方法。該方法利用差分時移(DTS)面板拾取時移來對速度模型進行擾動與判定。利用CFP技術(shù)進行速度反演主要有同步和異步兩種方法[12]。同步反演是聚焦算子與速度模型相結(jié)合，同步更新聚焦算子與速度模型。異步反演是將聚焦算子的計算與速度模型反演分離，在不更新速度模型的情況下，從反射數(shù)據(jù)中獲得聚焦算子。然后通過逆時聚焦算子的單向旅行時進行速度反演，更新算子的位置和速度模型。

在異步反演中，算子更新與速度模型擾動分離。首先利用CFP技術(shù)求取逆時聚焦算子，在求取的過程中需要拾取DTS面板上的時移。一般只需要迭代四次就可獲得逆時聚焦算子，當初始模型接近實際模型時，迭代兩次就可以求解出逆時聚焦算子。在之后的速度反演階段，不再需要拾取時移。所以較傳統(tǒng)層析中每次迭代都需要拾取殘差相比，異步反演大大減少了工作量。在速度反演階段，不需要重復(fù)偏移計算，只在初始階段對逆時聚焦算子進行單炮偏移選取初始的CFP。并且，逆時聚焦算子代表了一組從CFP到地表的單向旅行時，在速度反演階段沒有復(fù)雜界面的反射問題，降低了正演計算的難度。

在速度反演中，需要進行大量的射線追蹤，同時需要求解一個大型矩陣的逆，因而基于CFP技術(shù)的反演方法計算量極大。為此，本文分別從網(wǎng)格劃分、擾動方法和模型更新方法三個方面進行改進。在網(wǎng)格劃分方面可以減少運算時間的方法有二，一是Santos等[13]提出的分層網(wǎng)格化，第二個是Cox[14]提出的不規(guī)則網(wǎng)格化。結(jié)合CFP技術(shù)的特性，本文采用分層網(wǎng)格化。將同一反射界面上的聚焦算子分為一組，并相互連接起來。這樣的劃分方法不用顯式地劃分層位與定義斷層，在計算的過程中利用這些分層進行約束，加速收斂。并且在計算中，算子的互不干擾特性讓并行計算非常容易實現(xiàn)。在擾動方法方面，利用旅行時殘差的分布確定擾動方向，利用旅行時殘差的大小確定擾動量。同時，在擾動后射線旅行時的計算中，利用射線的出射角度和擾動量估計旅行時的變化量，避免大量重復(fù)的射線追蹤計算，大幅提升了計算效率。為了避免擾動過大，難以收斂，本文采取了分步擾動策略，先進行深度擾動再進行水平擾動。最后，利用收斂速度快、收斂性能較好的代數(shù)迭代法(Algebraic Reconstruction Technique，ART)[15]進行速度模型的更新。

1 共聚焦技術(shù)原理與反演方法

1.1 CFP理論

Berkhout[11]根據(jù)地震觀測的幾何關(guān)系提出了WRW模型，該模型將地震波的傳播過程分為震源的激發(fā)響應(yīng)、下行波場傳播算子、反射算子、上行波場傳播算子、檢波器的接收響應(yīng)，即

P(z0)=D(z0)×

(1)

式中：z0為激發(fā)、接收點深度；zm為第m個反射界面深度；P(z0)表示地震記錄；S(z0)和D(z0)分別為激發(fā)響應(yīng)和接受響應(yīng)；W(z0，zm)和W(zm，z0)分別是上行和下行波傳播效應(yīng)算子；R(zm)為反射算子。

共聚焦理論的聚焦過程分為兩步，激發(fā)聚焦與接收聚焦，聚焦過程通過逆時聚焦算子實現(xiàn)。逆時聚焦算子可以看作地下反射點激發(fā)、地表接收的單程旅行時數(shù)據(jù)。逆時聚焦算子表達式為

(2)

根據(jù)上述理論，CFP道集求取方法如下：首先，利用聚焦算子與共炮點道集對應(yīng)的道進行卷積，這就是第一步聚焦。其中，聚焦算子為逆時聚焦算子的按時間軸倒置。然后，將該道集中所有的響應(yīng)結(jié)果進行疊加，疊加結(jié)果為CFP道集的一道。對所有共炮點道集進行相同的過程，就生成了CFP道集。CFP道集Pj(z0，zm)(j=1，2，…，J，J為道數(shù))反映了地下反射點激發(fā)、地表接收的地震響應(yīng)，其表達式為[5]

(3)

比較CFP道集與逆時聚焦算子可以看出，速度模型隱含在傳播矩陣中，如果速度模型正確，那么逆時聚焦算子與CFP道集有相同的旅行時，即等時原理。用聚焦算子對CFP道集進行互相關(guān)運算，就可以得到DTS面板，這是第二步聚焦。在兩步聚焦過程中，消除了下行波與上行波的傳播效應(yīng)，所以差分時移面板就相當于在聚焦點處激發(fā)、聚焦點處接收所產(chǎn)生的響應(yīng)。因此在DTS面板中，所有的響應(yīng)都在零時刻發(fā)生，若偏離了零時刻，則表示逆時聚焦算子不準。逆時聚焦算子的求取主要利用了上述性質(zhì)，具體方法是從DTS面板上拾取偏離零時刻的響應(yīng)時間，將這個時間的一半作用在聚焦算子上，完成一次迭代，直到DTS面板位于零線為止。具體流程如圖1所示。

圖1 逆時聚焦算子求取流程

逆時聚焦算子是在不知道具體速度模型的情況下，通過數(shù)據(jù)驅(qū)動更新的方法從一組共炮點道集中獲得[16]。在本文的速度反演方法中，逆時聚焦算子是重要的輸入數(shù)據(jù)。

1.2 速度反演方法

1.2.1 速度模型更新方法

本文采用代數(shù)迭代法更新速度模型。首先將模型離散為網(wǎng)格，每個網(wǎng)格有一個對應(yīng)的參數(shù)值，然后用大型的線性方程組描述速度模型與旅行時的關(guān)系

(4)

將上式歸納可寫為

(5)

式中：sj為第j個網(wǎng)格的慢度；ti為第i條射線的旅行時；li,j表示第i條射線在第j個網(wǎng)格內(nèi)的長度；M為射線條數(shù)；N為網(wǎng)格總數(shù)。利用代數(shù)迭代法求解上述方程組。對于第i條射線，第j個網(wǎng)格的ART迭代公式為

(6)

式中：k為迭代次數(shù)；λ是松弛系數(shù)，取值區(qū)間為(0，2)。在計算完所有的焦點與射線后，即為完成一次迭代，每次迭代都以上次迭代的結(jié)果作為初值，直到達到收斂條件。具體迭代步驟如下：

(1)給定初始速度模型；

(7)

(4)計算第j個網(wǎng)格的慢度修正量，并轉(zhuǎn)化為速度修正量

(8)

(9)

(5)修正第j個網(wǎng)格的速度

(10)

(6)將上一輪的結(jié)果作為初值，重復(fù)步驟(2)～步驟(5)，直到達到收斂條件。

1.2.2 聚焦點位置估計

除了速度模型，聚焦算子的空間位置也是反演重要的一部分。若焦點位置不正確，由代數(shù)迭代算法反演出的速度模型也會不準確。為了找到最正確的聚焦算子的空間分布，需要分析模擬旅行時與逆時聚焦算子的旅行時的殘差分布，確定擾動量和擾動方向。

首先需要選取初始的共聚焦點位置，共聚焦點初始位置的選取決定著下一步反演的參數(shù)選取以及約束策略，所以相對科學的選取辦法可以減少反演的計算量，提高反演的質(zhì)量。本文利用逆時聚焦算子做單炮偏移，根據(jù)偏移后的能量分布選取初始的共聚焦點位置，即用初始速度模型對逆時聚焦算子進行偏移，能量會匯聚到一個區(qū)域，然后求取該區(qū)域的能量重心

(11)

式中：L為能量區(qū)樣點數(shù)；EK為樣點(xK，zK)的能量。

提取初始共聚焦點后，需要對其位置進行分析。圖2為當聚焦點有深度誤差時正確算子與估計算子計算的旅行時信息對比，可以看出，當聚焦點的深度比實際深度小時，整體射線路徑變短，計算旅行時比實際旅行時小。在進行代數(shù)迭代時，旅行時殘差數(shù)變大。在利用旅行時殘差進行多次代數(shù)迭代后，模型速度與焦點位置會達到一個耦合關(guān)系，反映的是與當前焦點位置最匹配的模型速度。在圖2例子中，因為估計算子的聚焦點深度比準確聚焦點的深度小，所以與估計算子相匹配的模型速度值也比準確的模型速度值小，從而才能達到耦合關(guān)系。圖3為在達到這種耦合關(guān)系時正確聚焦點與估計算子的聚焦點計算旅行時的相對關(guān)系。在圖3中，在分界線以左，殘差為正，在迭代過程中使模型速度變小；在分界線以右，殘差為負，反映在迭代過程中使速度變大。當左右的殘差對速度的效應(yīng)達到平衡時，速度就會收斂。同時，因為整體的射線路徑變短，導(dǎo)致最終模型速度比實際更小。相反，當聚焦點深度比實際深度大時，整體的射線路徑變長，導(dǎo)致最終模型速度比實際更大。

圖2 聚焦點深度誤差示意圖

圖3 殘差與速度耦合示意圖

圖4是當聚焦點有水平位置誤差時，正確算子與估計算子計算的旅行時信息對比。因為每次的擾動量很小，旅行時的變化量也較小。在算子只有水平誤差時，若估計算子在正確算子的右邊，左邊的射線路徑變長，右邊的射線路徑變短。估計算子計算的左邊旅行時變大，右邊旅行時變小，兩條旅行時曲線會相交。相反的情況下，若估計算子在正確算子的左邊，右邊的射線路徑變長，左邊的射線路徑變短，則算子左邊的旅行時變小，右邊的旅行時變大。

圖4 聚焦點水平誤差示意圖

當聚焦點只有水平誤差時，垂直向上的射線路徑與正確的算子路徑是相當?shù)?假設(shè)地面水平)。同時，在層狀模型中，算子上方的速度橫向變化較小，豎直向上的射線的旅行時與準確算子的旅行時基本相等，所以二者旅行時信息的交點在估計算子的橫向位置附近(圖5)，即在只有水平誤差時，估計算子橫向位置的旅行時殘差約等于零。

圖5 算子橫向誤差放大示意圖

根據(jù)以上分析，算子的擾動方法如下。

(1)因為深度的誤差對速度影響更大，為了避免速度不準確對橫向位置的影響，所以先擾動深度位置與速度模型。深度的擾動量由焦點上方的旅行時殘差決定，即

Δz=Δt×v

(12)

式中：Δz為深度擾動量；Δt為焦點位置上方的旅行時殘差；v為焦點附近的速度值。

(2)深度與速度模型擾動后，進行水平位置的擾動。因為水平位置的誤差與旅行時殘差的關(guān)系不明顯，所以將橫向位置的擾動給定一個范圍，在這個范圍內(nèi)尋找最匹配逆時聚焦算子的坐標。

(3)利用更新后的速度模型與擾動后的焦點進行正演，得到估計旅行時。利用估計旅行時與準確旅行時做差得到旅行時殘差，進行下一次迭代。

在迭代過程中，需要在一定范圍內(nèi)尋找最佳匹配逆時聚焦算子的坐標，本文利用疊加振幅進行判定。根據(jù)射線追蹤得到的旅行時，沿著旅行時曲線軌跡對逆時聚焦算子做曲線積分，就可以得到疊加振幅。疊加振幅越大，說明估計旅行時與逆時聚焦算子更匹配。

在橫向位置的擾動中，需要知道聚焦點在擾動范圍內(nèi)的所有旅行時，需要反復(fù)進行正演計算，計算量很大。本文采用近似的方法求取每次擾動后的旅行時(圖6)，近似公式為

圖6 水平擾動量計算示意圖

(13)

式中θ為射線出射角度。

在沒有約束的情況下，由于代數(shù)迭代算法的特性，旅行時殘差的分配對速度變化非常敏感，會陷入局部極值，并且在新一輪的殘差分配中會加劇這種狀況。由于算子反映的是反射波的單向旅行時，所以可以將同一反射面的算子分為一組，在同一組反射算子中利用插值，構(gòu)成一個反射面。這樣則不用顯式劃分層位或者斷層，只需要利用幾個具有幾何關(guān)系的算子，就可以劃分一個反射界面。通過反射界面定義的地質(zhì)體內(nèi)的速度是不變的，數(shù)值為層中所有網(wǎng)格速度的均值。在模型復(fù)雜的情況下，在分層約束達到收斂后，可以取消約束，以求得模型的細節(jié)特征。這種方法給定了算子的擾動方向與擾動大小，比一般的遍歷尋找更高效。圖7為速度反演的具體實施流程。

圖7 共聚焦點技術(shù)速度反演流程

在利用代數(shù)迭代反演速度模型的過程中，需要給定一個松弛因子。松弛因子如果太大，反演算法可能出現(xiàn)不收斂情況；如果太小，反演算法的收斂會變慢。同時，在擾動共聚焦點位置之前，需要利用殘差與代數(shù)迭代法更新一次速度模型。目的是避免初始速度模型與正確速度模型的差值過大，引起算子位置的過量擾動。并且，在每次更新速度模型后，都需對速度模型進行約束。

2 模型試驗

應(yīng)用圖8a所示低速體速度模型驗證本文方法。圖8b為初始速度模型，由準確速度模型平滑之后獲得。

首先利用共聚焦點技術(shù)求取逆時聚焦算子。以橫向位置為4500m處的逆時聚焦算子為例，初始的共聚焦點如圖8b星號所示，坐標為(4500m，1600m)，初始逆時聚焦算子為在初始模型上以初始共聚焦點為震源進行正演得到。首先利用初始的逆時聚焦算子求取第一個反射界面的逆時聚焦算子，迭代過程如圖9所示，DTS面板的迭代過程如圖10所示圖。因為第一個反射界面的同相軸在DTS面板上與零時很近，迭代三次就收斂了。

圖8 低速體速度模型

圖9 低速體模型第一個界面的聚焦算子迭代過程

圖10 低速體模型第一個界面的DTS面板迭代過程

第三個反射界面的迭代過程如圖11所示，可以看出，迭代四次后DTS面板上的響應(yīng)基本拉平，算子形態(tài)基本固定。考慮到初始算子與第三個反射界面的距離太遠，速度模型的誤差較大，這種計算效率完全可以接受。

圖11 低速體模型第三個反射界面迭代過程

在計算第四層的逆時聚焦算子前，可以充分利用已知信息。第三個反射界面與第四個反射界面距離很近，可以利用第三層的逆時聚焦算子作為第四層的初始逆時聚焦算子，計算結(jié)果如圖12所示。由圖12可以看出，因為初始聚焦算子選取十分接近正確的算子，所以收斂很快，并且在DTS面板上拾取時移也更加方便。

圖12 低速體模型第四個反射界面迭代過程

將求得的聚焦算子作為輸入數(shù)據(jù)，反演速度模型。首先，應(yīng)用初始速度模型對逆時聚焦算子進行單炮偏移，結(jié)果如圖13a所示。橫向位置為4500m處第四個反射界面的逆時聚集算子偏移結(jié)果如圖13b所示，可以看出偏移能量較為發(fā)散。在能量區(qū)利用式(11)提取初始的共聚焦點，如圖14所示。

圖13 低速體模型初始速度逆時聚焦算子偏移結(jié)果

圖14 低速體模型初始共聚焦點分布

利用上述初始共聚焦點進行迭代反演，速度反演的松弛因子設(shè)置為0.5，迭代20次后的結(jié)果如圖15所示，可見焦點位置基本歸位，低速體形狀基本還原，速度與層位反演效果較好。

圖15 低速體模型迭代20次后反演結(jié)果

圖16為橫向4500m處的反演速度曲線與真值的對比，可見低速體的邊界和速度值反演結(jié)果較準確，同時低速體下方的層位與速度已基本復(fù)原。

圖16 低速體模型x=4500m處反演與真實速度曲線的對比

在上述模型試驗中，初始的聚焦算子選取較精確，擾動量較小，為了驗證反演算法的適用性，選取較復(fù)雜的起伏界面速度模型(圖17a)進行試驗。模型共有五層，第二層底面起伏較大，速度分別為3500、3900、4300、4700、5000m/s，模型網(wǎng)格數(shù)為300×600，網(wǎng)格間距為5m。圖17b為初始的速度模型，初始速度模型為水平層狀速度模型平滑后產(chǎn)生。初始的共聚焦點每層設(shè)置29個焦點(圖17中的黑點所示)，均勻分布在界面上，共116個。

圖17 界面起伏速度模型

在求取逆時聚焦算子時，利用橫向位置(750m，800m)處的初始共聚焦點，在初始速度模型上做正演得到初始的逆時聚焦算子。利用初始聚焦算子求取橫向位置為750m處第一個反射界面的逆時聚焦算子，迭代過程如圖18、圖19所示。由于第一個界面反射界面較淺、同相軸清晰，迭代兩次后DTS面板的響應(yīng)被拉平，完成第一界面逆時聚焦算子的求取。

圖18 界面起伏速度模型第一個界面聚焦算子迭代過程

圖19 界面起伏速度模型第一個界面DTS面板迭代過程

利用求取的第一個界面的逆時聚焦算子作為第二界面的初始值，求取第二個界面的逆時聚焦算子，迭代過程如圖20、圖21所示。

圖20 界面起伏速度模型第二個界面聚焦算子迭代過程

圖21 界面起伏速度模型第二個界面DTS面板迭代過程

在求取所有的逆時聚焦算子后，利用初始速度模型與初始共聚焦點，進行反演計算，松弛因子設(shè)為0.5。圖22為迭代35次后的反演速度模型，圖23為橫向750m處的反演與真實速度曲線對比。

圖22 界面起伏速度模型反演結(jié)果

圖23 界面起伏模型x=750m處反演與真實速度曲線對比

可以看出，速度和層位基本還原，最大速度誤差約為10m/s，考慮到初始速度模型與真實模型差距較大，該誤差在可接受范圍內(nèi)。

總之，由上述模型試驗可知，即使在初始速度模型與真實模型差距較大情況下，本文方法求取的逆時聚焦算子穩(wěn)定性高、效率高，利用逆時聚焦算子反演的速度也較精確。

3 結(jié)論

總結(jié)上述成果，可以得出以下結(jié)論：

(1)利用共聚焦點技術(shù)，拾取DTS面板時移求取逆時聚焦算子的方法是可行的。即使初始的逆時聚焦算子與準確的算子相差較大的情況下，也能準確地求取逆時聚焦算子。但是在每次迭代中，需要人工拾取時移，有較大的工作量。

(2)利用逆時聚焦算子與ART反演方法能夠反演出一個可靠的速度模型，并且即使初始速度模型較差的情況下，也能得到較好的結(jié)果。

(3)通過ART的方法反演速度，不需要重復(fù)的偏移計算，并且在反演過程中只考慮反射界面上的多個幾何點，計算量大大降低。

(4)與其他反射層析成像方法相同，本文方法同樣要求反射波信號具有較高的信噪比，可準確拾取共反射點。