董良國， 黃超， 遲本鑫， 劉玉柱

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

?

基于地震數(shù)據(jù)子集的波形反演思路、方法與應(yīng)用

董良國， 黃超， 遲本鑫， 劉玉柱

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

地震數(shù)據(jù)與地下介質(zhì)物性參數(shù)之間的復(fù)雜關(guān)系，決定了地震全波形反演在理論方法上面臨著強烈的非線性難題.地下不同物性參數(shù)的不同分量在地震數(shù)據(jù)上具有不同的表現(xiàn)，勘探的不同階段對地下介質(zhì)模型的精度也具有不同的要求，這就決定了在地震全波形反演過程中不必時刻追求地震數(shù)據(jù)全部信息的匹配，部分信息的匹配就有可能解決現(xiàn)階段的某些問題，還可以一定程度上規(guī)避匹配全部地震信息所遇到的強烈非線性難題.基于這樣的考慮，我們提出了利用地震數(shù)據(jù)子集進行波形反演的思路，給出了統(tǒng)一的反演方法，并通過基于包絡(luò)數(shù)據(jù)子集以及反射波數(shù)據(jù)子集的波形反演的理論模型與實際資料反演試驗，證明了所提出的波形反演思路和方法的正確性.

全波形反演； 波形反演； 地震數(shù)據(jù)子集； 非線性； 目標(biāo)函數(shù)； 核函數(shù)

1 引言

全波形反演(FWI)是通過使模擬數(shù)據(jù)和觀測數(shù)據(jù)的信息達(dá)到最佳匹配來推測地下介質(zhì)物性參數(shù)的(Tarantola,1984).從理論上講，F(xiàn)WI可以考慮介質(zhì)的諸多屬性、考慮各種波的傳播現(xiàn)象、利用各種地震信息，在理論上是目前利用地震數(shù)據(jù)解決地下介質(zhì)性質(zhì)最徹底的方法.在近10年期間，由于計算能力的提升、實際觀測數(shù)據(jù)質(zhì)量的提高，尤其是逆時偏移技術(shù)的發(fā)展和成功應(yīng)用，推動了FWI的迅猛發(fā)展，掀起了FWI理論方法和實際應(yīng)用研究的熱潮.

然而，可以綜合利用各種地震信息的FWI在理論方法上還面臨著強烈的非線性難題(Virieux and Operto, 2009)，其根本原因在于地震波傳播的復(fù)雜性，即地震數(shù)據(jù)與地下介質(zhì)物性參數(shù)之間的復(fù)雜變化關(guān)系(Jannane et al.,1989；董良國等，2013).另外，目前FWI在實際應(yīng)用中也面臨著諸多實際問題，如觀測孔徑不夠?qū)挕⑷狈Φ皖l信息、震源子波難以確定、初始模型精度不高、復(fù)雜介質(zhì)中地震波傳播的準(zhǔn)確描述困難、信噪比較低，等等.由于上述問題的復(fù)雜性以及諸多客觀條件的限制，決定了FWI目前在實際中應(yīng)用并不普遍，在地震學(xué)領(lǐng)域FWI還沒有替代傳統(tǒng)的走時層析和有限頻層析技術(shù)，在地震勘探領(lǐng)域FWI也沒有替代傳統(tǒng)的地震數(shù)據(jù)處理和解釋流程.

面對目前復(fù)雜的地震勘探問題，在傳統(tǒng)的地震數(shù)據(jù)處理與解釋中，將地下模型變化的高低波數(shù)分量分步解決的流程是很有道理的，這也是將一個復(fù)雜問題分步解決的基本思路和策略.面對理論上的困難以及實際地震資料的復(fù)雜性，再考慮到目前地震勘探要解決的實際問題，目前正在發(fā)展的FWI技術(shù)也應(yīng)該放低身段，面向?qū)嶋H.也就是說，在全波形反演中沒有必要時刻追求全部地震信息的匹配，因為匹配部分信息就有可能解決目前地震成像中的宏觀速度模型建立問題，還可以一定程度上規(guī)避匹配全部地震信息極易遇到的強烈非線性難題.當(dāng)然，匹配部分信息在理論上肯定會降低反演的分辨率，但不切實際的“過高的”模型分辨率并不利于地震成像質(zhì)量的提高.

基于上述考慮，我們提出利用地震數(shù)據(jù)子集進行波形反演的思路與方法.下面，我們給出利用地震數(shù)據(jù)子集進行波形反演的基本思路，并給出基于不同地震數(shù)據(jù)子集的FWI統(tǒng)一的梯度計算公式.最后通過幾種不同的地震數(shù)據(jù)子集進行理論模型與實際數(shù)據(jù)的波形反演試驗，以證明基于不同地震數(shù)據(jù)子集的波形反演思路和方法的正確性.

2 基本思路與方法

2.1 基本思路

常規(guī)FWI是通過使模擬數(shù)據(jù)和觀測數(shù)據(jù)的信息達(dá)到最佳匹配來推斷地下介質(zhì)參數(shù)的.設(shè)xs、xr和x分別為炮點、檢波點和空間任意點的空間位置，時間域FWI最常用的最小平方目標(biāo)函數(shù)為

J(m) =h1(u(m),d,m)

-d(xr,t;xs)]2dtd3x，

(1)

其中，h1是在參數(shù)為m時衡量全部模擬數(shù)據(jù)u和觀測數(shù)據(jù)d之間匹配程度的函數(shù).T為觀測的最大時間，U為物理量u所在的空間(由于本文在時間域定義目標(biāo)函數(shù)和計算地震波場，這時U為實數(shù)空間).物理參數(shù)m(x)∈M，其中M為模型參數(shù)空間，一般為實空間.顯然，h1是從空間U×M到實數(shù)空間R的映射，J是從模型空間M到實數(shù)空間R的映射.

與此相對應(yīng)的描述地震數(shù)據(jù)攝動與模型參數(shù)攝動之間的關(guān)系為

Δd=KΔm，

(2)

其中，Δd、K、Δm分別為數(shù)據(jù)殘差、描述數(shù)據(jù)與模型參數(shù)之間關(guān)系的Fréchet導(dǎo)數(shù)或核函數(shù)、以及模型參數(shù)修改量.

在數(shù)據(jù)與模型參數(shù)之間的關(guān)系比較簡單的情況下，像(2)式這樣將地震數(shù)據(jù)和地下模型參數(shù)均看作一個整體的做法還是可行的.然而，地下模型具有性質(zhì)不同的多個分量(DifferentComponents)，例如高、低波數(shù)分量，而實際觀測地震數(shù)據(jù)中也有性質(zhì)不同的多個成分(DifferentContents)或數(shù)據(jù)子集(DifferentSubsets)，例如折射波、反射波、面波等不同震相，走時、振幅、相位等不同信息，等等.這些不同數(shù)據(jù)子集與模型參數(shù)不同成分之間的非線性程度是不同的(董良國等，2013)：模型的不同成分在不同震相上的體現(xiàn)不同，即使在同一震相上所體現(xiàn)出的信息也不同，例如速度低波數(shù)攝動主要體現(xiàn)在反射波走時上，而速度高波數(shù)攝動主要體現(xiàn)在反射波振幅上；初至波主要受淺中層結(jié)構(gòu)的影響，而淺、中、深層介質(zhì)性質(zhì)都會影響來自深層的反射波.不同數(shù)據(jù)子集與不同模型參數(shù)成分之間的這種復(fù)雜關(guān)系，決定了不同數(shù)據(jù)子集具有不同的反演能力，這就要求我們對它們之間的關(guān)系進行具體分析，從而降低FWI的非線性程度.具體地說，就是將傳統(tǒng)FWI中描述地震數(shù)據(jù)攝動與模型參數(shù)攝動之間的關(guān)系修改為

(3)

其中，Δdi、Δmj、Kij分別為數(shù)據(jù)子集殘差、模型參數(shù)分量的修改量、以及描述不同數(shù)據(jù)子集與不同模型參數(shù)成分之間關(guān)系的子核函數(shù)(Sub-Kernels).

通過上述關(guān)系發(fā)現(xiàn)，不同成分的模型參數(shù)修改量主要是由子核函數(shù)和數(shù)據(jù)子集殘差決定的.目前模型參數(shù)修改量的求取主要有兩種方法：伴隨狀態(tài)法(Tarantola, 1984;Trompetal., 2005;Plessix2006)將數(shù)據(jù)殘差沿核函數(shù)(波路徑)反投影，散射積分法(Zhao, 2005;Chenetal., 2007；Liuetal., 2015)則顯式計算和存儲核函數(shù)并求解反演(層析)方程組.而我們這里提出的基于數(shù)據(jù)子集的波形反演思路，就是根據(jù)勘探的不同階段對地下介質(zhì)模型精度要求的不同，通過分析地震數(shù)據(jù)子集和相應(yīng)的子核函數(shù)的性質(zhì)，以便在反演中將適當(dāng)?shù)牡卣饠?shù)據(jù)子集殘差沿著合理的子核函數(shù)進行反投影，從而決定模型參數(shù)的哪種波數(shù)成分在空間何處進行更新.

核函數(shù)的計算方法和特征分析是有限頻層析(Dahlenetal., 2000;Trompetal., 2005;Liuetal., 2009)、波動方程層析(Woodward, 1992)和全波形反演(Tarantola, 1984;Prattetal., 1998;VirieuxandOperto, 2009)中的核心問題(實際上，三者本質(zhì)上一致的).不同地震數(shù)據(jù)子集所包含的信息差別很大，體現(xiàn)不同地震數(shù)據(jù)子集與模型參數(shù)不同成分之間關(guān)系的子核函數(shù)也具有顯著差異，對于不同的地震信息在反演中當(dāng)然采用不同的子核函數(shù)，進而可以對模型中的高低波數(shù)成分進行解耦(Chietal., 2015).折射波對中淺層模型的背景速度較為敏感，將折射波的數(shù)據(jù)子集殘差沿折射波核函數(shù)反投影，可以較好地更新折射波波路徑覆蓋區(qū)域的背景速度.反射波的波形信息(Xuetal., 2012)或走時信息(Chietal., 2015)可以用來反演中深層的背景速度，前提是將反射波數(shù)據(jù)子集的殘差沿反射波子核函數(shù)(Rabbit-ear)進行反投影.否則，反射波通常只能更新模型深部的高波數(shù)成分.

在對核函數(shù)進行分解以及對不同子核函數(shù)的物理含義進行具體分析的基礎(chǔ)上，我們就可以利用部分地震信息進行波形反演，以降低常規(guī)FWI中匹配全部信息所引起的高度非線性問題.因為匹配整個地震道的標(biāo)準(zhǔn)實在太高，在初始模型還達(dá)不到要求情況下，匹配整個地震道所有信息的難度實在太大，有時也沒有必要(因為我們經(jīng)常要解決的僅僅是宏觀速度模型)，能把地震道的某個主要特征匹配好也許就可以達(dá)到某些要求了(如構(gòu)造成像).也就是說，要拋棄粗放式的對整個地震記錄的匹配方式，而是要仔細(xì)分析不同參數(shù)變化在不同地震數(shù)據(jù)子集上的具體體現(xiàn)，在精細(xì)分析觀測地震數(shù)據(jù)具體特點的基礎(chǔ)上，在FWI的不同階段采用不同的地震數(shù)據(jù)子集.

2.2 基于地震數(shù)據(jù)子集的波形反演方法

在上述思路指導(dǎo)下，下面我們給出基于地震數(shù)據(jù)子集的波形反演的具體方法.

為簡單起見，這里只是考慮最小平方目標(biāo)函數(shù).為此，我們把基于地震數(shù)據(jù)子集波形反演的最小平方目標(biāo)函數(shù)定義為

J(m) =h(u(m),d,m)

-R(d(xr,t;xs))]2dtd3x,

(4)

其中，R(u)和R(d)分別表示對數(shù)據(jù)u和d進行某種數(shù)學(xué)操作，R就是為提取某個地震數(shù)據(jù)子集的操作算子.需要說明的是，算子R是一個抽象的形式，在實際反演中，在不同的反演階段、根據(jù)不同的反演目標(biāo)，R要采用不同的具體形式.顯然，h是衡量模擬數(shù)據(jù)子集R(u(xr,t,m;xs))和觀測數(shù)據(jù)子集R(d(xr,t;xs))之間匹配程度的函數(shù).

設(shè)地震波場u(x,t,m;xs)所滿足的波動方程為

F(u(m),m)=0.

(5)

使(4)式的目標(biāo)函數(shù)最小化的求解過程通常采用局部優(yōu)化方法，其模型更新公式為

mk+1=mk+αkpk，

(6)

由(4)式可得

(7)

(8)

其中，計算伴隨場λ=λ(x,t,m;xs)的伴隨方程為

(9)

其中，星號*表示共軛，方程的右端項為伴隨震源.

可見，利用不同的地震數(shù)據(jù)子集進行波形反演時，梯度計算具有統(tǒng)一的形式，即(8)和(9)式，而且形式與常規(guī)FWI完全一致，伴隨方程(9)的形式也相同.只是面對不同的正問題、不同的目標(biāo)函數(shù)形式、以及不同的地震數(shù)據(jù)子集，計算伴隨場的伴隨震源不同而已.對于如(4)式的最小平方目標(biāo)函數(shù)，其伴隨震源不但取決于接收點處模擬數(shù)據(jù)和觀測數(shù)據(jù)的子集R(u(xr,t,m;xs))和R(d(xr,t;xs))之差，還與選取數(shù)據(jù)子集的方式(即算子R)有關(guān).

至于選取什么樣的地震數(shù)據(jù)子集，取決于實際地震數(shù)據(jù)特征和不同反演階段.通過利用具體的體現(xiàn)數(shù)據(jù)子集和不同模型參數(shù)分量之間關(guān)系的子核函數(shù)來確定反演的梯度，就可以實現(xiàn)分步驟、分尺度的地震波形反演.

3 基于地震數(shù)據(jù)子集的波形反演應(yīng)用舉例

從理論上講，只要能夠從地震數(shù)據(jù)中抽取具有一定意義的信息，均可以利用這些信息構(gòu)成的地震數(shù)據(jù)子集按照上述思路進行波形反演.體現(xiàn)地震數(shù)據(jù)某個特征的地震數(shù)據(jù)子集，它需要具備以下條件：(1)能夠體現(xiàn)地震數(shù)據(jù)的某個或某些變化特點；(2)在物理上和(或)從信號角度具有一定的含義；(3)可以通過一定的數(shù)學(xué)手段進行提取.

實際上，以前的一些學(xué)者提出的一些FWI反演方法和策略，也是不同程度地體現(xiàn)了基于地震數(shù)據(jù)子集進行波形反演的基本想法，這些方法和策略的目的都是通過利用部分信息來降低波形反演的非線性程度.這些數(shù)據(jù)子集可以是：

(1) 地震數(shù)據(jù)中某種特定震相構(gòu)成的數(shù)據(jù)子集.例如，若算子R是選取初至波的波形，相應(yīng)的反演則退化為初至波波形反演(Shenget.al.,2006)、高斯束初至波波形反演(劉玉柱等，2014)；若算子R是選取反射波的波形，相應(yīng)的反演則退化為反射波全波形反演(Xuetal., 2012).

(2) 地震數(shù)據(jù)中特定時間或特定空間觀測孔徑的信息構(gòu)成的數(shù)據(jù)子集.例如，為了降低反演的非線性程度，算子R可以是選取地震記錄某一個時窗波形的操作，這時相應(yīng)的反演就退化為層剝離全波形反演(Wangetal., 2009;Yangetal.,2013)；若算子R是選取不同偏移距的地震數(shù)據(jù)，相應(yīng)的反演就是采用分偏移距反演策略的全波形反演(Singh，1989；董良國等，2013).

(3) 地震數(shù)據(jù)經(jīng)過特定數(shù)學(xué)變換后構(gòu)成的數(shù)據(jù)子集.例如，若算子R是Laplace變換算子，這時上述基于數(shù)據(jù)子集的FWI就退化為Laplace域全波形反演(ShinandCha, 2008)，實際上就是通過Laplace變換算子、主要利用初至波的走時信息，更穩(wěn)健地為下一步的反演提供一個低波數(shù)初始速度模型.若算子R是對地震道進行積分的算子，這時相應(yīng)的反演就退化為基于道積分的全波形反演(Liuetal.，2011；Donnoetal.，2013)，這種方法可以在數(shù)據(jù)缺失低頻或者初始模型不好時為常規(guī)FWI提供一個較好的初始模型.若算子R是通過Hilbert變換提取地震道的包絡(luò)，相應(yīng)的反演就變?yōu)榛诎j(luò)的全波形反演(Chietal.,2013,2014;Wuetal.,2014;Huangetal.,2015)，在地震數(shù)據(jù)缺失低頻信息時，可以利用這種波形反演方法建立一個相對可靠的初始速度模型.如果算子R是把觀測數(shù)據(jù)和模擬數(shù)據(jù)之間的褶積運算以及通過Hilbert變換提取地震道包絡(luò)的運算綜合在一起的一種操作，這時上述基于數(shù)據(jù)子集的FWI就退化為不依賴子波、基于包絡(luò)的全波形反演方法(敖瑞德等，2015)，在數(shù)據(jù)缺失低頻信息、子波未知、先驗?zāi)Ｐ托畔⑷狈r，這種FWI方法可以提供一個較好的初始速度模型.而司空見慣的多尺度分頻反演策略(Bunksetal.，1995；SirgueandPratt，2004；劉國峰等，2012)，實際上也是利用低頻地震數(shù)據(jù)子集構(gòu)造的目標(biāo)函數(shù)非線性程度低的特點，來更穩(wěn)健地建立一個相對可靠的低波數(shù)初始速度模型，以此為基礎(chǔ)再逐步利用高頻成分的地震數(shù)據(jù)子集，從而在反演中逐漸提高反演的分辨率.

另外，地震數(shù)據(jù)中某種特定信息也可以構(gòu)成特定的地震數(shù)據(jù)子集，利用這樣的地震數(shù)據(jù)子集當(dāng)然也可以進行反演.例如，若算子R是利用相關(guān)確定初至波的走時，相應(yīng)的反演則退化為初至波波動方程走時反演(LuoandSchuster，1991)；若算子R是利用相關(guān)確定反射波的走時，相應(yīng)的反演則退化為反射波波形反演(Chietal.，2015).可以預(yù)測，若算子R是選取反射波振幅的某種操作，相應(yīng)的反演則退化為類似AVO技術(shù)的反射波振幅反演.基于這樣的地震波振幅等特定信息進行反演的梯度計算方法能否用(8)式和(9)式統(tǒng)一表達(dá)，目前我們還不能確定，需要作進一步的分析.

至于選取何種地震數(shù)據(jù)子集，取決于反演所處的階段和反演目的，還取決于實際地震資料的特點.如果是為了反演淺層模型，F(xiàn)WI過程中可以只選取初至波數(shù)據(jù)子集；如果是利用FWI為偏移成像提供宏觀速度模型，可以選取體現(xiàn)地震數(shù)據(jù)宏觀信息的數(shù)據(jù)子集，例如走時信息、地震道包絡(luò)、Laplace變換、積分變換、相位信息等；如果成像問題已經(jīng)解決，需要利用FWI進行精細(xì)儲層反演與描述，就需要考慮利用體現(xiàn)數(shù)據(jù)精細(xì)變化特征的地震數(shù)據(jù)子集，因為儲層的變化主要體現(xiàn)在地震數(shù)據(jù)的細(xì)微變化特征上.

我們提出的波形反演是一種概括性的思路，具體可以表現(xiàn)為許多種不同的波形反演方法.根據(jù)所采用的地震數(shù)據(jù)子集的不同，其中的一些方法可以反演模型的長波長分量，而另一些方法可以在初始模型較好情況下反演模型的精細(xì)結(jié)構(gòu).

如果是想通過FWI為成像提供一個具有更高精度的宏觀速度模型，在全波形反演中只需要匹配部分信息就有可能解決這樣的建模問題，同時可以降低反演的非線性程度.由于要反演的是速度變化的低波數(shù)成分，F(xiàn)WI中所選擇的部分信息主要應(yīng)該體現(xiàn)地震數(shù)據(jù)的宏觀變化特征.因此，下面我們通過體現(xiàn)地震數(shù)據(jù)宏觀特征的三個數(shù)據(jù)子集的例子，來說明基于地震數(shù)據(jù)子集的波形反演的思路與方法的有效性.至于不同數(shù)據(jù)子集反演時具體的梯度公式，均可以利用梯度的通式(8)和(9)得到，本文不再贅述.

需要說明的是，以下三種方法主要是利用了地震數(shù)據(jù)的宏觀特性，這些方法的反演結(jié)果當(dāng)然也就體現(xiàn)了模型的大尺度變化特征.從這個意義上說，下面所列舉的三個基于地震數(shù)據(jù)子集的波形反演方法，可以劃歸為速度建模方法的范疇，從三個方法的反演結(jié)果上也可以得出這樣的結(jié)論.但是，與目前通用的建模方法相比(Al-Yahya, 1989;MacKayandAbma, 1992;SymesandCarazzone,1991;LiuandBleistein,1995;SavaandFomel, 2003)，本文列舉的三個方法不需要人工干預(yù)，自動化程度更高.

3.1 基于包絡(luò)地震數(shù)據(jù)子集的波形反演應(yīng)用舉例

在實際地震資料波形反演中，低頻成分缺失與初始模型較差是目前高度非線性FWI面臨的兩個實際難題(二者均會引起Cycle-skipping現(xiàn)象)，同時，這兩個問題又是緊密聯(lián)系在一起的：地震數(shù)據(jù)中低頻成分缺失越嚴(yán)重，F(xiàn)WI就需要更高精度的初始模型.也就是說，如果能夠通過其他途徑建立一個高精度的初始模型，即使地震數(shù)據(jù)中缺失一定的低頻成分，F(xiàn)WI也可能得到一個較好的反演結(jié)果.包絡(luò)信息主要反映了地震道的宏觀特征，可以通過Hilbert變換提取地震道包絡(luò)來構(gòu)造目標(biāo)函數(shù)，從而反演速度模型的低波數(shù)成分.包絡(luò)是地震數(shù)據(jù)子集之一，當(dāng)然可以利用(8)和(9)式計算基于包絡(luò)地震數(shù)據(jù)子集的波形反演的梯度.

圖1和圖2是Marmousi模型試驗結(jié)果.由于地震數(shù)據(jù)中濾掉了5 Hz以下的頻率成分，沿深度呈梯度變化的速度初始模型(圖1a)精度又不高，常規(guī)FWI的反演結(jié)果比較差(圖1b)，反演中局部極值問題嚴(yán)重.而利用同樣的梯度初始模型，基于地震道包絡(luò)的波形反演結(jié)果(圖2a)在反映模型的宏觀變化趨勢方面明顯優(yōu)于圖1b，較好地反映了實際模型的低波數(shù)成分.以圖2a作為初始模型再進行常規(guī)的全波形反演，即使在地震數(shù)據(jù)缺失低頻成分、初始模型也很差的情況下，反演結(jié)果也非常理想(圖2b)，反演結(jié)果的高低波數(shù)成分都比較好.

圖1 梯度初始模型(a)以及傳統(tǒng)FWI反演結(jié)果(b)Fig.1 (a) The initial gradient model and (b) the inverted results of traditional FWI

圖2 利用圖1a的梯度初始模型、基于包絡(luò)數(shù)據(jù)子集的FWI反演結(jié)果(a)，以及利用圖2a作為初始模型的傳統(tǒng)FWI反演結(jié)果(b)Fig.2 (a) The inverted results by FWI of envelope data subset-based using the gradient initial model of Fig.1a, and (b) the inverted results of traditional FWI using the initial model of Fig.2a

圖3 初始速度模型(a)及RTM結(jié)果(a)Fig.3 (a) The initial velocity model and (b) the RTM imaging result

圖4 基于包絡(luò)FWI反演的速度(a)及RTM結(jié)果(b)Fig.4 (a) The inverted velocity model by FWI of envelope data subset-based and (b) the RTM imaging result

圖3是南海某二維測線的初始層速度模型以及逆時偏移結(jié)果，圖4是利用地震道包絡(luò)數(shù)據(jù)子集的FWI反演結(jié)果及相應(yīng)的RTM偏移結(jié)果.可以看出，通過包絡(luò)FWI反演的模型分辨率大幅度提高，呈現(xiàn)出比較精細(xì)的速度橫向變化(圖4a)，成像質(zhì)量顯著改善(圖4b).

3.2 低頻缺失時基于彈性波包絡(luò)數(shù)據(jù)子集的波形反演應(yīng)用舉例

包絡(luò)信息主要反映了地震道的宏觀特征，可以利用包絡(luò)來反演速度模型的低波數(shù)成分.利用聲波方程、基于單分量地震數(shù)據(jù)的包絡(luò)進行縱波速度的單參數(shù)反演，最近取得了較好的效果(Chi et al., 2013,2014; Wu et al.,2014).

在利用多分量地震數(shù)據(jù)進行彈性波全波形反演(EFWI)過程中，同樣存在低頻成分缺失與初始模型較差的難題，對于波長較小的橫波速度而言這個問題更加突出(Brossier et al., 2009；Baeten，2013).為此，我們利用彈性波多分量數(shù)據(jù)的包絡(luò)這個多分量地震數(shù)據(jù)子集，來反演地下介質(zhì)縱橫波速度的長波長分量，為以后的地下參數(shù)的精細(xì)反演提供較好的縱橫波速度初始模型.

圖5是Marmousi2縱波速度的真實模型以及用于EFWI的初始模型.由于兩個分量的地震數(shù)據(jù)中5Hz以下的頻率成分被濾掉，且縱橫波速度均采用沿深度呈常梯度變化的初始模型，所以，常規(guī)的EFWI反演結(jié)果出現(xiàn)嚴(yán)重的局部極值問題(圖6).而利用多分量地震道包絡(luò)信息的反演結(jié)果(圖7)卻較好地恢復(fù)了真實模型的長波長成分，以圖7結(jié)果作為初始模型再進行常規(guī)的EFWI，縱橫波速度均取得了比較好的反演結(jié)果(圖8).可見，盡管多分量地震數(shù)據(jù)中低頻缺失、初始模型精度也很差，但這種兩步法的EFWI反演方法卻得到了非常理想的反演結(jié)果.

從X=7.5 km處反演的縱橫波速度剖面上也可以看到(見圖9)，利用基于彈性波包絡(luò)數(shù)據(jù)子集的EFWI結(jié)果作為初始模型的常規(guī)EFWI反演結(jié)果(圖9中紅線)，在2.2 km深度以上幾乎與真實縱橫波速度(圖9中黑線)完全一致，而常規(guī)EFWI反演結(jié)果(圖9中藍(lán)線)與真實縱橫波速度差距很大，橫波速度尤其明顯，反演過程明顯陷入了局部極值.可見，在多分量地震數(shù)據(jù)中缺失低頻信息以及初始模型不好情況下，利用彈性波包絡(luò)這樣的數(shù)據(jù)子集，可以來反演地下介質(zhì)縱橫波速度的長波長分量，從而為下一步的精細(xì)反演提供較好的縱橫波速度初始模型.

3.3 基于反射波地震數(shù)據(jù)子集的波形反演應(yīng)用舉例

在實際地震資料全波形反演中，要反演中深層的速度結(jié)構(gòu)，必須利用反射波數(shù)據(jù).為了克服匹配全部地震信息的困難，可以通過波形的相關(guān)來確定模擬與觀測的反射波的時差，然后利用各地震道時差定義反演的目標(biāo)函數(shù)，這樣可以有效地反演中深部的宏觀速度變化.需要說明的是，這種基于反射波地震數(shù)據(jù)子集的波形反演方法(RFWI)與常規(guī)的偏移速度分析不同，本文是在數(shù)據(jù)域進行的，并且是一種不需人工干預(yù)的全自動速度反演方法.

圖5 Marmousi-2縱波速度真實模型(a)及常梯度初始模型(b)Fig.5 2D Marmousi-2 true model (a) and linear gradient initial model (b) for P-wave velocity

圖6 多分量地震數(shù)據(jù)中不含5 Hz以下頻率成分時的常規(guī)EFWI反演結(jié)果：(a)縱波速度; (b)橫波速度Fig.6 P-wave (a) and S-wave (b) velocity models obtained by conventional EFWI without low frequency below 5 Hz

圖7 多分量地震數(shù)據(jù)中不含5Hz以下頻率成分時，基于彈性波包絡(luò)數(shù)據(jù)子集的EFWI反演結(jié)果：(a) 縱波速度;(b) 橫波速度Fig.7 P-wave (a) and S-wave (b) velocity models obtained by FWI of elastic envelope data subset-based without low frequency below 5 Hz

圖8 多分量地震數(shù)據(jù)中不含5Hz以下頻率成分時，以圖7作為初始模型的常規(guī)EFWI反演結(jié)果：(a) 縱波速度;(b) 橫波速度Fig.8 P-wave (a) and S-wave (b) velocity models obtained by conventional EFWI using Fig.7 as the initial models without low frequency below 5 Hz

圖9 不同方法反演的縱橫波速度剖面對比(X=7.5 km處). 真實模型(黑實線)、常梯度初始模型(黑虛線)、常規(guī)EFWI反演結(jié)果(藍(lán)線)、利用基于彈性波包絡(luò)數(shù)據(jù)子集的EFWI反演結(jié)果作為初始模型的常規(guī)EFWI反演結(jié)果(紅線)Fig.9 Comparison between velocity profiles extracted from models recovered by conventional EFWI (blue line) and our two-step EFWI method (red line) for Vp and Vs at the location of X=7.5 km. The starting vertical gradient model and the true model are depicted with black dashed and solid lines, respectively

反射波相關(guān)目標(biāo)函數(shù)的構(gòu)建、以及體現(xiàn)反射波與模型長波長分量之間關(guān)系的子核函數(shù)的確定，是這種FWI反演方法的兩個關(guān)鍵.該方法的具體細(xì)節(jié)請見文獻(xiàn)(Chi et al., 2015)，不再贅述，這里只是想通過展示該方法在實際資料反演中的效果，來說明基于反射波地震數(shù)據(jù)子集的波形反演的有效性.

對中國東海一條二維拖纜地震數(shù)據(jù)，平滑后的疊前時間偏移速度分析結(jié)果見圖10a，對應(yīng)的高斯束疊前深度偏移結(jié)果見圖11.以圖10a作為初始模型，通過基于反射波地震數(shù)據(jù)子集進行波形反演，速度模型的更新量見圖10b.利用基于反射波地震數(shù)據(jù)子集波形反演結(jié)果，得到的高斯束疊前深度偏移結(jié)果見圖12.可以看到，成像質(zhì)量較RFWI反演前的圖11得到了大幅度提高，同相軸連續(xù)性顯著改善，許多斷層更加清晰，尤其是2 km深度處的上超以及削截等沉積現(xiàn)象在圖12上更加清晰.在基于反射波地震數(shù)據(jù)子集波形反演后的角度域共成像點道集上(圖13)，同相軸基本拉平，說明了波形反演得到的速度模型較初始模型更加準(zhǔn)確.

圖10 基于反射波地震數(shù)據(jù)子集的波形反演初始模型(a)及反演前后的速度修改量(b)Fig.10 (a) The initial velocity model for FWI of reflection data subset-based, and (b) the difference between the RFWI result and the initial velocity model

4 討論

根據(jù)實際觀測地震數(shù)據(jù)的具體特點，在不同的反演階段，選擇不同的地震數(shù)據(jù)子集，有針對性地構(gòu)造非線性程度更低的目標(biāo)函數(shù)，以實現(xiàn)分步驟、分尺度的全波形反演，其中每個步驟的反演都可以稱之為基于地震數(shù)據(jù)子集的波形反演方法.這是地震數(shù)據(jù)和地下參數(shù)之間復(fù)雜關(guān)系的客觀要求，也是地震數(shù)據(jù)處理與解釋中的不同階段不同反演目標(biāo)的客觀要求.

圖11 基于初始模型的疊前深度偏移結(jié)果Fig.11 The pre-depth migration result using the initial velocity model

圖12 以RFWI反演結(jié)果作為偏移速度模型的疊前深度偏移結(jié)果Fig.12 The pre-depth migration result using the RFWI result

圖13 RFWI前(a)后(b)的角度域共成像點道集Fig.13 The common image gathers with (a) before and (b) after RFWI result

當(dāng)然，有人會問：“這樣做還是全波形反演嗎？”.到底什么是FWI？ 我們現(xiàn)在都承認(rèn)，Tarantola(1984)和Lailly(1983)建立了FWI的理論框架，但他們當(dāng)時都沒有把他們的方法稱作FWI，他們更多地使用數(shù)據(jù)擬合(data fitting)的概念，全波形反演(FWI)一詞只是在1986年才出現(xiàn)(Pan and Phinney ，1986).在Tarantola以第一作者身份撰寫的所有文章中，我們也難以明確見到FWI的影子，但這并不影響他在FWI發(fā)展中的開創(chuàng)性貢獻(xiàn).而Tarantola在地震反演中的主要貢獻(xiàn)就是在數(shù)據(jù)擬合的框架下，搭起了參數(shù)反演和地震偏移成像之間的橋梁.

對地震記錄全部信息的匹配是地震反演的最高理想，但在不同的反演階段，面對諸多難題，明智而現(xiàn)實的作法是根據(jù)實際情況適當(dāng)降低反演標(biāo)準(zhǔn)，不再刻意追求對地震記錄全部信息的匹配，而是把對地震數(shù)據(jù)的部分信息進行匹配看作反演的階段目標(biāo)，這也就是基于地震數(shù)據(jù)子集的波形反演.

因此，將FWI人為定義為對地震記錄全部信息的匹配，必然會極大限制FWI的發(fā)展.用對地震數(shù)據(jù)的部分信息進行匹配來反演介質(zhì)參數(shù)的思路方法來定義FWI，這樣會給FWI更大的空間，會更好地解決實際問題.

5 結(jié)論

(1) 面對FWI面臨的強烈非線性問題以及實際資料中的諸多現(xiàn)實難題，本文提出了基于地震數(shù)據(jù)子集的波形反演思路，并給出了統(tǒng)一的反演方法.對不同的地震數(shù)據(jù)子集，梯度具有統(tǒng)一的計算形式，而且均可以通過地震波的正反向兩次傳播模擬進行計算.只是對不同的正問題、不同的地震數(shù)據(jù)子集以及不同形式的目標(biāo)函數(shù)，反向傳播的伴隨震源不同而已.

(2) 所謂基于地震數(shù)據(jù)子集的波形反演，就是要拋棄FWI中對整個地震記錄的匹配方式，而仔細(xì)研究不同參數(shù)在不同地震數(shù)據(jù)子集上的具體體現(xiàn)，在精細(xì)分析資料的基礎(chǔ)上，在FWI的不同階段采用不同的地震數(shù)據(jù)子集.并將該數(shù)據(jù)子集的殘差沿著合理的子核函數(shù)進行反投影，從而決定模型中的哪種波數(shù)成分在空間何處進行更新.

(3) 即使在初始模型精度不高、地震數(shù)據(jù)低頻缺失的情況下，利用地震道包絡(luò)以及反射波地震數(shù)據(jù)子集，通過構(gòu)建合理的目標(biāo)函數(shù)，可以較好地反演地下速度模型的長波長分量.理論模型以及實際資料的波形反演及成像結(jié)果，充分證明了基于地震數(shù)據(jù)子集的波形反演思路與方法的正確性.

致謝 感謝同濟大學(xué)程玖兵教授和王騰飛博士在高斯束疊前深度偏移方面提供的幫助.

Al Yahya K M. 1989. Velocity analysis by iterative profile migration.Geophysics, 54(6): 718-729. Ao R D, Dong L G, Chi B X. 2015. Source-independent envelope-based FWI to build an initial model.ChineseJ.Geophys. (in Chinese), 58(6): 1998-2010, doi: 10.6038/cjg20150615.

Baeten G, de Maag J W, Plessix R E, et al. 2013. The use of low frequencies in a full-waveform inversion and impedance inversion land seismic case study.GeophysicalProspecting, 61(4): 701-711.

Brossier R, Operto S, Virieux J. 2009. Seismic imaging of complex onshore structures by 2D elastic frequency-domain full-waveform inversion.Geophysics, 74(6): WCC105-WCC118.

Bunks C, Saleck F M, Zaleski S, et al. 1995. Multiscale seismic waveform inversion.Geophysics, 60(5): 1457-1473.

Chen P, Zhao L, Jordan T H. 2007. Full 3D tomography for the crustal structure of the Los Angeles region.BulletinoftheSeismologicalSocietyofAmerica, 97(4): 1094-1120, doi: 10.1785/0120060222.

Chi B X, Dong L G, Liu Y Z. 2013. Full waveform inversion based on envelope objective function. ∥75th Annual International Conference and Exhibition, EAGE, Extended Abstracts.

Chi B X, Dong L G, Liu Y Z. 2014. Full waveform inversion method using envelope objective function without low frequency data.JournalofAppliedGeophysics, 109: 36-46.

Chi B X, Dong L G, Liu Y Z. 2015. Correlation-based reflection full waveform inversion.Geophysics, 80(4): R189-R202, doi: 10.1190/GEO2014-0345.1.

Dahlen F A, Hung S H, Nolet G. 2000. Fréchet kernels for finite-frequency traveltimes—I. Theory.GeophysicalJournalInternational, 141(1):157-174, doi: 10.1046/j.1365-246X.2000.00070.x.Dong L G, Chi B X, Tao J X, et al. 2013. Objective function behavior in acoustic wave full-waveform inversion.ChineseJ.Geophys. (in Chinese), 56(10): 3445-3460, doi: 10.6038/cjg20131020.Donno D, Chauris H, Calandra H. 2013. Estimating the Background Velocity Model with the Normalized Integration Method. ∥ 75th Annual International Conference and Exhibition, EAGE, Extended Abstracts. Huang C, Dong L G, Chi B X, et al. 2015. Elastic envelope inversion using multicomponent seismic data with filtered-out low frequency.AppliedGeophysics,12(3):362-377.

Jannane M, Beydoun W B, Crase E, et al. 1989. Wavelengths of earth structures that can be resolved from seismic reflection data.Geophysics, 54(7): 906-910.

Lailly P. 1983. The seismic inverse problems as a sequence of before stack migration. ∥ Conference on Inverse Scattering Theory and Application, Society of Industrial and Applied Mathematics, Proceedings, 206-220.

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, doi: 10.6038/j.issn.0001-5733.2012.04.030.

Liu J, Chauris H, Calandra H. 2011. The Normalized Integration Method—an alternative to Full Waveform Inversion?. ∥ 17th European Meeting of Environmental and Engineering Geophysics of the Near Surface Geoscience Division of EAGE, Leicester, UK.Liu Y Z, Dong L G, Wang Y W, et al. 2009. Sensitivity kernels for seismic Fresnel volume tomography.Geophysics, 74(5): U35-U46, doi: 10.1190/1.3169600.

Liu Y Z, Xie C, Yang J Z. 2014. Gaussian beam first-arrival waveform inversion based on Born wavepath.ChineseJ.Geophys. (in Chinese), 57(9): 2900-2909, doi: 10.6038/cjg20140915.

Liu Y Z, Yang J Z, Chi B X, et al. 2015. An improved scattering-integral approach for frequency-domain full waveform inversion.GeophysicalJournalInternational, 202(3): 1827-1842.

Liu Z Y, Bleistein N. 1995. Migration velocity analysis: Theory and an iterative algorithm.Geophysics, 60(1): 142-153.

Luo Y, Schuster G T. 1991. Wave-equation traveltime inversion.Geophysics, 56(5): 645-653.

MacKay S, Abma R. 1992. Imaging and velocity estimation with depth-focusing analysis.Geophysics, 57(12): 1608-1622.

Pan G S, Phinney R A. 1986. Full waveform inversion ofp-τseismogram. ∥ 56th Ann. Internat. Mtg. Soc. Expl. Geophys., Expanded Abstracts. http:∥dx.doi.org/10.1190/SEGEAB.5.Plessix R E. 2006. A review of the adjoint-state method for computing the gradient of a functional with geophysical applications.Geophys.J.Int., 167(2): 495-503, doi: 10.1111/j.1365-246X.2006.02978.x.

Pratt R G, Shin C, Hick G. 1998. Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion.GeophysicalJournalInternational, 133(2): 341-362, doi: 10.1046/j.1365-246X.1998.00498.x.

Sava P, Fomel S. 2003. Angle-domain common-image gathers by wavefield continuation methods.Geophysics, 68(3): 1065-1074.

Sheng J M, Leeds A, Buddensiek M, et al. 2006. Early arrival waveform tomography on near-surface refraction data.Geophysics, 71(4): U47-U57. Shin C, Cha Y H. 2008. Waveform inversion in the Laplace domain.Geophys.J.Int., 173(3): 922-931.

Singh S C, West G F, Bregman N D, et al. 1989. Full waveform inversion of reflection data.J.Geophys.Res., 94(B2): 1777-1794.

Sirgue L, Pratt R G. 2004. Efficient waveform inversion and imaging: A strategy for selecting temporal frequencies.Geophysics, 69(1): 231-248.

Symes W W, Carazzone J J. 1991. Velocity inversion by differential semblance optimization.Geophysics, 56(5): 654-663.

Tarantola A. 1984. Inversion of seismic reflection data in the acoustic approximation.Geophysics, 49(8): 1259-1266, doi: 10.1190/1.1441754.

Tromp J, Tape C, Liu Q Y. 2005. Seismic tomography, adjoint methods, time reversal and banana-doughnut kernels.GeophysicalJournalInternational, 160(1): 195-216, doi: 10.1111/j.1365-246X.2004.02453.x.Virieux J, Operto S. 2009. An overview of full-waveform inversion in exploration geophysics.Geophysics, 74(6):WCC1-WCC26, doi: 10.1190/1.3238367.

Wang Y H, Rao Y. 2009. Reflection seismic waveform tomography.J.Geophys.Res., 114(B3): B03304, doi:10.1029/2008JB005916.

Woodward M J. 1992. Wave-equation tomography.Geophysics, 57(1):15-26, doi: 10.1190/1.1443179.

Wu R S, Luo J R, Wu B Y. 2014. Seismic envelope inversion and modulation signal model.Geophysics, 79(3): WA13-WA24.

Xu S, Wang D, Chen F, et al. 2012. Full waveform inversion for reflected seismic data. ∥ 74th Annual International Conference and Exhibition, EAGE, Extended Abstracts. Yang J Z, Liu Y Z, Dong L G. 2013. Time-windowed frequency domain full waveform inversion using phase-encoded simultaneous sources. ∥ 75th Annual International Conference and Exhibition, EAGE, Extended Abstracts.Zhao L. 2005. Fréchet kernels for imaging regional earth structure based on three-dimensional reference models.BulletinoftheSeismologicalSocietyofAmerica, 95(6): 2066-2080, doi: 10.1785/0120050081.

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

敖瑞德, 董良國, 遲本鑫. 2015. 不依賴子波、基于包絡(luò)的FWI初始模型建立方法研究. 地球物理學(xué)報, 58(6): 1998-2010, doi: 10.6038/cjg20150615.

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

劉國峰, 劉洪, 孟小紅等. 2012. 頻率域波形反演中與頻率相關(guān)的影響因素分析. 地球物理學(xué)報, 55(4): 1345-1353, doi: 10.6038/j.issn.0001-5733.2012.04.030.

劉玉柱, 謝春, 楊積忠. 2014. 基于Born波路徑的高斯束初至波波形反演. 地球物理學(xué)報, 57(9): 2900-2909, doi: 10.6038/cjg20140915.

(本文編輯 胡素芳)

Strategy and application of waveform inversion based on seismic data subset

DONG Liang-Guo， HUANG Chao， CHI Ben-Xin， LIU Yu-Zhu

StateKeyLaboratoryofMarineGeology,TongjiUniversity,Shanghai200092,China

The main difficulty in seismic full waveform inversion (FWI) is the strong nonlinearity, which is caused by the complexity of seismic wave propagation. Different components of elastic parameters result in different characteristics of seismic data. Meanwhile, different inversion accuracy is required in different stages of exploration or exploitation. So, it is really not necessary to pursue matching all the seismic information during the inversion. Some problems can be solved by matching part of seismic information and the strong nonlinearity can also be alleviated by this way. According to this consideration, a generalized FWI strategy and method based on seismic data subsets is presented. For different seismic data subsets, the gradients of the misfits have the same form and can be calculated by two modeling applications just as traditional FWI, and the only difference in calculating the gradients is the adjoint sources.During the waveform inversion based on seismic data subsets, the responses of seismic data to different scales of perturbation on different media parameters should be analyzed intensively. Different seismic data subsets should be used in different stages of full waveform inversion. And then the residual of this seismic data subset is back-projected along the reasonable sub-kernels to decide where and which component of the medium parameter need to be updated. As examples, envelope and reflection data subsets are used in FWI with synthetic and real data to prove the validity and effectiveness of our presented FWI strategy and method. Especially, in the absence of low frequency contents, reasonable misfits can be used to recover the background compressional and shear velocity using these FWI methods.

Full waveform inversion; Waveform inversion; Seismic data subset; Nonlinearity; Misfit function; Sensitivity kernel

10.6038/cjg20151024.

Dong L G, Huang C, Chi B X，et al. 2015. Strategy and application of waveform inversion based on seismic data subset.ChineseJ.Geophys. (in Chinese),58(10):3735-3745,doi:10.6038/cjg20151024.

國家油氣重大專項(2011ZX05005-005-007HZ)和國家自然科學(xué)基金(41274116和41474034)項目資助.

董良國，男，同濟大學(xué)海洋與地球科學(xué)學(xué)院教授，博士生導(dǎo)師，主要從事地震波傳播理論與數(shù)值模擬、地震波反演等方面的研究. E-mail:dlg@#edu.cn

10.6038/cjg20151024

P631

2015-01-15，2015-09-06收修定稿

董良國， 黃超， 遲本鑫等. 2015. 基于地震數(shù)據(jù)子集的波形反演思路、方法與應(yīng)用.地球物理學(xué)報，58(10):3735-3745,