史才旺 何兵壽*

(①中國(guó)海洋大學(xué)，山東青島 266100； ②青島海洋科學(xué)與技術(shù)國(guó)家實(shí)驗(yàn)室海洋礦產(chǎn)資源評(píng)價(jià)與探測(cè)技術(shù)功能實(shí)驗(yàn)室,山東青島 266071； ③海底科學(xué)與探測(cè)技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室，山東青島 266100)

1 引言

全波形反演(FWI)利用了地震記錄中的振幅、相位等信息，是一種分辨率較高的速度建模方法。Tarantola[1]最早提出了FWI的基本方法，該方法基于廣義最小二乘反演理論，利用理論波場(chǎng)與反傳播的殘差波場(chǎng)的互相關(guān)構(gòu)造梯度方向。此后，F(xiàn)WI方法得到不斷完善，Bunks等[2]首先提出了時(shí)間域多尺度反演的理念，從低頻逐步反演到高頻，顯著降低了FWI的非線性性。Pratt[3]首先提出了頻率域的FWI思想，其優(yōu)勢(shì)在于僅利用幾個(gè)離散的頻率就能對(duì)模型實(shí)現(xiàn)高精度重建。此外，Shin等[4]的Laplace域反演和Sirgue等[5]的混合域反演方法進(jìn)一步豐富了FWI理論。

FWI技術(shù)在過(guò)去的30年間取得了長(zhǎng)足的進(jìn)步，現(xiàn)在已經(jīng)開(kāi)始應(yīng)用于實(shí)際資料。FWI本身是一套完美的數(shù)學(xué)體系，然而在實(shí)際應(yīng)用中，由于地震記錄總是會(huì)受到不同程度的破壞，如低頻能量的缺失、各種噪聲的影響等，F(xiàn)WI 很難取得理想的效果。實(shí)際地震記錄包含各種隨機(jī)噪聲和相干噪聲，在噪聲較強(qiáng)時(shí)，能量較弱的反射同相軸甚至可能被完全掩蓋。在模擬記錄中不含噪聲，而由殘差記錄計(jì)算出的梯度會(huì)包含噪聲，因此適當(dāng)?shù)娜ピ胩幚韺?duì)FWI十分必要。不同的相干噪聲有不同的針對(duì)性壓制方法。對(duì)于隨機(jī)噪聲，常規(guī)的策略是對(duì)信號(hào)進(jìn)行濾波處理，以犧牲部分有效信號(hào)予以去除。對(duì)于頻帶與有效信號(hào)重合的噪聲，則一般通過(guò)有效信號(hào)之間的相關(guān)性予以壓制，如SVD方法[6,7]、小波變換[8]和近年興起的Curvelet變換[9,10]等。若將這些方法直接作用于地震記錄，當(dāng)信噪比較低時(shí)對(duì)波形傷害很大，會(huì)導(dǎo)致FWI數(shù)據(jù)的不匹配。

對(duì)于含噪地震數(shù)據(jù)，目前一般使用幾種特殊的目標(biāo)函數(shù)改善反演效果。常規(guī)FWI基于L2范數(shù)，其實(shí)用性較強(qiáng)，但抗噪能力較弱。Tarantola[11]和Hugues等[12]的實(shí)驗(yàn)和實(shí)際資料處理結(jié)果表明，L1范數(shù)對(duì)噪聲并不敏感。Huber范數(shù)和混合范數(shù)都是L1范數(shù)與L2范數(shù)的結(jié)合，由于其帶有L1范數(shù)的特性，因而抗噪性較強(qiáng)，但是這兩種范數(shù)都需要指定具體閾值用于切換L1/L2范數(shù)，若閾值不合適可能導(dǎo)致反演效果較差[13,14]。Brossier等[15]的彈性波實(shí)驗(yàn)結(jié)果也表明，L1范數(shù)比L2范數(shù)抗噪性強(qiáng)，對(duì)存在明顯異常值的地震記錄也能得到更好的反演結(jié)果。但是L1范數(shù)本身也具有一定缺陷，如當(dāng)數(shù)據(jù)誤差接近于零時(shí)，基于L1范數(shù)求取的目標(biāo)函數(shù)梯度會(huì)出現(xiàn)不穩(wěn)定，而L2范數(shù)則沒(méi)有這個(gè)問(wèn)題。當(dāng)?shù)卣鹩涗浿性肼晣?yán)重時(shí)，L1范數(shù)也難以得到正確的結(jié)果。此外，成景旺等[16]使用了基于柯西分布的目標(biāo)函數(shù)，其對(duì)數(shù)據(jù)異常、線性噪聲具有良好的壓制能力。

目前，大多數(shù)基于L1范數(shù)的目標(biāo)函數(shù)對(duì)數(shù)據(jù)異常及一些相干噪聲壓制效果良好，但對(duì)隨機(jī)噪聲壓制能力有限[17]，前人在L2范數(shù)框架下提出了多種改進(jìn)方法。Brossier等[18]指出，在進(jìn)行含隨機(jī)噪聲的頻率域反演時(shí)，通過(guò)增加頻率個(gè)數(shù)可以在一定程度上改善反演效果。Oh等[19]在單尺度、全頻域彈性波反演時(shí)構(gòu)建了一個(gè)降噪函數(shù)，針對(duì)不同頻率估算信噪比，以增大高信噪比頻段在梯度中的占比。這種方法需要針對(duì)實(shí)際情況指定一個(gè)合適的降噪因子；此外，該方法只適合單尺度反演，當(dāng)高頻段信噪比較低時(shí)，會(huì)降低反演分辨率。Kwon等[20]利用平面波編碼方式提高編碼記錄的信噪比，改善了含噪彈性波記錄的反演效果。

為進(jìn)一步降低隨機(jī)噪聲對(duì)反演結(jié)果的影響，本文結(jié)合主成分分析，通過(guò)對(duì)梯度進(jìn)行額外處理，實(shí)現(xiàn)了一種抗噪性更強(qiáng)的L2范數(shù)FWI方法。模型實(shí)驗(yàn)表明，該方法在信噪比較低時(shí)仍能得到正確的反演結(jié)果。

2 FWI基本理論

頻率域聲波方程可以視為一個(gè)關(guān)于波場(chǎng)值U的方程組

AU=S

(1)

式中：A為阻抗矩陣；U為由所有網(wǎng)格點(diǎn)波場(chǎng)值組成的列向量；S為震源列向量。

在最小二乘意義下，F(xiàn)WI的目標(biāo)函數(shù)為

(2)

式中：E為目標(biāo)函數(shù)；m為待反演參數(shù)，文中m代表速度；Ns為炮數(shù)； δdi為第i炮各檢波點(diǎn)殘差，它是模擬地震記錄與實(shí)際觀測(cè)地震記錄之差，i為炮點(diǎn)序號(hào)； *代表取共軛。由此推導(dǎo)出目標(biāo)函數(shù)關(guān)于速度參數(shù)的梯度

(3)

為加快誤差函數(shù)的收斂速度，本文利用L-BFGS法對(duì)梯度方向進(jìn)行優(yōu)化[21]，隨后通過(guò)拋物線擬合法求取一個(gè)合理的步長(zhǎng)，從而實(shí)現(xiàn)模型更新。

3 隨機(jī)噪聲對(duì)FWI的影響

為了詳細(xì)說(shuō)明隨機(jī)噪聲對(duì)FWI的影響，使用Marmousi 2模型進(jìn)行含噪數(shù)據(jù)單頻反演實(shí)驗(yàn)。首先對(duì)Marmousi 2部分模型(圖1a)采用1000m×1000m的高斯濾波器進(jìn)行平滑(默認(rèn)海水層速度為1500m/s)，得到初始模型(圖1b)。對(duì)Marmousi 2部分模型及初始模型(圖1)進(jìn)行正演得到頻率域2.0Hz地震記錄的實(shí)部(圖2)，對(duì)其進(jìn)行2.0Hz的單頻FWI，得到2.0Hz單頻FWI結(jié)果(圖3)，采用以擬Hessian矩陣[22]對(duì)角元素為初始矩陣的L-BFGS法對(duì)梯度進(jìn)行了優(yōu)化。由FWI反演結(jié)果可見(jiàn)，含噪數(shù)據(jù)反演結(jié)果(圖3b)與不含噪數(shù)據(jù)反演結(jié)果(圖3a)存在較大差距，且前者的迭代次數(shù)更多。

FWI追求的是地震記錄的誤差最小，但并不能識(shí)別數(shù)據(jù)中的噪聲，而噪聲與地質(zhì)模型無(wú)關(guān)，因此噪聲的存在可能導(dǎo)致模型誤差不斷擴(kuò)大。圖4為2.0Hz單頻FWI相對(duì)誤差曲線。由圖可見(jiàn)，其變化趨勢(shì)為先減小后增大。這是因?yàn)樵诜囱莩跗?，初始模型與真實(shí)模型相差較大，導(dǎo)致地震記錄差距明顯，因而殘

圖1 Marmousi 2部分模型(a)及初始模型(b)

模型網(wǎng)格數(shù)為175×450，網(wǎng)格間距為20m。震源為主頻8Hz的雷克子波，原始記錄共有112炮數(shù)據(jù)，炮點(diǎn)深度為20m，炮間距為80m，共有450道接收，檢波器均勻分布于模型表面

圖2 頻率域2.0Hz地震記錄的實(shí)部 (a)不含噪聲； (b)含高斯白噪聲(信噪比為1)

圖3 2.0Hz單頻FWI結(jié)果 (a)不含噪聲； (b)含高斯白噪聲

圖4 2.0Hz單頻FWI相對(duì)誤差曲線 模型誤差errm根據(jù)標(biāo)準(zhǔn)模型mstd、FWI結(jié)果mFWI 按errm=計(jì)算[23]

差記錄中有效信號(hào)能量強(qiáng)，噪聲能量弱，使反演中的梯度較準(zhǔn)確。當(dāng)反演繼續(xù)進(jìn)行，反演模型與真實(shí)模型較接近時(shí)，殘差中有效信號(hào)能量弱，噪聲能量占主導(dǎo)，此時(shí)梯度受噪聲破壞嚴(yán)重，導(dǎo)致模型朝著錯(cuò)誤方向更新。若能在模型誤差發(fā)散前及時(shí)使反演過(guò)程收斂，則反演結(jié)果的可靠性將大為提高。

基于以上分析，本文引入一種新的梯度處理方式——主成分分析(PCA)，在反演中能達(dá)到以下目的：①在信噪比較高時(shí)能獲得準(zhǔn)確的梯度，不影響其收斂速度；②在信噪比較低時(shí)對(duì)梯度進(jìn)行優(yōu)化，在提取有效信息的同時(shí)阻止模型朝著錯(cuò)誤方向更新。

4 基于PCA的梯度重構(gòu)

PCA是一種多元統(tǒng)計(jì)分析方法，其目的是使用少數(shù)幾個(gè)特征變量描述一組具有相關(guān)性的多個(gè)變量，這些特征變量互不相關(guān)，為原始的多個(gè)變量的線性組合[24]。PCA最初是一種降維方法，Liu等[25]利用這種方法對(duì)地震資料進(jìn)行壓縮從而大大降低了FWI的計(jì)算量； 段超然等[26]在此基礎(chǔ)上利用PCA實(shí)現(xiàn)了高效的波場(chǎng)重構(gòu)反演，為常規(guī)FWI提供了可靠的初始模型。在圖像處理方面，PCA也是一種良好的去噪方法，在電磁數(shù)據(jù)去噪中也有較多應(yīng)用[27,28]。

在常規(guī)FWI中，總梯度是由單炮梯度直接疊加得到的。在本文中，每一炮的梯度都被視為描述總梯度的一個(gè)變量，而總梯度是所有炮的梯度的一種線性組合。假設(shè)某次反演迭代獲取了N炮的梯度信息，每一炮梯度含有M個(gè)空間點(diǎn)的速度梯度。首先將每一炮的梯度拉伸為一個(gè)列優(yōu)先存放的列向量，并將這些列向量組合為一個(gè)M×N階的矩陣A，稱為設(shè)計(jì)矩陣。矩陣A的第i(i=1，2，…，N)列代表第i炮記錄得到的所有空間點(diǎn)的梯度；第j(j=1，2，…，M)行代表第j個(gè)空間點(diǎn)每一炮對(duì)應(yīng)的梯度值。對(duì)A進(jìn)行PCA，返回三個(gè)矩陣變量，分別為主成分系數(shù)矩陣P、主成分變量矩陣A′及其對(duì)應(yīng)的特征值矩陣λ，其中A′的每一列稱為一個(gè)主成分，它們滿足

i=1,2,…,M;j=1,2,…,N

(4)

特征值矩陣λ描述了每一個(gè)主成分對(duì)矩陣A的貢獻(xiàn)值，λ由大到小排列，值越大則此主成分越重要，其對(duì)應(yīng)的主成分排位也越靠前。上式表示的過(guò)程可以改寫為

A=A′×P-1

(5)

即原始梯度可以通過(guò)主成分和主成分系數(shù)進(jìn)行重構(gòu)。由于A的主要信息集中于特征值較大的前幾個(gè)主成分中，因而只要選取幾個(gè)重要的主成分進(jìn)行重構(gòu)就可以還原矩陣A的大部分信息，其他主成分則可以舍棄,該過(guò)程可記為

(6)

式中：B′也為M×N階矩陣，它是由A′中只保留前k個(gè)主成分得到的；B為重構(gòu)得到的新梯度，將B按行進(jìn)行疊加便可得到最終的梯度場(chǎng)。

地震記錄中的噪聲會(huì)通過(guò)波場(chǎng)的互相關(guān)運(yùn)算映射到梯度中，當(dāng)噪聲具有隨機(jī)性時(shí)，可認(rèn)為各炮梯度中包含的噪聲也互不相關(guān)，與此同時(shí)，各炮的有效梯度則具有很強(qiáng)的相關(guān)性。當(dāng)有效梯度能量總體占優(yōu)時(shí)，主成分重構(gòu)能得到高信噪比的梯度信息。隨著反演的進(jìn)行，有效梯度能量減弱，常規(guī)FWI得到的梯度中噪聲能量占主導(dǎo)，依據(jù)此梯度更新得到的新模型則會(huì)引入大量噪聲，導(dǎo)致模型誤差擴(kuò)散。在各炮噪聲不相關(guān)的情況下，各主成分貢獻(xiàn)值差距不大。由于只選取了少量主成分參與重構(gòu)，此時(shí)重構(gòu)的梯度與原始梯度相差較大，可以使目標(biāo)函數(shù)不降低或者極緩慢地降低。這種特性可以使反演過(guò)程快速收斂，不會(huì)引起誤差擴(kuò)散。

主成分重構(gòu)需要選定使用的主成分個(gè)數(shù)，Kaiser準(zhǔn)則是一種常用的選取主成分個(gè)數(shù)的經(jīng)驗(yàn)方法。Kaiser準(zhǔn)則建議只保留特征值大于1的主成分，其他主成分則舍棄。本文即采用這種方法，不過(guò)應(yīng)用Kaiser準(zhǔn)則需要先對(duì)設(shè)計(jì)矩陣A進(jìn)行標(biāo)準(zhǔn)化，將其每一列轉(zhuǎn)換為均值為0、方差為1的標(biāo)準(zhǔn)化向量，而重構(gòu)得到的梯度也要去除標(biāo)準(zhǔn)化影響。事實(shí)上，由于FWI需要不斷迭代，因此其對(duì)主成分個(gè)數(shù)沒(méi)有較高要求，直接利用第一主成分進(jìn)行反演依然能得到正確結(jié)果，只是其收斂速度遜于Kaiser準(zhǔn)則方法(附錄A)。

5 模型實(shí)驗(yàn)

在頻率域反演中，噪聲的影響與自身強(qiáng)度和頻譜分布都有關(guān)，因而有必要選取不同強(qiáng)度、不同性質(zhì)的隨機(jī)噪聲分別進(jìn)行實(shí)驗(yàn)。在加入白噪聲的情形下，由于白噪聲的能量在頻譜上均勻分布，低頻段受噪聲影響嚴(yán)重，震源主頻處信噪比較高。因此，在反演含白噪聲的地震數(shù)據(jù)時(shí)，往往低頻段的反演效果很差，直接導(dǎo)致高頻段反演失敗。為更好地突出文中方法的優(yōu)勢(shì)，文中又利用頻譜能量分布不均衡的高斯有色噪聲進(jìn)行了實(shí)驗(yàn)。

5.1 含高斯白噪聲的FWI

使用Marmousi 2部分模型及初始模型(圖1)進(jìn)行實(shí)驗(yàn)，觀測(cè)系統(tǒng)保持不變。在頻率域采用單頻連續(xù)反演策略，頻率點(diǎn)由低到高依次為2.0、2.5、3.75、5.625、8.125、11.25Hz，頻率迭代的終止條件是前后兩次迭代的反演結(jié)果相對(duì)誤差小于0.001。頻率選擇是頻率域FWI反演精度的重要影響因素，頻率選取越密集、頻率越高，反演精度越高。就本文的實(shí)際情況而言，11.25Hz的反演結(jié)果已經(jīng)接近原模型。對(duì)于不含噪聲的情況，常規(guī)FWI給出了準(zhǔn)確的反演結(jié)果(圖5)。

5.1.1 常規(guī)FWI

為說(shuō)明本文方法應(yīng)用于不同信噪比地震數(shù)據(jù)時(shí)的反演效果，分別在原始記錄中加入不同強(qiáng)度的高斯白噪聲。圖6、圖7分別為信噪比為1.0、0.5時(shí)不同頻率點(diǎn)原始記錄傅里葉變換的實(shí)部，用以說(shuō)明各頻率的信噪比。利用常規(guī)FWI對(duì)上述地震記錄進(jìn)行反演，圖8為常規(guī)FWI在5.625Hz頻點(diǎn)的第一次迭代獲取的梯度，可見(jiàn)梯度形態(tài)或多或少都偏離了真實(shí)模型。由常規(guī)FWI對(duì)含高斯白噪聲地震數(shù)據(jù)的反演結(jié)果(圖9)可見(jiàn)：由于信噪比的影響，反演效果不理想，其中淺層反演較為準(zhǔn)確，深層受噪聲破壞嚴(yán)重。如精確模型的最大速度為4700m/s，反演結(jié)果則分別達(dá)到了5580、6685m/s，嚴(yán)重偏離了真實(shí)值，導(dǎo)致反演結(jié)果誤差較大。

圖5 不含噪聲地震數(shù)據(jù)的FWI結(jié)果

圖7 信噪比為0.5時(shí)不同頻率點(diǎn)的原始記錄傅里葉變換的實(shí)部 (a)2.0Hz； (b)3.75Hz； (c)11.25Hz

5.1.2 基于單炮的主成分梯度重構(gòu)FWI

首先，將每一炮的梯度單獨(dú)轉(zhuǎn)換為設(shè)計(jì)矩陣A的一列，直接利用主成分梯度重構(gòu)方法進(jìn)行反演。此時(shí)A是一個(gè)78750×112階的矩陣，矩陣標(biāo)準(zhǔn)化后利用Kaiser準(zhǔn)則選取特征值大于1的主成分參與重構(gòu)。

圖10為基于單炮的主成分梯度重構(gòu)FWI在5.625Hz頻點(diǎn)的第一次迭代獲取的梯度，與常規(guī)FWI在5.625Hz頻點(diǎn)的第一次迭代獲取的梯度(圖8)進(jìn)行對(duì)比可知，前者的精度更高。圖11為基于單炮的主成分梯度重構(gòu)FWI對(duì)含高斯白噪聲地震數(shù)據(jù)的反演結(jié)果，與常規(guī)FWI對(duì)含高斯白噪聲地震數(shù)據(jù)的反演結(jié)果(圖9)相比，前者更加貼合真實(shí)模型。圖12為基于單炮的主成分梯度重構(gòu)FWI的2.0Hz單頻反演相對(duì)誤差曲線。由圖可見(jiàn)，基于單炮的主成分梯度重構(gòu)FWI的反演過(guò)程穩(wěn)健，即使在信噪比最低的第一個(gè)頻點(diǎn)，反演也沒(méi)有出現(xiàn)誤差擴(kuò)散現(xiàn)象，具體表現(xiàn)為：當(dāng)信噪比為1.0時(shí)，反演結(jié)果幾乎完全還原了精確模型的特征(圖12a)；當(dāng)信噪比為0.5時(shí)，模型淺層得到了較好的還原，深層雖然稍差，但是也基本體現(xiàn)了精確模型的輪廓，并且最大速度約為5000m/s，速度誤差在合理范圍內(nèi)(圖12b)。

圖8 常規(guī)FWI在5.625Hz頻點(diǎn)的第一次迭代獲取的梯度 (a)信噪比為1.0； (b)信噪比為0.5

5.1.3 基于組合炮的主成分梯度重構(gòu)FWI

在本文方法的理論假設(shè)中，每一炮的梯度被視為對(duì)總梯度的一個(gè)描述。在主成分分析和重構(gòu)中，這些梯度組成了一個(gè)訓(xùn)練集，在理想情況下該訓(xùn)練集中的梯度應(yīng)當(dāng)都是由某個(gè)公共部分加上隨機(jī)噪聲組成。實(shí)際上每一炮的梯度只反映了炮點(diǎn)附近的梯度信息，頻率較高時(shí)這一點(diǎn)尤為明顯，這必定會(huì)對(duì)梯度重構(gòu)產(chǎn)生一定影響。為此，本文提出了基于組合炮的主成分梯度重構(gòu)方法。顧名思義，這種方法首先將所有炮梯度分成L個(gè)組，每個(gè)組內(nèi)的梯度先累加得到一個(gè)組合梯度，將其變形為設(shè)計(jì)矩陣A的一列，之后在再對(duì)A進(jìn)行主成分分析和重構(gòu)。由于主要目的是使訓(xùn)練集的公共成分盡量相似，所以本文采用等間隔采樣的方法進(jìn)行炮分組(圖13)，盡量使每個(gè)組合里的炮點(diǎn)能夠均勻地覆蓋整個(gè)模型。這里將所有炮點(diǎn)按照空間位置均勻分成10個(gè)組，并對(duì)10個(gè)組合梯度進(jìn)行主成分分析和重構(gòu)。

圖10 基于單炮的主成分梯度重構(gòu)FWI在5.625Hz 頻點(diǎn)的第一次迭代獲取的梯度 (a)信噪比為1.0； (b)信噪比為0.5

圖11 基于單炮的主成分梯度重構(gòu)FWI對(duì)含高斯白噪聲 地震數(shù)據(jù)的反演結(jié)果 (a)信噪比為1.0； (b)信噪比為0.5

圖12 基于單炮的主成分梯度重構(gòu)FWI的2.0Hz 單頻反演相對(duì)誤差曲線 (a)信噪比為1.0； (b)信噪比為0.5

圖13 炮分組方法示意

圖14 基于組合炮的主成分梯度重構(gòu)FWI對(duì) 含高斯白噪聲地震數(shù)據(jù)的反演結(jié)果 (a)信噪比為1.0； (b)信噪比為0.5

利用此方法再進(jìn)行反演實(shí)驗(yàn)，其他參數(shù)不變，得到基于組合炮的主成分梯度重構(gòu)FWI對(duì)含高斯白噪聲地震數(shù)據(jù)的反演結(jié)果(圖14)及其2.0Hz時(shí)的模型相對(duì)誤差曲線(圖15)。表1為含高斯白噪聲時(shí)不同方法的反演精度對(duì)比，可見(jiàn)改進(jìn)后的FWI方法得到的結(jié)果更加接近真實(shí)模型。

圖15 基于組合炮的主成分梯度重構(gòu)FWI的 2.0Hz單頻反演相對(duì)誤差曲線 (a)信噪比為1.0； (b) 信噪比為0.5表1 含高斯白噪聲時(shí)不同方法的反演精度對(duì)比

方 法迭代次數(shù)相對(duì)誤差信噪比1.0信噪比0.5信噪比1.0信噪比0.5常規(guī)FWI1892040.1280.222單炮主成分梯度重構(gòu)FWI91900.0990.123組合炮主成分梯度重構(gòu)FWI91930.0970.115

5.2 含高斯有色噪聲的FWI

頻率域FWI從低頻開(kāi)始逐步提高頻率，低頻段的反演結(jié)果對(duì)高頻段影響嚴(yán)重。為進(jìn)一步說(shuō)明主成分梯度重構(gòu)方法對(duì)多尺度反演的普適性，本文在地震記錄中引入有色噪聲。有色噪聲的能量隨頻率變化，實(shí)驗(yàn)中加入的有色噪聲能量的頻率分布與地震記錄一致，即任意頻率的信號(hào)能量和噪聲能量都滿足

(7)

式中：Sf為頻率f處的單頻信號(hào)能量；Nf為單頻噪聲能量。

利用常規(guī)FWI和基于組合炮的主成分梯度重構(gòu)FWI進(jìn)行反演，得到含高斯有色噪聲地震數(shù)據(jù)反演結(jié)果(圖16)。表2為含有色噪聲時(shí)不同方法的反演相對(duì)誤差對(duì)比。可見(jiàn)：組合炮主成分梯度重構(gòu)FWI的反演相對(duì)誤差明顯小于常規(guī)FWI；對(duì)于常規(guī)FWI，由于信噪比較高，低頻段的反演結(jié)果相對(duì)準(zhǔn)確，但是后續(xù)頻段的誤差累積依然使反演結(jié)果的相對(duì)誤差較大。

圖16 含高斯有色噪聲地震數(shù)據(jù)反演結(jié)果 (a)常規(guī)FWI； (b)基于組合炮的主成分梯度重構(gòu)FWI表2 含有色噪聲時(shí)不同方法的反演相對(duì)誤差對(duì)比

頻率/Hz2.02.52.755.6258.12511.25常規(guī)FWI0.1360.1270.1240.1270.1280.137組合炮主成分梯度重構(gòu)FWI0.1300.1250.1160.1090.1060.105

5.3 計(jì)算量對(duì)比

主成分分析和重構(gòu)過(guò)程的運(yùn)算量都很小，基本不會(huì)使反演過(guò)程變慢。同時(shí)，由于常規(guī)FWI的誤差函數(shù)常常無(wú)法收斂(圖4)，因此迭代次數(shù)會(huì)達(dá)到預(yù)設(shè)的最大值，從而造成常規(guī)FWI的耗時(shí)很長(zhǎng)且效果差。表3列出了基于組合炮的主成分梯度重構(gòu)FWI與常規(guī)FWI的計(jì)算量，表中所有實(shí)驗(yàn)基于同一平臺(tái)，所用的CPU型號(hào)為Intel Core i5-4460，主頻為3.20GHz，內(nèi)存為8G，利用UMFPACK軟件包提供的LU分解法進(jìn)行頻率域正演。由表3可見(jiàn)：在含白噪聲的數(shù)據(jù)反演中，改進(jìn)FWI的運(yùn)算時(shí)間不到常規(guī)FWI的一半；在含有色噪聲的數(shù)據(jù)反演中，運(yùn)算時(shí)間也大幅減少。

表3 基于組合炮的主成分梯度重構(gòu)FWI 與常規(guī)FWI的計(jì)算量

綜上所述，當(dāng)?shù)卣鹩涗浿泻性肼晻r(shí)，主成分梯度重構(gòu)FWI較常規(guī)FWI的抗噪性更強(qiáng)，模型誤差收斂速度也較快，計(jì)算量更小。此外，基于組合炮的梯度重構(gòu)FWI比基于單炮的梯度重構(gòu)FWI的效果更好，反演的相對(duì)誤差更小。

6 結(jié)束語(yǔ)

隨機(jī)噪聲會(huì)對(duì)FWI的梯度產(chǎn)生破壞，甚至導(dǎo)致反演失敗。在研究隨機(jī)噪聲對(duì)FWI影響的基礎(chǔ)上，本文引入了基于PCA的FWI方法。利用PCA將各炮梯度重新組合為一組不相關(guān)的主成分，選取部分主成分重構(gòu)梯度以壓制噪聲的擴(kuò)散，此方法在低信噪比情況下實(shí)現(xiàn)了FWI。

模型實(shí)驗(yàn)表明，基于PCA和梯度重構(gòu)的FWI方法具有較高的抗噪性。此外，本文對(duì)比了基于單炮和組合炮的主成分梯度重構(gòu)FWI方法，結(jié)果表明后者的效果更佳。

尚需指出，本文所添加的噪聲均服從高斯分布，且所用的經(jīng)典的PCA的魯棒性仍然不足，有關(guān)其他類型的噪聲對(duì)波形反演的影響及其壓制策略、引入改進(jìn)的高魯棒性PCA進(jìn)一步提高對(duì)噪聲和數(shù)據(jù)異常的適應(yīng)性等問(wèn)題是今后的研究方向。

附錄A 只選取第一主成分進(jìn)行梯度重構(gòu)的反演效果

本文中使用Kaiser準(zhǔn)則的目的是選定特定數(shù)量的主成分參與重構(gòu)，這是一種經(jīng)驗(yàn)做法。一般可以認(rèn)為特征值大于1的主成分中有效信號(hào)居多，反之則噪聲能量占優(yōu)?？紤]到FWI的過(guò)程需要不斷迭代，因此，即使某次重構(gòu)的梯度缺失了部分有效信息，在之后的迭代中也可以得到彌補(bǔ)，但相應(yīng)的反演過(guò)程的收斂速度可能會(huì)受到影響。針對(duì)文中所述的基于單炮的主成分梯度重構(gòu)FWI，為與之對(duì)比，在其他各參數(shù)不變的情況下，僅利用第一主成分進(jìn)行梯度重構(gòu)。

圖A1為含高斯白噪聲地震數(shù)據(jù)的反演結(jié)果，圖A2為含高斯白噪聲地震數(shù)據(jù)的反演相對(duì)誤差曲線。由圖可見(jiàn)，應(yīng)用Kaiser準(zhǔn)則的梯度重構(gòu)方法的反演結(jié)果更為精確。表A1為含高斯白噪聲地震數(shù)據(jù)的反演相對(duì)誤差對(duì)比。由表A1可見(jiàn)，Kaiser準(zhǔn)則方法最終模型的相對(duì)誤差達(dá)到0.099，只利用第一主成分的相對(duì)誤差為0.101。在收斂速度方面，應(yīng)用Kaiser準(zhǔn)則的梯度重構(gòu)方法顯然更快。值得注意的是，Kaiser準(zhǔn)則的作用與參與PCA 的梯度數(shù)量有關(guān)，即梯度數(shù)量越多，兩種方法差距越明顯。而當(dāng)使用基于組合炮的主成分梯度重構(gòu)方法時(shí)，Kaiser準(zhǔn)則本身選取的主成分?jǐn)?shù)量就較少，因而兩種方法差距并不大。

圖A1 含高斯白噪聲地震數(shù)據(jù)的反演結(jié)果 (a)應(yīng)用Kaiser準(zhǔn)則； (b)只應(yīng)用第一主成分

圖A2 含高斯白噪聲地震數(shù)據(jù)的反演相對(duì)誤差曲線 (a)2.0Hz； (b)8.125Hz表A1 含高斯白噪聲地震數(shù)據(jù)的反演相對(duì)誤差對(duì)比

頻率/Hz2.02.52.755.6258.12511.25應(yīng)用Kaiser準(zhǔn)則0.1310.1270.1180.1080.1010.099只應(yīng)用第一主成分0.1310.1270.1180.1100.1030.101

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

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

[3] Pratt R G.Seismic waveform inversion in the frequency domain,Part 1:Theory and verification in a physical scale model.Geophysics,1999,64(3):888-901.

[4] Shin C,Cha Y H.Waveform inversion in the Laplace domain.Geophysical Journal International,2008,173(3):922-931.

[5] Sirgue L,Etgen J T,Albertin U.3D frequency domain waveform inversion using time domain finite difference methods.EAGE 70th Conference and Technical Exhibition,2008.

[6] Bekara M,Van der Baan M.Local SVD/ICA for signal enhancement of pre-stack seismic data.EAGE 68th Conference and Technical Exhibition,2006.

[7] 沈鴻雁,李慶春.頻域奇異值分解 (SVD) 地震波場(chǎng)去噪.石油地球物理勘探,2010,45(2):185-189. Shen Hongyan,Li Qingchun.SVD(singular value decomposition) seismic wave field noise elimination in frequency domain.OGP,2010,45(2):185-189.

[8] 高靜懷,毛劍,滿蔚仕等.疊前地震資料噪聲衰減的小波域方法研究.地球物理學(xué)報(bào),2006,49(4):1155-1163. Gao Jinghuai,Mao Jian,Man Weishi et al.On the denoising method of prestack seismic data in wavelet domain.Chinese Journal of Geophysics,2006,49(4):1155-1163.

[9] 張恒磊,劉天佑,張?jiān)拼?基于高階相關(guān)的 Curvelet 域和空間域的傾角掃描噪聲壓制方法.石油地球物理勘探,2010,45(2):208-214. Zhang Henglei,Liu Tianyou,Zhang Yuncui.High order correlation based dip angle scanning noise elimination method in Curvelet domain and space domain.OGP,2010,45(2):208-214.

[10] 張華,陳小宏,李紅星等.曲波變換三維地震數(shù)據(jù)去噪技術(shù).石油地球物理勘探,2017,52(2):226-232. Zhang Hua,Chen Xiaohong,Li Hongxing et al.3D seismic data de-noising approach based on Curvelet transform.OGP,2017,52(2):226-232.

[11] Tarantola A.Inverse Problem Theory:Methods for Data Fitting and Model Parameter Estimation.Elsevier,Amsterdam and New York,1987.

[12] Hugues A D,Tarantola A.Multiparameter L1norm waveform fitting: interpretation of Gulf of Mexico reflection seismograms.Geophysics,1999,64(4):1023-1035.

[13] Ha T,Chung W,Shin C.Waveform inversion using a back-propagation algorithm and a Huber function norm.Geophysics,2009,74(3):R15-R24.

[14] Bube K P.Hybrid L1/L2minimization with applica-tions to tomography.Geophysics,1997,62(4):1183-1195.

[15] Brossier R,Operto S,Virieux J.Robust elastic frequency-domain full-waveform inversion using the L1norm.Geophysical Research Letters,2009,36(20):L20310.

[16] 成景旺,呂曉春,顧漢明等.基于柯西分布的頻率域全波形反演.石油地球物理勘探,2014,49(5):940-945. Cheng Jingwang,Lü Xiaochun,Gu Hanming et al.Full waveform inversion with Cauchy distribution in the frequency domain.OGP,2014,49(5):940-945.

[17] Oh J W,Min D J.Robust scaling strategy for frequency-domain acoustic full waveform inversion.ASEG 23rd Geophysical Conference Extended Abstracts,2013,1-4.

[18] Brossier R,Operto S,Virieux J.Which data residual norm for robust elastic frequency-domain full waveform inversion?.Geophysics,2010,75(3):R37-R46.

[19] Oh J W,Min D J.Spectral filtering of gradient for L2-norm frequency-domain elastic waveform inversion.Geophysical Journal International,2013,193(2):820-840.

[20] Kwon T,Seol S J,Byun J.Reduction of noise effect on elastic full-waveform inversion by using plane-wave data.SEG Technical Program Expanded Abstracts,2016,35:1131-1135.

[21] 苗永康.基于L-BFGS算法的時(shí)間域全波形反演.石油地球物理勘探,2015,50(3):469-474. Miao Yongkang.Full waveform inversion in time domain based on limited-memory BFGS algorithm.OGP,2015,50(3):469-474.

[22] Shin C,Jang S,Min D J.Improved amplitude preservation for prestack depth migration by inverse.Geophysical Prospecting,2008,49(5):592-606.

[23] Ben-Hadj-Ali H,Operto S,Virieux J.An efficient frequency-domain full waveform inversion method using simultaneous encoded sources.Geophysics,2011,76(4):R109-R124.

[24] 王學(xué)民.應(yīng)用多元分析.上海：上海財(cái)經(jīng)大學(xué)出版社,2014,251-256.

[25] Liu C,Han M,Han L et al.Application of principal component analysis for frequency-domain full waveform inversion.SEG Technical Program Expanded Abstracts,2012,31:1-5.

[26] 段超然,韓立國(guó).主成分分析波場(chǎng)重構(gòu)反演與全波形反演聯(lián)合速度重構(gòu).石油地球物理勘探,2016,51(6):1134-1140. Duan Chaoran,Han Liguo.A joint velocity reconstruction method: principal component analysis based wavefield-reconstruction inversion combined with full waveform inversion.OGP,2016,51(6):1134-1140.

[27] 王凌群,李冰冰,林君等.航空電磁數(shù)據(jù)主成分濾波重構(gòu)的噪聲去除方法.地球物理學(xué)報(bào),2015,58(8):2803-2811. Wang Lingqun,Li Bingbing,Lin Jun et al.Noise removal based on reconstruction of filtered principal components.Chinese Journal of Geophysics,2015,58(8):2803-2811.

[28] Zhu K G,Ma M Y,Che H W et al.PC-based artificial neural network inversion for airborne time-domain electromagnetic data.Applied Geophysics,2012,9(1):1-8.