史才旺 何兵壽*
(①中國海洋大學(xué),山東青島 266100; ②青島海洋科學(xué)與技術(shù)國家實(shí)驗(yàn)室海洋礦產(chǎn)資源評價(jià)與探測技術(shù)功能實(shí)驗(yàn)室,山東青島 266071; ③海底科學(xué)與探測技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室,山東青島 266100)
全波形反演(FWI)利用了地震記錄中的振幅、相位等信息,是一種分辨率較高的速度建模方法。Tarantola[1]最早提出了FWI的基本方法,該方法基于廣義最小二乘反演理論,利用理論波場與反傳播的殘差波場的互相關(guān)構(gòu)造梯度方向。此后,F(xiàn)WI方法得到不斷完善,Bunks等[2]首先提出了時(shí)間域多尺度反演的理念,從低頻逐步反演到高頻,顯著降低了FWI的非線性性。Pratt[3]首先提出了頻率域的FWI思想,其優(yōu)勢在于僅利用幾個(gè)離散的頻率就能對模型實(shí)現(xiàn)高精度重建。此外,Shin等[4]的Laplace域反演和Sirgue等[5]的混合域反演方法進(jìn)一步豐富了FWI理論。
FWI技術(shù)在過去的30年間取得了長足的進(jìn)步,現(xiàn)在已經(jīng)開始應(yīng)用于實(shí)際資料。FWI本身是一套完美的數(shù)學(xué)體系,然而在實(shí)際應(yīng)用中,由于地震記錄總是會受到不同程度的破壞,如低頻能量的缺失、各種噪聲的影響等,F(xiàn)WI 很難取得理想的效果。實(shí)際地震記錄包含各種隨機(jī)噪聲和相干噪聲,在噪聲較強(qiáng)時(shí),能量較弱的反射同相軸甚至可能被完全掩蓋。在模擬記錄中不含噪聲,而由殘差記錄計(jì)算出的梯度會包含噪聲,因此適當(dāng)?shù)娜ピ胩幚韺WI十分必要。不同的相干噪聲有不同的針對性壓制方法。對于隨機(jī)噪聲,常規(guī)的策略是對信號進(jìn)行濾波處理,以犧牲部分有效信號予以去除。對于頻帶與有效信號重合的噪聲,則一般通過有效信號之間的相關(guān)性予以壓制,如SVD方法[6,7]、小波變換[8]和近年興起的Curvelet變換[9,10]等。若將這些方法直接作用于地震記錄,當(dāng)信噪比較低時(shí)對波形傷害很大,會導(dǎo)致FWI數(shù)據(jù)的不匹配。
對于含噪地震數(shù)據(jù),目前一般使用幾種特殊的目標(biāo)函數(shù)改善反演效果。常規(guī)FWI基于L2范數(shù),其實(shí)用性較強(qiáng),但抗噪能力較弱。Tarantola[11]和Hugues等[12]的實(shí)驗(yàn)和實(shí)際資料處理結(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]的彈性波實(shí)驗(yàn)結(jié)果也表明,L1范數(shù)比L2范數(shù)抗噪性強(qiáng),對存在明顯異常值的地震記錄也能得到更好的反演結(jié)果。但是L1范數(shù)本身也具有一定缺陷,如當(dāng)數(shù)據(jù)誤差接近于零時(shí),基于L1范數(shù)求取的目標(biāo)函數(shù)梯度會出現(xiàn)不穩(wěn)定,而L2范數(shù)則沒有這個(gè)問題。當(dāng)?shù)卣鹩涗浿性肼晣?yán)重時(shí),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í),通過增加頻率個(gè)數(shù)可以在一定程度上改善反演效果。Oh等[19]在單尺度、全頻域彈性波反演時(shí)構(gòu)建了一個(gè)降噪函數(shù),針對不同頻率估算信噪比,以增大高信噪比頻段在梯度中的占比。這種方法需要針對實(shí)際情況指定一個(gè)合適的降噪因子;此外,該方法只適合單尺度反演,當(dāng)高頻段信噪比較低時(shí),會降低反演分辨率。Kwon等[20]利用平面波編碼方式提高編碼記錄的信噪比,改善了含噪彈性波記錄的反演效果。
為進(jìn)一步降低隨機(jī)噪聲對反演結(jié)果的影響,本文結(jié)合主成分分析,通過對梯度進(jìn)行額外處理,實(shí)現(xiàn)了一種抗噪性更強(qiáng)的L2范數(shù)FWI方法。模型實(shí)驗(yàn)表明,該方法在信噪比較低時(shí)仍能得到正確的反演結(jié)果。
頻率域聲波方程可以視為一個(gè)關(guān)于波場值U的方程組
AU=S
(1)
式中:A為阻抗矩陣;U為由所有網(wǎng)格點(diǎn)波場值組成的列向量;S為震源列向量。
在最小二乘意義下,F(xiàn)WI的目標(biāo)函數(shù)為
(2)
式中:E為目標(biāo)函數(shù);m為待反演參數(shù),文中m代表速度;Ns為炮數(shù); δdi為第i炮各檢波點(diǎn)殘差,它是模擬地震記錄與實(shí)際觀測地震記錄之差,i為炮點(diǎn)序號; *代表取共軛。由此推導(dǎo)出目標(biāo)函數(shù)關(guān)于速度參數(shù)的梯度
(3)
為加快誤差函數(shù)的收斂速度,本文利用L-BFGS法對梯度方向進(jìn)行優(yōu)化[21],隨后通過拋物線擬合法求取一個(gè)合理的步長,從而實(shí)現(xiàn)模型更新。
為了詳細(xì)說明隨機(jī)噪聲對FWI的影響,使用Marmousi 2模型進(jìn)行含噪數(shù)據(jù)單頻反演實(shí)驗(yàn)。首先對Marmousi 2部分模型(圖1a)采用1000m×1000m的高斯濾波器進(jìn)行平滑(默認(rèn)海水層速度為1500m/s),得到初始模型(圖1b)。對Marmousi 2部分模型及初始模型(圖1)進(jìn)行正演得到頻率域2.0Hz地震記錄的實(shí)部(圖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相對誤差曲線。由圖可見,其變化趨勢為先減小后增大。這是因?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相對誤差曲線 模型誤差errm根據(jù)標(biāo)準(zhǔn)模型mstd、FWI結(jié)果mFWI 按errm=計(jì)算[23]
差記錄中有效信號能量強(qiáng),噪聲能量弱,使反演中的梯度較準(zhǔn)確。當(dāng)反演繼續(xù)進(jìn)行,反演模型與真實(shí)模型較接近時(shí),殘差中有效信號能量弱,噪聲能量占主導(dǎo),此時(shí)梯度受噪聲破壞嚴(yán)重,導(dǎo)致模型朝著錯(cuò)誤方向更新。若能在模型誤差發(fā)散前及時(shí)使反演過程收斂,則反演結(jié)果的可靠性將大為提高。
基于以上分析,本文引入一種新的梯度處理方式——主成分分析(PCA),在反演中能達(dá)到以下目的:①在信噪比較高時(shí)能獲得準(zhǔn)確的梯度,不影響其收斂速度;②在信噪比較低時(shí)對梯度進(jìn)行優(yōu)化,在提取有效信息的同時(shí)阻止模型朝著錯(cuò)誤方向更新。
PCA是一種多元統(tǒng)計(jì)分析方法,其目的是使用少數(shù)幾個(gè)特征變量描述一組具有相關(guān)性的多個(gè)變量,這些特征變量互不相關(guān),為原始的多個(gè)變量的線性組合[24]。PCA最初是一種降維方法,Liu等[25]利用這種方法對地震資料進(jìn)行壓縮從而大大降低了FWI的計(jì)算量; 段超然等[26]在此基礎(chǔ)上利用PCA實(shí)現(xiàn)了高效的波場重構(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)每一炮對應(yīng)的梯度值。對A進(jìn)行PCA,返回三個(gè)矩陣變量,分別為主成分系數(shù)矩陣P、主成分變量矩陣A′及其對應(yīng)的特征值矩陣λ,其中A′的每一列稱為一個(gè)主成分,它們滿足
i=1,2,…,M;j=1,2,…,N
(4)
特征值矩陣λ描述了每一個(gè)主成分對矩陣A的貢獻(xiàn)值,λ由大到小排列,值越大則此主成分越重要,其對應(yīng)的主成分排位也越靠前。上式表示的過程可以改寫為
A=A′×P-1
(5)
即原始梯度可以通過主成分和主成分系數(shù)進(jìn)行重構(gòu)。由于A的主要信息集中于特征值較大的前幾個(gè)主成分中,因而只要選取幾個(gè)重要的主成分進(jìn)行重構(gòu)就可以還原矩陣A的大部分信息,其他主成分則可以舍棄,該過程可記為
(6)
式中:B′也為M×N階矩陣,它是由A′中只保留前k個(gè)主成分得到的;B為重構(gòu)得到的新梯度,將B按行進(jìn)行疊加便可得到最終的梯度場。
地震記錄中的噪聲會通過波場的互相關(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ù)此梯度更新得到的新模型則會引入大量噪聲,導(dǎo)致模型誤差擴(kuò)散。在各炮噪聲不相關(guān)的情況下,各主成分貢獻(xiàn)值差距不大。由于只選取了少量主成分參與重構(gòu),此時(shí)重構(gòu)的梯度與原始梯度相差較大,可以使目標(biāo)函數(shù)不降低或者極緩慢地降低。這種特性可以使反演過程快速收斂,不會引起誤差擴(kuò)散。
主成分重構(gòu)需要選定使用的主成分個(gè)數(shù),Kaiser準(zhǔn)則是一種常用的選取主成分個(gè)數(shù)的經(jīng)驗(yàn)方法。Kaiser準(zhǔn)則建議只保留特征值大于1的主成分,其他主成分則舍棄。本文即采用這種方法,不過應(yīng)用Kaiser準(zhǔn)則需要先對設(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需要不斷迭代,因此其對主成分個(gè)數(shù)沒有較高要求,直接利用第一主成分進(jìn)行反演依然能得到正確結(jié)果,只是其收斂速度遜于Kaiser準(zhǔn)則方法(附錄A)。
在頻率域反演中,噪聲的影響與自身強(qiáng)度和頻譜分布都有關(guān),因而有必要選取不同強(qiáng)度、不同性質(zhì)的隨機(jī)噪聲分別進(jìn)行實(shí)驗(yàn)。在加入白噪聲的情形下,由于白噪聲的能量在頻譜上均勻分布,低頻段受噪聲影響嚴(yán)重,震源主頻處信噪比較高。因此,在反演含白噪聲的地震數(shù)據(jù)時(shí),往往低頻段的反演效果很差,直接導(dǎo)致高頻段反演失敗。為更好地突出文中方法的優(yōu)勢,文中又利用頻譜能量分布不均衡的高斯有色噪聲進(jìn)行了實(shí)驗(yàn)。
使用Marmousi 2部分模型及初始模型(圖1)進(jìn)行實(shí)驗(yàn),觀測系統(tǒng)保持不變。在頻率域采用單頻連續(xù)反演策略,頻率點(diǎn)由低到高依次為2.0、2.5、3.75、5.625、8.125、11.25Hz,頻率迭代的終止條件是前后兩次迭代的反演結(jié)果相對誤差小于0.001。頻率選擇是頻率域FWI反演精度的重要影響因素,頻率選取越密集、頻率越高,反演精度越高。就本文的實(shí)際情況而言,11.25Hz的反演結(jié)果已經(jīng)接近原模型。對于不含噪聲的情況,常規(guī)FWI給出了準(zhǔn)確的反演結(jié)果(圖5)。
5.1.1 常規(guī)FWI
為說明本文方法應(yīng)用于不同信噪比地震數(shù)據(jù)時(shí)的反演效果,分別在原始記錄中加入不同強(qiáng)度的高斯白噪聲。圖6、圖7分別為信噪比為1.0、0.5時(shí)不同頻率點(diǎn)原始記錄傅里葉變換的實(shí)部,用以說明各頻率的信噪比。利用常規(guī)FWI對上述地震記錄進(jìn)行反演,圖8為常規(guī)FWI在5.625Hz頻點(diǎn)的第一次迭代獲取的梯度,可見梯度形態(tài)或多或少都偏離了真實(shí)模型。由常規(guī)FWI對含高斯白噪聲地震數(shù)據(jù)的反演結(jié)果(圖9)可見:由于信噪比的影響,反演效果不理想,其中淺層反演較為準(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)行對比可知,前者的精度更高。圖11為基于單炮的主成分梯度重構(gòu)FWI對含高斯白噪聲地震數(shù)據(jù)的反演結(jié)果,與常規(guī)FWI對含高斯白噪聲地震數(shù)據(jù)的反演結(jié)果(圖9)相比,前者更加貼合真實(shí)模型。圖12為基于單炮的主成分梯度重構(gòu)FWI的2.0Hz單頻反演相對誤差曲線。由圖可見,基于單炮的主成分梯度重構(gòu)FWI的反演過程穩(wěn)健,即使在信噪比最低的第一個(gè)頻點(diǎn),反演也沒有出現(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è)中,每一炮的梯度被視為對總梯度的一個(gè)描述。在主成分分析和重構(gòu)中,這些梯度組成了一個(gè)訓(xùn)練集,在理想情況下該訓(xùn)練集中的梯度應(yīng)當(dāng)都是由某個(gè)公共部分加上隨機(jī)噪聲組成。實(shí)際上每一炮的梯度只反映了炮點(diǎn)附近的梯度信息,頻率較高時(shí)這一點(diǎn)尤為明顯,這必定會對梯度重構(gòu)產(chǎn)生一定影響。為此,本文提出了基于組合炮的主成分梯度重構(gòu)方法。顧名思義,這種方法首先將所有炮梯度分成L個(gè)組,每個(gè)組內(nèi)的梯度先累加得到一個(gè)組合梯度,將其變形為設(shè)計(jì)矩陣A的一列,之后在再對A進(jìn)行主成分分析和重構(gòu)。由于主要目的是使訓(xùn)練集的公共成分盡量相似,所以本文采用等間隔采樣的方法進(jìn)行炮分組(圖13),盡量使每個(gè)組合里的炮點(diǎn)能夠均勻地覆蓋整個(gè)模型。這里將所有炮點(diǎn)按照空間位置均勻分成10個(gè)組,并對10個(gè)組合梯度進(jìn)行主成分分析和重構(gòu)。
圖10 基于單炮的主成分梯度重構(gòu)FWI在5.625Hz 頻點(diǎn)的第一次迭代獲取的梯度 (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í)驗(yàn),其他參數(shù)不變,得到基于組合炮的主成分梯度重構(gòu)FWI對含高斯白噪聲地震數(shù)據(jù)的反演結(jié)果(圖14)及其2.0Hz時(shí)的模型相對誤差曲線(圖15)。表1為含高斯白噪聲時(shí)不同方法的反演精度對比,可見改進(jìn)后的FWI方法得到的結(jié)果更加接近真實(shí)模型。
圖15 基于組合炮的主成分梯度重構(gòu)FWI的 2.0Hz單頻反演相對誤差曲線 (a)信噪比為1.0; (b) 信噪比為0.5表1 含高斯白噪聲時(shí)不同方法的反演精度對比
方 法迭代次數(shù)相對誤差信噪比1.0信噪比0.5信噪比1.0信噪比0.5常規(guī)FWI1892040.1280.222單炮主成分梯度重構(gòu)FWI91900.0990.123組合炮主成分梯度重構(gòu)FWI91930.0970.115
頻率域FWI從低頻開始逐步提高頻率,低頻段的反演結(jié)果對高頻段影響嚴(yán)重。為進(jìn)一步說明主成分梯度重構(gòu)方法對多尺度反演的普適性,本文在地震記錄中引入有色噪聲。有色噪聲的能量隨頻率變化,實(shí)驗(yàn)中加入的有色噪聲能量的頻率分布與地震記錄一致,即任意頻率的信號能量和噪聲能量都滿足
(7)
式中:Sf為頻率f處的單頻信號能量;Nf為單頻噪聲能量。
利用常規(guī)FWI和基于組合炮的主成分梯度重構(gòu)FWI進(jìn)行反演,得到含高斯有色噪聲地震數(shù)據(jù)反演結(jié)果(圖16)。表2為含有色噪聲時(shí)不同方法的反演相對誤差對比??梢姡航M合炮主成分梯度重構(gòu)FWI的反演相對誤差明顯小于常規(guī)FWI;對于常規(guī)FWI,由于信噪比較高,低頻段的反演結(jié)果相對準(zhǔn)確,但是后續(xù)頻段的誤差累積依然使反演結(jié)果的相對誤差較大。
圖16 含高斯有色噪聲地震數(shù)據(jù)反演結(jié)果 (a)常規(guī)FWI; (b)基于組合炮的主成分梯度重構(gòu)FWI表2 含有色噪聲時(shí)不同方法的反演相對誤差對比
頻率/Hz2.02.52.755.6258.12511.25常規(guī)FWI0.1360.1270.1240.1270.1280.137組合炮主成分梯度重構(gòu)FWI0.1300.1250.1160.1090.1060.105
主成分分析和重構(gòu)過程的運(yùn)算量都很小,基本不會使反演過程變慢。同時(shí),由于常規(guī)FWI的誤差函數(shù)常常無法收斂(圖4),因此迭代次數(shù)會達(dá)到預(yù)設(shè)的最大值,從而造成常規(guī)FWI的耗時(shí)很長且效果差。表3列出了基于組合炮的主成分梯度重構(gòu)FWI與常規(guī)FWI的計(jì)算量,表中所有實(shí)驗(yàn)基于同一平臺,所用的CPU型號為Intel Core i5-4460,主頻為3.20GHz,內(nèi)存為8G,利用UMFPACK軟件包提供的LU分解法進(jìn)行頻率域正演。由表3可見:在含白噪聲的數(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的效果更好,反演的相對誤差更小。
隨機(jī)噪聲會對FWI的梯度產(chǎn)生破壞,甚至導(dǎo)致反演失敗。在研究隨機(jī)噪聲對FWI影響的基礎(chǔ)上,本文引入了基于PCA的FWI方法。利用PCA將各炮梯度重新組合為一組不相關(guān)的主成分,選取部分主成分重構(gòu)梯度以壓制噪聲的擴(kuò)散,此方法在低信噪比情況下實(shí)現(xiàn)了FWI。
模型實(shí)驗(yàn)表明,基于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)驗(yàn)做法。一般可以認(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)方法時(shí),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é)報(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] 段超然,韓立國.主成分分析波場重構(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.