一種基于動態(tài)目標泛函的全波形反演多解性評估方法

黃建平, 崔 超

(1.中國石油大學地球科學與技術學院,山東青島 266580;2.海洋國家實驗室海洋礦產資源評價與探測技術功能實驗室,山東青島 266071)

一種基于動態(tài)目標泛函的全波形反演多解性評估方法

黃建平1,2, 崔 超1,2

(1.中國石油大學地球科學與技術學院,山東青島 266580;2.海洋國家實驗室海洋礦產資源評價與探測技術功能實驗室,山東青島 266071)

全波形反演是一種精度較高的地下速度建模方法,但其存在計算量巨大、反演多解性、初始模型依賴性強、適用條件苛刻等困難,生產實用化仍然較為困難。將基于動態(tài)目標泛函的全波形反演方法應用于反演解的討論中,利用已有的反演解作為先驗信息,通過動態(tài)調整目標泛函,在保證地下介質模型能夠準確解釋地震數據的前提下得到盡可能不同的反演結果。進一步根據數據約束強弱差異選取不同的反演策略,新方法能夠對數據覆蓋較弱區(qū)域的反演結果進行有效評估。在實現算法和反演流程的基礎上,通過典型海底模型進行試算。結果表明，新方法能夠在保證解的準確性前提下對全波形反演的特征給出合理分析,相對于傳統(tǒng)方法計算效率也有了明顯提高。

全波形反演; 多解性; 目標泛函; 可靠性評估

全波形反演作為目前最有發(fā)展?jié)摿Φ乃俣葓鼋Ｊ侄沃籟1-2],在理論發(fā)展和實際應用中受到了廣泛關注[1,3-9]。然而由于地球物理理論不完善[10-14]、野外數據質量差[15-17]、初始模型構建困難[18-20]等多方面因素的影響,全波形反演結果往往不夠理想。對于全波形反演解的討論和評估成為近年來的研究重點之一。前人提出了全局尋優(yōu)法[21-26]、概率分析法[27]、“Bootstraping”法[28]等反演解的分析方法。由于地球物理問題具有數據量大、計算量大等特征[29],以及波形反演問題本身所具有的復雜性[30-31],傳統(tǒng)方法往往需要極大的計算、存儲量,并且適用條件苛刻[30],目前缺少一種快速高效的全波形反演問題解的評估方法。為此,筆者引入一種基于動態(tài)目標泛函的全波形反演解的評估方法[29],根據數據對模型敏感程度不同給出不同的求解策略,并對該方法的適用性及優(yōu)缺點進行討論。

1 方法原理

1.1 傳統(tǒng)的全波形反演理論

全波形反演通過對目標泛函的優(yōu)化達到求解地下最優(yōu)模型的目的。本文中將目標泛函表示為χ(m),目標泛函可以表示為空間目標泛函核函數的積分形式,即

χ(m)=?L(m)dxdt=〈L(m)〉.

(1)

式中,L為χ(m)的核函數;m為模型參數;x為空間坐標;t為時間采樣點。利用鏈式法則,目標泛函關于模型的一階導數可以轉化為其對于數據的一階導數與數據關于模型的一階導數的乘積,即

(2)

式中,u為正演數據;δm為模型擾動量;δu=muδm為數據對于模型的敏感核函數。對于該核函數的求解,可利用有限差分法對模型的不同參數進行擾動,然而由于全波形反演問題中的參數數目眾多,并且正演的計算量巨大,現階段上式的有限差分解法應用較為困難。利用伴隨狀態(tài)法[32]求解該核函數是一種較為有效的策略,將波動方程表示為

B(u)=s.

(3)

式中,B為正演算子;s為震源項。對式(3)關于模型參數m求導有

(4)

(5)

(6)

(7)

(8)

(9)

(10)

同時,正演算子的初值條件轉化為伴隨算子的終值條件。式(10)中震源波場的伴隨場為觀測數據與正演數據殘差的反傳波場。利用上述結論,公式(2)可以簡化為

(11)

(12)

由式(12)可見,目標泛函關于模型參數的導數項可以表示為震源正傳波場與地震數據殘差反傳波場的零偏移互相關。因此目標泛函關于模型參數導數項的求解可以簡化為兩次正演過程,一次為震源波場的正傳,另一次為檢波器波場殘差的反傳。

全波形反演通過對模型參數的迭代更新求得最優(yōu)解,其迭代公式可以表示為

mk+1=mk-αkHkmχ.

(13)

式中,k為迭代次數;Hk為對梯度的修正項,將模型參數m表示為列向量，Hk為一大規(guī)模矩陣,可采用共軛梯度法、L-BFGS方法等;αk為步長,可由拋物插值方法獲得[7]。

1.2 基于動態(tài)目標泛函的全波形反演

傳統(tǒng)的全波形反演理論往往只能給出一個反演解,即通過對目標泛函的逐步尋優(yōu)獲得能夠解釋數據的最優(yōu)模型,但是由于地球物理反演問題本身所具有的特征,能有效解釋同一觀測數據的介質模型往往有多種,主要原因包括:

(1)地球物理反演問題本身是一個欠定問題,由于觀測數據的有限性和對地下介質先驗知識的局限性,能夠解釋同一觀測數據的模型往往有多個。

(2)全波形反演問題是一個強的非線性問題,目標泛函的性態(tài)較為復雜[20],存在較多的局部極小,對于全局極小的查找極為困難。因此全波形反演的反演結果往往是眾多局部極值之一。

(3)由于觀測系統(tǒng)照明的限制,地下介質存在對于觀測數據沒有影響或者影響極小的區(qū)域,因此目標泛函關于模型的導數矩陣在這些位置為0或者接近于0,即模型參數在該區(qū)域取任意值都可以對數據做較好地解釋,全波形反演方法對此模型參數不具有反演能力或者反演能力很微弱。

Rawlinson等[29]通過修改目標泛函,將已有的反演結果作為先驗信息,在傳統(tǒng)目標泛函的基礎上添加動態(tài)正則化項,在保證模型能夠有效解釋觀測數據的基礎上,盡量保證每次反演所得到的反演解不同。通過求解最后多組最優(yōu)解的均值與均方差,從而對全波形反演多解性問題進行分析,并對反演結果的可靠性進行定量評估。修改后的目標泛函可以表示為

χ1=χ+S(m).

(14)

其中

式中,S(m)為動態(tài)正則化項;Mj為第j個已有的反演解;λ、p、ε為正則化項中人為控制參數。對于上述目標泛函關于模型參數求導可得

(15)

其中

由公式(15)可見,修改后的目標泛函在已有的反演結果位置處生成局部極大值,能夠保證后續(xù)的反演結果與已有反演解不同。

1.3 目標泛函的性態(tài)分析與反演策略

圖1給出了一個一維不同參數情況下的動態(tài)正則化項性態(tài),令M=50,m在1～100變化。分析可見,λ主要決定動態(tài)正則化項的分布范圍,取值過大則動態(tài)正則化項只能在已有解很小的鄰域內進行約束,導致不同反演解之間十分接近;而取值過小會使動態(tài)正則化項的作用范圍過大,影響目標泛函整體性態(tài),甚至導致反演失敗。p的作用與λ相似,在實際應用中往往取值為1。ε主要決定動態(tài)正則化項在已有解位置處的峰值大小,其取值越大對應的峰值越小。本文中每次反演均從初始模型開始,并令ε=0?？偟膩碚f,選取參數的原則為:正則化項的作用既要具有一定區(qū)域性,又能夠有效改變目標泛函在已有解位置處的性態(tài)。

圖1 不同參數取值情況下的動態(tài)正則化項的分布性態(tài)Fig.1 Regularization character of dynamic misfit function with different parameters

需要著重指出的是:在數據對模型變化不敏感的區(qū)域,即反演前后模型變化極小的區(qū)域,正則化項在初始模型位置生成一極大值,修正后的目標泛函在該位置處關于模型參數的導數始終為0,因此不能生成多個最優(yōu)解。然而實際上模型參數在該區(qū)域取任意值對觀測數據都可以作較好地解釋,為了能夠有效解決這一問題,在生成第一個最優(yōu)解時查找出反演前后模型參數變化小于某一閾值的區(qū)域,認為數據在該區(qū)域對于模型的敏感性極小。為保證修正后的目標泛函在該區(qū)域關于模型參數的導數項不為0,可在反演過程中對Mj在該區(qū)域添加一定量的隨機擾動,以保證每次所得到的最優(yōu)解不同。

1.4 基于動態(tài)目標泛函的全波形反演流程

本文中利用已有的反演解作為先驗信息，通過動態(tài)地修改目標泛函以獲得盡可能不同的反演解。其反演的流程如圖2所示。其中k表示內部循環(huán)的迭代次數(式(13))，j表示不同的外部循環(huán)次數。將終止條件(1)表示為目標泛函小于給定閾值或反演迭代次數k達到給定上限。實際上，每次內部循環(huán)都是一次常規(guī)意義的全波形反演過程，在每次內部循環(huán)結束時，返回一個反演結果Mj，將其作為已知反演解代入動態(tài)目標泛函，然后根據判斷條件(2)決定是否繼續(xù)生成反演解。由于每次反演都利用初始模型m0作為輸入，因此內部循環(huán)開始時令k=0。將終止條件(2)表示為所得解的個數達到給定上限。

圖2 基于動態(tài)目標泛函的全波形反演方法流程Fig.2 Flowchart of full waveform inversion by using dynamic function

2 模型測試

為了測試本文中方法的適用性,采用如圖3(a)所示的海底地下模型。其中模型縱橫向尺寸為4 km×8 km, 每100 m做一次正演模擬,共80炮,在海水表面從0開始每20 m放置一個檢波器。地震子波采用主頻為15 Hz的雷克子波,反演過程中假定震源子波和水層速度已知。對水層以下采用三角平滑方式得到反演的初始模型(圖3(b))。

圖4為采用傳統(tǒng)的全波形反演方法得到的第50、100、250次迭代的反演,結果可見:由于給出的初始模型較為準確,隨著迭代次數的增加,海底模型細節(jié)的刻畫逐漸趨于明顯,高波數成分逐漸被恢復。圖5抽取了位于2 km位置處的縱向速度剖面對比,可見隨著迭代次數的增加,地下介質分界面準確歸位,尤其是位于2～2.2 km位置處的低速體被有效反演。但是由于觀測系統(tǒng)的限制,全波形反演對模型的右下側區(qū)域反演能力不足,反演結果較為模糊。而在模型的左側,由于深層傾斜高速體的強反射作用受到較強照明,因而反演結果較好。在深層高速基地區(qū)域由于不能接收到有效的反射信息,模型更新較弱。

圖3 反演所用的準確速度場與初始速度場Fig.3 True model and initial model used by inversion

圖4 常規(guī)全波形反演不同迭代次數后的反演結果Fig.4 Inversion results after traditional full waveform inversion for different iteration times

圖5 由圖3和圖4中模型中抽取2 km位置處的單道對比Fig.5 Single vertical velocity profile comparision from Fig.3 and Fig.4 at 2 km

為了能夠定量地衡量全波形反演對不同模型區(qū)域的反演能力,并對反演模型可靠性進行評估,本文中采用基于動態(tài)目標泛函的全波形反演方法進行反演。共得到10個最優(yōu)解,其中參數的設置為:λ=7×10-7,p=1,ε=0。由于動態(tài)目標泛函反演的目的在于對反演的多解性和結果的可靠性做評估,并非得到十分準確的反演結果,且基于較少次數的迭代就可以有效地反映出全波形反演對模型的更新能力。為了有效的降低計算量,每次反演過程只迭代20次。對反演結果進行加權平均得到的平均模型如圖6所示,可見加權平均后的結果對模型仍具有一定的恢復能力,淺層層位更為清晰。對于10次基于動態(tài)目標泛函所得到的反演結果做均方差分析,如圖7所示。淺層數據約束較強的中間區(qū)域(區(qū)域A)多次反演結果的均方差較小。表明全波形反演在數據覆蓋強的區(qū)域具有穩(wěn)定的反演能力,所得到的反演結果可靠性也比較高,而在深層基底(區(qū)域B)及模型右下側數據照明較弱區(qū)域(區(qū)域C)均方差較大,說明不同反演次數得到的反演結果相差較大。觀察圖8可見不同反演結果(Mj)對應的整體數據殘差變化較小,說明深層基底(B)以及數據照明較弱區(qū)域(C)模型參數的變化對數據影響不大,全波形反演在該區(qū)域反演能力不足且反演結果不可靠。對于這一現象的具體解釋為:對于深層基底(B),由于不能接收穿過該區(qū)域的反射波,因此該區(qū)域的速度場變化對于數據幾乎不產生影響;對于模型右下側照明較弱區(qū)域(C),由于缺少較強的反射界面以及觀測系統(tǒng)的限制,該區(qū)域速度的變化對于目標泛函的影響較小。

圖6 利用動態(tài)的目標泛函得到的10個反演結果的平均值Fig.6 Average of 10 inversion results using full waveform inversion with dynamic misfit function

圖7 利用動態(tài)的目標泛函得到的10個反演結果的歸一化后的均方差Fig.7 Mean square deviation of 10 inversion results using full waveform inversion with dynamic misfit function

圖8 利用動態(tài)的目標泛函得到的10個反演結果對應的歸一化后的目標泛函Fig.8 Normalized value of misfit function corresponding to 10 inversion results using dynamic misfit function

3 結 論

(1)在數據約束較強的區(qū)域,通過修改目標泛函在已知解位置處的性態(tài)能夠保證多次反演所得到的反演結果不同。

(2)在數據約束較弱的區(qū)域,通過對已有反演解的強制擾動,保證每次所得到的解均不同,從而對該區(qū)域的反演結果可信度做出評估。

(3)該方法不需要對后驗協(xié)方差矩陣求解,也避免了“Bootstraping”方法只能適用于超定問題的局限性,只需要在不同的求解過程中修改目標泛函,因此本文方法具有較強的適用性。

[1] TARANTOLA A. Inversion of seismic reflection data in the acoustic approximation [J]. Geophysics, 1984,49(8):1259-1266.

[2] PLESSIX R E. A review of the adjoint-state method for computing the gradient of a functional with geophysical applications [J]. Geophys J Int, 2006,167(2):495-503.

[3] PRATT R, SHIN C, HICKS G. Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion [J]. Geophys J Int, 1998,133(2):341-362.

[4] SHIN C, YOUNG H C. Waveform inversion in the Laplace domain [J]. Geophys J Int, 2008,173(3):922-931.

[5] SHIN C, YOUNG H C. Waveform inversion in the Laplace-Fourier domain [J]. Geophys J Int, 2009,177(3):1067-1079

[6] 楊積忠,劉玉柱,董良國.變密度聲波方程多參數全波形反演策略[J].地球物理學報,2014,57(2):628-643.

YANG Jizhong, LIU Yuzhu, DONG Liangguo. A multi-parameter full waveform inversion strategy for acoustic media with variable density [J]. Chinese Journal Geophysics, 2014,57(2):628-643.

[7] K?HN D. Time domain 2D elastic full waveform tomography [D]. Kiel: Kiel University, 2011.

[8] MORA P. Nonlinear two-dimensional elastic inversion of multi-offset seismic data [J]. Geophysics, 1987,52(9):1211-1228.

[9] VIRIEUX J, OPERTO S. An overview of full-waveform inversion in exploration geophysics [J]. Geophysics, 2009,74(6):WCC1-WCC26.

[10] ALAN R L. Fourth-order finite-difference P-SV seismograms [J]. Geophysics, 1988,53(11):1425-1436.

[11] 李慶洋,李振春,黃建平,等. 基于貼體全交錯網格的起伏地表正演模擬影響因素 [J]. 吉林大學學報(地球科學版),2016,46(3):920-929.

LI Qingyang, LI Zhenchun, HUANG Jianping, et al. Factor analysis of seismic modeling with topography based on a fully staggered body-fitted grids [J]. Journal of Jilin University (Earth Science Edition),2016,46(3):920-929.

[12] 李娜,李振春,黃建平,等.Lebedev網格與標準交錯網格耦合機制下的復雜各向異性正演模擬[J].石油地球物理勘探,2014,49(1):121-131.

LI Na, LI Zhenchun, HUANG Jianping, et al. Numerical simulation with coupling Lebedev and standard staggered grid schemes for complex anisotropic media [J]. Oil Geophysical Prospecting, 2014,49(1):121-131.

[13] 李振春,李慶洋,黃建平,等.一種穩(wěn)定的高精度雙變網格正演模擬與逆時偏移方法[J]. 石油物探, 2014,53(2):127-136.

LI Zhenchun, LI Qingyang, HUANG Jianping, et al. A stable and high-precision dual-variable grid forward modeling and reverse time migration method [J]. Geophysical Prospecting for Petroleum, 2014,53(2):127-136.

[14] 黃建平,楊宇,李振春,等. 基于 M-PML 邊界的 Lebedev 網格起伏地表正演模擬方法及穩(wěn)定性分析[J].中國石油大學學報(自然科學版), 2016,40(4):47-56.

HUANG Jianping, YANG Yu, LI Zhenchun, et al. Lebedev grid finite difference modeling for irregular free surface and stability analysis based on M-PML boundary condition [J]. Journal of China University of Petroleum (Edition of Natural Science), 2016,40(4):47-56.

[15] BIONDI B, ALMOMIN A. Tomographic full-waveform inversion (TFWI) by combining FWI and wave-equation migration velocity analysis [J]. The Leading Edge, 2013,32(9):1074-1080.

[16] BIONDI B, ALMOMIN A. Simultaneous inversion of full data bandwidth by tomographic full-waveform inversion [J]. Geophysics, 2014,79(3):WA129-WA140.

[17] WANG H, SINGH S C, AUDEBERT F, et al. Inversion of seismic refraction and reflection data for building long-wavelength velocity models [J]. Geophysics, 2015,80(2):R81-R93.

[18] FICHTNER A, KENNETT B L N, IGEL H, et al. Theoretical background for continental-and global-scale full-waveform inversion in the time-frequency domain[J]. Geophys J Int, 2008,175(2):665-685.

[19] MORA P. Inversion=migration+tomography[J]. Geophysics, 1989,54(12):1575-1586.

[20] 董良國,遲本鑫,陶紀霞,等.聲波全波形反演目標函數性態(tài)[J].地球物理學報,2013,56(10):3445-3460.

DONG Liangguo, CHI Benxin, TAO Jixia, et al. Objective function behavior in acoustic full-waveform inversion [J]. Chinese Journal Geophysics, 2013,56(10):3445-3460.

[21] LOMAX A, SNIEDER R. Finding sets of acceptable solutions with a genetic algorithm with application to surface wave group dispersion in Europe [J]. Geophys Res Lett, 1994,21(24):2617-2620.

[22] SAMBRIDGE M. Exploring multidimensional landscapes without a map [J]. Inverse Problems, 1998,14(3):427-440.

[23] SAMBRIDGE M. Geophysical inversion with a neighborhood algorithm-I. searching parameter space [J]. Geophys J Int, 1999a,138(2):479-494.

[24] SAMBRIDGE M. Geophysical inversion with a neighborhood algorithm-II. appraising the ensemble [J]. Geophys J Int, 1999b,138(3):727-746.

[25] KOPER K D, WYSESSION M E, WIENS D A. Multimodal function optimization with a niching genetic algorithm: a seismological example [J].Bull Seism Soc Am, 1999,89(4):978-988.

[26] MOSEGAARD K, SAMBRIDGE M. Monte Carlo analysis of inverse problems [J]. Inverse Problems, 2002,18(3):R29-R54.

[27] TARANTOLA A. Inverse problem theory [M]. Amsterdam:Elsevier, 1987:44-57.

[28] EFRON B, TIBSHIRANI, ROBERT J. An introduction to the bootstrap [M]. New York: Chapman and Hall/CRC, 1993.

[29] RAWLINSON N, M SAMBRIDGE, SAYGIN E. A dynamic objective function technique for generating multiple solution models in seismic tomography [J]. Geophys J Int, 2008,174(1):295-308.

[31] 李振春,雍鵬,黃建平,等. 基于矢量波場分離彈性波逆時偏移成像[J]. 中國石油大學學報(自然科學版),2016,40(1):42-48.

LI Zhenchun, YONG Peng, HUANG Jianping, et al. Elastic wave reverse time migration based on vector wavefield seperation [J]. Journal of China University of Petroleum (Edition of Natural Science), 2016,40(1):42-48.

[32] FICHTNER A, TRAMPERT J. Hessian kernels of seismic data functionals based upon adjoint techniques [J]. Geophys J Int, 2011,185(2):775-798.

Multi-solutionanalysisoffullwaveforminversionbasedondynamicmisfitfunction

HUANG Jianping1,2， CUI Chao1,2

(1.SchoolofGeosciencesinChinaUniversityofPetroleum,Qingdao266580,China;2.LaboratoryforMarineMineralResources,QingdaoNationalLaboratoryforMarineScienceandTechnology,Qingdao266071,China)

Full wave-form inversion (FWI) is one of the most promising velocity model building tools for its high precision and resolution which theoretically can utilize all the information carried by seismic record. However, because of the multi-solution property of the geophysical inversion problem and the limitation of the quality of seismic data, FWI can only generate one of the local optimize solutions which partly explain the seismic data. The analysis of the inversion result and the multi-solution property is one of the key problems in geophysical inversion. In order to mitigate the computational burden and strict application conditions in conventional result analysis methods, we introduce a dynamic function to analyze the multi-solution property of FWI. The technique designed for generating multiple solutions is based on the modification of objective function using the information of the previous results. The searching for a new model is achieved by adding a feedback term which creates a local maximum at each point in parameter space filled with previously computed models. For the model where the objective function defines a local maximum, the gradient is randomly perturbed to create a family of distinct solutions which helps evaluate the model where poorly illuminated by seismic data. The application on the sea model verifies that the method is efficient in terms of analyzing the inversion results of conventional FWI while keeping the accuracy. The merit of the proposed method is that it only needs to produce a relatively small ensemble of solutions, since each model will substantially differ from all others to the extent permitted by the data and the computational burden compared with the traditional method is significantly reduced.

full waveform inversion; multi-solution; objective function; reliability analysis

2017-02-18

山東省重大創(chuàng)新工程(2017CXGC1608，2017CXGC1602)；中國科學院戰(zhàn)略性先導科技專項(XDA14010303)；國際合作重點項目(41720104006)；國家油氣重大專項課題(2016ZX05006-004，2016ZX05014001，2016ZX05002)

黃建平(1982-)，男，教授，博士，研究方向為地震波傳播與成像技術。 E-mail:jphuang@upc.edu.cn。

1673-5005(2017)06-0064-07

10.3969/j.issn.1673-5005.2017.06.007

P 631.4

A

黃建平,崔超. 一種基于動態(tài)目標泛函的全波形反演多解性評估方法[J].中國石油大學學報(自然科學版),2017,41(6):64-70.

HUANG Jianping, CUI Chao. Multi-solution analysis of full waveform inversion based on dynamic misfit function [J].Journal of China University of Petroleum (Edition of Natural Science), 2017,41(6):64-70.

(編輯 修榮榮)