基于主成分分析和梯度重構(gòu)的全波形反演

史才旺 何兵壽*

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

1 引言

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

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

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

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

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

2 FWI基本理論

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

AU=S

(1)

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

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

(2)

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

(3)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

(4)

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

A=A′×P-1

(5)

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

(6)

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

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

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

5 模型實驗

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

5.1 含高斯白噪聲的FWI

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

5.1.1 常規(guī)FWI

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

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

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

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

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

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

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

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

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

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

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

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

圖13 炮分組方法示意

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

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

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

方 法迭代次數(shù)相對誤差信噪比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從低頻開始逐步提高頻率，低頻段的反演結(jié)果對高頻段影響嚴(yán)重。為進(jìn)一步說明主成分梯度重構(gòu)方法對多尺度反演的普適性，本文在地震記錄中引入有色噪聲。有色噪聲的能量隨頻率變化，實驗中加入的有色噪聲能量的頻率分布與地震記錄一致，即任意頻率的信號能量和噪聲能量都滿足

(7)

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

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

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

頻率/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 計算量對比

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

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

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

6 結(jié)束語

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

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

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

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

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

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

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

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

頻率/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) 地震波場去噪.石油地球物理勘探,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é)報,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] 張恒磊,劉天佑,張云翠.基于高階相關(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算法的時間域全波形反演.石油地球物理勘探,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)用多元分析.上海：上海財經(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] 段超然,韓立國.主成分分析波場重構(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é)報,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.