王彥飛

中國科學(xué)院地質(zhì)與地球物理研究所中國科學(xué)院油氣資源研究重點實驗室，北京 100029

1 引 言

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

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

本文研究干涉偏移方法和偏移反演成像方法對于地震成像效果的影響，探討二者在提高成像分辨率上的異同.對于偏移反演，基于不適定性考慮以及最小二乘偏移反演計算的局限性，提出應(yīng)用預(yù)條件正則化的偏移反演手段進(jìn)行計算及數(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所示.因此需要對Born算子成像方程做改進(jìn).下面簡單敘述一下兩種改進(jìn)方式：干涉偏移成像和偏移反演成像.

2.2 干涉偏移成像

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

對于震源在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方程，對于VSP幾何成像，在遠(yuǎn)場假設(shè)條件下，根據(jù)互易原理，He等[4］導(dǎo)出了干涉偏移成像公式：

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

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

2.3 偏移反演成像

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

（1）給定接收器和源以及點擴散模型m；

（2）計算點擴散響應(yīng)d=Lm；

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

（4）反演計算：

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

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

或

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

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

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

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

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

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

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

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

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

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

4.1 正則化模型

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

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

則地層反射模型m可以通過

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

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

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

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

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

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

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

直接計算可得：

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

其中，gk為第k次迭代的梯度，即gk=LT（Lmk-d）+αDTDmk；hk為搜索方向；參數(shù)βk為共軛方向參數(shù)因子；νk為步長因子.在本文的計算中，預(yù)條件矩陣P取為對稱矩陣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時，則算法產(chǎn)生的步長可能非常小，導(dǎo)致算法表現(xiàn)變差.為此，我們考慮PRP方法，即選取參數(shù)因子βk為

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

5 數(shù)值模擬

5.1 繞射點模型直接偏移成像及偏移反演成像對比

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

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

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

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

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

6 結(jié) 語

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

致 謝 感謝兩位審稿人對論文的內(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é)報，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ù)值反演問題[全球變化與地球系統(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]崔巖，王彥飛，楊長春.帶先驗知識的波阻抗反演正則化方法研究.地球物理學(xué)報，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]李振華，王彥飛，楊長春.正則化偏移成像的全局優(yōu)化快速算法.地球物理學(xué)報，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]王彥飛.反演問題的計算方法及其應(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]肖庭延，于慎根，王彥飛.反問題的數(shù)值解法[計算方法叢書］.北京：科學(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.