二維聲波頻率域正演

二維聲波頻率域正演方法研究目錄TOC\o"1-5"\h\z\o"CurrentDocument"1引言 1\o"CurrentDocument"2二維聲波頻率域正演原理 1\o"CurrentDocument"2.1頻率-空間域聲波方程的離散 1\o"CurrentDocument"PML邊界條件的引入 2\o"CurrentDocument"2.3大型稀疏矩陣方程組的求解 4\o"CurrentDocument"3頻率域聲波正演模擬實(shí)驗(yàn) 4\o"CurrentDocument"3.1均勻介質(zhì) 4\o"CurrentDocument"corner-edge模型 6\o"CurrentDocument"3結(jié)論與認(rèn)識(shí) 7\o"CurrentDocument"參考文獻(xiàn) 81引言隨著目前國內(nèi)外勘探開發(fā)的不斷深入，人們對(duì)高精度成像，精確地震反演等技術(shù)提出了更高的要求。全波形反演是近些年發(fā)展起來的一種反演方法，它能夠充分利用地震資料中的時(shí)間，振幅，相位等信息，具有高精度多參數(shù)建模的能力⑴。全波形反演既可以在時(shí)間域?qū)崿F(xiàn)，也可以在頻率域?qū)崿F(xiàn)，頻率域反演因?yàn)槠涮焐亩喑叨忍匦远艿綇V泛關(guān)注。頻率域反演的核心問題是高效的頻率域正演模擬方法，它決定了反演的精度和效率，因而近年來頻率域正演方法發(fā)展迅速。就數(shù)值模擬本身來說，頻率域數(shù)值模擬相對(duì)于時(shí)間域模擬也有著許多優(yōu)點(diǎn)。在頻率域計(jì)算波場的衰減效應(yīng)比時(shí)間域更方便。而且對(duì)于二維多炮數(shù)值模擬，頻率域模擬方法比時(shí)間域模擬方法更有效。另外，頻率域模擬方法是基于單個(gè)頻率對(duì)所有空間網(wǎng)格點(diǎn)進(jìn)行求解，誤差將分配到各個(gè)網(wǎng)格點(diǎn)。由于各個(gè)頻率間是獨(dú)立計(jì)算，所以不存在誤差的累計(jì)，不同頻率波場間可以方便地進(jìn)行并行計(jì)算，也可以考慮進(jìn)行長時(shí)間地震波數(shù)值模擬。2二維聲波頻率域正演原理2.1頻率-空間域聲波方程的離散時(shí)間域常密度二維聲波方程可寫為：(2-1)1d2u_d2ud2u(2-1)v2dt2dx2+dz2+S對(duì)上式作傅里葉變換得到：(2-2)co2d2Ud(2-2)Pu=kks其中U(X,Z,3)是u(x,Z,t)的傅立葉變換，S(x,Z,3)是震源函數(shù)s(x,z,t)的傅立葉變換，3為角頻率，數(shù)值上滿足3=2nf。可以看到，式(2-2)其實(shí)是一個(gè)關(guān)于波場值U(X,Z,3)的方程組，它可以寫成如下形式：AU=-S "一3)A稱為阻抗矩陣，它依賴于地下介質(zhì)的參數(shù)以及特定頻率。U是一個(gè)列向量，它是把特定頻率的波場值寫成了列向量的形式，容易知道，當(dāng)有限差分的網(wǎng)格規(guī)

模為nxXnz時(shí),U和S都是有nxXnz個(gè)元素的列向量，阻抗矩陣A是一個(gè)nxXnz階的方形矩陣。這時(shí)，我們就把求取波場值的問題轉(zhuǎn)化為了一個(gè)方程組的求解問題，不過這個(gè)方程組階數(shù)巨大，需要用特殊方法求解。頻率域求解的一個(gè)重要問題是如何將式（2-2）寫成式（2-3）的形式，這其中涉及到差分階數(shù)和差分方法的問題。目前使用較多的有五點(diǎn)差分法，九點(diǎn)差分法⑵和二十五點(diǎn)差分法【3】，對(duì)于本文的聲波方程，九點(diǎn)差分法具有足夠的精度，下面介紹九點(diǎn)差分法的實(shí)現(xiàn)。常規(guī)的十字型差分法只使用離散點(diǎn)上下及左右方向的點(diǎn)，九點(diǎn)差分法用到了離散點(diǎn)周圍所有點(diǎn)，包含左上，左下，右上，右下等方位（圖2-1）。它其實(shí)是將常規(guī)正交網(wǎng)格和45度旋轉(zhuǎn)網(wǎng)格結(jié)合的產(chǎn)物。圖2-1圖2-1九點(diǎn)差分格式示意圖混合網(wǎng)格(i,j)處的波場值U對(duì)空間的二階導(dǎo)數(shù)可以寫為：(2-4)*U,Ui,j+i-2Ui,j+Uq-iX 'Ui+ij+i-2Ui,j+Uj-ij-i(2-4)彩=a( )+b( 而波場值Uj,j則使用其周圍八個(gè)節(jié)點(diǎn)的波場值及其自身的值作加權(quán)平均:(2-5)㈤=cUi,j+d(Ui,j+i(2-5)+e（Ui+i,j+i+Ui+i,j_i+Ui—ij+i+Uj-Lj-i）其中，最優(yōu)化系數(shù)a=0.5461,b=l-a,c=0.6248,d=0.09381,e=（l-c-4d）/4。依據(jù)上面的思路，可以將頻率域聲波方程進(jìn)行離散，進(jìn)而寫成式（2-3）的形式。2.2PML邊界條件的引入正演時(shí)如果不對(duì)邊界進(jìn)行處理，則模擬得到的波場包含有大量反射波，無法

得到正確的結(jié)果。在模擬時(shí)間域的聲波傳播時(shí)，如果不添加吸收邊界條件，在波傳到邊界之前不產(chǎn)生邊界反射，只有在初至波到達(dá)邊界才產(chǎn)生反射，然后與有效波加。但頻率域聲波模擬不同，它的每一個(gè)單頻切片包含所有時(shí)間，所以如果不加吸收邊界條件，即使在零時(shí)刻也得不到正確的解⑷。PML邊界是目前使用最為廣泛的邊界條件，且PML邊界本身就是在頻率域進(jìn)行推導(dǎo)和證明的，因而頻率域PML邊界的引入較時(shí)間域更為簡單。頻率域聲波方程帶有PML邊界的形式可以寫為：(2-6)0）2 d2U du d2U du(2-6) 2U= Df+S

vz oxz ox ozz ozko)2ico+d⑵)ko)2ico+d⑵)心Y=ia）+d（x）7'一一舫Ui,j+(dk2+君)一舫Ui,j+(dk2+君)Uj+Lj+(dk2+M)Ui,j(dk2+黑)Ui，j+i+(dk?+急)Ui,j—i+(ek2(2-7)Ab+2Ax2（22AaAb2Cavk—彩—彩—彩(2-7)Ab+2Ax2+Iek2+Iek2+U+Iek2+Ui+Lj-1在所有網(wǎng)格點(diǎn)處進(jìn)行空間離散，并將其寫成式（2-3）所示的方程組形式。式（2-7）的系數(shù)構(gòu)成了阻抗矩陣A?？梢灾?，矩陣A是一個(gè)大型的稀疏矩陣,頻率域模擬的關(guān)鍵問題就是此稀疏矩陣方程組的求解。矩陣A的結(jié)構(gòu)如圖（2-2）所示。2.3大型稀疏矩陣方程組的求解在地震波場的正演模擬中，網(wǎng)格節(jié)點(diǎn)的數(shù)目往往達(dá)到數(shù)十萬，因而阻抗矩陣A將會(huì)是一個(gè)數(shù)十萬階的矩陣，常規(guī)的消元法是無法求解如此大規(guī)模的方程的。在數(shù)值計(jì)算領(lǐng)域，求解大型矩陣方程的方法往往有兩種，一種是迭代法，給定一個(gè)初始解，通過迭代更新使其不斷逼近真實(shí)解。另一種是直接法，其代表方法是LU分解法，將系數(shù)矩陣分解為一個(gè)上三角矩陣和一個(gè)下三角矩陣的乘積，從而便于求解。LU分解法使得頻率域波場模擬相對(duì)時(shí)間域有著巨大優(yōu)勢，在觀測系統(tǒng)和頻率不變時(shí)，LU分解的中間結(jié)果是可以重復(fù)利用的，當(dāng)完成了第一炮的求解后，其余炮可以調(diào)用己經(jīng)分解好的中間結(jié)果，從而在極短的時(shí)間內(nèi)完成求解。LU分解一般可以調(diào)用第三方的數(shù)學(xué)函數(shù)庫完成，這一類的函數(shù)庫有UMFPACK,SUPERLU等。但是這一類函數(shù)包都是在Linux/Unix系統(tǒng)下編寫的，Windows系統(tǒng)調(diào)用起來較為麻煩，可喜的是UMFPACK包可以在MATLAB中調(diào)用，而MATLAB本身就有強(qiáng)大的矩陣運(yùn)算功能以及完備的函數(shù)庫，因而本文選擇了MATLAB環(huán)境搭配UMFPACK進(jìn)行模擬。3頻率域聲波正演模擬實(shí)驗(yàn)3.1均勻介質(zhì)設(shè)i^一個(gè)2000m*2000m均勻介質(zhì)模型，空間網(wǎng)格步長Ax=Az=10mo縱波速度為2000m/So時(shí)間步長At=Imso震源為雷克子波，主頻為10Hz,道間距10m,震源位于模型中央。PML層厚度50層，理論反射系數(shù)IO*。分別對(duì)頻率為7.32Hz,12.2Hz,15.87Hz的分量進(jìn)行模擬，得到它們的實(shí)部和虛部如下圖：

虛部:05C01C0015002CG00 600 1000 1600 20000 500 1033 1E00 2000圖3-17.32Hz單頻波場圖3-212.21Hz單頻波場實(shí)部:圖虛部:05C01C0015002CG00 600 1000 1600 20000 500 1033 1E00 2000圖3-17.32Hz單頻波場圖3-212.21Hz單頻波場實(shí)部:圖3-315.87Hz單頻波場20000虛含8:o50010001600ECO1033 1500 2333可以看到，隨著頻率的增加，波峰和波谷交替越來越密集，波長越來越小。在0-36.63HZ內(nèi)以0.49Hz的頻率間隔進(jìn)行各頻率的波場模擬，對(duì)得到的單頻波場進(jìn)行傅立葉反變換進(jìn)而得到時(shí)間域的波場。選取其中第406ms的波前快照如圖(3.4)。利用時(shí)間域正演得到相同時(shí)刻的波前快照如圖(3-5),兩者的波場基本是一致的，但是也存在著可以接受的微小誤差。

02004006008001100012001400160018002000100015OTX/m圖3-5時(shí)間域406ms波前快照1CO31400160018D02000020040060080010DD100015OTX/m圖3-5時(shí)間域406ms波前快照1CO31400160018D02000020040060080010DD1200圖3-4頻率域經(jīng)過反變換的406ms波前快職圖3-6corner-edge模型Corner-edge模型的網(wǎng)格規(guī)模為200x200,網(wǎng)格間距10mo介質(zhì)速度為1500m/s和2500m/s。用雷克子波激發(fā)，震源位于模型頂層的中間。接收點(diǎn)位于地表，道間距10m,共200道。通過正演獲取12.2Hz的單頻波場如下圖:實(shí)部:1000廢部:150020000 500 1000 1933 2000實(shí)部:1000廢部:150020000 500 1000 1933 2000圖3-7corner-edge模型12.2Hz單頻波場在0-36.63HZ之間以間隔0.12Hz對(duì)各頻率進(jìn)行正演，對(duì)得到的各單頻波場進(jìn)行傅立葉反變換，得到820ms時(shí)的波前快照如圖(3-8)o將接收點(diǎn)接收到的各單頻波場做傅立葉反變換可以得到時(shí)間域的地震記錄，如圖(3-9)。從時(shí)間域的波前快照和地震記錄中可以明顯看到來自地下的反射波和尖滅點(diǎn)的繞射波，也證明了本文的頻率域波場模擬的結(jié)果是正確的。0 500 1000 1500 20000500100015000 500 1000 1500 20000500100015002000圖3-8820ms波前快照?qǐng)D3-9corner-edge模型地震記錄3結(jié)論與認(rèn)識(shí)二維頻率域正演相比時(shí)間域正演在效率上有著巨大優(yōu)勢，在保證模擬精度的前提下，頻率域方法可以利用LU分解的中間結(jié)果方便快速的計(jì)算多炮的單頻波場，并且它能夠非常有效的模擬各項(xiàng)異性介質(zhì)中波的傳播以及其衰減效應(yīng)，可以7方便地模擬粘彈性介質(zhì)中不同頻率波場能量的衰減。但是頻率域正演也有自身的缺點(diǎn)，最重要的缺點(diǎn)就是其內(nèi)存需求比時(shí)間域大得多，尤其是在3D模型中⑸。另一個(gè)缺陷就是，在高性能計(jì)算快速發(fā)展的今天，單炮情形下的基于LU分解的頻率域正演無法適應(yīng)MPI的并行算法以及GPU加速,但如果采用易于并行的迭代法進(jìn)行求解，則正演效率會(huì)大大降低。頻率域正演若要大規(guī)模應(yīng)用于3D模型的全波形反演，需要解決其內(nèi)存需求過大的問題以及難以并行的問題。參考文獻(xiàn)胡光輝，方伍寶等.全波形反演方法及應(yīng)用[M].北京：石油工業(yè)出版社，2014JoC,ShinC,SuhJH.Anoptimal9-point,finite-difference,frequency-space,2-Dscalarwaveextrapolator[J].Geophysics,1996,61(2):529-537ChangsooShin,HeejeungSohn.Afrequency-space2-Dscalarwaveextrapolatorusingextended25-pointfinite-differenceoperator_J]?Geophysics,1998,63(1):28