變密度聲波方程多參數(shù)全波形反演策略

楊積忠，劉玉柱，董良國

同濟(jì)大學(xué)海洋地質(zhì)國家重點實驗室，上海 200092

1 引言

地震波全波形反演（FWI）是近30年來發(fā)展起來的高分辨率地震數(shù)據(jù)處理方法.它利用地震觀測記錄的波形信息提取地下介質(zhì)的彈性參數(shù)（如速度、密度和彈性模量等）.全波形反演既可以為偏移成像提供精確的背景速度（劉國峰等，2012），也可以應(yīng)用到儲層定量描述等領(lǐng)域（單蕊等，2011；李國平等，2011），在油氣勘探開發(fā)中具有重要作用.

由于全波形反演利用了地震波場的全部信息，理論上來講，應(yīng)該盡可能多的考慮影響波傳播的因素，特別是將大孔徑、寬方位數(shù)據(jù)用于全波形反演時，AVO／AVA效應(yīng)和方位各向異性效應(yīng)等需要特別考慮，這就要求將聲波近似單參數(shù)全波形反演擴(kuò)展到多參數(shù)全波形反演，包括縱橫波速度、密度、吸收衰減效應(yīng)、各向異性參數(shù)等.

然而，各參數(shù)之間的相互耦合，增加了多參數(shù)全波形反演的非線性程度.分析不同參數(shù)之間的相互關(guān)系，研究不同參數(shù)化方式對反演的影響程度，是建立多參數(shù)反演策略的有效途徑.Tarantola（1986）和Mora（1987）指出，在P波震源反射地震波形反演中，彈性模量是一種比較差的參數(shù)化方式；對近偏移距數(shù)據(jù)（短波長成分）來說，阻抗是相對較好的參數(shù)化方式，而對遠(yuǎn)偏移距數(shù)據(jù)（長波長成分）來說，速度是比較好的參數(shù)化方式，但是，不管利用哪種數(shù)據(jù)信息，密度都不能被很好地估計.Kohn等（2012）指出，在彈性波多參數(shù)波形反演中，不管利用哪種參數(shù)化方式，彈性模量、縱橫波速度、阻抗都可以被較好地反演出來，而密度反演結(jié)果較差；但總的來說，利用速度、阻抗參數(shù)化方式得到的反演結(jié)果要優(yōu)于模量參數(shù)化方式得到的反演結(jié)果.Forgues和Lambaré（1997）通過輻射模式分析指出，密度很難利用全波形反演得到.因此，大部分學(xué)者在反演過程中，將密度固定成某一常數(shù)或利用速度通過經(jīng)驗關(guān)系換算得到密度.

密度是巖石物性參數(shù)的重要組成部分，如果能通過反演同時得到比較可靠的縱橫波速度、密度資料，這對儲層評價、巖性解釋和油藏描述等都具有重要作用.彈性波多參數(shù)全波形反演的研究實例證明，雖然可以得到相對可靠的縱橫波速度，但密度并不能被很好地估計（Choi et al.，2008）.因此，需要對密度反演進(jìn)行深入研究，提出有利于密度反演的多參數(shù)全波形反演策略.Jeong等（2012）提出了一種頻率域密度反演策略，首先將密度固定為任意常數(shù)，通過單參數(shù)反演得到彈性模量，進(jìn)而將換算得到的縱橫波速度作為初始模型，再同時反演速度和密度，有效提高了密度的反演精度.Prieux等（2013）先利用大偏移距數(shù)據(jù)單參數(shù)反演速度參數(shù)，再利用全偏移距數(shù)據(jù)同時反演速度、密度參數(shù)，有效提高了密度反演的穩(wěn)定性，但密度反演結(jié)果仍然不夠理想.他們同時指出，利用速度－密度參數(shù)化方式反演得到的速度、密度和阻抗結(jié)果要好于利用速度－阻抗參數(shù)化方式得到的反演結(jié)果.

本文利用速度－密度參數(shù)化方式，通過輻射模式及目標(biāo)函數(shù)敏感度分析，研究了速度、密度在反演過程中的相互影響，提出了一種有利于速度、密度同時反演的策略，并通過球狀異常模型和 Marmousi－2模型的數(shù)值實驗驗證了該策略的有效性.

2 反演方法簡介

全波形反演通過擬合模擬數(shù)據(jù)umod和觀測數(shù)據(jù)uobs，使波場殘差δu達(dá)到最小，進(jìn)而獲取地下介質(zhì)的物性參數(shù)m.目標(biāo)函數(shù)可以表達(dá)為

可以利用擬牛頓法更新模型參數(shù)mn，使目標(biāo)函數(shù)E達(dá)到最?。?/p>

選定某個初始矩陣B0k，對（3）式做m次遞推可得：

因而，只需記錄m（一般取3～20，本文取m＝5）個向量對，就可以按公式（4）構(gòu)造出近似矩陣，具體算法可參考Nocedal and Wright（2006）.

頻率域變密度聲波方程可以表達(dá)為

其中，κ為模量，ρ為密度，ω為角頻率，P（x，ω）為頻率域壓力波場，S（ω）是炮點xs處的震源函數(shù).目標(biāo)函數(shù)對模量和密度的梯度可以表達(dá)為（Tarantola，1984；Pratt and Worthington,1900)

其中，Pf為正傳波場，Pb為反傳殘差波場.梯度可以看成是正傳波場和反傳殘差波場零延遲的互相關(guān)，而各參數(shù)之間輻射模式的不同導(dǎo)致了梯度在具體表達(dá)形式上的差異.

由于模量－密度并不是一種好的參數(shù)化方式（Tarantola，1986；Mora，1987；Kohn et al.，2012），本文選取速度－密度為模型參數(shù).目標(biāo)函數(shù)對其它參數(shù)的梯度可以利用鏈?zhǔn)椒▌t（Mora，1987）求得：

利用模量、速度、密度之間的關(guān)系：

可以得到速度、密度參數(shù)化時對應(yīng)的梯度表達(dá)式：

本文的研究都基于公式（9）進(jìn)行，在具體實驗過程中，選取慢度（速度的倒數(shù)）和密度的倒數(shù)作為模型參數(shù)來更新模型（Mulder and Plessix，2004）.

3 速度、密度同時反演策略

Forgues和Lambaré（1997）在Ray＋Born近似條件下導(dǎo)出了變密度聲波方程各種參數(shù)化方式對應(yīng)的輻射模式的數(shù)學(xué)表達(dá)式，并通過對Hessian矩陣的特征值分解分析指出，密度反演比較困難.在Ray＋Born近似條件下，速度、密度參數(shù)化方式對應(yīng)的輻射模式如圖1所示.

圖1 速度－密度參數(shù)化方式對應(yīng)的輻射模式Fig.1 Radiation pattern as a function of diffraction angle for velocity and density

從圖1可以看出：（1）速度對應(yīng)的輻射模式是各向同性的，密度只有在中小角度情況下散射能量，在大角度情況下幾乎不散射能量，這說明速度對各種角度信息的數(shù)據(jù)都比較敏感，但密度只對中小角度的數(shù)據(jù)敏感；（2）在中小角度情況下，速度、密度的輻射模式基本一致，在這種情況下，速度、密度是耦合在一起的.因此當(dāng)觀測數(shù)據(jù)中只含有中小角度信息時，要將速度和密度完全區(qū)分開來是比較困難的（盡管此時可以得到較好的阻抗信息）.

Jannane等（1989）以反射地震數(shù)據(jù)的L2范數(shù)作為目標(biāo)函數(shù)，通過地震波傳播正演模擬，得到了目標(biāo)函數(shù)隨速度和波阻抗攝動尺度之間的變化關(guān)系.結(jié)果發(fā)現(xiàn)，反射地震數(shù)據(jù)對速度的長波長成分和阻抗的短波長成分敏感，難以利用地震數(shù)據(jù)來反演中等攝動尺度的速度信息.董良國等（2013）在此基礎(chǔ)上，詳細(xì)分析了目標(biāo)函數(shù)隨不同物性參數(shù)（主要為速度和密度）擾動的變化關(guān)系，并指出：（1）地震反射波對密度的較小尺度攝動具有很強(qiáng)的敏感度，而中長波長的密度攝動對反射波基本沒有影響.根本原因在于密度的攝動不改變地震波的走時，而小尺度的密度攝動會引起反射波振幅的變化，但中長尺度的密度攝動不但不能改變地震波走時，對反射波的振幅也基本沒有影響.因此，利用反射波無法反演較大尺度的密度攝動成份；（2）小偏移距反射波可以較好地反演小尺度的波阻抗攝動成份，但要同時反演速度和密度的小尺度攝動可能比較困難，因為小偏移距反射波振幅的變化是小尺度速度攝動和小尺度密度攝動二者共同作用的結(jié)果，兩者的作用是耦合在一起的；（3）由于密度的中長尺度攝動對反射波基本沒有影響，因此，在速度和密度同步擾動的中長尺度攝動區(qū)域，目標(biāo)函數(shù)的變化主要是速度的中長尺度攝動引起的；（4）密度攝動的高頻成份（小尺度攝動）在不同偏移距反射波上都有很好的反映，增大偏移距，反射波對密度的反演能力有稍許提高；（5）隨偏移距的增大，反射波對速度攝動的小尺度成份的反演能力逐漸降低.而大偏移距反射波數(shù)據(jù)中蘊含著豐富的速度的中長波長信息，有利于反演背景速度（盡管非線性程度較高）；（6）大偏移距反射波數(shù)據(jù)中，由密度的小尺度攝動（高頻成份）所引起的目標(biāo)函數(shù)值的變化要遠(yuǎn)遠(yuǎn)小于由速度的中長波長攝動所引起的目標(biāo)函數(shù)值的變化.

從上面的分析可以看出，輻射模式分析和目標(biāo)函數(shù)敏感度分析的結(jié)論具有較好的一致性.如果直接利用全偏移距、大孔徑反射波數(shù)據(jù)對速度、密度進(jìn)行同時反演，可以得到速度的各種尺度的攝動成份和密度的高頻攝動成份，但由速度的中長波長攝動所引起的目標(biāo)函數(shù)的強(qiáng)烈非線性會嚴(yán)重影響密度的反演結(jié)果.前人的研究結(jié)果也表明，直接對反射地震數(shù)據(jù)進(jìn)行速度、密度同時反演，雖然可以得到比較可靠的速度反演結(jié)果，但密度反演結(jié)果較差（Choi et al.，2008；Prieux et al.，2013）.因此，在進(jìn)行速度、密度反演的過程中，可以先利用大偏移距數(shù)據(jù)單參數(shù)反演得到速度的中長波長成份，然后再同時對速度、密度進(jìn)行反演，這是一種比較合理的反演思路.但是，反演過程中無法消除密度的高頻攝動對速度反演造成的影響.而且，在具體的實現(xiàn)過程中，速度、密度波長變化的復(fù)雜度無法與偏移距的大小建立明確的對應(yīng)方式，并且反射地震數(shù)據(jù)中偏移距的大小本來就是一個相對概念，如何有效地將大、小偏移距的數(shù)據(jù)區(qū)分開來成為影響反演結(jié)果的重要因素.

因此，本文提出如下的分步反演策略.第一步，利用速度－密度方程同時反演速度和密度，這樣可以得到比較可靠的速度反演結(jié)果，而密度反演結(jié)果較差.第二步，用第一步反演得到的速度作為初始速度模型，舍棄第一步反演得到的密度，利用最初給定的初始密度作為第二步的密度初始模型，繼續(xù)進(jìn)行雙參數(shù)同時反演，這樣可以同時得到相對可靠的速度、密度反演結(jié)果.為了進(jìn)一步提高反演精度，可以將第二步反演得到的速度、密度結(jié)果作為初始模型，進(jìn)行下一輪的速度、密度聯(lián)合反演.

4 數(shù)值實驗

為了進(jìn)一步說明速度、密度在反演過程中的相互影響，驗證速度、密度分步聯(lián)合反演策略的有效性，進(jìn)行如下數(shù)值實驗.

4.1 速度、密度相互影響

4.1.1 球狀異常模型

本節(jié)數(shù)值實驗基于球狀異常模型.球狀異常體位于均勻背景模型中央，異常體半徑100m.速度模型均勻背景值為3000m·s－1，球狀異常值為3300m·s－1，密度模型均勻背景值為2000kg·m－3，球狀異常值為2200kg·m－3.采用四周觀測方式，使得目標(biāo)體在各個方向的照明均勻，進(jìn)而消除觀測系統(tǒng)的影響.模型縱橫向大小為2km×2km，網(wǎng)格離散間距10m.數(shù)值模擬152炮數(shù)據(jù)，記錄長度1.5s，時間采樣間隔1ms，震源函數(shù)選用主頻7Hz的Ricker子波.

首先，進(jìn)行速度（v）、密度（ρ）單參數(shù)全波形反演.合成觀測數(shù)據(jù)時，需要反演的模型參數(shù)中含有球狀異常體，另外一個模型參數(shù)為均勻背景模型，初始模型為均勻背景模型.反演結(jié)果如圖2所示，圖3顯示的是穿過球狀異常體中心的縱向剖面.從圖中可以看出，在此種觀測方式下，模型各個方向照明均勻，速度、密度反演結(jié)果與真實模型非常接近.

圖2 單參數(shù)全波形反演結(jié)果（a）速度；（b）密度.Fig.2 Mono－parameter full－waveform inversion results for velocity（a）and density（b）.

接下來，進(jìn)行速度、密度同時反演，以考察反演過程中速度、密度之間的相互影響.合成觀測數(shù)據(jù)時，其中一個模型參數(shù)中含有球狀異常體，另外一個模型參數(shù)為均勻背景模型.反演過程中所用的初始速度、密度模型為均勻背景模型.

圖4顯示的是速度含有球狀異常體，密度為均勻背景模型時雙參數(shù)聯(lián)合反演的結(jié)果，圖5為穿過球狀異常體中心的縱向剖面.很顯然，在反演過程中，由速度異常造成的數(shù)據(jù)殘差同時反投影到了速度、密度反演結(jié)果當(dāng)中.從圖4和圖5中可以看出，速度反演結(jié)果與真實模型非常接近，但密度反演結(jié)果卻遠(yuǎn)遠(yuǎn)偏離于真實模型.

圖3 穿過球狀異常體中心的單參數(shù)反演結(jié)果縱向剖面（a）速度；（b）密度.Fig.3 Vertical profiles across the center of the inclusion for（a）P－wave velocity and（b）density

圖4 雙參數(shù)全波形反演結(jié)果（a）速度；（b）密度.Fig.4 Multi－parameter full－waveform inversion results for velocity（a）and density（b）

圖5 穿過球狀異常體中心的雙參數(shù)反演結(jié)果縱向剖面（a）速度；（b）密度.Fig.5 Vertical profiles across the center of the inclusion for（a）P－wave velocity and（b）density

圖6 顯示的是密度含有球狀異常體、速度為均勻背景模型時雙參數(shù)聯(lián)合反演的結(jié)果，圖7為穿過球狀異常體中心的縱向剖面.可以發(fā)現(xiàn)，在反演過程中，由密度異常造成的數(shù)據(jù)殘差也同時反投影到了速度、密度反演結(jié)果當(dāng)中.從圖6和圖7中可以看出，速度、密度反演結(jié)果都與真實模型比較接近，但密度反演結(jié)果分辨率較低，還需進(jìn)一步的反演.將圖6、7與圖4、5進(jìn)行對比，可知速度對密度反演的影響較大，而密度對速度反演的影響相對較小.密度反演依賴于比較可靠的速度信息，而速度的準(zhǔn)確反演也離不開有效的密度信息.

4.1.2 Marmousi－2模型

本節(jié)數(shù)值實驗基于 Marmousi－2模型（Martin et al.，2006），橫波速度為0，僅反演縱波速度與密度.模型縱橫向大小為8km×3.5km，網(wǎng)格離散間距20m.數(shù)值模擬79炮數(shù)據(jù)，第一炮位于地表水平方向100m處，炮間距100m，檢波點分布于地表所有網(wǎng)格點上，數(shù)據(jù)記錄長度5.7s，時間采樣間隔2ms，震源函數(shù)選用主頻7Hz的Ricker子波.為了提高計算效率，實驗過程中利用了相位編碼同步震源技術(shù)（Krebs et al.，2009）.為了降低同步震源相位編碼技術(shù)造成的干涉噪聲，需要改變每個頻率組循環(huán)中相位編碼的隨機(jī)數(shù)（Krebs et al.，2009），以盡量減少同步震源波場的相干性.

（1）密度對速度反演的影響

首先將密度固定，研究密度在速度反演中所起的作用.真實速度模型和高斯平滑初始速度模型如圖8所示.

分別利用以下四種不同密度模型進(jìn)行單參數(shù)（速度）全波形反演：

（1）常密度，密度值為水層密度1010kg·m－3，如圖9a所示；

（2）水層密度1010kg·m－3，水層以下密度值為2000kg·m－3，如圖9b所示；

（3）高斯平滑初始密度，如圖9c所示；

（4）真實密度，如圖9d所示；

速度反演結(jié)果如圖10所示.由圖10a可以看出，當(dāng)密度信息完全不準(zhǔn)時，速度反演結(jié)果中有許多高頻假象，這是因為進(jìn)行單參數(shù)反演時，由密度變化引起的波場殘差全部投影到速度反演結(jié)果當(dāng)中，導(dǎo)致速度反演結(jié)果較差，特別是有強(qiáng)阻抗界面存在時，這種影響更加明顯.由圖10（b—d）可以看出，考慮強(qiáng)阻抗界面的影響時，由密度變化引起的波場殘差對速度反演的影響并不是特別明顯，但這種影響會導(dǎo)致速度反演結(jié)果在某些局部地方過高或過低，進(jìn)而影響速度的后續(xù)使用.圖11顯示了圖10（b—d）的縱向速度剖面，可以更清楚的看到這種影響.由此我們可以得知，利用速度、密度參數(shù)化方式單參數(shù)反演得到的速度結(jié)果并不能滿足需要，特別是將其用于后續(xù)的多參數(shù)全波形反演.

圖6 雙參數(shù)全波形反演結(jié)果（a）速度；（b）密度.Fig.6 Multi－parameter full－waveform inversion results for velocity（a）and density（b）

圖7 穿過球狀異常體中心的雙參數(shù)反演結(jié)果縱向剖面（a）速度；（b）密度.Fig.7 Vertical profiles across the center of the inclusion for（a）P－wave velocity and（b）density

圖8 （a）真實速度；（b）高斯平滑初始速度Fig.8 （a）The true velocity model；（b）the Gauss－smoothed initial velocity

圖9 密度模型（a）常密度，密度值為水層密度1010kg·m－3；（b）水層密度1010kg·m－3，水層以下密度值為2000kg·m－3；（c）高斯平滑初始密度模型；（d）真實密度模型.Fig.9 The density model（a）Constant density model which equals 1010kg·m－3；（b）In sediments，density equals 2000kg·m－3，in water，density equals 1010kg·m－3；（c）The Gauss－smoothed initial density model；（d）True density model.

（2）速度對密度反演的影響

將速度固定，研究速度對密度反演的影響程度.真實密度模型和高斯平滑初始密度模型如圖12所示.

分別利用圖10（b—d）速度模型與圖8a所示真實速度模型進(jìn)行單參數(shù)（密度）全波形反演，密度反演結(jié)果如圖13所示.從圖中可以看出，當(dāng)速度信息準(zhǔn)確時，可以得到比較準(zhǔn)確的密度反演結(jié)果.但是，速度模型的微小變化會引起密度反演結(jié)果的巨大變化，并且有可能導(dǎo)致密度反演結(jié)果遠(yuǎn)遠(yuǎn)偏離真實情況，這說明密度反演依賴于相對可靠的速度信息.圖14顯示了圖13（b—d）的縱向剖面.

4.2 速度、密度分步聯(lián)合反演策略

由上面的分析可知，密度反演依賴于比較可靠的速度信息，而速度的準(zhǔn)確反演也離不開有效的密度信息，單參數(shù)反演得到的速度結(jié)果并不能用于重建密度模型，需要進(jìn)行速度、密度分步聯(lián)合反演.接下來通過如下數(shù)值實驗，證實本文第3節(jié)提出的反演策略的有效性.

圖10 單參數(shù)速度反演結(jié)果（a）利用圖9a密度模型單參數(shù)反演得到的速度結(jié)果；（b）利用圖9b密度模型單參數(shù)反演得到的速度結(jié)果；（c）利用圖9c密度模型單參數(shù)反演得到的密度結(jié)果；（d）利用真實密度模型單參數(shù)反演得到的速度結(jié)果.Fig.10 Mono－parameter full－waveform inversion results for velocity（a）The inverted velocity using the density model shown in Fig.9a；（b）The inverted velocity using the density model shown in Fig.9b；（c）The inverted density using the density model shown in Fig.9c；（d）The inverted velocity using the true density model as shown in Fig.9d.

圖11 水平方向3km（a）、4km（b）、5km（c）處的縱向速度剖面黑線代表真實速度模型，紅線代表初始速度模型，藍(lán)線代表用圖9b密度模型單參數(shù)反演得到的速度結(jié)果，綠線代表用圖9c密度模型單參數(shù)反演得到的速度結(jié)果，粉線代表用真實密度模型單參數(shù)反演得到的速度結(jié)果.Fig.11 Vertical profiles at horizontal positions of 3km （a），4km （b）and 5km （c）of the inverted results shown in Fig.10The black line indicates the true velocity，the red line denotes the initial velocity，the blue line indicates the reconstructed velocity using the density shown in Fig.9b，the green line denotes the reconstructed velocity using the density shown in Fig.9c，the pink line denotes the reconstructed velocity using the true density.

4.2.1 球狀異常模型

本節(jié)數(shù)值實驗基于球狀異常模型，模型參數(shù)和觀測系統(tǒng)設(shè)置與4.1.1節(jié)完全相同.合成觀測數(shù)據(jù)時，真實速度、密度模型中均含有球狀異常體.

圖12 （a）真實密度；（b）高斯平滑密度Fig.12 （a）The true density model；（b）the Gauss－smoothed initial density

圖13 單參數(shù)密度反演結(jié)果（a）利用圖10b速度模型單參數(shù)反演得到的密度結(jié)果；（b）利用圖10c速度模型單參數(shù)反演得到的密度結(jié)果；（c）利用圖10d速度模型單參數(shù)反演得到的密度結(jié)果；（d）利用真實速度模型單參數(shù)反演得到的密度結(jié)果.Fig.13 Mono－parameter full－waveform inversion results for density（a）The inverted density using the velocity model shown in Fig.10b；（b）The inverted density using the velocity model shown in Fig.10c；（c）The inverted density using the velocity model shown in Fig.10d；（d）The inverted density using the true velocity model.

首先，利用均勻背景模型作為初始速度、密度模型，對速度、密度進(jìn)行同時反演.反演結(jié)果如圖15所示，圖16顯示了反演結(jié)果的縱向剖面.從圖中可以看到，速度反演結(jié)果與真實模型具有較好的一致性，但密度反演結(jié)果卻遠(yuǎn)偏離于真實模型，誤差較大.

接下來，將第一步反演得到的速度作為第二階段的初始速度模型，初始密度模型依然利用均勻背景模型，再一次進(jìn)行速度、密度聯(lián)合反演.反演結(jié)果和縱向剖面如圖17和18所示.從圖中可以看到，經(jīng)過第二階段的聯(lián)合反演，可以得到相對比較準(zhǔn)確的速度、密度反演結(jié)果，但密度反演結(jié)果的分辨率不高.

為了進(jìn)一步提高密度反演精度，利用第二階段反演得到的速度、密度作為初始模型，進(jìn)行第三階段的雙參數(shù)聯(lián)合反演.反演結(jié)果和縱向剖面如圖19和20所示，不難發(fā)現(xiàn)，密度反演結(jié)果的分辨率較第二階段有了進(jìn)一步的提升.

4.2.2 Marmousi－2模型

本節(jié)數(shù)值實驗基于 Marmousi－2模型（Martin et al.，2006），模型參數(shù)和觀測系統(tǒng)設(shè)置與4.1.2節(jié)完全相同.為了提高計算效率，實驗過程中利用了相位編碼同步震源技術(shù)（Krebs et al.，2009）.為了降低同步震源相位編碼技術(shù)造成的干涉噪聲，需要改變每個頻率組循環(huán)中相位編碼的隨機(jī)數(shù)（Krebs et al.，2009），以盡量減少同步震源波場的相干性.

圖14 水平方向3km（a）、4km（b）、5km（c）處的縱向密度剖面黑線代表真實密度模型，紅線代表初始密度模型，綠線代表用圖10c速度模型單參數(shù)反演得到的密度結(jié)果，藍(lán)線代表用圖10d速度模型單參數(shù)反演得到的密度結(jié)果，粉線代表用真實速度模型單參數(shù)反演得到的密度結(jié)果.Fig.14 Vertical profiles at horizontal positions of 3km （a），4km （b）and 5km （c）of the models shown in Fig.13The black line indicates the true density，the red line denotes the initial density，the green line denotes the reconstructed density using the velocity shown in Fig.10c，the blue line indicates the reconstructed density using the velocity shown in Fig.10d，the pink line denotes the reconstructed density using the true velocity.

圖15 第一階段雙參數(shù)全波形反演結(jié)果（a）速度；（b）密度.Fig.15 The inverted results in the first stage for velocity（a）and density（b）

圖16 穿過球狀異常體中心的第一階段雙參數(shù)反演結(jié)果縱向剖面（a）速度；（b）密度.Fig.16 Vertical profiles across the center of the inclusion for（a）P－wave velocity and（b）density

圖17 第二階段雙參數(shù)全波形反演結(jié)果（a）速度；（b）密度.Fig.17 The inverted results in the second stage for velocity（a）and density（b）

圖18 穿過球狀異常體中心的第二階段雙參數(shù)反演結(jié)果縱向剖面（a）速度；（b）密度.Fig.18 Vertical profiles across the center of the inclusion for（a）P－wave velocity and（b）density

圖19 第三階段雙參數(shù)全波形反演結(jié)果（a）速度；（b）密度.Fig.19 The inverted results in the third stage for velocity（a）and density（b）

真實速度、密度模型如圖8a和圖12a所示，高斯平滑初始速度、密度模型如圖8b和圖12b所示.首先，利用給定的初始速度、密度模型，對速度、密度進(jìn)行同時反演.反演結(jié)果如圖21所示，圖22顯示了反演結(jié)果的縱向剖面.從圖中可以看到，雖然密度結(jié)果遠(yuǎn)偏離于真實模型，但速度反演結(jié)果與真實模型具有較好的一致性.

接下來，將第一步反演得到的速度作為第二階段的初始速度模型，初始密度模型依然利用最初給定的初始模型，再一次進(jìn)行速度、密度聯(lián)合反演.反演結(jié)果和縱向剖面如圖23和24所示.從圖中可以看到，經(jīng)過第二階段的聯(lián)合反演，可以得到比較準(zhǔn)確的速度、密度反演結(jié)果，并且第二階段速度反演結(jié)果的分辨率有進(jìn)一步的提升.

圖20 穿過球狀異常體中心的第三階段雙參數(shù)反演結(jié)果縱向剖面（a）速度；（b）密度.Fig.20 Vertical profiles across the center of the inclusion for（a）P－wave velocity and（b）density

圖21 第一階段雙參數(shù)反演結(jié)果（a）速度；（b）密度.Fig.21 The inverted results in the first stage：（a）P－wave and（b）density

圖22 第一階段雙參數(shù)反演結(jié)果縱向剖面，分別位于水平方向3km（a1，b1）、4km（a2，b2）和5km（a3，b3）處的速度反演結(jié)果（A）與密度反演結(jié)果（B）黑線代表真實模型，紅線代表初始模型，綠線代表反演結(jié)果.Fig.22 Vertical profiles at horizontal positions of 3km （a1，b1），4km （a2，b2）and 5km （a3，b3）of the models shown in Fig.21for P－wave velocity（A）and density（B）The black line indicates the true model，the red line denotes the initial model and the green line denotes the reconstructed model.

圖23 第二階段雙參數(shù)反演結(jié)果（a）速度；（b）密度.Fig.23 The inverted results in the second stage：（a）velocity；（b）density

圖24 第二階段雙參數(shù)反演結(jié)果縱向剖面，分別位于水平方向3km（a1，b1）、4km（a2，b2）和5km（a3，b3）處的（A）速度反演結(jié)果與（B）密度反演結(jié)果黑線代表真實模型，紅線代表初始模型，綠線代表反演結(jié)果.Fig.24 Vertical profiles at horizontal positions of 3km （a1，b1），4km （a2，b2）and 5km （a3，b3）of the models shown in Fig.23for（A）P－wave velocity and（B）densityThe black line indicates the true model，the red line denotes the initial model and the green line denotes the reconstructed model.

圖25 第三階段雙參數(shù)反演結(jié)果（a）速度；（b）密度.Fig.25 The inverted results in the third stage：（a）velocity；（b）density

圖26 第三階段雙參數(shù)反演結(jié)果縱向剖面，分別位于水平方向3km（a1，b1）、4km（a2，b2）和5km（a3，b3）處的（A）速度反演結(jié)果與（B）密度反演結(jié)果黑線代表真實模型，紅線代表初始模型，綠線代表反演結(jié)果.Fig.26 Vertical profiles at horizontal positions of 3km （a1，b1），4km （a2，b2）and 5km （a3，b3）of the models shown in Fig.25for（A）P－wave velocity and（B）densityThe black line indicates the true model，the red line denotes the initial model and the green line denotes the reconstructed model.

從圖21—24可以看出，進(jìn)行分步多參數(shù)全波形反演，可以得到比較好的速度、密度反演結(jié)果，接下來進(jìn)行第三階段的多參數(shù)聯(lián)合反演，看是否可以進(jìn)一步提高反演精度.此時，利用第二階段反演得到的速度、密度結(jié)果作為第三階段的初始模型.第三階段的反演結(jié)果和縱向剖面如圖25和26所示.不難發(fā)現(xiàn)，密度反演結(jié)果的分辨率較第二階段有了進(jìn)一步的提升，但提升空間有限，主要是由于第二階段密度反演結(jié)果已非常接近真實模型.

5 結(jié)論與討論

反演策略對多參數(shù)全波形反演至關(guān)重要.本文通過數(shù)值實驗證實，密度反演依賴于比較可靠的速度信息，而速度的準(zhǔn)確反演也離不開有效的密度信息.固定密度，單參數(shù)反演得到的速度往往并不能滿足后續(xù)對速度的進(jìn)一步要求.因此，本文提出了有利于速度、密度同時反演的策略.該策略充分考慮了多參數(shù)全波形反演中目標(biāo)函數(shù)對速度、密度的敏感性差異，分步實現(xiàn)對速度、密度的反演，在一定程度上降低了多參數(shù)全波形反演的非線性程度.第一步，利用給定的初始模型對速度、密度進(jìn)行同時反演.此時，可以得到相對可靠的速度反演結(jié)果，而密度反演結(jié)果較差.第二步，將第一步反演得到的速度作為第二步的初始速度模型，舍棄第一步反演得到的密度，利用最初給定的密度作為第二步的初始密度模型，進(jìn)行雙參數(shù)同時反演，這樣可以同時得到比較可靠的速度、密度反演結(jié)果.為了進(jìn)一步提高反演精度，可以將第二步反演得到的速度、密度作為初始模型，繼續(xù)進(jìn)行下一輪的雙參數(shù)同時反演.數(shù)值實驗證實了該策略的有效性.盡管本文是基于變密度聲波方程研究多參數(shù)全波形反演策略，但該策略同樣可以應(yīng)用到彈性波或各向異性的多參數(shù)全波形反演中.

Choi Y，Min D J，Shin C.2008.Two－dimensional waveform inversion of multi－component data in acoustic－elastic coupled media.GeophysicalProspecting，56（6）：863－881.

Dong L G，Chi B X，Tao J X，et al.2013.Objective function behavior in acoustic full－waveform inversion.ChineseJ.Geophys.（in Chinese），56（10）：3445－3460.

Forgues E，LambaréG.1997.Parameterization study for acoustic and elastic ray＋ Born inversion.JournalofSeismic Exploration，6：253－277.

Jannane M，Beydount W，Crase E，et al.1989.Wavelengths of Earth structures that can be resolved from seismic reflection data.Geophysics，54（7）：906－910.

Jeong W，Lee H Y，Min D J.2012.Full waveform inversion strategy for density in the frequency domain.Geophys.J.Int.，188（3）：1221－1242.

Kohn D，De Nil D，Kurzmann A，et al.2012.On the influence of model parametrization in elastic full waveform tomography.Geophys.J.Int.，191（1）：325－345.

Krebs J，Anderson J，Hinkley D，et al.2009.Fast full－wavefield seismic inversion using encoded sources.Geophysics，74（6）：WCC177－WCC188.

Li G P，Yao F C，Shi Y M，et al.2011.Pre－stack wave equation inversion and its application in frequency domain.OGP，46（3）：411－416.

Liu G F，Liu H，Meng X H，et al.2012.Frequency－related factors analysis in frequency domain waveform inversion.ChineseJ.Geophys.（in Chinese），55（4）：1345－1353.

Mart in G S，Wiley R，Marfurt K J.2006.Marmousi2—an elastic upgrade for Marmousi.TheLeadingEdge，25（2）：156－166.Mora P R.1987.Nonlinear two－dimensional elastic inversion of multi－offset seismic data.Geophysics，52（9）：1211－1228.

Mulder W A，Plessix R E.2004.A comparison between one－way and two－way wave－equation migration.Geophysics，69（6）：1491－1504.

Nocedal J，Wright S J.2006.Numerical optimization.Springer Science＋Business Media，LLC.

Prat t R G，Worthington M H.1990.Inverse theory applied to multisource cross－h(huán)ole tomography，Part 1：Acoustic waveequation method.GeophysicalProspecting，38（3）：287－310.

Prieux V，Brossier R，Operto S，et al.2013.Multiparameter full waveform inversion of multicomponent ocean－bottom－cable data from the Valhall field.Part 1：imaging compressional wave speed，density and attenuation.Geophys.J.Int.，194（3）：1640－1664.

Shan R，Bian A F，Yu W H，et al.2011.Pre－stack full waveform inversion for reservoir prediction.OilGeophysical Prospecting，46（4）：629－633.

Tarantola A.1984.Inversion of seismic reflection data in the acoustic approximation.Geophysics，1984，49（8）：1259－1266.

Tarantola A.1986.A strategy for nonlinear inversion of seismic reflection data.Geophysics，51（10）：1893－1903.

附中文參考文獻(xiàn)

董良國，遲本鑫，陶紀(jì)霞等.2013.聲波全波形反演目標(biāo)函數(shù)性態(tài).地球物理學(xué)報，56（10）：3445－3460.

李國平，姚逢昌，石玉梅等.2011.頻率域疊前波動方程反演及其應(yīng)用.石油地球物理勘探，46（3）：411－416.

劉國峰，劉洪，孟小紅等.2012.頻率域波形反演中與頻率相關(guān)的影響因素分析.地球物理學(xué)報，55（4）：1345－1353.

單蕊，卞愛飛，於文輝等.2011.利用疊前全波形反演進(jìn)行儲層預(yù)測.石油地球物理勘探，46（4）：629－633.