孫宏宇， 韓立國(guó),2*， 韓淼， 王志強(qiáng)

1 吉林大學(xué) 地球探測(cè)科學(xué)與技術(shù)學(xué)院， 長(zhǎng)春 1300262 國(guó)土資源部應(yīng)用地球物理重點(diǎn)實(shí)驗(yàn)室， 長(zhǎng)春 1300263 中國(guó)地質(zhì)調(diào)查局 油氣資源調(diào)查中心， 北京 100029

?

基于可視性分析與能量補(bǔ)償?shù)慕饘俚V彈性波全波形反演

孫宏宇1， 韓立國(guó)1,2*， 韓淼3， 王志強(qiáng)1

1 吉林大學(xué) 地球探測(cè)科學(xué)與技術(shù)學(xué)院， 長(zhǎng)春 1300262 國(guó)土資源部應(yīng)用地球物理重點(diǎn)實(shí)驗(yàn)室， 長(zhǎng)春 1300263 中國(guó)地質(zhì)調(diào)查局 油氣資源調(diào)查中心， 北京 100029

采用彈性波全波形反演方法精確重建深部金屬礦多參數(shù)模型，建模過程采用基于地震照明的反演策略.首先給出基于照明理論的觀測(cè)系統(tǒng)可視性定義，利用可視性分析構(gòu)建新的目標(biāo)函數(shù)，對(duì)反演目標(biāo)可視性較高的炮檢對(duì)接收到的地震記錄在波場(chǎng)匹配時(shí)占有更高的權(quán)重，確保了參與反演計(jì)算中的地震數(shù)據(jù)的有效性；其次將給定觀測(cè)系統(tǒng)對(duì)地下介質(zhì)的彈性波場(chǎng)照明強(qiáng)度作為優(yōu)化因子，根據(jù)地震波在波阻抗界面處的能量分配特點(diǎn)，自適應(yīng)補(bǔ)償波場(chǎng)能量分布和優(yōu)化速度梯度，以提高彈性波全波形反演過程的穩(wěn)定性和反演結(jié)果的精度.理論模型和金屬礦模型反演試驗(yàn)結(jié)果表明，基于可視性分析和能量補(bǔ)償?shù)姆囱莶呗钥梢允箯椥圆ㄈㄐ畏囱莞斓厥諗康侥繕?biāo)函數(shù)的全局極小值，獲得適用于金屬礦高分辨率地震偏移成像的多參數(shù)模型.

金屬礦地震； 全波形反演； 彈性波； 可視性分析； 能量補(bǔ)償

1 引言

近年來(lái)，國(guó)內(nèi)外金屬礦地震勘探技術(shù)發(fā)展較為迅速，已經(jīng)成為尋找深部礦產(chǎn)的重要方法之一.我國(guó)在金屬礦地震理論和應(yīng)用方面也開展了一些探索，取得了一定進(jìn)展（呂慶田等，2010，2014；廉玉廣等，2011）.金屬礦區(qū)人文環(huán)境與地質(zhì)構(gòu)造復(fù)雜、目標(biāo)地質(zhì)體與圍巖波阻抗差較小，造成了金屬礦反射地震信號(hào)能量弱和呈不連續(xù)的非層狀反射，或者散射特征，往往還會(huì)受到強(qiáng)背景噪聲的干擾，導(dǎo)致利用常規(guī)的速度建模技術(shù)構(gòu)建的金屬礦地震速度模型的精度較低，偏移成像效果差.此外，金屬礦目標(biāo)地質(zhì)體與圍巖之間波阻抗差異較小、或者存在漸變等特征，常規(guī)的縱波勘探難以取得預(yù)期效果.而基于彈性波多參數(shù)的多波多分量勘探方法有望克服傳統(tǒng)的縱波地震勘探反射信息弱的問題（李慶春等，2007），但是常規(guī)的速度建模和成像方法的精度和分辨率對(duì)于大傾角等復(fù)雜介質(zhì)都會(huì)有很大的誤差（韓建光等，2013），為此需要發(fā)展高精度地震波速度和密度等多參數(shù)模型重建和偏移成像的方法.

地震全波形反演（FWI）理論于20世紀(jì)80年代建立（Tarantola,1984; Mora,1987）并逐步得到完善和改進(jìn)，近年來(lái)更是取得了迅速發(fā)展（Virieux and Operto，2009；韓淼，2014）.同時(shí)，彈性波全波形反演（EFWI）的計(jì)算速度和精度也有了很大的提高（Castellanos et al.,2011; Wang et al., 2012; Xiong et al., 2013; Son et al., 2014）.繼Pratt等（1998）將聲波近似的全波形反演方法引入到頻率域后，Brossier等（2009）在頻率域?qū)崿F(xiàn)了彈性波FWI，并采用多尺度策略提高反演精度.而時(shí)間域的多尺度反演方法能靈活選取時(shí)窗分離波場(chǎng)，便于降低反演問題的非線性性（Boonyasiriwat et al.,2009; Fichtner et al.,2013），更適用于處理波場(chǎng)復(fù)雜的實(shí)際資料.由于波場(chǎng)中各屬性參數(shù)相互耦合，因此利用彈性波全波形反演方法重建地下介質(zhì)波場(chǎng)參數(shù)分布更加困難.K?hn等（2012）詳細(xì)討論了時(shí)域EFWI選用不同的參數(shù)組合，如拉梅系數(shù)，地震波速度或波阻抗作為反演參數(shù)時(shí)，參數(shù)間的相互影響.Xu和McMechan（2014）綜合考慮各個(gè)參數(shù)的梯度變化，在得到多參數(shù)相互獨(dú)立的搜索方向和步長(zhǎng)后，沿著多參數(shù)合成的搜索方向?qū)ふ乙粋€(gè)最優(yōu)步長(zhǎng)，以較快的收斂速度實(shí)現(xiàn)任意參數(shù)組合的反演；Groos等（2014）引入黏彈性介質(zhì)的品質(zhì)因子同時(shí)反演淺部體波和瑞雷波；Ren和Liu （2015）提出聯(lián)合變差分階數(shù)正演和基于子波的多尺度反演策略.這些方法的提出使得彈性波全波形反演更加成熟，已經(jīng)成功應(yīng)用到了海洋（Prieux et al., 2013；Raknes et al., 2015）和陸地地震勘探實(shí)際數(shù)據(jù)的多參數(shù)建模中（Plessix and Pérez, 2015；Borisov and Singh, 2015）.

地震照明技術(shù)通過建立地質(zhì)模型人工模擬野外采集方式，分析地震波在介質(zhì)中傳播過程中能量的分布，評(píng)價(jià)觀測(cè)系統(tǒng)參數(shù)（謝小碧等，2013），提高AVO/AVA分析精度（Askim et al., 2010）和優(yōu)化偏移成像剖面（Zhou et al.,2014）.Yang等（2013）將照明能量作為加權(quán)因子引入到成像域聲波全波形反演中，消除建立速度模型過程中照明不均的影響.Bian等（2015）利用照明技術(shù)估計(jì)全波形反演中觀測(cè)系統(tǒng)參數(shù)對(duì)反演結(jié)果的作用.陳永芮等（2013）利用雙向照明優(yōu)化觀測(cè)系統(tǒng)并對(duì)聲波全波形反演梯度進(jìn)行能量補(bǔ)償.應(yīng)用全波形反演處理實(shí)際數(shù)據(jù)時(shí)，可以根據(jù)先驗(yàn)?zāi)Ｐ瓦M(jìn)行可視性分析找到最優(yōu)炮點(diǎn)分布和最佳接收范圍，并在對(duì)應(yīng)位置加密震源和檢波器，使觀測(cè)數(shù)據(jù)包含更多的有效信息，提高反演精度和穩(wěn)定性.由于時(shí)間域全波形反演的計(jì)算復(fù)雜度與模型大小及觀測(cè)系統(tǒng)有直接關(guān)系，隨著震源點(diǎn)數(shù)的增加，計(jì)算量成正比增長(zhǎng).Wang和Rao（2009）在聲波全波形反演中采用層剝離的反演方法，解決了有限偏移距的波場(chǎng)效應(yīng)和反演誤差傳遞等問題，逐層反演以重構(gòu)深部縱波速度模型，但這種方法的計(jì)算量較大，正比于對(duì)底層的人為分層數(shù)目.本文提出一種基于可視性分析與自適應(yīng)能量補(bǔ)償?shù)膹椥圆ㄈㄐ畏囱莘椒?，并?yīng)用于復(fù)雜的金屬礦地震勘探多參數(shù)建模處理.首先將觀測(cè)系統(tǒng)中單個(gè)震源和檢波器對(duì)目標(biāo)體的總照明強(qiáng)度定義為該炮檢對(duì)對(duì)目標(biāo)體的可視性，對(duì)彈性波全波形反演中的觀測(cè)系統(tǒng)進(jìn)行可視性分析并計(jì)算出可視性分析平面圖，將可視性作為加權(quán)因子引入到全波形反演中目標(biāo)函數(shù)的計(jì)算，給出基于可視性的目標(biāo)函數(shù)，以增加觀測(cè)波場(chǎng)和模擬波場(chǎng)相互匹配過程中與目的層相關(guān)性較高的炮檢對(duì)的殘差權(quán)重.其次，對(duì)于金屬礦復(fù)雜地區(qū)存在的照明陰影區(qū)而導(dǎo)致彈性波場(chǎng)能量不均勻分布、成像效果差的問題，根據(jù)給定的觀測(cè)系統(tǒng)計(jì)算波場(chǎng)對(duì)地下介質(zhì)的雙向照明強(qiáng)度，構(gòu)建自適應(yīng)的加權(quán)梯度函數(shù)，使得地震波場(chǎng)能量近似均衡分布，從而改善全波形反演效果.

2 方法原理

2.1 彈性波波形反演

全波形反演中的目標(biāo)函數(shù)E可以表示為

-uobs（xs,xr,t）]2dt,

（1）目標(biāo)函數(shù)E對(duì)模型參數(shù)m求導(dǎo)可得到梯度方向?yàn)?/p>

（2）

其中，ucal和uobs分別是模擬波場(chǎng)和觀測(cè)波場(chǎng)，δd是兩者間的殘差.若雅克比矩陣?u/?m已知，則對(duì)模型空間中所有擾動(dòng)進(jìn)行積分可計(jì)算數(shù)據(jù)空間波場(chǎng)的總變化值. 同理，對(duì)數(shù)據(jù)空間的波場(chǎng)擾動(dòng)求積分可以得到模型空間的變化量為（Tarantola,2005）

（3）

（4）

其中，*代表相應(yīng)的共軛矩陣.由于這種反問題的解并非唯一，即δu≠δu′，δm≠δm′，但由于[?u（m）/?m]*=?u（m）/?m，因此若將數(shù)據(jù)空間的擾動(dòng)看作殘差，則從數(shù)據(jù)空間到模型空間的映射就等同于目標(biāo)函數(shù)的梯度：

（5）

因此，根據(jù)一階擾動(dòng)形式下的彈性波方程和格林函數(shù)，可以得到縱波速度vp、橫波速度vs和密度ρ的梯度表達(dá)式為（Mora,1987）

（6）

（8）

其中，ux和uz分別是在震源點(diǎn)激發(fā)的彈性波場(chǎng)的水平分量和垂直分量，φx和φz分別是在檢波點(diǎn)位置激發(fā)的反傳殘差波場(chǎng)的水平分量和垂直分量.選用合適的優(yōu)化方法，沿著目標(biāo)函數(shù)的負(fù)梯度方向，利用迭代公式（9）對(duì)初始模型m0迭代更新，即可得到目標(biāo)函數(shù)的全局最優(yōu)解，公式（9）為

mn+1=mn-αnPδmn，

（9）

其中，αn為第n次迭代的步長(zhǎng)，P為預(yù)處理算子.若選用牛頓法優(yōu)化參數(shù)模型，則P為海森矩陣的逆，若選用高斯牛頓法優(yōu)化參數(shù)模型，則P為對(duì)角海森矩陣的逆.2.2 基于可視性的目標(biāo)函數(shù)

金屬礦地震勘探由于地表?xiàng)l件的復(fù)雜和地下介質(zhì)的不均勻分布等因素的干擾，只有部分位置的炮點(diǎn)和檢波點(diǎn)能激發(fā)并接收到與探測(cè)地質(zhì)體有關(guān)的地震波信息，通過波動(dòng)方程照明分析，可以得到觀測(cè)系統(tǒng)中單個(gè)震源和檢波點(diǎn)對(duì)目的層的歸一化總照明強(qiáng)度即觀測(cè)系統(tǒng)的可視性，它能直觀地描述觀測(cè)系統(tǒng)中任意炮撿對(duì)的觀測(cè)數(shù)據(jù)對(duì)于成像或反演的有效作用程度.計(jì)算基于可視性的目標(biāo)函數(shù)來(lái)優(yōu)化全波形反演中對(duì)應(yīng)炮檢對(duì)的模擬地震記錄與觀測(cè)記錄之間的殘差，提高了與目的層有關(guān)的地震記錄在目標(biāo)函數(shù)中的權(quán)重，重點(diǎn)反演深部目標(biāo)體的構(gòu)造參數(shù)信息.對(duì)于二維觀測(cè)系統(tǒng)中的任意炮檢對(duì)（xs,xr），其對(duì)地下介質(zhì)的照明強(qiáng)度為

（10）

（11）

（12）

（13）其中，vis表示地震波由xs處的震源激發(fā)xr處接收時(shí)成像目標(biāo)范圍內(nèi)整體的照明強(qiáng)度，即炮檢對(duì)（xs,xr）的可視性.

為便于估計(jì)任意一炮檢對(duì)所得到的地震記錄在波形反演計(jì)算時(shí)波場(chǎng)匹配過程中的作用，我們用可視性的平均值作為判別標(biāo)準(zhǔn)，大于平均值說明該炮檢對(duì)接收到的地震記錄與目的層的相關(guān)性較高，因此經(jīng)預(yù)處理后可以得到該觀測(cè)系統(tǒng)對(duì)目的層的相對(duì)可視性vis′為

vis′（xs,xr）=φ×vis（xs,xr），

（14）

（15）

其中，ns和nr分別為炮點(diǎn)和檢波點(diǎn)的個(gè)數(shù).φ為vis的平均值的倒數(shù).當(dāng)相對(duì)可視性的值大于1時(shí)，即該炮檢對(duì)的可視性大于觀測(cè)系統(tǒng)對(duì)目的層的可視性平均值，將其引入到波形反演中目標(biāo)函數(shù)的計(jì)算，可以得到基于可視性的彈性波全波形反演的目標(biāo)函數(shù)為

[ucal（xs,xr,m,t）-uobs（xs,xr,t）]2dt,

（16）

此時(shí)，與目標(biāo)體相關(guān)性較高的地震數(shù)據(jù)在波場(chǎng)匹配時(shí)將會(huì)占有更多的權(quán)重.

圖1a和c分別為觀測(cè)系統(tǒng)中某炮在MarmousiII模型及反演初始模型上正演的水平分量地震記錄，圖1b和d為相應(yīng)的垂直分量地震記錄.該炮對(duì)應(yīng)的所有檢波點(diǎn)對(duì)目標(biāo)體的相對(duì)可視性曲線如圖2所示，該觀測(cè)系統(tǒng)中共有400個(gè)檢波器.可以看出，位于第120到190號(hào)和第280到400號(hào)之間的檢波器與目標(biāo)體的相關(guān)性較高，可以接收到更多來(lái)自于反演目標(biāo)體的波場(chǎng)信息.以此可視性作用于彈性波全波形反演中波場(chǎng)水平分量和垂直分量的殘差（圖3a和c，對(duì)應(yīng)的位置殘差能量增強(qiáng)（圖3b和d）可以使彈性波全波形反演計(jì)算時(shí)充分利用更多的有效數(shù)據(jù)，提高反演精度和穩(wěn)定性.

圖1 MarmousiII模型上正演的彈性波場(chǎng)（a）x分量、（b）z分量和初始模型上正演的彈性波場(chǎng)（c）x分量、（d）z分量

圖2 觀測(cè)系統(tǒng)中某炮及其對(duì)應(yīng)所有檢波點(diǎn)的相對(duì)可視性曲線

圖3 x分量觀測(cè)與模擬記錄的（a）原始?xì)埐詈停╞）經(jīng)能量加權(quán)后的殘差以及z分量觀測(cè)與模擬記錄的（c）原始?xì)埐詈停╠）經(jīng)能量加權(quán)后的殘差

2.3 自適應(yīng)能量補(bǔ)償

正是因?yàn)榄h(huán)保的自覺遵守，以及理念超前，所以無(wú)錫工廠自2012年起連年被無(wú)錫新區(qū)評(píng)為綠標(biāo)企業(yè)。愛克發(fā)用行動(dòng)證明了印刷業(yè)完全可以成為綠色環(huán)保產(chǎn)業(yè)。

反射能量較強(qiáng)的地震記錄在波形反演數(shù)據(jù)擬合時(shí)占有更大的比重，而能量較弱的波場(chǎng)對(duì)應(yīng)的速度結(jié)構(gòu)在反演過程中得不到較好的恢復(fù)，即使其速度結(jié)構(gòu)與真實(shí)模型相差很大，因此，這種波場(chǎng)能量的不均勻分布會(huì)影響反演結(jié)果的穩(wěn)定性.兩種原因會(huì)導(dǎo)致地震波能量在地下介質(zhì)中不均勻分布，一是由于地震波的幾何擴(kuò)散現(xiàn)象，波場(chǎng)能量隨著傳播深度的增加而減小；二是由于地下介質(zhì)的復(fù)雜性，使得波在其傳播方向上存在照明陰影區(qū).利用波動(dòng)方程理論得到的雙向照明強(qiáng)度自適應(yīng)加權(quán)梯度，可以同時(shí)補(bǔ)償因觀測(cè)系統(tǒng)中偏移距不足導(dǎo)致深部波場(chǎng)能量弱的問題和復(fù)雜介質(zhì)中存在照明盲區(qū)而引起的能量損失，加權(quán)因子取決于在地震波的傳播方向?qū)?yīng)位置處的照明強(qiáng)度，這在數(shù)學(xué)上類似于利用海森矩陣優(yōu)化梯度，得到牛頓類優(yōu)化方法中的模型更新量.地震波在彈性介質(zhì)分界面處的能量分配情況可以用Zoeppritz方程描述.波場(chǎng)由高速介質(zhì)向低速介質(zhì)傳播時(shí)，透射系數(shù)T>1，地震波的透射能量強(qiáng)于反射能量，界面下部的照明強(qiáng)度高于界面上部，此時(shí)可以用雙向照明強(qiáng)度反比例加權(quán)于相應(yīng)參數(shù)的梯度；相反，地震波場(chǎng)由低速向高速介質(zhì)中傳播時(shí)，透射系數(shù)T<1，則地震波的反射能量強(qiáng)于透射能量，界面上部的照明強(qiáng)度高于界面下部，此時(shí)可以用歸一化后的照明強(qiáng)度正比例加權(quán)于相應(yīng)參數(shù)的梯度.由此可得到彈性波反演中經(jīng)自適應(yīng)能量補(bǔ)償后的模型更新量δm′為

（17）

其中，I是彈性波場(chǎng)照明強(qiáng)度.

3 模型試驗(yàn)與效果分析

3.1 MarmousiII模型

采用MarmousiII模型（圖4）來(lái)測(cè)試本文的反演方法.模型網(wǎng)格大小為174×500，網(wǎng)格間距為20 m，橫向長(zhǎng)10 km，縱向深3.48 km.觀測(cè)系統(tǒng)模擬水中激發(fā)縱波，以海底多波觀測(cè)方式采集接收縱波和轉(zhuǎn)換波，100個(gè)震源水平分布于海水面以下40 m處，400個(gè)三分量的檢波器位于海底.為了分析觀測(cè)系統(tǒng)的可視性，設(shè)目的層為位于深部斷層下的含油藏的圈閉構(gòu)造，如圖4中黑框所示.采樣率為2 ms，地震記錄長(zhǎng)度為6s.震源為主頻10 Hz的雷克子波.反演選用LBFGS優(yōu)化方法，最大迭代次數(shù)為100次.

對(duì)于給定的目標(biāo)范圍，經(jīng)照明統(tǒng)計(jì)可以得到該觀測(cè)系統(tǒng)對(duì)地下介質(zhì)的相對(duì)可視性如圖5a所示，圖中橫縱坐標(biāo)分別代表該觀測(cè)系統(tǒng)中檢波器和震源的序號(hào)，各點(diǎn)的數(shù)值代表當(dāng)?shù)卣鸩ㄓ煽v坐標(biāo)對(duì)應(yīng)位置

處的震源激發(fā)，橫坐標(biāo)對(duì)應(yīng)位置處的檢波器接收時(shí)，該炮檢對(duì)記錄對(duì)目的層的相對(duì)可視性.可以看出該觀測(cè)系統(tǒng)中，第40到第80個(gè)震源，第100到第230個(gè)檢波器，以及第280到第400個(gè)檢波器所對(duì)應(yīng)的相對(duì)可視性均大于1，說明這個(gè)范圍內(nèi)的炮檢對(duì)接收到的波場(chǎng)對(duì)目的層的成像起到主要作用.圖5b為該觀測(cè)系統(tǒng)對(duì)地下介質(zhì)的彈性波場(chǎng)照明強(qiáng)度，由于觀測(cè)系統(tǒng)中震源和檢波器主要分布在模型中部，因此模型兩側(cè)的照明較弱，中上部照明較強(qiáng).此外由于地震波的幾何擴(kuò)散現(xiàn)象以及復(fù)雜介質(zhì)引起的照明不均，模型深部以及復(fù)雜斷裂帶下方的照明能量較弱.我們應(yīng)用可視性分析方法提高對(duì)目的層成像作用較大的炮檢點(diǎn)對(duì)應(yīng)的水平分量和垂直分量的殘差在波場(chǎng)匹配中的權(quán)重，同時(shí)采用自適應(yīng)算法進(jìn)行能量補(bǔ)償.

圖6為兩種反演方法均迭代100次后的反演結(jié)果.縱波速度、橫波速度和密度三種參數(shù)均得到了較好的恢復(fù)，可以看出，基于可視性分析和能量補(bǔ)償?shù)腅FWI的反演結(jié)果要明顯優(yōu)于常規(guī)EFWI.從反演結(jié)果、真實(shí)模型和初始模型中提取模型中水平距離為5.2 km處的參數(shù)變化曲線，如圖7所示.相比于常規(guī)EFWI，本文提出的EFWI方法反演的縱波速度、橫波速度和密度曲線更加趨近于真實(shí)模型，對(duì)復(fù)雜構(gòu)造的有效重構(gòu)同時(shí)加快了反演的收斂速度，使本文方法在相同的迭代次數(shù)條件下獲得更加精確的反演結(jié)果.3.2 金屬礦模型彈性波全波形反演對(duì)安徽省廬樅盆地某礦床的復(fù)雜地球物理模型（圖8）采用彈性波全波形反演方法重建速度和密度多參數(shù)模型.該模型總體呈復(fù)雜的似層狀、平緩?fù)哥R狀，空間上表現(xiàn)為穹隆狀，中心以侵染狀磁鐵礦為主，富、厚礦多環(huán)于四周.礦床由多個(gè)礦體組成，其中規(guī)模較大的只有3個(gè)礦體，其余均為小礦體（呂慶田等，2010；廉玉廣等，2011）.對(duì)該模型進(jìn)行重采樣后得到的縱波速度模型如圖9a所示.模型網(wǎng)格點(diǎn)數(shù)為120×400，網(wǎng)格間距為20 m，模型橫向長(zhǎng)8 km，縱向深2.4 km.主礦體埋深距地表1.3～1.8 km間，且厚度變化較大.橫波速度和密度可以在縱波速度基礎(chǔ)上由經(jīng)驗(yàn)公式得出（Brocher, 2005; Xu and McMechan,2014），如圖9b和c所示，公式為

圖4 MarmousiII模型

圖5 （a）觀測(cè)系統(tǒng)對(duì)指定目的層的可視性和（b）彈性波場(chǎng)照明強(qiáng)度

圖6 反演結(jié)果

圖7 不同方法反演的參數(shù)剖面對(duì)比（x=5.2 km處）：真實(shí)模型（黑實(shí)線）、初始模型（黑虛線）、常規(guī)EFWI反演結(jié)果（藍(lán)線）、基于可視性分析和能量補(bǔ)償?shù)腅FWI反演結(jié)果（紅線）

圖8 S1線過羅河礦體地球物理模型（據(jù)廉玉廣等，2011）

圖9 金屬礦多參數(shù)模型

圖10 （a）觀測(cè)系統(tǒng)對(duì)指定礦體的可視性和（b）彈性波場(chǎng)照明強(qiáng)度

圖11 金屬礦多參數(shù)模型反演結(jié)果

（18）

（19）

應(yīng)用于彈性波全波形反演的初始模型如圖9d、e和f所示.采樣率為2ms，地震記錄長(zhǎng)度為2s.震源子波為主頻10Hz的雷克子波.反演選用LBFGS優(yōu)化方法，最大迭代次數(shù)為100次.觀測(cè)系統(tǒng)中100個(gè)震源水平分布于地表下40m處，起始震源位于水平距離40m處，震源間距為80m.400個(gè)檢波器深度均為20m，起始坐標(biāo)和間距分別為0m和20m.

選定的目標(biāo)范圍即為礦體所在的位置，對(duì)該觀測(cè)系統(tǒng)進(jìn)行可視性分析的結(jié)果如圖10a所示，彈性波場(chǎng)照明強(qiáng)度如圖10b所示，彈性波的能量除了在地表較強(qiáng)，在礦體部位也有所聚集.將照明分析的結(jié)果引入到全波形反演中，得到基于地震照明的彈性波全波形反演結(jié)果如圖11d、e和f所示.與常規(guī)彈性波全波形反演的結(jié)果（圖11a、b和c）對(duì)比可以看出，礦體輪廓的反演更加清晰，說明在反演過程中引入可視性分析和照明補(bǔ)償，可以有效地改善反演效果.

4 結(jié)論

（1）基于可視性分析的彈性波全波形反演策略可以適用于金屬礦等復(fù)雜介質(zhì)構(gòu)造的多參數(shù)建模.通過照明統(tǒng)計(jì)分析可以得到給定觀測(cè)系統(tǒng)對(duì)反演目的層的可視性，將其應(yīng)用于波場(chǎng)匹配過程中對(duì)應(yīng)炮檢點(diǎn)的殘差能量計(jì)算.構(gòu)建基于可視性的反演目標(biāo)函數(shù)，可以充分利用來(lái)自于目的層的有效地震信息，提高反演精度.

（2）由于觀測(cè)系統(tǒng)對(duì)地下介質(zhì)的照明強(qiáng)度可以直觀地描述地震波能量的分布情況，利用Zoeppritz方程將彈性波場(chǎng)照明強(qiáng)度引入到全波形反演中，使得地震波場(chǎng)的能量自適應(yīng)地均衡分布和優(yōu)化梯度，改善反演質(zhì)量.

（3）彈性波全波形反演方法可以高精度地重構(gòu)深部金屬礦地質(zhì)體的多參數(shù)模型，包括縱橫波速度和密度，可以實(shí)現(xiàn)金屬礦區(qū)地下介質(zhì)構(gòu)造的高精度偏移成像.由于彈性波全波形反演可以利用多波多分量全波場(chǎng)信息，克服了單一的縱波地震勘探時(shí)由于弱阻抗差可能帶來(lái)的不利影響，為有效地應(yīng)用地震勘探技術(shù)開展深部礦產(chǎn)資源探測(cè)提供了新途徑.

致謝 感謝國(guó)家自然科學(xué)基金（41374115）和國(guó)家科技重大專項(xiàng)（2011ZX0502500104）的支持，感謝評(píng)審專家的寶貴意見，在論文研究過程得到了王者江教授、鞏向博副教授、張鳳蛟博士的有益幫助，在此一并致謝.

Askim O J, Albertin U, Matson K, et al. 2010. Wave-equation-based illumination determination for Amplitude versus Angle risk assessment. ∥80th SEG Annual Meeting Expanded Abstracts. Expanded Abstracts, 3263-3267.

Bian A F, Zou Z H, Zhou H W, et al. 2015. Evaluation of multi-scale full waveform inversion with marine vertical cable data.JournalofEarthScience, 26（4）: 481-486.

Boonyasiriwat C, Valasek P, Routh P, et al. 2009. An efficient multiscale method for time-domain waveform tomography.Geophysics, 74（6）: WCC59-WCC68. Borisov D, Singh S C. 2015. Three-dimensional elastic full waveform inversion in a marine environment using multicomponent ocean-bottom cables: a synthetic study.GeophysicalJournalInternational, 201（3）: 1215-1234.

Brocher T M. 2005. Empirical relations between elastic wavespeeds and density in the earth′s crust.BulletinoftheSeismologicalSocietyofAmerica, 95（6）: 2081-2092.

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.

Castellanos C, Etienne V, Hu G H, et al. 2011. Algorithmic and methodological developments towards full waveform inversion in 3D elastic media. ∥81st SEG Annual Meeting Expanded Abstracts. Expanded Abstracts, 2793-2798.

Chen Y R, Li Z C, Qin N, et al. 2013. Full waveform inversion with wave equation bi-directional illumination optimization.ProgressinGeophysics（in Chinese）, 28（6）: 3015-3021, doi: 10.6038/pg20130624.

Fichtner A, Trampert J, Cupillard P, et al. 2013. Multiscale full waveform inversion.GeophysicalJournalInternational, 194（1）: 534-556.

Groos L, Sch?fer M, Forbriger T, et al. 2014. The role of attenuation in 2D full-waveform inversion of shallow-seismic body and Rayleigh waves.Geophysics, 79（6）: R247-R261.

Han J G, Wang Y, Lu J, et al. 2013. Elastic wave Kirchhoff pre-stack depth migration velocity analysis.OilGeophysicalProspecting（in Chinese）, 48（5）: 694-699.

Han M. 2014. Methods and application of full waveform inversion for abyssal seismic data [Phd′s thesis]. （in Chinese）. Changchun: Jilin University.

K?hn D, De Nil D, Kurzmann A, et al. 2012. On the influence of model parametrization in elastic full waveform tomography.GeophysicalJournalInternational, 191（1）: 325-345.

Li Q C, Liu G H, Li Q. 2007. Multi-component seismic explorations-directions of seismic exploration for metal deposits. （in Chinese）. Xian: Shanxi Scientific and Technique Press.

Lian Y G, Lu Q T, Han L G, et al. 2011. Research of seismic modeling of complex metal ore body.ActaGeologicaSinica（in Chinese）, 85（5）:887-898.

Lü Q T, DONG S W, SHI D N, et al. 2014. Lithosphere architecture and geodynamic model of middle and lower reaches of Yangtze metallogenic belt: A review from SinoProbe.ActaGeologicaSinica（in Chinese）, 30（4）: 889-906.

Lü Q T, HAN L G, Yan J Y et al. 2010. Seismic imaging of volcanic hydrothermal iron-sulfur deposits and its hosting structure in Luzong ore district.ActaGeologicaSinica（in Chinese）, 26（9）: 2598-2612.

Mora P. 1987. Nonlinear two-dimensional elastic inversion of multioffset seismic data.Geophysics, 52（9）: 1211-1228.

Plessix R é, Pérez Solano C A. 2015. Modified surface boundary conditions for elastic waveform inversion of low-frequency wide-angle active land seismic data.GeophysicalJournalInternational, 201（3）: 1324-1334.

Pratt R G, Shin C, Hicks G J. 1998. Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion.Geophysics, 133（2）: 341-362.

Prieux V, Brossier R, Operto S, et al. 2013. Multiparameter full waveform inversion of multicomponent ocean-bottom-cable data from the Valhall field. Part 2: imaging compressive-wave and shear-wave velocities.GeophysicalJournalInternational, 194（3）: 1665-1681.

Raknes E B, Arntsen B, Weibull W. 2015. Three-dimensional elastic full waveform inversion using seismic data from the Sleipner area.GeophysicalJournalInternational, 202（3）:

1877-1894.

Ren Z M, Liu Y. 2015. Elastic full-waveform inversion using the second-generation wavelet and an adaptive-operator-length scheme.Geophysics, 80（4）: R155-R173.

Son W, Pyun S, Shin C, et al. 2014. Laplace-domain wave equation modeling and full waveform inversion in 3D isotropic elastic media.JournalofAppliedGeophysics, 105: 120-132.

Tarantola A. 1984. Inversion of seismic reflection data in the acoustic approximation.Geophysics, 49（8）: 1259-1266.

Tarantola A. 2005. Inverse problem theory and methods for model parameter estimation. Paris, France: SIAM, 128-133.

Virieux J, Operto S. 2009. An overview of full-waveform inversion in exploration geophysics.Geophysics, 74（6）: WCC127-WCC152.

Wang J, Zhou H, Tian Y K, et al. 2012. A new scheme for elastic full waveform inversion based on velocity-stress wave equations in time domain. ∥82nd SEG Annual Meeting Expanded Abstracts. Expanded Abstracts, 1-5.

Wang Y H, Rao Y. 2009. Reflection seismic waveform tomography.JournalofGeophysicalResearch:SolidEarth, 114（B3）: B03304.

Xiong J L, Lin Y, Abubakar A, et al. 2013. 2.5-D forward and inverse modelling of full-waveform elastic seismic survey.GeophysicalJournalInternational, 193（2）: 938-948.

Xu K, McMechan G A. 2014. 2D frequency-domain elastic full-waveform inversion using time-domain modeling and a multistep-length gradient approach.Geophysics, 79（2）: R41-R53.

Yang T N, Shragge J, Sava P. 2013. Illumination compensation for image-domain wavefield tomography.Geophysics, 78（5）: U65-U76.

Zhou H M, Chen S C, Ren H R, et al. 2014. One-way wave equation least-squares migration based on illumination compensation.ChineseJ.Geophys., 57（5）: 726-738.

Xie X B, He Y Q, Li P M. 2013. Seismic illumination analysis and its applications in seismic survey design.ChineseJ.Geophys. （in Chinese）, 56（5）: 1568-1581, doi: 10.6038/cjg20130515.

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

陳永芮, 李振春, 秦寧等. 2013. 波動(dòng)方程雙向照明優(yōu)化的全波形反演. 地球物理學(xué)進(jìn)展, 28（6）: 3015-3021, doi: 10.6038/pg20130624.韓建光, 王赟, 蘆俊等. 2013. 彈性波Kirchhoff疊前深度偏移速度分析. 石油地球物理勘探, 48（5）: 694-699.

韓淼. 2014. 深水區(qū)地震全波形反演策略與應(yīng)用[博士論文]. 長(zhǎng)春: 吉林大學(xué).

李慶春, 劉國(guó)華, 李勤. 2007. 多波多分量地震-金屬礦地震勘探的方向. 西安: 陜西科學(xué)技術(shù)出版社.

廉玉廣, 呂慶田, 韓立國(guó)等. 2011. 復(fù)雜金屬礦體地震波正演模擬研究——以廬樅盆地羅河、泥河和大包莊礦床為例. 地質(zhì)學(xué)報(bào), 85（5）: 887-899.

呂慶田, 董樹文, 史大年, 等. 2014. 長(zhǎng)江中下游成礦帶巖石圈結(jié)構(gòu)與成礦動(dòng)力學(xué)模型——深部探測(cè)（SinoProbe）綜述. 巖石學(xué)報(bào), 30（4）: 889-9064.

呂慶田, 韓立國(guó), 嚴(yán)加永等. 2010. 廬樅礦集區(qū)火山氣液型鐵、硫礦床及控礦構(gòu)造的反射地震成像. 巖石學(xué)報(bào), 26（9）: 2598-2612.

謝小碧, 何永清, 李培明. 2013. 地震照明分析及其在地震采集設(shè)計(jì)中的應(yīng)用. 地球物理學(xué)報(bào), 56（5）: 1568-1581, doi: 10.6038/cjg20130515.

（本文編輯 張正鋒）

Elastic full waveform inversion based on visibility analysis and energy compensation for metallic deposit exploration

SUN Hong-Yu1，HAN Li-Guo1,2*，HAN Miao3，WANG Zhi-Qiang1

1CollegeofGeo-explorationScienceandTechnology,JilinUniversity,Changchun130026,China2KeyLaboratoryofAppliedGeophysics,TheMinistryofLandandResources，Changchun130026,China3OilandGasSurveyofChinaGeologicalSurvey,Beijing100029,China

The weak and discontinuous reflections and scattering features in complex environments make it difficult to construct velocity models suitable for imaging in conventional P- wave seismic exploration for metallic deposits. Multi-wave and multi-component seismic exploration technology may have the advantages in effort to solve this problem. Elastic full waveform inversion （EFWI） exploiting the multi-component and pre-stack seismic data can reconstruct highly accurate multi-parameter models for imaging complex deposit structure in the deep subsurface. The EFWI method based on visibility analysis and adaptive energy compensation is introduced in this paper which can utilize prior information more adequately relevant to vector seismic fields from the target.For a specific geometry in seismic acquisition, only signals from sources and receivers within a certain extent have effect on the inversion and imaging about the objective body. According to the contribution for a single shot-geophone pair to the imaging of the objective body, the statistical analysis method is used to obtain the whole illumination intensity of the target zone which is defined as the visibility of the single source-receiver pair. To improve the inversion result by utilizing the seismic information related to the target zone as much as possible without increasing computation and acquisition costs, we can construct the objective function of EFWI based on geometry visibility to increase the proportion of residuals related to the target zone during wavefield matching. FWI is a data-fitting procedure between observed and calculated data. Wavefields with weaker energy make less contribution to the misfit of objective function even though the velocity contrast is bigger. The uneven distribution of energy due to several cases can be expressed by seismic illumination. So bi-directional illumination intensity of the elastic wavefield is used as a weighting factor to adaptively optimize and balance the gradients of EFWI. Zoeppritz equations control the process of automatic compensation because they describe the energy distribution of seismic waves across the interface of impedance.The EFWI algorithm proposed in this paper is implemented by the finite-difference method to calculate the elastic wavefields in the time domain and LBFGS optimization method to update the initial models. The multi-parameter gradients （density, P-wave and S-wave velocity） are calculated by the cross-correlation between forward and backward wavefields in the time domain. First, the validity of our inversion method is verified on the MarmousiII model. The accuracy of inversion results with visibility analysis and energy compensation is improved compared with conventional EFWI. The model test results based on one metallic model of Luzong Basin also illustrate that our new method makes the inversion converge to global minimum of objective function faster and provides accurate and multi-parameter models which can be applied to high-resolution seismic migration imaging of metallic deposits.The numerical results show that the elastic full waveform inversion algorithm based on visibility analysis utilizes the effective information from the target zone adequately. Introducing the illumination intensity of the elastic wavefield to EFWI by transmission coefficients between two layers can also balance the energy distribution of wavefields adaptively. Both of the two inversion strategies with illumination can improve the accuracy of imaging. The EFWI method can accurately reconstruct deep and multi-parameter models for metallic deposits, including density, P and S wave velocity, which makes it possible to obtain high-resolution migration imaging about deposit structures.

Metallic deposits； Full waveform inversion； Elastic wave； Visibility analysis； Energy compensation

國(guó)家高技術(shù)研究發(fā)展計(jì)劃（863計(jì)劃）重大項(xiàng)目課題（2014AA06A605）和國(guó)家深部探測(cè)專項(xiàng)第3項(xiàng)目（SinoProbe-03）聯(lián)合資助.

孫宏宇，女，1992年生，碩士研究生，主要從事地震波場(chǎng)正反演研究. E-mail: sunhongyu2014@foxmail.com

*通訊作者 韓立國(guó)，男，1961年生，教授，博士生導(dǎo)師，主要從事地震數(shù)據(jù)處理解釋工作. E-mail: hanliguo@jlu.edu.cn

10.6038/cjg20151222.

10.6038/cjg20151222

P631

2015-11-03，2015-11-26收修定稿

孫宏宇， 韓立國(guó)， 韓淼等. 2015. 基于可視性分析與能量補(bǔ)償?shù)慕饘俚V彈性波全波形反演.地球物理學(xué)報(bào)，58（12）:4605-4616,

Sun H Y, Han L G, Han M, et al. 2015. Elastic full waveform inversion based on visibility analysis and energy compensation for metallic deposit exploration.ChineseJ.Geophys. （in Chinese）,58（12）:4605-4616,doi:10.6038/cjg20151222.