霍元媛,楊睿,潘紀(jì)順,胡金虎
1. 華北水利水電大學(xué)地球科學(xué)與工程學(xué)院,鄭州 450046
2. 自然資源部天然氣水合物重點(diǎn)實(shí)驗(yàn)室,青島海洋地質(zhì)研究所,青島 266237
3. 青島海洋科學(xué)與技術(shù)國(guó)家實(shí)驗(yàn)室海洋礦產(chǎn)資源評(píng)價(jià)與探測(cè)技術(shù)功能實(shí)驗(yàn)室,青島 266237
全波形反演是一種綜合運(yùn)用地震波場(chǎng)的運(yùn)動(dòng)學(xué)和動(dòng)力學(xué)信息,通過不斷匹配正演波場(chǎng)與實(shí)測(cè)波場(chǎng)來調(diào)整模型,以獲得地下準(zhǔn)確速度場(chǎng)的反演方法。20世紀(jì)80年代, Tarantola[1]最早提出了基于廣義最小二乘法的時(shí)域波動(dòng)方程波形反演;隨后Pratt等[2-5]將全波形反演理論推廣到頻率域,形成了頻率域全波形反演方法;Shin等[6]針對(duì)常規(guī)時(shí)域和頻域方法無法高精度恢復(fù)速度場(chǎng)長(zhǎng)波長(zhǎng)分量等問題,提出了拉普拉斯域波形反演,并將其改進(jìn)到拉普拉斯-傅里葉域。
作為目前理論最先進(jìn)且精度最高的建模、成像手段,全波形反演被國(guó)內(nèi)外眾多學(xué)者廣泛應(yīng)用于海陸地震資料處理中。20世紀(jì)80年代,Gauthier等[7]和 Mora[8]最早實(shí)現(xiàn)了二維模擬地震數(shù)據(jù)的全波形反演,驗(yàn)證了理論可行性;90年代,Pratt[2]和Song等[9]利用寬角透射波信息重構(gòu)速度模型,進(jìn)一步發(fā)展了全波形反演方法;2004 年,Ravaut 等[10]將全波形反演應(yīng)用于深部地層,實(shí)現(xiàn)了整個(gè)巖石圈的反演成像;2009 年,Smithyman 等[11]提出了利用大偏移距、折射信息來進(jìn)行波形反演,獲得了優(yōu)于旅行時(shí)反演的結(jié)果;2010 年,Sirgue等[12]對(duì)挪威北海油田Valhall 地區(qū)的OBC數(shù)據(jù)進(jìn)行了全波形反演,對(duì)該區(qū)的淺層低速氣云進(jìn)行了準(zhǔn)確描述并對(duì)周邊充氣斷裂構(gòu)造進(jìn)行了精細(xì)刻畫;隨后全波形反演陸續(xù)應(yīng)用于海底電纜(OBC)、常規(guī)拖纜、變深度拖纜等不同采集類型的海上地震資料[13-16]。
相對(duì)于海上資料波形反演開展的如火如荼,針對(duì)陸上地震資料的波形反演研究工作比較少。2012年,Plessix等[17]利用低頻、大偏移距觀測(cè)系統(tǒng)采集了一條二維地震數(shù)據(jù),從1.5 Hz的低頻開始進(jìn)行全波形反演,取得了較好的反演結(jié)果,驗(yàn)證了陸上資料開展全波形反演的可行性,但由于這種觀測(cè)系統(tǒng)在實(shí)際采集中推廣難度較大,因此不具有普遍性。楊勤勇等[18]分析認(rèn)為陸上資料全波形反演的瓶頸在于:①全波場(chǎng)信息的獲取嚴(yán)重不足。全波形反演要求全方位角、大偏移距的觀測(cè)系統(tǒng),而陸上地震勘探多是基于反射波勘探的滾動(dòng)采集,偏移距較小、地震波傳播路徑正交程度差,限制了全波場(chǎng)信息的獲?。虎诘皖l信息嚴(yán)重缺失。海上數(shù)據(jù)去除海洋噪聲干擾后的最低可用頻帶為3.0~3.5 Hz,傳統(tǒng)陸上檢波器的頻帶有效范圍為6 Hz以上;③陸地資料信噪比偏低。波形反演對(duì)噪聲敏感度高,而陸地資料通常干擾波復(fù)雜,信噪比明顯偏低,給數(shù)據(jù)預(yù)處理提出了較高要求,制約了全波形反演的適用性;④誤差函數(shù)收斂效果差。陸上震源、檢波器等不確定因素導(dǎo)致波動(dòng)方程很難完全模擬全部波形,影響了誤差函數(shù)的收斂方向。鑒于以上情況,開展陸上全波形反演相關(guān)研究的難度是較大的,需要克服從理論到算法的多重困難。
綜合考慮全波形反演的適用性可以發(fā)現(xiàn),其與多個(gè)方面的因素有關(guān),如高信噪比(寬方位、大偏移距)的原始地震數(shù)據(jù)、精確的初始速度模型、準(zhǔn)確的震源子波、精細(xì)的正演模擬算法以及高性能的計(jì)算設(shè)備等[19]。這些因素涵蓋了從采集到處理、從數(shù)據(jù)到模型、從算法到硬件等多個(gè)方面。這充分說明全波形反演是一個(gè)復(fù)雜的系統(tǒng)工程,但其特性又展現(xiàn)出良好的研究潛力,如果能獲得突破,全波形反演將為地球物理勘探揭開新的篇章。
天然氣水合物是一種新型潔凈能源,其作為一種潛在的煤炭石油替代品得到了廣泛的研究。似海底反射(BSR)是識(shí)別天然氣水合物的最直觀證據(jù)之一,為了正確識(shí)別BSR,許多學(xué)者利用BSR上方(可能的水合物層)、下方(可能的游離氣層)沉積層存在特殊彈性參數(shù)結(jié)構(gòu)的特征,對(duì)天然氣水合物及游離氣的分布特征和形成機(jī)制進(jìn)行了分析,同時(shí)研究了BSR上、下方的地震速度結(jié)構(gòu)[20-21]。全波形反演作為一種高精度的速度建模及成像手段,在水合物的研究中發(fā)揮了重要作用。
Singh等[22]于1993年首次提出了全波形反演的方法,并將這一方法應(yīng)用于溫哥華岸外的天然氣水合物反射層速度結(jié)構(gòu)的研究中。他的研究結(jié)果也被大洋鉆探146航次證實(shí),此后全波形反演的方法得到了廣泛的應(yīng)用和發(fā)展。Yuan等[23]在溫哥華北部的Cascadia大陸邊緣對(duì)多道地震數(shù)據(jù)進(jìn)行了全波形反演,得到的速度結(jié)構(gòu)表明:含有水合物地層的速度比不含水合物或游離氣地層速度下降的更快,并估算出了水合物的濃度。Minshull等[24]提出的波形反演主要分三步,首先由均方根速度計(jì)算出主要反射層之間的層速度,然后用蒙特卡洛最大能量法獲得精確的層速度,最后在頻率-慢度域用共軛梯度法進(jìn)行波形的擬合,取得了良好的效果。Pecher等[25]對(duì)秘魯岸外的BSR精細(xì)速度結(jié)構(gòu)進(jìn)行了波形反演,并嘗試了氣體飽和度的估算方法。Sain等[26]對(duì)Makran大陸邊緣資料進(jìn)行的波形反演,揭示了海底500~800 m深度存在明顯的BSR,其下的低速帶表明下方沉積物中存在大量氣體(圖1),且此處的BSR與布萊克海臺(tái)處的地震特征類似,并分析這些游離氣是由于溫壓條件的改變?cè)斐伤衔锓€(wěn)定帶的基底上移,從而導(dǎo)致水合物分解而產(chǎn)生。Westbrook[27]對(duì)挪威Storegga 地區(qū)大陸邊緣的地震資料進(jìn)行了三維層析反演、二維射線追蹤反演以及一維波形反演,BSR上方的速度異常揭示了水合物的存在;2012年,Jaiswal[28]等在印度Krishna-Godavari盆地進(jìn)行了二維多道地震的頻率域波形反演,速度模型顯示有斑片狀的水合物存在。
國(guó)內(nèi)的天然氣水合物波形反演起步較晚。宋海斌最早針對(duì)日本東南海海槽增生楔的雙BSR速度結(jié)構(gòu)開展了研究,提出了與Singh類似的水平層狀彈性介質(zhì)的截距時(shí)間-水平慢度域波形反演思路[20],即首先根據(jù)τ-p域走時(shí)數(shù)據(jù)應(yīng)用非??焖倌M退火算法求得速度結(jié)構(gòu)的長(zhǎng)波長(zhǎng)部分,然后將第一步的速度模型作為初始模型,利用共軛梯度法求得速度的短波長(zhǎng)擾動(dòng)部分。唐勇等通過地震似海底反射層的識(shí)別,以BSR作為天然氣水合物穩(wěn)定帶的底,通過波形反演、AVO技術(shù)來確定水合物穩(wěn)定帶的上界[29],其研究的亮點(diǎn)在于他們還考慮了水合物穩(wěn)定帶的飽和度、孔隙度等一系列影響因素,來估算水合物穩(wěn)定帶的厚度。徐寧從理論上推導(dǎo)了基于遺傳算法的疊前全波形反演算法,全面系統(tǒng)地描述了計(jì)算流程,探討了反演中不同參數(shù)對(duì)反演收斂速度和算法穩(wěn)定性的影響[30],并基于前述算法及流程,以泥底辟型水合物為樣本予以驗(yàn)證,同時(shí)還研究了含水合物地層的速度、流體勢(shì)能、密度和波阻抗等屬性特征。 霍元媛等提出了一種基于遺傳算法的全波形反演方法[21](圖2),將縱波速度作為反演參數(shù),通過模型參數(shù)編碼,隨機(jī)生成初始模型種群,經(jīng)過遺傳算法的迭代,獲得了高分辨率的地層速度結(jié)構(gòu),并在中國(guó)南海的天然氣水合物研究中得以應(yīng)用,得到了分辨率高于常規(guī)速度分析的似海底反射層速度結(jié)構(gòu),并識(shí)別出似海底反射層的速度異常。Sun通過混合反演方法(疊前AVA反演、基于遺傳算法的全波形反演及疊后反演),獲取了可靠的地球物理屬性參數(shù),并基于此對(duì)水合物分解引起的海底滑坡這一現(xiàn)象進(jìn)行了數(shù)值模擬,隨后,他又從模擬結(jié)果出發(fā)研究了水合物分解形成氣煙囪的過程及地震響應(yīng)特征[31]。王吉亮利用 OBS 和多道地震數(shù)據(jù),進(jìn)行走時(shí)反演和全波形反演,并基于反演結(jié)果建立了墨西哥灣過 GC955 站位剖面的精細(xì)速度結(jié)構(gòu),解釋了水合物和游離氣的分布范圍,同時(shí)他還將獲得的速度結(jié)構(gòu)應(yīng)用于深度偏移,得到了過 GC955 站位的深度剖面,利用水道-天然堤沉積體系的空間展布,結(jié)合速度模型、鉆井測(cè)井資料,刻畫了水合物穩(wěn)定帶在區(qū)域內(nèi)的平面展布[32]。劉斌對(duì)稀疏OBS數(shù)據(jù)進(jìn)行了全波形反演,獲得了穩(wěn)定的高分辨率反演結(jié)果,能夠清晰地顯示出含水合物沉積層的速度異常[33]。
圖 1 BSR及其下方氣體在波形反演結(jié)果中的顯示[26]Fig.1 Display of BSR and its underlying gas in waveform inversion results[26]
圖 2 全波形反演流程圖[21]Fig.2 Workflow for prestack waveform inversion (PSWI) [21]
綜合分析前人研究成果,波形反演在海域地震數(shù)據(jù)方面的研究最早,也最深入。隨著海域天然氣水合物研究的不斷深入,波形反演方法也自然而然的與其結(jié)合,并取得了良好的發(fā)展。從波形反演的適用性和水合物的賦存特性等角度考慮,可從正演模擬與子波、反演模型及算法、預(yù)處理等幾個(gè)重點(diǎn)方面入手,進(jìn)一步提升水合物波形反演的精確度和有效性,加強(qiáng)波形反演方法在水合物探測(cè)及成藏模式等方面的研究。
地震正演模擬是對(duì)特定的地質(zhì)、地球物理問題做適當(dāng)?shù)暮?jiǎn)化形成數(shù)學(xué)模型,并采用數(shù)值計(jì)算的方法獲取地震響應(yīng)的過程。全波形反演是根據(jù)在檢波器處收集到的地震數(shù)據(jù)去模擬出地下介質(zhì)模型[34],因此,地震波正演模擬是全波形反演必不可少的前期步驟。要想兼顧地震采集方式、資料特點(diǎn)、算法穩(wěn)定性、計(jì)算效率等多個(gè)方面,并取得高分辨率的正演模擬效果,就要選擇合適的方法。
現(xiàn)有的正演模擬方法大致可以分為兩類:以射線理論為基礎(chǔ)的方法和以波動(dòng)方程理論為基礎(chǔ)的方法。波動(dòng)方程法又可以分為波動(dòng)方程的數(shù)值解法和解析算法兩類。以射線理論為基礎(chǔ)的方法主要有射線追蹤法、辛幾何法等;波動(dòng)方程的數(shù)值解法主要有有限差分法、有限元法、傅里葉偽譜法、有限體積法、廣義屏法等;解析算法主要有反射率法等。
以射線理論為基礎(chǔ)的方法計(jì)算速度快、易于實(shí)現(xiàn),但在處理橫向速度變化劇烈的介質(zhì)(如含水合物地層)時(shí),穩(wěn)定性明顯降低,且缺少地震波的動(dòng)力學(xué)信息,沒有考慮全波場(chǎng)特征,因此不能滿足全波形反演的需要。波動(dòng)方程數(shù)值解法不僅能保持地震波的運(yùn)動(dòng)學(xué)特征,而且還能保持地震波的動(dòng)力學(xué)特征,能夠較為精確地模擬任意復(fù)雜介質(zhì)的地震波波場(chǎng),但是計(jì)算量大,計(jì)算效率較低。
Kennett提出了計(jì)算分層均勻介質(zhì)所有響應(yīng)的廣義反射透射系數(shù)矩陣方法[35]。該方法的基本思想是首先利用Fourier-Hankel變換,將時(shí)間-空間域內(nèi)的波動(dòng)方程轉(zhuǎn)化到頻率-慢度域,然后根據(jù)實(shí)際問題求解,得到震源引起的響應(yīng),最后對(duì)該響應(yīng)進(jìn)行Fourier-Hankel積分反變換完成時(shí)間-空間域內(nèi)介質(zhì)響應(yīng)的計(jì)算。隨后,Kennett等進(jìn)行了一系列的改進(jìn)使得廣義反射透射系數(shù)矩陣法的基本理論已經(jīng)比較完善[36]。廣義反射透射系數(shù)矩陣方法能夠合成分層均勻介質(zhì)的完全響應(yīng),允許反演過程中考慮調(diào)諧效應(yīng),來自多次波和轉(zhuǎn)換波的干擾,以及在存在速度梯度的情況下的反射和透射效應(yīng)[28],且計(jì)算效率高,因此大部分的水合物波形反演都采用該方法進(jìn)行正演計(jì)算[20-21,37]。
但是廣義反射透射系數(shù)矩陣法只適用于垂向水平層狀非均勻介質(zhì),當(dāng)水合物賦存區(qū)域構(gòu)造復(fù)雜時(shí),該方法就無法準(zhǔn)確描述波場(chǎng)。隨著計(jì)算機(jī)性能的逐漸提高,不少學(xué)者開展了波動(dòng)方程法的正演模擬相關(guān)研究。Marfurt對(duì)有限差分法和有限元法的計(jì)算結(jié)果分別在時(shí)間域和頻率域進(jìn)行了對(duì)比,認(rèn)為在頻率域有限元法的解最為精確,但是最耗機(jī)時(shí)[38]。何彥鋒等從正演模擬的計(jì)算精度、計(jì)算效率、頻散特性以及適用性等方面對(duì)邊界元法和有限差分法進(jìn)行了比較研究,認(rèn)為在處理內(nèi)部非均質(zhì)和高頻計(jì)算時(shí),有限差分法更有效[39]。Virieux首次使用二階交錯(cuò)網(wǎng)格有限差分方法進(jìn)行波場(chǎng)正演模擬[40],隨后這一方法得到不斷發(fā)展,目前,高階交錯(cuò)網(wǎng)格有限差分法已經(jīng)廣泛應(yīng)用于全波形反演。Barnes等使用軸向?qū)ΨQ的有限差分法來正演模擬各項(xiàng)異性介質(zhì)中的彈性波傳播,對(duì)日本南海海槽的天然氣水合物開展了研究[41];Pratt等采用有限差分法對(duì)加拿大Mallik地區(qū)的水合物進(jìn)行了波形反演,獲得了高分辨地震速度結(jié)構(gòu)[42];Jaiswal等使用了頻率域全波形反演,在印度Krishna-Godavari盆地進(jìn)行了縱波速度與衰減系數(shù)反演,從而確定了天然氣水合物分布范圍[28]。除此以外,還有不少有代表性的正演方法涌現(xiàn),如基于三維27點(diǎn)加權(quán)平均算子的粘聲波正演法[43]、省資源有限體積法[44]、不連續(xù)的 Galer-kin方法[45],張偉等[46]將GPU/CPU 并行計(jì)算技術(shù)應(yīng)用到全波形反演中,極大地提高了計(jì)算效率。Kim等采用聲波-彈性波耦合波動(dòng)方程進(jìn)行有限元法的正演,對(duì)日本海Ulleung盆地的天然氣水合物帶的特征和分布進(jìn)行了調(diào)查[47]。Guasch對(duì)比了聲波和彈性波模擬波場(chǎng)信號(hào),認(rèn)為聲波全波形反演在參數(shù)合適的時(shí)候仍然可以獲得較好的結(jié)果[48]。從對(duì)比圖(圖3)中可以明顯看出彈性波波場(chǎng)記錄比聲波波場(chǎng)記錄更復(fù)雜。盡管震源是純炸藥震源,接收器是水下檢波器,但雙重轉(zhuǎn)換的P波和其他彈性效應(yīng)為聲波增加了顯著的尾波,并改變了它們的 AVO 及其他振幅特性,但它們并不會(huì)影響主要縱波的旅行時(shí)。Guasch的研究也說明,在參數(shù)化合適的情況下,聲波全波形反演仍然可以獲得理想的反演結(jié)果。
能否獲得準(zhǔn)確的震源子波是全波形反演較為關(guān)鍵的步驟。何兵紅等通過數(shù)值示例表明子波相位變化引起波形變化,會(huì)影響全波形反演的梯度計(jì)算[49]。同時(shí)子波主頻也是影響全波形反演的重要因素,主頻過大或過小都會(huì)引起反演結(jié)果的不準(zhǔn)確。全波形反演對(duì)子波振幅不是特別敏感,但當(dāng)振幅過大時(shí)會(huì)干擾梯度的求取。為了使振幅匹配,需要將子波和原始數(shù)據(jù)進(jìn)行組合校正。 Crutchley的做法是將選定位置處一定范圍內(nèi)的CMP道集進(jìn)行時(shí)移,使得海底波峰振幅對(duì)齊,然后將這些道集疊加;再對(duì)首個(gè)海底多次波應(yīng)用同樣的處理,并將其極性反轉(zhuǎn),與海底極性相同;然后將所有這些疊加后的海底子波與海底多次子波求平均,以獲得最終的震源子波[50]。Xia等則認(rèn)為反演過程中應(yīng)同時(shí)考慮子波,而不是將子波固定[51]。還有的學(xué)者為了克服震源子波的影響,提出了不依賴于震源子波的全波形反演方法,如頻率域的卷積波場(chǎng)和反卷積波場(chǎng)方法[52-58]。
敖瑞德等通過對(duì)Marmousi速度模型的數(shù)值實(shí)驗(yàn),發(fā)現(xiàn)即使采用相對(duì)精確的高斯平滑初始模型,但使用錯(cuò)誤的子波,常規(guī)的全波形反演仍然無法得到好的反演結(jié)果[59](圖4)。秦寧總結(jié)了獲得震源子波的兩類方法:一類是利用原始地震數(shù)據(jù)進(jìn)行子波估計(jì);另一類是直接從原始地震數(shù)據(jù)中反演子波,或者進(jìn)行速度與子波的多參數(shù)全波形反演[60]。
圖 3 Marmousi 模型的聲波方程正演(左)和彈性波方程正演(右)[48]Fig.3 A gather generated by a single source in acoustic (left) and elastic (right) versions of the 3D Marmousi model [48]
圖 4 Marmousi 準(zhǔn)確模型(上)以及子波錯(cuò)誤時(shí)的全波形反演結(jié)果(下)[59]Fig.4 The Marmousi model (up) and FWI result with an incorrect wavelet (down) [59]
圖 5 波形反演對(duì)初始模型的靈敏度對(duì)比a. 平滑后的真實(shí)模型,b. 利用FATT得到的模型,c. 立體層析得到的模型[64]。Fig.5 Initial models for FWIa. Smoothed version of the true model, b. FATT model, and c stereotomography model [64].
子波估計(jì)的準(zhǔn)確性,會(huì)嚴(yán)重影響反演過程中數(shù)據(jù)的擬合,直接決定反演成功與否。在天然氣水合物的波形反演中,大多數(shù)學(xué)者都假設(shè)海底子波是震源子波和反射系數(shù)序列的褶積,并從海底反射及其首次多次波中提取震源子波[24-26],以期獲得較好的擬合結(jié)果。
全波形反演由Tarantola[1]最早提出,經(jīng)歷了時(shí)間域的多尺度反演方法、頻率域全波形反演、頻率迭代的多尺度全波形反演、拉普拉斯域的全波形反演等發(fā)展階段,由時(shí)間域到頻率域,再到拉普拉斯域,由海洋到陸地,由2D到3D不斷進(jìn)步并實(shí)用化。在這個(gè)過程中,全波形反演方法的發(fā)展曾遭遇不少障礙,其中最大的瓶頸就是實(shí)際資料中低頻分量與初始速度模型的耦合問題[18]。早在1986年,Gauthier等最早對(duì)二維地震資料進(jìn)行全波形反演時(shí),就指出了全波形反演對(duì)于初始模型的嚴(yán)重依賴[7]。Virieux等認(rèn)為受初始模型精度及低頻缺失的影響,局部?jī)?yōu)化算法無法阻止目標(biāo)函數(shù)向局部最小值收斂[61]。Nikhil等通過模型實(shí)例揭示了一個(gè)成功的全波形反演需要足夠的低頻成分和精確的初始速度模型的耦合[62]。由于全波形反演常受周波躍遷問題困擾,一個(gè)好的初始模型能夠幫助波形反演跳出局部最優(yōu)解的陷阱。秦寧認(rèn)為初始模型的準(zhǔn)確程度決定了全波形反演的成功與否,通過對(duì)不同初始模型進(jìn)行反演,發(fā)現(xiàn)平滑速度模型的反演結(jié)果優(yōu)于常梯度速度模型,分析其原因是速度場(chǎng)的低頻成分需要在初始速度模型中給予準(zhǔn)確的表達(dá),平滑后的速度模型低頻成分保留的較好,高頻成分在反演迭代的過程中可以得到有效恢復(fù),而常梯度速度模型低頻成分的準(zhǔn)確度低會(huì)使得反演無法收斂到全局最優(yōu)解[60]。因此,在實(shí)際資料的應(yīng)用中,一個(gè)好的初始模型是全波形反演成功的前提條件。
Virieux等總結(jié)了三種全波形反演初始模型的建立方式:初至旅行時(shí)層析成像、立體層析成像和拉普拉斯域反演[61]。
初至旅行時(shí)層析成像是指利用各種類型地震波的初至旅行時(shí)信息來反演地下介質(zhì)速度。它只利用地震波的旅行時(shí)信息,數(shù)據(jù)計(jì)算量小,并且旅行時(shí)反演的非線性程度相對(duì)較弱,正演模擬和反演迭代相對(duì)穩(wěn)定,因此,在實(shí)際生產(chǎn)中的應(yīng)用最為廣泛。Brenders等實(shí)現(xiàn)了旅行時(shí)層析與全波形的聯(lián)合反演,認(rèn)為旅行時(shí)層析所獲得的初始模型的精度、豐富的低頻成分以及大偏移距信息為全波形反演提供了重要的基礎(chǔ)[63]。Prieux等評(píng)估了二維頻率域全波形反演的敏感性,包括對(duì)采集到的最小偏移距、對(duì)自由表面效應(yīng)以及對(duì)初至旅行時(shí)反演得到的初始模型和立體層析成像得到的初始模型精度等各方面影響因素[64](圖5)。評(píng)估結(jié)果表明,初至層析能夠幫助波形反演得到比較好的淺層模型,而立體層析能夠有效恢復(fù)速度的長(zhǎng)波長(zhǎng)分量。
但是當(dāng)實(shí)際資料中存在低速帶時(shí),如何拾取可靠的初至旅行時(shí)是一個(gè)難題,因此,Billette等提出了立體層析成像技術(shù)[65-69]。立體層析成像是一種基于反演數(shù)據(jù)體的旅行時(shí)間和局部相干同相軸斜率的層析成像方法,與常規(guī)旅行時(shí)反射層析成像相比,立體層析成像的最大特點(diǎn)是數(shù)據(jù)空間的分布可以是稀疏的、不連續(xù)的[65],它可以半自動(dòng)拾取局部相干同相軸,且更簡(jiǎn)便,拾取密度更高。這種方法除了利用地震波走時(shí),還使用炮、檢點(diǎn)位置與炮、檢點(diǎn)處射線的局部傳播方向來約束速度模型[66],比傳統(tǒng)層析成像算法更穩(wěn)定。劉定進(jìn)等采用高斯束層析與全波形反演進(jìn)行聯(lián)合速度建模,取得了較好的效果[70]。Prieux等聯(lián)合折射和反射立體層析成像建立了頻率域全波形反演初始模型,并應(yīng)用于Valhall油田資料的全波形反演,為研究低速氣層和儲(chǔ)層下方深部構(gòu)造提供了更可靠的初始模型[64]。
拉普拉斯域波形反演算法穩(wěn)健并且不依賴于初始模型,能夠給頻率域波形反演提供一個(gè)光滑的初始速度模型,Shin等將這種波形反演與直接在頻率域進(jìn)行的波形反演在多種模型上進(jìn)行了對(duì)比,認(rèn)為將拉普拉斯域波形反演結(jié)果作為初始模型的反演結(jié)果更好[71]。雖然拉普拉斯域波形反演能夠恢復(fù)模型的低頻成分,但是反演深度容易受偏移距和衰減的影響,Shin等利用阻尼波場(chǎng)的低頻分量,提出了 Laplace-Fourier波形反演算法,并將其應(yīng)用于頻率小于地震帶寬的最小頻率的資料[71]。劉張聚建立了一種基于時(shí)域加權(quán)的拉普拉斯-頻率域彈性波全波形反演方法,他通過引入拉普拉斯衰減因子降低了全波形反演對(duì)于低頻成分的依賴性,形成了兼具三種域優(yōu)勢(shì)的時(shí)間-拉普拉斯-頻率域彈性波全波形反演方法[72]。
據(jù)王連坤等[73]總結(jié),除走時(shí)層析、拉普拉斯域反演以外,偏移速度分析也是主要的初始模型建立方法。偏移速度分析是利用波場(chǎng)傳播算子產(chǎn)生擴(kuò)展的成像道集,并通過一定的判定準(zhǔn)則來獲取最佳偏移速度的方法[19]。其中,波場(chǎng)傳播算子可以是動(dòng)校正疊加算子、傾斜校正疊加算子,也可以是克?;舴蚱扑阕?、疊前時(shí)間偏移算子等;成像道集可以是角度域,也可以是地下偏移距域等;判定準(zhǔn)則可以是道集拉平、能量聚焦或能量疊加最大化。這種方法需要在每次迭代更新速度模型后進(jìn)行疊前深度偏移,因此計(jì)算成本較高,國(guó)內(nèi)外學(xué)者很少利用偏移速度分析方法來構(gòu)建全波形反演的初始模型,其主要還是用于偏移成像領(lǐng)域。
針對(duì)天然氣水合物的波形反演,出于計(jì)算效率和保證算法收斂的考慮,一般都使用旅行時(shí)層析算法先獲得速度模型的長(zhǎng)波長(zhǎng)分量,然后再用波形反演重建該模型的細(xì)節(jié)部分,得到模型的短波長(zhǎng)分量,這樣就可以得到一個(gè)比較精確的速度模型[20-21,24-25,37,50]。
Brossier等通過海上、陸上含噪模型數(shù)據(jù)的計(jì)算,對(duì)比了最小平方L2范數(shù)、最小絕對(duì)值L1范數(shù)、Huber 準(zhǔn)則和混合準(zhǔn)則等4種最小化函數(shù)在頻率域波形反演中的表現(xiàn)[74](圖6),其中Huber 準(zhǔn)則和混合準(zhǔn)則是L1范數(shù)和L2范數(shù)二者的聯(lián)合函數(shù)。計(jì)算結(jié)果表明,L1范數(shù)對(duì)噪音不太敏感,能提供最可靠的速度模型;L2范數(shù)在存在均勻白噪的情況下,如果能犧牲計(jì)算效率而減小頻率采樣間隔也能夠得到可靠的結(jié)果;而Huber 準(zhǔn)則和混合準(zhǔn)則在轉(zhuǎn)換門閥選擇合理的情況下,也能夠提供令人滿意的結(jié)果。Liu等通過對(duì)互相關(guān)函數(shù)與L2范數(shù)等不同目標(biāo)函數(shù)的對(duì)比,闡明了互相關(guān)函數(shù)作為目標(biāo)函數(shù)的優(yōu)勢(shì)[75]。Guasch等將維納濾波器引入目標(biāo)函數(shù),并將該方法推廣到三維,在缺失低頻信息的情況下獲得了很好的反演結(jié)果[48]。Shin等提出利用對(duì)數(shù)波場(chǎng)來構(gòu)建目標(biāo)函數(shù),其允許考慮三種類型的目標(biāo)函數(shù),即只有振幅、只有相位或二者皆有的目標(biāo)函數(shù),在對(duì)不同數(shù)據(jù)進(jìn)行常規(guī)目標(biāo)函數(shù)和對(duì)數(shù)波場(chǎng)目標(biāo)函數(shù)的波形反演后發(fā)現(xiàn),對(duì)于合成數(shù)據(jù)來說,新方法獲得的速度模型要優(yōu)于傳統(tǒng)波形反演;對(duì)實(shí)際地震數(shù)據(jù)應(yīng)用新方法進(jìn)行反演,則能更好地揭示近地表的高頻特征[71]。此后,Shin等又使用對(duì)數(shù)目標(biāo)函數(shù)方法的變體即全對(duì)數(shù)方法進(jìn)行全波形反演,并將結(jié)果與常規(guī)最小二乘法進(jìn)行比較,發(fā)現(xiàn)對(duì)數(shù)反演更加穩(wěn)健,分析原因可能是對(duì)數(shù)方法對(duì)殘余波場(chǎng)的振幅產(chǎn)生了自然調(diào)節(jié),其使計(jì)算穩(wěn)定,從而改善反演結(jié)果[76]。Bednar等的研究認(rèn)為,基于相位的對(duì)數(shù)方法比傳統(tǒng)方法更實(shí)用[77]。Pyun等分別使用僅考慮振幅的對(duì)數(shù)波場(chǎng)目標(biāo)函數(shù)和常規(guī)波場(chǎng)目標(biāo)函數(shù)進(jìn)行波形反演,認(rèn)為后者的反演結(jié)果較好,但仍比對(duì)數(shù)波場(chǎng)目標(biāo)函數(shù)和基于相位的對(duì)數(shù)波場(chǎng)目標(biāo)函數(shù)的結(jié)果要差,因此不建議反演中單獨(dú)使用基于振幅的對(duì)數(shù)波場(chǎng)目標(biāo)函數(shù)[78]。廖文等提出改進(jìn)算法中固定不變的目標(biāo)函數(shù),在反演過程中有針對(duì)性地調(diào)整目標(biāo)函數(shù),以達(dá)到在不增加計(jì)算量的前提下降低頻率域全波形反演算法對(duì)初始模型的依賴,提高算法的適用性的目的[79]。
圖 6 不同目標(biāo)函數(shù)對(duì)含噪數(shù)據(jù)的縱波速度全波形反演結(jié)果a、b: L2范數(shù), c、d: L1范數(shù),e、f: Huber 準(zhǔn)則,g、h: 混合準(zhǔn)則[74]。Fig.6 Reconstructed left VP and right VS models for the first Valhall test with the noisy data after the two FWI steps.a, b: L2-norm; c, d: L1-norm; e, f: Huber criterion; g, h: hybrid criterion[74] .
時(shí)間域的波形反演主要有兩類目標(biāo)函數(shù):一類是將目標(biāo)函數(shù)定義為最小化預(yù)測(cè)波場(chǎng)與觀測(cè)波場(chǎng)之間的誤差[51],另一類是將目標(biāo)函數(shù)定義為規(guī)則化的互相關(guān)函數(shù)[76]。最早Singh進(jìn)行的全波形反演是最小化采樣點(diǎn)間的實(shí)際地震記錄和合成地震記錄的差值,獲得了BSR層的反射速度結(jié)構(gòu)和水合物層厚度,其結(jié)果由ODP驗(yàn)證[22]。張寶金等引入了參數(shù)和觀測(cè)數(shù)據(jù)的先驗(yàn)概率信息,以實(shí)際數(shù)據(jù)和正演計(jì)算數(shù)據(jù)殘差的加權(quán)能量為目標(biāo)函數(shù),對(duì)中國(guó)南海北部陸坡水合物鉆探區(qū)的資料進(jìn)行了波形反演[80]。徐寧則認(rèn)為傳統(tǒng)的誤差目標(biāo)函數(shù)收斂速度較慢,而相關(guān)形式的目標(biāo)函數(shù)對(duì)波形的絕對(duì)旅行時(shí)要求不是很嚴(yán),并且具有一定的抗噪能力[30],因此選擇和Sen等相同的目標(biāo)函數(shù)[26],即觀測(cè)數(shù)據(jù)與理論數(shù)據(jù)的相關(guān)函數(shù)。宋海斌等采用τ-p域走時(shí)計(jì)算值與觀測(cè)值的均方差作為目標(biāo)函數(shù),反演了日本東南海海槽雙BSR的速度結(jié)構(gòu)[20]。
現(xiàn)有的地球物理反演方法主要有兩大類:一類是線性反演方法,另一類是非線性反演方法。線性反演主要是指Born近似和Rytov近似這類反演方法,這種方法僅取前向小角度輻射作近似處理,忽略了前后向波長(zhǎng)耦合及多次反射及折射效應(yīng)等構(gòu)造間的相互作用[81],因此不適用于反演疊前地震數(shù)據(jù)。非線性反演方法主要用于解決不確定性高、非線性程度高、計(jì)算規(guī)模大的優(yōu)化問題。按是否具備全局搜索能力,可以進(jìn)一步將非線性反演方法分為兩類[82]:①基于目標(biāo)函數(shù)局部微分特性的局部方法;②基于隨機(jī)搜索及兼有智能性的全局方法。其中局部?jī)?yōu)化方法主要有梯度類方法,包括最速下降法、共軛梯度法等。另外,常見的還有牛頓類方法,主要包括牛頓法、擬牛頓法、高斯-牛頓法等。全局優(yōu)化方法主要有蒙特卡洛法、遺傳算法、模擬退火、神經(jīng)網(wǎng)絡(luò)等。
最速下降法等梯度類方法首先對(duì)梯度求解,獲得目標(biāo)函數(shù)的最快上升方向,然后沿負(fù)梯度方向?qū)Τ跏寄P瓦M(jìn)行迭代更新。這種方法具有良好的全局收斂性,所需要的存儲(chǔ)量比較小,結(jié)構(gòu)簡(jiǎn)單,易于實(shí)現(xiàn),Pratt等將最速下降法應(yīng)用于頻率域聲波方程和彈性波方程的全波形反演,并對(duì)該方法進(jìn)行了可靠性評(píng)估[4]。由于梯度類算法利用的是目標(biāo)函數(shù)的一階偏導(dǎo),所以其搜索過程收斂較慢。牛頓類算法不僅利用了目標(biāo)函數(shù)的二階導(dǎo)數(shù),還利用了其梯度信息,因此在收斂速度上比梯度類算法更有優(yōu)勢(shì)[83-85],Tarantola在將牛頓法應(yīng)用到時(shí)間域全波形反演算法的過程中發(fā)現(xiàn),牛頓法不僅需要計(jì)算目標(biāo)函數(shù)的二階導(dǎo)數(shù),還要求海森矩陣必須滿足非奇異性[1],這也對(duì)計(jì)算提出了較高要求。高斯-牛頓法不需要對(duì)海森矩陣精確求解,因此比牛頓法的計(jì)算效率高,Hu等在進(jìn)行二維多頻率全波形反演時(shí)應(yīng)用了此方法[86],不但證明了方法的可靠性,還提供了可參考的工作方案。擬牛頓法則只需要計(jì)算目標(biāo)函數(shù)的一階導(dǎo)數(shù),具有良好的全局收斂性和超線性收斂速度。劉璐等將該方法應(yīng)用于二維頻率域全波形反演,雖然成功求取了速度參數(shù)[87],但同時(shí)也暴露出擬牛頓法的劣勢(shì),即在每次迭代過程中都需要構(gòu)建逆海森矩陣,而且需要求解線性方程組來確定下降方向,存儲(chǔ)量巨大,因此擬牛頓法只適用于中小規(guī)模問題,受此限制,直接將該方法應(yīng)用于全波形反演的研究較少。共軛梯度法是利用共軛性對(duì)最速下降法進(jìn)行修正來構(gòu)造更好的搜索方向,既不用計(jì)算海森矩陣及其逆矩陣,又不需要目標(biāo)函數(shù)的二階導(dǎo)數(shù)信息,因此大大提高了收斂速度,是全波形反演最常用的方法之一。高鳳霞等的研究很具有代表性,他們將幾種局部?jī)?yōu)化方法應(yīng)用于同一個(gè)數(shù)據(jù)的頻率域二維聲波方程速度參數(shù)全波形反演中(圖7),從收斂速度、反演精度以及計(jì)算時(shí)間等方面進(jìn)行了對(duì)比分析,認(rèn)為高斯-牛頓法收斂速度最快,但穩(wěn)定性稍差,最速下降法的收斂速度最慢;共軛梯度方法反演精度較高,梯度法和高斯-牛頓法的精度較低,最速下降法精度最低;最速下降方法和共軛梯度法只需計(jì)算梯度,計(jì)算用時(shí)最少[88]。
圖 7 不同優(yōu)化算法得到的全波形反演速度模型a. 最速下降法,b. 對(duì)角海森尺度化方法,c. 共軛梯度法,d. 對(duì)角海森+共軛梯度法,e. 高斯-牛頓法,f. L-BFGS方法[88]。Fig.7 FWI results of (a) steepest-descent method ,(b) diagonal Hessian matrix scaled method, (c) conjugate gradient method, (d) diagonal Hessian matrix scaled gradient + conjugate gradient method, (e) Gauss-Newton method and, (f) L-BFGS method [88]
圖 8 預(yù)處理流程圖Fig.8 Workflow for precondition
局部?jī)?yōu)化方法收斂速度快,計(jì)算用時(shí)少,但反演結(jié)果強(qiáng)烈依賴于初始模型,對(duì)于海量的地球物理數(shù)據(jù)而言,選取一個(gè)好的初始模型絕非易事。全局優(yōu)化方法不需要求導(dǎo)等運(yùn)算,不依賴于初始模型,在解空間內(nèi)隨機(jī)搜索以獲得全局最優(yōu)解,但是其隨機(jī)搜索的特征也導(dǎo)致算法收斂慢、時(shí)間代價(jià)高。因此如何將兩者的優(yōu)點(diǎn)結(jié)合起來是目前反演問題的一個(gè)重要研究方向。
Singh等的工作就很好地印證了這一點(diǎn),他們用蒙特卡洛法來尋找全局最優(yōu)解,以避免陷入局部最優(yōu)解的陷阱,但其搜索方向完全隨機(jī),單參數(shù)反演時(shí)尚可,多參數(shù)同時(shí)反演時(shí)計(jì)算效率極其低下[22]。Cary等針對(duì)上述問題,提出先用蒙特卡洛法確定全局最優(yōu)解的范圍,然后用局部?jī)?yōu)化方法在這一范圍內(nèi)進(jìn)一步搜索獲得全局最優(yōu)解[89]的思路,從反演結(jié)果來看,在指明全局優(yōu)化算法的收斂方向上,該方法有不錯(cuò)的表現(xiàn),但由于蒙特卡洛法依然是一種高算力需求的方法,在降低全局優(yōu)化算法時(shí)間代價(jià)的效果方面仍比較有限。
遺傳算法是模擬生物在自然環(huán)境中的遺傳和進(jìn)化過程而形成的一種自適應(yīng)全局搜索優(yōu)化算法[90]。Mallick提出一種基于遺傳算法的全波形反演算法,應(yīng)用于美國(guó)含氣砂巖的研究[91],隨后Mallick等將疊前基于遺傳算法的全波形反演與疊后反演相結(jié)合提出一種混合反演方法,成功用于印度安達(dá)曼海的水合物勘探[92]。徐寧建立了基于遺傳算法的疊前全波形反演流程,并對(duì)算法的性能進(jìn)行了分析[30]。Huo等提出了基于遺傳算法和共軛梯度法的全波形反演,用于對(duì)中國(guó)南海北部海域的天然氣水合物研究[21]。另外,模擬退火法和神經(jīng)網(wǎng)絡(luò)法在全波形反演的尋優(yōu)過程中也取得了很好的效果[93]。比如宋海斌等采用快速模擬退火和共軛梯度相結(jié)合的思路,先求取速度場(chǎng)的低頻背景,再求取高頻分量,對(duì)日本東南海海槽的水合物速度結(jié)構(gòu)開展了研究[20]。而劉鵬程、Fallat等學(xué)者將模擬退火、神經(jīng)網(wǎng)絡(luò)與單純形法、均勻設(shè)計(jì)等方法結(jié)合起來,有效地提高了算法的效率[94-95]。
主流的優(yōu)化方法有很多,每一個(gè)方法都有其自身的優(yōu)缺點(diǎn)及不同的適用條件。隨著算法的發(fā)展,這些優(yōu)化方法發(fā)展出了很多分支,也都各具特點(diǎn),并不存在絕對(duì)的優(yōu)劣之分。實(shí)際使用的時(shí)候,需要根據(jù)具體情況,靈活處理。在天然氣水合物的全波形反演中,無論在時(shí)間域還是頻率域,由于實(shí)際地震資料常常缺少低頻成分,采用常規(guī)方法難以得到理想的初始模型,因此采用全局優(yōu)化方法與局部?jī)?yōu)化方法相結(jié)合的混合算法,兼顧計(jì)算效率與反演精度,往往能取得高質(zhì)量的反演結(jié)果。大量的實(shí)踐證明,這種折衷的辦法能更有效地促進(jìn)全波形反演的實(shí)際應(yīng)用。
全波形反演的預(yù)處理主要指對(duì)原始地震數(shù)據(jù)的處理,主要包括靜校正、多域保幅去噪、振幅補(bǔ)償、歸一化處理等,其目的是為了使原始數(shù)據(jù)與正演模擬數(shù)據(jù)能夠更好地匹配(圖8)。由于全波形反演對(duì)數(shù)據(jù)中存在的噪聲非常敏感,因此在進(jìn)行反演之前做好預(yù)處理工作至關(guān)重要,既要做到最大程度地消除噪音的影響,又不能破壞波的動(dòng)力學(xué)特征,實(shí)施難度很大。
海洋地震數(shù)據(jù)通常受多次波的影響最大,同時(shí),受采集條件的限制,目前常規(guī)開展二維地震資料反演的最高頻率通常在30 Hz以下,三維理論模型反演的最高頻率在20 Hz以下,也就是說全波形反演無法恢復(fù)到全頻帶,從而影響了對(duì)模型的真實(shí)重現(xiàn)。此外,受野外采集觀測(cè)系統(tǒng)設(shè)計(jì)的限制,原始地震資料是有限觀測(cè)孔徑的數(shù)據(jù)[96],這給波形反演帶來了很大困難。鑒于上述情況,如研究目的包含波形反演相關(guān)工作,在施工設(shè)計(jì)中,應(yīng)酌情增加長(zhǎng)排列、大偏移距的采集方式,將為波形反演提供豐富的低頻成分,有助于反演得到較好的背景速度場(chǎng),另外,寬方位采集還可以增加數(shù)據(jù)的完備性,有利于全波形反演。
近年來,隨著原始地震數(shù)據(jù)采集質(zhì)量以及超算能力的大幅提升,全波形反演逐漸由理論走向?qū)嵺`、從二維走向三維,尤其在天然氣水合物的研究中發(fā)揮了越來越重要的作用。從大量的實(shí)踐案例中可以得出結(jié)論,全波形反演作為一種高精度的速度建模及成像方法,在水合物的研究中,特別是在揭示BSR附近沉積層的地震速度、彈性參數(shù)等方面,展現(xiàn)出很好的適用性和優(yōu)越性。同時(shí),由于地震波傳播的復(fù)雜性,全波形反演仍舊面臨許多難題,如原始資料缺乏低頻成分、復(fù)雜介質(zhì)中地震波傳播難以準(zhǔn)確描述、震源子波估計(jì)不準(zhǔn)、初始模型建立精度偏低、目標(biāo)函數(shù)適用性不廣、尋優(yōu)算法性能仍待提高等。這些都嚴(yán)重制約了全波形反演的快速發(fā)展,因此全波形反演目前仍無法替代走時(shí)層析等傳統(tǒng)技術(shù)。
為了解決這些難題,許多學(xué)者提出了一些很有意義的思路、方法。如董良國(guó)等認(rèn)為,考慮到實(shí)際地震資料的復(fù)雜性,不一定要強(qiáng)求全部地震信息的匹配,可以根據(jù)不同勘探階段對(duì)反演結(jié)果精度要求的不同,構(gòu)建不同的目標(biāo)函數(shù),利用不同的地震數(shù)據(jù)子集進(jìn)行波形反演,以解決全波形反演的強(qiáng)非線性問題[97]。至于如何構(gòu)建目標(biāo)函數(shù)、選擇地震數(shù)據(jù)子集,取決于反演的階段和目的,以及資料的品質(zhì)、特點(diǎn)。王毓瑋等從彈性波有限差分正演模擬出發(fā),分析了多波多分量地震數(shù)據(jù)及其不同的數(shù)據(jù)子集隨縱、橫波速度,縱、橫波阻抗,密度以及拉梅常數(shù)等4種彈性參數(shù)的不同空間攝動(dòng)尺度的變化關(guān)系,并系統(tǒng)地分析了這些變化關(guān)系在不同目標(biāo)函數(shù)下的變化性態(tài),為彈性波全波形反演中目標(biāo)函數(shù)的選取與反演策略的制定提供了理論指導(dǎo)[98]。敖瑞德等首先從理論上分析了全波形反演不依賴震源子波的可行性,提出了不依賴于子波,基于包絡(luò)的全波形反演初始模型構(gòu)建方法,并通過數(shù)值試驗(yàn)驗(yàn)證了其正確性[59]。胡勇等在前人的研究基礎(chǔ)上[99-106],提出時(shí)頻域包絡(luò)反演方法,從而為全波形反演的初始速度模型構(gòu)建提供了一個(gè)有效途徑[107]。
在天然氣水合物的研究中,最初的一些全波形反演都是以縱波速度作為唯一的參數(shù),通過經(jīng)驗(yàn)公式求取橫波速度、密度等其他參數(shù)的。這種做法雖然具有一定的理論基礎(chǔ),但由于缺乏對(duì)理論值和實(shí)際值的差異評(píng)估,給反演結(jié)果帶來了很大的不確定性,因此利用實(shí)測(cè)數(shù)據(jù)開展多參數(shù)聯(lián)合反演必然成為未來的發(fā)展趨勢(shì)。Pratt等對(duì)加拿大Mallik地區(qū)的水合物層速度和Q值進(jìn)行了波形反演,獲得了高分辨地震速度結(jié)構(gòu)[42];Jaiswal等在印度Krishna-Godavari盆地進(jìn)行了縱波速度與衰減系數(shù)的聯(lián)合反演,從而確定了天然氣水合物分布范圍[28]。Nikhil等對(duì)挪威北海的Tommeliten Alpha野外資料實(shí)現(xiàn)了各向異性三維聲波全波形反演,這種方法與Ratcliffe[108]的方法原則上相似,可對(duì)空間可變、傾斜的橫向各向同性及各向異性介質(zhì)進(jìn)行反演,并且這種方法使用五參數(shù)對(duì)各向異性進(jìn)行參數(shù)化[62]。這些多參數(shù)聯(lián)合反演的研究為我們提供了全新的思路,可能成為未來全波形反演更廣泛應(yīng)用的突破口。當(dāng)然值得注意的是,多參數(shù)聯(lián)合反演對(duì)于反演算法的精度和穩(wěn)定性提出了更高的要求。全波形反演問題本身的高度非線性化,疊加海量地震數(shù)據(jù)的超大計(jì)算規(guī)模,對(duì)于計(jì)算機(jī)硬件的超高計(jì)算性能都是巨大的挑戰(zhàn),需要聯(lián)合多學(xué)科研究人員共同努力。