地震波干涉偏移及預(yù)條件正則化最小二乘偏移成像方法對(duì)比

王彥飛

中國(guó)科學(xué)院地質(zhì)與地球物理研究所中國(guó)科學(xué)院油氣資源研究重點(diǎn)實(shí)驗(yàn)室，北京 100029

1 引 言

通常的偏移手段對(duì)于垂向地震剖面（VSP）成像效果不是太好，因此人們發(fā)展了很多新的方法進(jìn)行研究，比如近期發(fā)展起來(lái)的干涉偏移方法和偏移反演方法.干涉偏移（Interferometric migration），與光學(xué)中的干涉相關(guān).干涉是光學(xué)物理的概念，通常指的是兩列或兩列以上的波在空間中重疊時(shí)發(fā)生疊加從而形成新波形的現(xiàn)象.Claerbout[1］指出對(duì)自由地表接收到的從底部來(lái)的透射地震記錄進(jìn)行自相關(guān)，等價(jià)于其自激自收模擬記錄（包含負(fù)時(shí)間記錄以及零時(shí)刻的脈沖響應(yīng)）.在忽略儀器響應(yīng)的情況下，地震檢波器記錄到的地震信號(hào)等效于格林函數(shù)和地震子波的卷積.地震波干涉法的核心思想就是對(duì)記錄的地震數(shù)據(jù)進(jìn)行重新校準(zhǔn)，得到以其中一個(gè)檢波器為虛擬震源的新的地震記錄.然后對(duì)虛擬震源求取格林函數(shù)并對(duì)新的地震記錄進(jìn)行偏移運(yùn)算.地震波干涉偏移的目的是為了改善傳統(tǒng)直接偏移方法的成像效果：如展寬成像區(qū)域并提高成像的分辨率.由于干涉思想的采用，該方法可以把弱信號(hào)（低信噪比、多次波）變?yōu)橛杏玫男盘?hào)，也就是所謂的被動(dòng)源成像（Passive seismic imaging）.該方法在最近的10年來(lái)得到了極大的發(fā)展[2-13］.與之相對(duì)應(yīng)，為了改善傳統(tǒng)直接偏移方法的成像效果的另一種方法是地震波偏移及反演成像方法，即在傳統(tǒng)直接偏移方法的基礎(chǔ)上，通過(guò)增加反演迭代的次數(shù)來(lái)達(dá)到高分辨率成像的目的[14-18］.

地震波偏移成像是基于對(duì)Born近似成像算子方程表達(dá)的深化.通?？梢詫懗傻谝活愃阕颖磉_(dá)的形式[19-21］：Lm=d.其中，L為正向模型算子（由核函數(shù)表征）；m為地層反射模型；d為觀測(cè)到的地震信號(hào).在地震勘探領(lǐng)域，試圖通過(guò)求mmig=L＊d=L＊Lm來(lái)得到對(duì)地層的偏移成像，其中L＊為L(zhǎng)的伴隨算子，定義作 （Lx，y）= （x，L＊y）；mmig表示模糊的地層反射偏移圖像；L＊L表示模糊積分核算子，也稱作分辨率函數(shù)或點(diǎn)擴(kuò)展函數(shù).該函數(shù)對(duì)于優(yōu)化勘測(cè)規(guī)劃，比如說(shuō)選取最優(yōu)偏移孔徑，合適的采樣以及選擇偏移參數(shù)起著十分重要的作用.顯然，對(duì)偏移成像作逆運(yùn)算可以得到真實(shí)的反射模型m=（L＊L）-1mmig=（L＊L）-1L＊d.但是，這樣做的缺陷：（1）L＊L是比L 更為不適定的算子，因而 （L＊L）-1對(duì)數(shù)據(jù)誤差／噪音高度敏感；（2）求逆運(yùn)算的計(jì)算量巨大.對(duì)于偏移反演成像，發(fā)展起來(lái)的方法有基于PSF核函數(shù)的反演法[22］，基于CG迭代的最小二乘法[15-16］，非穩(wěn)態(tài)的匹配濾波法[23］，正則化反演成像法和非單調(diào)梯度迭代法[18，24-25］.此外，更多的地球物理反演計(jì)算方法以期用到偏移反演計(jì)算上，如衰減的最小二乘法[26-27］，帶偏差原則的重開始共軛梯度法[28］，全局收斂的信賴域算法和無(wú)記憶擬牛頓法[29-30］，基于Bayes的統(tǒng)計(jì)推理的反演計(jì)算[31-33］以及奇異值展開法[17，19，34］.這些方法有的已 經(jīng)用到 地震偏移反演成像上.但這些方法對(duì)于偏移成像的反演計(jì)算可能收斂的非常慢[35-38］.

本文研究干涉偏移方法和偏移反演成像方法對(duì)于地震成像效果的影響，探討二者在提高成像分辨率上的異同.對(duì)于偏移反演，基于不適定性考慮以及最小二乘偏移反演計(jì)算的局限性，提出應(yīng)用預(yù)條件正則化的偏移反演手段進(jìn)行計(jì)算及數(shù)值模擬.

2 偏移成像原理

2.1 Born算子成像方程

設(shè)Green函數(shù)G（rg｜rs）滿足 Helmholtz方程

其中，d（rg｜rs，ω）表示接收到的反射數(shù)據(jù)；w（ω）為頻率域小波，m（r0）為地層反射函數(shù)，Ω為成像區(qū)間.定義算子L為

則得到Born算子成像方程為

由Born算子成像方程可以得到偏移公式：

但由Born算子成像方程直接做偏移，成像效果較差，如采集腳印和偏移劃弧現(xiàn)象很難消除，如圖1b所示.因此需要對(duì)Born算子成像方程做改進(jìn).下面簡(jiǎn)單敘述一下兩種改進(jìn)方式：干涉偏移成像和偏移反演成像.

2.2 干涉偏移成像

在忽略儀器響應(yīng)的情況下，地震檢波器或地震儀記錄到的地震信號(hào)等效于格林函數(shù)和地震子波的卷積.地震波干涉法的核心思想就是對(duì)記錄的地震信號(hào)進(jìn)行一定的數(shù)學(xué)處理，得到以其中一個(gè)檢波器為震源的新的地震記錄.如果對(duì)地震波場(chǎng)的格林函數(shù)進(jìn)行運(yùn)算后，能產(chǎn)生一個(gè)虛震源記錄，那么在對(duì)地震子波做卷積之后，這種關(guān)系將仍然成立.因此，地震波干涉法的數(shù)學(xué)實(shí)現(xiàn)方式體現(xiàn)在虛擬震源格林函數(shù)的提取方法上[39］.圖2為關(guān)于VSP地震道干涉成像示意圖.

對(duì)于震源在rs、檢波器分別在rg1和rg2的數(shù)據(jù)d（rg1｜rs）和d（rg2｜rs），則可以定義它們的互相關(guān)公式為

其中，φ（·）表示互相關(guān).基于上述的互相關(guān)描述可以導(dǎo)出干涉偏移公式.設(shè)格林函數(shù)G（r0｜rA）和G（r0｜rB）分別滿足相應(yīng)的Helmholtz方程，對(duì)于VSP幾何成像，在遠(yuǎn)場(chǎng)假設(shè)條件下，根據(jù)互易原理，He等[4］導(dǎo)出了干涉偏移成像公式：

其中，Swell表示 VSP測(cè)井，G（r0｜rA）＊表示G（r0｜rA）的復(fù)共軛，rB，rA??S.按照互相關(guān)的定義，可定義

（8）式可以看做是虛擬震源在A點(diǎn)，由在B點(diǎn)的檢波器記錄得到的地震道記錄.因此，公式（7）表明干涉偏移其實(shí)是對(duì)互相關(guān)得到的新的地震記錄做偏移.

2.3 偏移反演成像

與干涉偏移不同的是，偏移反演成像的目的就是通過(guò)對(duì)偏移Green核函數(shù)和偏移圖像作反卷積獲得地球深部反射模型[15，17，22］.但二者的共同目的是提高地震成像的分辨率.由第2.1節(jié)的Born近似，可得記錄和反演合成地震數(shù)據(jù)按以下步驟實(shí)現(xiàn)：

（1）給定接收器和源以及點(diǎn)擴(kuò)散模型m；

（2）計(jì)算點(diǎn)擴(kuò)散響應(yīng)d=Lm；

（3）偏移成像mmig=L＊d=L＊Lm；

（4）反演計(jì)算：

其中，L＊定義作L的伴隨算子.

按照 Yilmaz的思想[40］，偏移反演需要計(jì)算（L＊L）-1L＊d.直接求逆不可取，因此計(jì)算實(shí)現(xiàn)上轉(zhuǎn)化為求解一個(gè)偏移改變量能量最低的最小二乘問(wèn)題

或

對(duì)于上述半正定的線性系統(tǒng)，有很多算法，比如利用改進(jìn)的Cholesky分解法[17］.但是，該計(jì)算方法只是在中等規(guī)模的問(wèn)題時(shí)可以利用，因?yàn)槠溆?jì)算量可以達(dá)到O（1／6n3），這里n為變量的長(zhǎng)度.對(duì)于大規(guī)模的偏移反演計(jì)算，現(xiàn)實(shí)可行的方法應(yīng)該是迭代法.比如對(duì)于最小二乘偏移反演的迭代實(shí)現(xiàn)，Nemeth等[15］和Sj?berg等[16］采用了共軛梯度方法，其中的計(jì)算量主要集中在矩陣-向量積的計(jì)算上.但該方法由于共軛下降算法線性收斂速度的限制以及二次終止特性，將會(huì)導(dǎo)致計(jì)算時(shí)間漫長(zhǎng)，從而影響了成像的分辨率和振幅的保真度，因此有必要研究一些改進(jìn)措施.有關(guān)近期偏移成像振幅補(bǔ)償和偏移反演成像算法的研究工作可見(jiàn)文獻(xiàn)[14-16，18，22-25，41-44］.

3 干涉偏移成像的數(shù)值實(shí)現(xiàn)

Yu和 Schuster以 及 He等[4，12，13］考 慮 了 干 涉偏移公式（7）的數(shù)值實(shí)現(xiàn)問(wèn)題.記Φ（rg1，rg2，ω）為φ（rg1，rg2，t）的頻率域表示，則對(duì)互相關(guān)的地震道φ（rg1，rg2，t）做偏移并對(duì)所有的位置rg1和rg2做疊加可得偏移公式：

其中，τAB指的是波從A點(diǎn)傳播到B點(diǎn)的走時(shí).從上式可以看出，干涉偏移并不需要計(jì)算從震源到檢波器的走時(shí)，但需要計(jì)算從虛擬震源到反射點(diǎn)的走時(shí)，對(duì)于三維偏移來(lái)說(shuō)計(jì)算量巨大.He等[4］考慮了（12）式的經(jīng)濟(jì)型實(shí)現(xiàn)方式.數(shù)值上，上述積分可以做如下實(shí)現(xiàn)

干涉偏移的好處是可以利用表面多次波成像，從而展寬成像區(qū)域提高分辨率.具體計(jì)算過(guò)程如下：

（1）分解垂直地震剖面數(shù)據(jù)為上行波dup和下行波ddown；

（2）通過(guò)濾波獲得直達(dá)波ddirect，下行波ddown減去直達(dá)波ddirect得到多次波dmultiple；

（3）對(duì)直達(dá)波和多次波做互相關(guān)得到虛擬的地表地震剖面數(shù)據(jù)Φ；

（4）對(duì)上述的虛擬地表地震剖面數(shù)據(jù)Φ做偏移.

4 偏移反演成像的迭代正則化

4.1 正則化模型

正則化方法指的是尋找一族正則算子Rα（L，d）使得理論上觀測(cè)數(shù)據(jù)d的噪聲或誤差水平趨于0時(shí)，Rα（L，d）收斂于能夠求出的真實(shí)的地層反射模型m.Tikhonov標(biāo)準(zhǔn)正則化指的是求解一個(gè)極小化的正則泛函[17，34，45-46］：

其中，Ω[·］為Tikhonov穩(wěn)定子；α∈ （0，1）為可調(diào)節(jié)的正則參數(shù)，用于平衡不穩(wěn)定性及光滑性.Ω[m］由用戶根據(jù)問(wèn)題的需要設(shè)定.如果Ω[m］為其中D為離散的對(duì)稱正半定算子，表征對(duì)模型m施加的先驗(yàn)約束，則Rα（L，·）具有如下形式：

則地層反射模型m可以通過(guò)

得到.對(duì)于上述的非線性規(guī)劃問(wèn)題（14）式，有很多算法可以用來(lái)求解[17，34］.但對(duì)于地震偏移反演成像來(lái)說(shuō)，直接求逆計(jì)算量巨大，有必要研究合適可行的優(yōu)化算法.我們研究迭代正則化方法，其迭代形式如下：

若f（ε）恒為 （ 0，1／‖LTL‖）中的常數(shù)，則對(duì)應(yīng)著Landweber-Fridman迭代法；若f（ε）=argminJα（mkεεrk），其中rk=LTLmk-LTd，則對(duì)應(yīng)著最速下降法[17，34］.很明顯，若令初始猜測(cè)值m0=0，并令k=1，f（ε）≡1，則一步梯度迭代法或一步Landweber-Fridman迭代法就是習(xí)慣上經(jīng)常采用的Kirchhoff時(shí)間偏移公式.熟知一步梯度迭代法或一步Landweber-Fridman迭代法遠(yuǎn)遠(yuǎn)沒(méi)有達(dá)到收斂的精度，因此直接的偏移會(huì)導(dǎo)致成像分辨率降低以及振幅恢復(fù)不準(zhǔn)的缺點(diǎn)[18］.

4.2 預(yù)條件共軛梯度迭代方法

下面給出選取正則算子Rα（L，·）的預(yù)條件形式.文獻(xiàn)[25］探討了一系列地震信號(hào)反演中梯度法的預(yù)條件問(wèn)題.這里導(dǎo)出簡(jiǎn)單梯度迭代法和共軛梯度迭代法的預(yù)條件正則化公式.設(shè)P是一個(gè)對(duì)稱正定矩陣，則條件數(shù)cond（P-1LTL）＜cond（LTL）且P-1LTL的譜分布將比LTL的譜分布更集中.令C為一個(gè)非奇異矩陣，并定義分解式P=CCT.注意到公式（14）等價(jià)于極小化下面的一個(gè)二次規(guī)劃問(wèn)題[47］

于是預(yù)條件問(wèn)題可以寫成

其中，A=C-1（LTL+αDTD）C-T，b=C-1LTd，z=CTm.Qα[m］和的梯度分別由下式算得

于是預(yù)條件的最速下降梯度迭代法可以寫成：

直接計(jì)算可得：

共軛梯度法具有二次終止性和較好的穩(wěn)定性，該方法比最速下降方法稍微復(fù)雜一點(diǎn).與最速下降方法不同的是，該方法在每個(gè)迭代步不必均采用最速下降方向，而考慮一個(gè)共軛方向hk-1.即如下的迭代公式：

其中，gk為第k次迭代的梯度，即gk=LT（Lmk-d）+αDTDmk；hk為搜索方向；參數(shù)βk為共軛方向參數(shù)因子；νk為步長(zhǎng)因子.在本文的計(jì)算中，預(yù)條件矩陣P取為對(duì)稱矩陣LTL的逆的逼近.

4.3 共軛方向參數(shù)因子的優(yōu)化選取

上述共軛方向參數(shù)因子βk的選取是基于Fletcher-Reeves（FR）方法，但FR方法并不是數(shù)值表現(xiàn)最好的.記ηk為搜索方向hk與負(fù)梯度方向-gk之間的夾角，并令預(yù)條件矩陣為單位陣，則cos（ηk）=當(dāng)搜索方向與負(fù)梯度方向-gk之間的夾角ηk接近π／2時(shí)，則算法產(chǎn)生的步長(zhǎng)可能非常小，導(dǎo)致算法表現(xiàn)變差.為此，我們考慮PRP方法，即選取參數(shù)因子βk為

當(dāng)搜索方向與負(fù)梯度方向之間的夾角ηk接近π／2并且迭代步長(zhǎng)非常小以及 ‖gk+1-gk‖ ＜ ‖gk‖時(shí)，則有可知下一次迭代的搜索方向?qū)⒖拷?fù)梯度方向-gk+1，從而避免了步長(zhǎng)非常小的情況.考慮到預(yù)條件的情形，參數(shù)因子βk可以取為由于共軛梯度法產(chǎn)生共軛單調(diào)下降方向，該方向可能偏離目標(biāo)值.因此，在多次迭代時(shí)，我們可以考慮重開始技巧，即對(duì)所有的k取β＊k=max｛βk，0｝，其中βk可以為FR的或PRP的.

5 數(shù)值模擬

5.1 繞射點(diǎn)模型直接偏移成像及偏移反演成像對(duì)比

為了驗(yàn)證偏移反演算法在提高成像分辨率和振幅恢復(fù)方面的能力，我們進(jìn)行數(shù)值試驗(yàn).由于d為觀測(cè)值，因此不可避免地帶有各種噪聲，假定該噪聲具有可加性，則

其中，δ為 （0，1）區(qū)間的噪聲水平，rand（size（dtrue））具有與dtrue維數(shù)一致的Gauss隨機(jī)噪聲.dtrue由計(jì)算機(jī)做正向模擬生成.在本次試驗(yàn)中，我們?nèi)≡肼曀溅臑?.05.在構(gòu)造正則算子Rα（L，d）時(shí)，需要選取正則參數(shù).注意到如果α遠(yuǎn)遠(yuǎn)大于1，則將導(dǎo)致逼近問(wèn)題對(duì)原問(wèn)題過(guò)于光滑，從而計(jì)算結(jié)果與原問(wèn)題的解相差甚遠(yuǎn)；相反，如果α遠(yuǎn)遠(yuǎn)小于1，則逼近問(wèn)題未能將算子的譜做很好的改良，因而增加了誤差／噪聲的傳播[34，48］.在本次試驗(yàn)中，我們先驗(yàn)的選取α，即于區(qū)間（0，1）中取定一個(gè)合適的α值，且α＞δ2，這樣算法是可以保證穩(wěn)定收斂的[48］.設(shè)地下埋深有12個(gè)繞射散射點(diǎn)，采用56炮疊加，地震記錄如圖1a所示.直接的偏移成像結(jié)果見(jiàn)圖1b.由圖1b可以看出，直接的偏移對(duì)噪聲控制能力差，成像分辨率質(zhì)量不高.接著用正則化偏移反演進(jìn)行計(jì)算.在預(yù)條件梯度迭代正則化偏移反演中，在噪聲水平為0.05的情況下，采用如下的初始參數(shù)設(shè)置：初始模型賦零向量，正則參數(shù)取為α＞δ2，正半定算子D取為單位矩陣，預(yù)條件矩陣取對(duì)稱矩陣LTL的逆的逼近.成像的結(jié)果見(jiàn)圖3a.同時(shí)畫出了預(yù)條件共軛梯度迭代的誤差分布，見(jiàn)圖3b.為了驗(yàn)證正則偏移反演方法的振幅保真度明顯好于標(biāo)準(zhǔn)的偏移算法，給出振幅對(duì)比如圖4a和圖4b.從圖3和圖4很明顯看出，給出的預(yù)條件正則偏移反演方法具有快速的收斂性，只用了兩次迭代則達(dá)到了很好的收斂效果.同時(shí)，很好地壓制了噪聲并且成像的分辨率和振幅的保真度要比標(biāo)準(zhǔn)的偏移成像結(jié)果好.

5.2 層速度模型干涉偏移成像及偏移反演成像對(duì)比

為了對(duì)比干涉偏移以及偏移反演算法在提高成像分辨率方面的能力，進(jìn)行層速度模型VSP數(shù)據(jù)數(shù)值模擬.設(shè)一個(gè)5層速度模型，各層的厚度和速度是逐漸變化的，原始速度模型見(jiàn)圖5a.在數(shù)據(jù)模擬中，取10個(gè)檢波器放置在豎直井中，并進(jìn)行80炮疊加；因?yàn)榈卣饠?shù)據(jù)通常是帶各種各樣干擾的，因此我們?nèi)≡肼曀溅臑?.05進(jìn)行隨機(jī)干擾，獲得圖5b所示的原始地震記錄.添加的噪聲如圖6a所示.

在預(yù)條件梯度迭代正則化偏移反演中，我們?cè)谠肼曀綖?.05的情況下，采用如下的初始參數(shù)設(shè)置：初始模型賦零向量，正則參數(shù)取為α＞δ2，正半定算子D取為單位矩陣，預(yù)條件矩陣取對(duì)稱矩陣LTL的逆的逼近.標(biāo)準(zhǔn)的偏移成像結(jié)果見(jiàn)圖6b.由圖6b可以看出，直接的偏移對(duì)噪聲控制能力差，成像分辨率質(zhì)量不高.接著分別用干涉偏移方法和預(yù)條件正則化偏移反演進(jìn)行計(jì)算，成像的結(jié)果分別見(jiàn)圖7a和圖7b.很明顯干涉偏移和正則偏移反演方法很好地壓制了噪聲，成像的分辨率要比標(biāo)準(zhǔn)的偏移成像結(jié)果好.對(duì)比圖7a和圖7b發(fā)現(xiàn)，盡管干涉偏移和正則化偏移反演均能展寬成像區(qū)域并提高成像的分辨率，但后者似乎能夠提供更寬的成像區(qū)域和更高的成像分辨率，而前者對(duì)噪聲的壓制要更高一些.

6 結(jié) 語(yǔ)

最近幾十年Kirchhoff偏移在工業(yè)界得到了廣泛的應(yīng)用.但隨著計(jì)算技術(shù)的進(jìn)步以及應(yīng)用的需求，更高質(zhì)量的偏移成像需求越來(lái)越高.本文比較了直接的偏移、正則化偏移反演和干涉偏移三種方法的成像能力對(duì)比.特別是提出了預(yù)條件的梯度迭代正則化偏移反演方法.數(shù)值試驗(yàn)中，采用了帶隨機(jī)噪聲的模擬地震剖面.數(shù)值模擬表明，干涉偏移和正則化偏移反演均能很好地控制噪聲傳播，成像效果較直接的偏移具有更好的分辨率，因此為這兩種方法更廣泛的應(yīng)用提供了有希望的前景.

致 謝 感謝兩位審稿人對(duì)論文的內(nèi)容和格式上提出的良好建議.

（References）

[1]Claerbout J.Synthesis of a layered medium from its acoustic transmission response.Geophysics，1968，33（2）：264-269.

[2]Bakulin A，Calvert R.Virtual source：New method for imaging and 4Dbelow complex overburden.74th Annual International Meeting，SEG，Expanded Abstracts，2004：2477-2480.

[3]常旭，劉伊克，王輝等.地震相干偏移與數(shù)據(jù)自參照偏移的關(guān)系.地球物理學(xué)報(bào)，2009，52（11）：2840-2845.Chang X，Liu Y K，Wang H，et al.Seismic interferometric migration and data-referenced-only migration.Chinese J.Geophys.（in Chinese），2009，52（11）：2840-2845.

[4]He R Q， Hornby B，Schuster G.3Dwave-equation interferometric migration of VSP free-surface multiples.Geophysics，2007，72（5）：S195-S203.

[5]Schuster G T，Yu J，Sheng J，Rickett J.Interferometric／daylight seismic imaging.Geophysical Journal International，2004，157（2）：838-852.

[6]Schuster G T，Katz L，F(xiàn)ollowill F，et al.Autocorrelogram migration：Theory.Geophysics，2003，68（5）：1685-1694.

[7]Snieder R，Grêt A，Douma H，et al.Coda wave interferometry for estimating nonlinear behavior in seismic velocity.Science，2002，295（5563）：2253-2255.

[8]Wang Y B，Luo Y，Schuster G T.Interferometric interpolation of missing seismic data.Geophysics，2009，74（3）：SI37-SI45.

[9]Wapenaar K.Retrieving the elastodynamic Green′s function of an arbitrary inhomogeneous medium by cross correlation.Physics Review Letters，2004，93（25）：254301-1-254301-4.

[10]Wapenaar K，Draganov D，Robertsson J O A，et al.Seismic Interferometry：History and Present Status.Geophysics Reprint Series No.26，Society of Exploration Geophysicists，2008.

[11]Wapenaar K，F(xiàn)okkema J.Green′s function representations for seismic interferometry.Geophysics，2006，71（4）：SI33-SI46.

[12]Yu J H，Katz L，F(xiàn)ollowill F，et al.Autocorrelogram migration：IVSPWD test.Geophysics，2003，68（1）：297-307.

[13]Yu J H，Schuster G T.Crosscorrelogram migration of inverse vertical seismic profile data.Geophysics，2006，71（1）：S1-S11.

[14]Hu J X，Schuster G T，Valasek P A.Poststack migration deconvolution.Geophysics，2001，66（3）：939-952.

[15]Nemeth T，Wu C J，Schuster G T.Least-squares migration of incomplete reflection data.Geophysics，1999，64（1）：208-221.

[16]Sj?berg T A，Gelius L J，Lecomte I.2Ddeconvolution of seismic image blur.Expanded Abstracts，SEG 73th Annual Meeting，Dallas，2003.

[17]王彥飛，斯捷潘諾娃I E，提塔連科V N等.地球物理數(shù)值反演問(wèn)題[全球變化與地球系統(tǒng)科學(xué)系列］.北京：高等教育出版社，2011.Wang Y F，Stepanova I E，Titarenko V N，et al.Inverse Problems in Geophysics and Solution Methods[Series in Global Change and Earth System Science］.Beijing：Higher Education Press，2011.

[18]Wang Y F，Yang C C.Accelerating migration deconvolution using a non-monotone gradient method.Geophysics，2010，75（4）：S131-S137.

[19]Aki K，Richards P G.Quantitative Seismology：Theory and Methods.San Francisco：W H Freeman and Company，1980.

[20]Claerbout J F.Imaging the Earth′s Interior.Oxford：Blackwell，1985.

[21]Stolt R H，Benson A K.Seismic Migration：Theory and Practice，Vol.5，Handbook of Geophysical Exploration，Section I.Seismic Exploration.London：Geophysical Press，1986.

[22]Lecomte I.Resolution and illumination analyses in PSDM：A ray-based approach.The Leading Edge，2008，27（5）：650-663.

[23]Guitton A.Amplitude and kinematic corrections of migrated images for nonunitary imaging operators.Geophysics，2004，69（4）：1017-1024.

[24]Sacchi M D，Wang J，Kuehl H.Regularized migration／inversion：new generation of imaging algorithms.CSEG Recorder，2006，31（S1）：54-59.

[25]Wang Y F.Preconditioning non-monotone gradient method for retrieval of seismic reflection signals.Advances in Computational Mathematics，2012，36（2）：353-376.

[26]Marquardt D W.An algorithm for least-squares estimation of nonlinear parameters.SIAM J.Appl.Math.，1963，11（2）：431-441.

[27]Levenberg K.A method for the solution of certain nonlinear problems in least squares.Quart.Appl.Math.，1944，2：164-166.

[28]崔巖，王彥飛，楊長(zhǎng)春.帶先驗(yàn)知識(shí)的波阻抗反演正則化方法研究.地球物理學(xué)報(bào)，2009，52（8）：2135-2141.Cui Y，Wang Y F，Yang C C.Regularizing method with a priori knowledge for seismic impedance inversion.Chinese J.Geophys.（in Chinese），2009，52（8）：2135-2141.

[29]李振華，王彥飛，楊長(zhǎng)春.正則化偏移成像的全局優(yōu)化快速算法.地球物理學(xué)報(bào)，2011，54（3）：828-834.Li Z H，Wang Y F，Yang C C.A fast global optimization algorithm for regularized migration imaging.Chinese J.Geophys.（in Chinese），2011，54（3）：828-834.

[30]Wang Y F，Yuan Y X.Convergence and regularity of trust region methods for nonlinear ill-posed inverse problems.Inverse Problems，2005，21（3）：821-838.

[31]Tarantola A.Inverse Problems Theory：Methods for Data Fitting and Model Parameter Estimation.Amsterdam：Elsevier，1987.

[32]Ulrych T J，Sacchi M D，Woodbury A.A Bayesian tour to inversion.Geophysics，2000，66：55-69.

[33]楊文采.地球物理反演和地震層析成像.北京：地質(zhì)出版社，1989.Yang W C.Geophysical Inversion and Seismic Tomography（in Chinese）.Beijing：Geological Publishing House，1989.

[34]王彥飛.反演問(wèn)題的計(jì)算方法及其應(yīng)用[當(dāng)代科學(xué)前沿論叢］.北京：高等教育出版社，2007.Wang Y F.Computational Methods for Inverse Problems and Their Applications[New Frontiers of Sciences］（in Chinese）.Beijing：Higher Education Press，2007.

[35]Nolet G.Solving or resolving inadequate and noisy tomographic systems.J.Comp.Phys.，1985，61（3）：463-482.

[36]Nolet G，Snieder R.Solving large linear inverse problems by projection.Geophys.J.Int.，1990，103（2）：565-568.

[37]VanDecar J C，Snieder R.Obtaining smooth solutions to large，linear，inverse problems.Geophysics，1994，59（5）：818-829.

[38]Trampert J，Leveque J J.Simultaneous iterative reconstruction technique：physical interpretation based on the generalized least squares solution.Journal of Geophysical Research，1990，95（B9）：12553-12559.

[39]陶毅，符力耘，孫偉家等.地震波干涉法研究進(jìn)展綜述.地球物理學(xué)進(jìn)展，2010，25（5）：1775-1784.Tao Y，F(xiàn)u L Y，Sun W J，et al.A review of seismic interferometry.Progress in Geophysics（in Chinese），2010，25（5）：1775-1784.

[40]Yilmaz O.Seismic Data Analysis.Investigations in Geophysics No.10.Society of Exploration Geophysicists，Tulsa，Okla，2001.

[41]Chen J，Schuster G T.Resolution limits of migrated images.Geophysics，1999，64（4）：1046-1053.

[42]Schuster G T.Green′s functions for migration.Abstracts of Expanded Abstracts 67th SEG meeting，1997：1754-1758.

[43]Schuster G T，Hu J X.Green′s function for migration：Continuous recording geometry.Geophysics，2000，65（1）：167-175.

[44]Yu J H，Hu J X，Schuster G T，Estill R.Prestack migration deconvolution.Geophysics，2006，71（2）：S53-S62.

[45]Tikhonov A N， Arsenin V Y.Solutions of Ill-posed Problems.New York：John Wiley and Sons，1977.

[46]Engl H W，Hanke M，Neubauer A.Regularization of Inverse Problems.Dorcrecht：Kluwer Academic Publishers，1996.

[47]袁亞湘.非線性規(guī)劃數(shù)值方法.上海：上?？茖W(xué)與技術(shù)出版社，1993.Yuan Y X.Numerical Methods for Nonlinear Programming（in Chinese）.Shanghai：Shanghai Science and Technology Press，1993.

[48]肖庭延，于慎根，王彥飛.反問(wèn)題的數(shù)值解法[計(jì)算方法叢書］.北京：科學(xué)出版社，2003.Xiao Y F，Yu S G，Wang Y F.Numerical Methods for the Solution of Inverse Problems[Series in Numerical Methods］（in Chinese）.Beijing：Science Press，2003.