孟祥翠，王尚旭，唐跟陽，劉志遠(yuǎn)

1 中國(guó)石油大學(xué)(北京)CNPC物探重點(diǎn)實(shí)驗(yàn)室, 北京 102249 2 中國(guó)石油大學(xué)(北京)油氣資源與探測(cè)國(guó)家重點(diǎn)實(shí)驗(yàn)室, 北京 102249 3 中國(guó)石化石油勘探開發(fā)研究院, 北京 100083

石油地球物理

二維隨機(jī)介質(zhì)中的聲波輻射傳輸蒙特卡洛數(shù)值模擬及應(yīng)用

孟祥翠1,2*，王尚旭1,2，唐跟陽1,2，劉志遠(yuǎn)3

1 中國(guó)石油大學(xué)(北京)CNPC物探重點(diǎn)實(shí)驗(yàn)室, 北京 102249 2 中國(guó)石油大學(xué)(北京)油氣資源與探測(cè)國(guó)家重點(diǎn)實(shí)驗(yàn)室, 北京 102249 3 中國(guó)石化石油勘探開發(fā)研究院, 北京 100083

為了模擬富含小尺度非均質(zhì)體的勘探區(qū)內(nèi)地震波能量的變化，本文將天然地震中的輻射傳輸理論引入到地震勘探頻段內(nèi)，給出了基于Born散射系數(shù)的二維聲波輻射傳輸方程，并介紹了蒙特卡洛數(shù)值模擬的思路。為了檢驗(yàn)基于Born散射系數(shù)的聲波輻射傳輸蒙特卡洛方法的正確性和實(shí)用性，將輻射傳輸蒙特卡洛模擬的結(jié)果與聲波方程有限差分模擬的結(jié)果進(jìn)行對(duì)比。結(jié)果表明，輻射傳輸蒙特卡洛模擬的結(jié)果與聲波有限差分模擬的結(jié)果基本吻合。同時(shí)，相比于聲波方程有限差分?jǐn)?shù)值模擬，輻射傳輸蒙特卡洛模擬的計(jì)算成本遠(yuǎn)小于聲波有限差分?jǐn)?shù)值模擬。因此，對(duì)于復(fù)雜的非均質(zhì)儲(chǔ)層，輻射傳輸蒙特卡洛模擬是一種有效的地震散射波能量建模方法。而且，對(duì)于物理模型實(shí)驗(yàn)觀測(cè)到的地震散射波的數(shù)據(jù)，我們采用基于輻射傳輸蒙特卡洛模擬的網(wǎng)格掃描方法來估算非均勻介質(zhì)的統(tǒng)計(jì)參數(shù)－相關(guān)長(zhǎng)度和擾動(dòng)強(qiáng)度。計(jì)算結(jié)果與物理模型的真實(shí)值基本吻合，表明輻射傳輸蒙特卡洛模擬可以被用來估算非均質(zhì)儲(chǔ)層的統(tǒng)計(jì)參數(shù)。

Born散射系數(shù)；隨機(jī)介質(zhì)；輻射傳輸理論；蒙特卡洛模擬

0 引言

長(zhǎng)期以來，油氣地震勘探都以層狀介質(zhì)模型為基礎(chǔ)，但是現(xiàn)在勘探開發(fā)的重點(diǎn)逐漸轉(zhuǎn)向復(fù)雜的非均質(zhì)儲(chǔ)層，例如復(fù)雜砂體構(gòu)造區(qū)、碳酸鹽巖孔洞區(qū)以及火山巖區(qū)。地質(zhì)上常用基于數(shù)學(xué)統(tǒng)計(jì)理論的隨機(jī)模型來描述該類非均質(zhì)儲(chǔ)層[1-3]。地震波在該類介質(zhì)中傳播時(shí)，會(huì)發(fā)生多次散射，產(chǎn)生Coda波(也叫地震尾波)。由于多次散射形成的Coda波的相位是隨機(jī)的，在傳統(tǒng)的地震數(shù)據(jù)處理中，Coda波一般被視為噪聲。但是，研究表明散射波的包絡(luò)(即散射波的能量)與隨機(jī)介質(zhì)屬性密切相關(guān)[4-5]?？梢岳蒙⑸洳òj(luò)的變化，來反演隨機(jī)介質(zhì)的屬性參量[6-9]，例如非均質(zhì)體的特征長(zhǎng)度(表征非均質(zhì)體大小)和擾動(dòng)強(qiáng)度(表征非均質(zhì)體速度變化程度)。因此，準(zhǔn)確地描述隨機(jī)介質(zhì)中的散射波能量的變化至關(guān)重要。

在天然地震中，有很多理論來描述地震波能量在隨機(jī)介質(zhì)中的傳播，例如：波動(dòng)方程有限差分?jǐn)?shù)值模擬法[10]，擴(kuò)散理論[11-12]，Markov近似法[13]等。目前，最常用的是輻射傳輸理論，輻射傳輸理論忽略散射波的相位信息，考慮地震波能量在空間時(shí)間上的變化。輻射傳輸方程可以準(zhǔn)確的從波動(dòng)方程中推導(dǎo)出來。早期，出于簡(jiǎn)便的考慮，一般假設(shè)散射是各向同性的[14-15]，在這種情況下，可以直接得到解析解。但是，各向同性散射的假設(shè)，只能解釋Coda的產(chǎn)生，不能準(zhǔn)確的描述地震波脈沖展寬的現(xiàn)象；隨后，引入了各向異性散射系數(shù)[16]。各向異性散射下的輻射傳輸方程無法得到解析解，一般采用蒙特卡洛數(shù)值模擬法來求解[17-19]。蒙特卡洛模擬法的引入，使輻射傳輸理論更加簡(jiǎn)便實(shí)用。

本文將天然地震中的輻射傳輸理論引入地震勘探頻段。首先，介紹了基于Born散射系數(shù)的聲波輻射傳輸方程以及蒙特卡洛數(shù)值模擬的基本思想。接著，給出地震勘探中常用的高斯隨機(jī)模型和指數(shù)隨機(jī)模型各自對(duì)應(yīng)的4個(gè)例子，用輻射傳輸蒙特卡洛數(shù)值方法模擬地震波能量的變化。為了檢驗(yàn)基于Born散射系數(shù)的聲波輻射傳輸蒙特卡洛方法的正確性和實(shí)用性，將輻射傳輸蒙特卡洛模擬結(jié)果與聲波方程有限差模擬結(jié)果進(jìn)行對(duì)比。最后，采用基于輻射傳輸蒙特卡洛數(shù)值模擬的網(wǎng)格掃描方法，從物理模型實(shí)驗(yàn)采集到的地震散射數(shù)據(jù)來反演模型樣品的統(tǒng)計(jì)參數(shù)。

1 聲波輻射傳輸理論

1.1 隨機(jī)介質(zhì)

隨機(jī)介質(zhì)模型由大、小兩種尺度的非均勻性組成。大尺度非均勻性描述介質(zhì)的背景特性，用速度或密度的平均值表示；小尺度非均勻性則是加在背景模型上的隨機(jī)擾動(dòng)，用密度或速度的擾動(dòng)表示。當(dāng)隨機(jī)過程是平穩(wěn)的，隨機(jī)介質(zhì)可以通過速度擾動(dòng)的自相關(guān)函數(shù)以及對(duì)應(yīng)的功率譜密度函數(shù)來表征。常用的兩種隨機(jī)介質(zhì)模型是高斯型隨機(jī)模型和指數(shù)型隨機(jī)模型。二維各向同性高斯型和指數(shù)型隨機(jī)模型對(duì)應(yīng)的自相關(guān)函數(shù)和功率譜密度函數(shù)分別是：

1)二維高斯型隨機(jī)模型

自相關(guān)函數(shù)：

功率譜密度函數(shù)：

其中，a是相關(guān)長(zhǎng)度,ε是擾動(dòng)強(qiáng)度，m是波數(shù)。

2)二維指數(shù)型隨機(jī)模型

自相關(guān)函數(shù)：

功率譜密度函數(shù)：

為了減少獨(dú)立的介質(zhì)參數(shù)的個(gè)數(shù)，根據(jù)Birch法則，速度和密度擾動(dòng)表示為：

v0和ρ0是平均速度和平均密度，δv和δρ是速度和密度擾動(dòng)，υ=0.8[20]。

1.2 聲波輻射傳輸方程

輻射傳輸方程描述了地震波在隨機(jī)介質(zhì)中傳播時(shí)其能量的變化規(guī)律，可以由波動(dòng)方程嚴(yán)格推導(dǎo)出來[21-24]。二維聲波輻射傳輸方程可以表示為[20]：

對(duì)于方程(6)中的散射系數(shù)，可選用從Born近似中導(dǎo)出的單散射系數(shù)[19], 二維情況下：

指數(shù)型隨機(jī)模型對(duì)應(yīng)的散射系數(shù)：

高斯型隨機(jī)模型對(duì)應(yīng)的散射系數(shù)：

全散射系數(shù)可以通過對(duì)散射系數(shù)進(jìn)行積分獲得：

高斯型隨機(jī)介質(zhì)全散射系數(shù)為：

B0是修正的第一類貝塞爾函數(shù)，0表示零階。

指數(shù)型隨機(jī)介質(zhì)全散射系數(shù)為：

1.3 蒙特卡洛模擬

一般采用蒙特卡洛模擬的方法來求解方程 (6)。Wegler和Yoshimoto詳細(xì)介紹了蒙特卡洛模擬的算法[18-19]。蒙特卡洛方法求解輻射傳輸方程的基本思想是：通過大量的粒子運(yùn)動(dòng)來描述能量的傳輸過程。首先，在源的位置，向任意方向發(fā)射大量的粒子，每個(gè)粒子都代表一個(gè)單位能量，粒子沿一定方向以速度0v運(yùn)動(dòng)。然后，當(dāng)粒子運(yùn)動(dòng)到一定距離時(shí)會(huì)發(fā)生散射，該運(yùn)動(dòng)距離是由散射自由程(全散射系數(shù)的倒數(shù))來控制；當(dāng)發(fā)生散射時(shí)，粒子的運(yùn)動(dòng)方向會(huì)發(fā)生改變，散射角度是由歸一化的Born散射系數(shù)來確定。如果介質(zhì)是多層的，粒子運(yùn)動(dòng)到界面時(shí)，在界面是發(fā)生反射還是透射是根據(jù)平均反射和透射能量系數(shù)來確定的。在粒子的運(yùn)動(dòng)過程中，記錄每個(gè)粒子的運(yùn)動(dòng)軌跡。最后，統(tǒng)計(jì)在t時(shí)刻經(jīng)過接收點(diǎn)R處的粒子數(shù)，即可求得能量強(qiáng)度。

2 輻射傳輸蒙特卡洛數(shù)值模擬和波動(dòng)方程有限差分模擬結(jié)果對(duì)比

圖1 歸一化的散射系數(shù)g(θ)/ g(0)隨散射角θ的變化Fig. 1 The normalized scattering coeff i cient g(θ)/ g(0) varies with scattering angle θ

本文分別選用兩種隨機(jī)介質(zhì)模型：高斯型隨機(jī)模型和指數(shù)型隨機(jī)模型[25-29]。隨機(jī)介質(zhì)建模時(shí)，參數(shù)如表1所示，模型的背景速度v0=3 000 m/s，ε分別選ε=5%和ε=15%，非均質(zhì)體特征長(zhǎng)度分別選20 m和100 m。隨機(jī)模型的縱橫尺寸分別為8 000 m和8 000 m，中心放炮，空間采樣間隔dx=dy=1 m，時(shí)間采樣間隔為dt=0.1 ms，傳播時(shí)間為3 s，震源子波選用雷克子波，主頻是30 Hz。采用時(shí)間二階和空間四階交錯(cuò)網(wǎng)格差分格式來模擬隨機(jī)非均勻介質(zhì)中的地震波場(chǎng)，在模型上下左右都采用PML吸收邊界條件。記錄距離震源分別為1 500 m和3 000 m的圓上的地震波場(chǎng)，在每個(gè)圓上都設(shè)置72個(gè)等間距的接收點(diǎn)，如圖2所示。

表1 隨機(jī)模型參數(shù)Table 1 Parameters for random models

波動(dòng)方程有限差分模擬的是應(yīng)力場(chǎng)在空間時(shí)間上的變化，輻射傳輸蒙特卡洛模擬是能量在時(shí)間空間上的變化，為了將有限差分結(jié)果和輻射傳輸蒙特卡洛結(jié)果相比較，需要對(duì)有限差分模擬的地震波場(chǎng)求均方(MS)包絡(luò)。求取均方包絡(luò)的方法如下：對(duì)共炮檢距道集的每一道求平方，再對(duì)所有道求和取平均，最后以某一時(shí)間窗口進(jìn)行圓滑得到包絡(luò)。另外，有限差分模擬時(shí)采用的震源是有一定延展時(shí)間的雷克子波，輻射傳輸蒙特卡洛模擬相當(dāng)于震源是點(diǎn)脈沖震源，考慮到這點(diǎn)，還需要對(duì)輻射傳輸蒙特卡洛的結(jié)果進(jìn)行處理：對(duì)雷克子波求其平方，然后將輻射傳輸蒙特卡洛的結(jié)果與雷克子波取平方的結(jié)果求卷積。最后，用均勻介質(zhì)1 500 m處接收得到的包絡(luò)對(duì)處理后的輻射傳輸蒙特卡洛結(jié)果和波動(dòng)方程有限差分包絡(luò)進(jìn)行歸一化。

圖3是ε為5%的高斯隨機(jī)介質(zhì)蒙特卡洛輻射傳輸蒙特卡洛的結(jié)果和波動(dòng)方程有限差分結(jié)果對(duì)比圖，左圖是a=100 m，分別在1 500 m和3 000 m處接收到的結(jié)果，右圖是a=20 m，分別在1 500 m和3 000 m處接收的結(jié)果。不論是寬角散射還是前向散射，輻射傳輸蒙特卡洛模擬結(jié)果和波動(dòng)方程有限差分計(jì)算的結(jié)果都吻合的比較好；而且，輻射傳輸方法可以比較準(zhǔn)確的描述初至?xí)r間、前面直達(dá)波包絡(luò)的展寬現(xiàn)象以及后面尾波能量的變化。圖4是ε為15%的高斯隨機(jī)介質(zhì)輻射傳輸蒙特卡洛模擬結(jié)果和波動(dòng)方程有限差分結(jié)果對(duì)比圖?？梢钥吹?，輻射傳輸方法和聲波有限差分方法得到的散射波的包絡(luò)也基本吻合，只是前面部分有些許偏差。圖5和圖6是指數(shù)隨機(jī)介質(zhì)時(shí)，蒙特卡洛輻射傳輸結(jié)果和波動(dòng)方程有限差分結(jié)果對(duì)比圖，比較結(jié)果與高斯隨機(jī)介質(zhì)模型類似。

圖2 有限差分模型，縱橫長(zhǎng)度8 000×8 000 m。源在中心位置，分別在距離震源1 500 m和3 000 m的圓上接收，每個(gè)圓上包含72個(gè)等間距的接收點(diǎn)Fig. 2 A fi nite difference model with 8 000 m times 8 000 m grid. The source is located in the center and receivers are distributed along circles with radius 1 500 m and 3 000 m

圖7是4個(gè)指數(shù)隨機(jī)模型(見表1)，分別在1 500 m處接收，采用輻射傳輸蒙特卡洛模擬和波動(dòng)方程有限差分模擬用時(shí)的比較?？梢钥吹?，相比于寬角散射(ka=1.26)，前向散射輻射(ka=6.28)輻射傳輸蒙特卡洛模擬用時(shí)增加很多，這主要是由于當(dāng)ka比較大時(shí)，散射發(fā)生時(shí)散射角度主要集中在0o，我們采用取舍法來選擇散射角度則會(huì)比較費(fèi)時(shí)[19]。但是，無論是寬角散射還是前向散射情況，聲波輻射傳輸蒙特卡洛模擬計(jì)算所需時(shí)間遠(yuǎn)小于隨機(jī)介質(zhì)聲波方程有限差分?jǐn)?shù)值模擬的計(jì)算所需時(shí)間。

圖3 高斯隨機(jī)介質(zhì)模型，不同參數(shù)取值時(shí)兩種方法結(jié)果的對(duì)比。左圖是ε=5%，a=100 m 分別在R=1 500 m(上圖)和=3000 m (下圖)接收時(shí)波動(dòng)方程有限差分包絡(luò)(黑線)和輻射傳輸結(jié)果(紅線)的對(duì)比，右圖是ε=5%，a=20 m 分別在R=1 500 m(上圖)和R=3 000 m(下圖)接收時(shí)波動(dòng)方程有限差分包絡(luò)(黑線)和輻射傳輸結(jié)果(紅線)的對(duì)比Fig. 3 Numerical solutions for 2-D equation of radiative transfer (red) in comparison to numerical solutions of the acoustic wave equation (black) for Gaussian random medium with standard deviation of velocity fl uctuation ε=5%. Correlation distances of the random medium are a=100 m (left) and a=20 m (right). Source receiver distance are R=1 500 m (top) and R=3 000 m (bottom)

圖4 高斯隨機(jī)介質(zhì)模型，ε=15%，其它與圖3類似Fig. 4 Same as Fig. 3 for Gaussian random medium but with ε=15%

圖5 指數(shù)介質(zhì)模型，不同參數(shù)取值時(shí)兩種方法結(jié)果的對(duì)比。左圖是ε=5%，a=100 m分別在R=1 500 m和R=3 000 m接收時(shí)波動(dòng)方程有限差分包絡(luò)(黑色)和輻射傳輸結(jié)果(紅色)的對(duì)比，右圖是ε=5%，a=20 m分別在R=1 500 m和R=3 000 m接收時(shí)波動(dòng)方程有限差分包絡(luò)(黑色)和輻射傳輸結(jié)果(紅色)的對(duì)比Fig. 5 Numerical solutions for 2-D equation of radiative transfer (red) in comparison to numerical solution of the acoustic wave equation (black) for exponent random medium with standard deviation of velocity fl uctuation ε=5%. Correlation distances of the random medium are a=100 m (left) and a=20 m (right). Source receiver distance are R=1 500 m (top) and R=3 000 m (bottom).

圖6 指數(shù)隨機(jī)介質(zhì)模型，ε=15%，其它與圖5類似Fig. 6 Same as Fig. 5 for exponent random medium but with ε=15%

圖7 指數(shù)介質(zhì)隨機(jī)模型：模型1、模型2、模型3以及模型4，分別在1 500 m處接收時(shí)，輻射傳輸模擬(RT)和波動(dòng)方程有限差分模擬(FDM)用時(shí)比較,其中蒙特卡洛模擬時(shí)發(fā)射的震源粒子數(shù)為1 000 000Fig. 7 Times for 2-D equation of radiative transfer simulation (RT) in comparison to that for the acoustic wave equation fi nite difference method (FDM) for exponent random models (model 1, model 2, model 3 and model 4). The number of particles used in RT is 1 000 000

3 隨機(jī)介質(zhì)統(tǒng)計(jì)參數(shù)的估計(jì)

從方程(6)可知，當(dāng)已知地震子波和隨機(jī)介質(zhì)的平均速度時(shí)，穿過隨機(jī)介質(zhì)的地震散射波的能量只依賴于介質(zhì)的統(tǒng)計(jì)參數(shù)：相關(guān)長(zhǎng)度a和擾動(dòng)強(qiáng)度ε。因此，對(duì)采集到的地震散射數(shù)據(jù)，我們可以采用基于輻射傳輸方法的網(wǎng)格掃描方法，來估算隨機(jī)介質(zhì)的相關(guān)長(zhǎng)度a和擾動(dòng)強(qiáng)度ε。首先，選定要用的隨機(jī)模型的類型，例如高斯隨機(jī)模型，給定不同的相關(guān)長(zhǎng)度a和擾動(dòng)強(qiáng)度ε的值(不同的隨機(jī)模型)；然后，根據(jù)輻射傳輸蒙特卡洛方法來計(jì)算出不同模型時(shí)的地震波能量；最后，將輻射傳輸模擬的結(jié)果與采集到的地震波場(chǎng)的包絡(luò)進(jìn)行對(duì)比，找出與觀測(cè)到的地震波場(chǎng)包絡(luò)吻合最好時(shí)對(duì)應(yīng)的模型，該模型的相關(guān)長(zhǎng)度和擾動(dòng)強(qiáng)度就是隨機(jī)介質(zhì)的統(tǒng)計(jì)參數(shù)的值。為了評(píng)價(jià)不同模型輻射傳輸蒙特卡洛模擬的結(jié)果與觀測(cè)得到的地震波場(chǎng)包絡(luò)的吻合程度，我們采用下面的公式：

其中，it是時(shí)間，N是時(shí)間窗口的個(gè)數(shù)，jr表示第j個(gè)接收點(diǎn)，RT是接收點(diǎn)個(gè)數(shù)。Data是輻射傳輸模擬的結(jié)果，是觀測(cè)得到的地震波場(chǎng)的包絡(luò)。

表2 觀測(cè)系統(tǒng)參數(shù)Table 2 Parameters of seismic observation system

為了分析輻射傳輸蒙特卡洛模擬估算隨機(jī)介質(zhì)統(tǒng)計(jì)參數(shù)的準(zhǔn)確與否，我們用從物理模型實(shí)驗(yàn)采集得到的地震數(shù)據(jù)來估算介質(zhì)的統(tǒng)計(jì)參數(shù)，因?yàn)槲锢砟Ｐ偷幕緟?shù)是確定的。隨機(jī)模型樣品是選用玻璃珠和環(huán)氧樹脂來制作，環(huán)氧樹脂作為背景介質(zhì)，玻璃珠作為非均質(zhì)體，制作時(shí)玻璃珠隨機(jī)的分布于環(huán)氧樹脂中。隨機(jī)模型樣品與野外地質(zhì)構(gòu)造尺寸比為1：5 000，即模型1 mm代表野外地質(zhì)構(gòu)造的5 m；速度比為1：1。玻璃珠的粒徑為6 mm(野外為30 m)，速度為4 300 m/s。隨機(jī)模型樣品是長(zhǎng)方體，大小為50×30×5 cm3，即野外隨機(jī)介質(zhì)大小為2 500×1 500×250 m3，圖8b是樣品切面示意圖。模型樣品的平均速度為3 426 m/s,平均密度為 1.9 g/cm3。在水槽中采集透射波場(chǎng)，如圖8a所示，只沿樣品中線采集一條測(cè)線，觀測(cè)系統(tǒng)參數(shù)見表2。震源子波主頻大約是250 kHz(野外為50 Hz)。圖9是第41炮原始地震記錄，可以看到，不僅有直達(dá)的P波，P波后面還有明顯的Coda波。我們抽取共偏移距道集，根據(jù)前面介紹的方法，求取不同偏移距對(duì)應(yīng)的地震波的包絡(luò)，然后做總能量歸一化處理，得到歸一化的地震波包絡(luò)，如圖9b所示。

圖8 物理模型實(shí)驗(yàn)示意圖，(a)透射波水槽采集照片；(b) 實(shí)驗(yàn)采集觀測(cè)系統(tǒng)示意圖Fig. 8 An illustration of physical model experiments. (a) A photograph of the transmitted wave acquisition modelling system，(b) the transmitted wave experimental conf i guration

圖9 (a)第41炮原始地震記錄，(b)歸一化的地震包絡(luò)隨偏移距的變化情況Fig. 9 (a) The raw common shot gather (shot 41), (b) the normalized envelopes at different offsets

圖10 三層模型以及觀測(cè)系統(tǒng)示意圖Fig. 10 The three-layer model and acquisition conf i guration

在進(jìn)行輻射傳輸蒙特卡洛模擬時(shí)，我們選用與物理模型實(shí)驗(yàn)一致的三層模型(圖10)，上層和下層都是水，中間層是非均勻介質(zhì)。在地球物理中，常用來描述非均勻介質(zhì)的兩種隨機(jī)模型是高斯隨機(jī)模型和指數(shù)隨機(jī)模型；相比而言，高斯隨機(jī)模型更加平滑，其非均質(zhì)體大小比較單一[20]。由于物理模型樣品中的非均質(zhì)體是玻璃珠，而且粒徑大小基本一樣，所以此處隨機(jī)介質(zhì)模型我們選用高斯隨機(jī)模型。高斯隨機(jī)模型從(a, ε)=(5, 0.05)變到 (a, ε)=(40, 0.4)，相關(guān)長(zhǎng)度a以2.5 m的步長(zhǎng)從5 m變到40 m，對(duì)于每一個(gè)相關(guān)長(zhǎng)度a，擾動(dòng)強(qiáng)度ε都以0.025的步長(zhǎng)從0.05變到0.4。對(duì)于這些不同的高斯隨機(jī)模型，分別采用輻射傳輸蒙特卡洛模擬計(jì)算不同偏移距處的地震波能量的情況，炮點(diǎn)和接收點(diǎn)的相對(duì)位置于物理模型實(shí)驗(yàn)一致，如圖10所示。根據(jù)公式(13)求取不同模型下輻射傳輸蒙特卡洛模擬的結(jié)果與實(shí)驗(yàn)采集得到的地震波包絡(luò)的偏差χ，圖11a是χ的等值線圖，可以看到 當(dāng)(a, ε)=(27.5, 20%) 時(shí)，偏差最小。在物理模型實(shí)驗(yàn)中，隨機(jī)模型玻璃珠粒徑是30 m，我們通過輻射傳輸方法估算得到相關(guān)長(zhǎng)度是27.5 m，考慮到相關(guān)長(zhǎng)度是一個(gè)統(tǒng)計(jì)值，可以說該結(jié)果是正確可信的；模型樣品的平均速度是3 426 m/s, 玻璃珠的速度是 4 300 m/s，速度變化強(qiáng)度大約是 25.5%，估算得到的擾動(dòng)強(qiáng)度為20%，兩者也基本吻合。

圖11 (a) χ隨(a, ε)變化的等值線圖，(b)(a, ε)=(27.5 , 20%) 對(duì)應(yīng)模型的輻射傳輸模擬結(jié)果(藍(lán)色)和實(shí)驗(yàn)采集地震波場(chǎng)包絡(luò)(黑色)對(duì)比圖Fig. 11 (a) The contour plot of χ as a function of (a, ε), (b) Comparisons between the RT results (blue line) with the envelopes computed from the observed data (black line)

4 結(jié)論

本文針對(duì)常用的高斯隨機(jī)模型和指數(shù)隨機(jī)模型，分別采用基于Born散射系數(shù)的輻射傳輸蒙特卡洛方法模擬了寬角散射(ka=1.26)和前向散射(ka=6.28)兩種情況下，速度擾動(dòng)標(biāo)準(zhǔn)差較小(ε=5%)和較大(ε=15%)時(shí)地震波能量的傳輸情況。將輻射傳輸蒙特卡洛計(jì)算的結(jié)果與聲波方程有限差模擬的結(jié)果進(jìn)行了對(duì)比，兩者吻合的比較好，表明輻射傳輸蒙特卡洛數(shù)值模擬可以比較準(zhǔn)確的描述隨機(jī)介質(zhì)中散射波的能量變化；另外，相比于聲波方程有限差分?jǐn)?shù)值模擬，聲波輻射傳輸蒙特卡洛模擬的計(jì)算成本遠(yuǎn)小于聲波有限差分?jǐn)?shù)值模擬的計(jì)算成本。因此，對(duì)于復(fù)雜的非均質(zhì)儲(chǔ)層，輻射傳輸蒙特卡洛模擬是一種有效的地震散射波能量建模方法。對(duì)物理模型實(shí)驗(yàn)觀測(cè)到的地震散射數(shù)據(jù)，采用基于輻射傳輸數(shù)模擬的網(wǎng)格掃描方法來估算隨機(jī)介質(zhì)的統(tǒng)計(jì)參數(shù)-相關(guān)長(zhǎng)度和擾動(dòng)強(qiáng)度。計(jì)算結(jié)果與物理模型的真實(shí)值吻合較好，說明輻射傳輸蒙特卡洛模擬方法可以用來估算隨機(jī)介質(zhì)的統(tǒng)計(jì)參數(shù)。

[1] 奚先, 姚姚, 顧漢明. 隨機(jī)溶洞介質(zhì)模型及其波場(chǎng)模擬[J]. 地球物理學(xué)進(jìn)展, 2005, 20(2), 365-369. [XI X, YAO Y, GU H M. Random cavity medium model and the wave fi eld simulation[J]. Progress In Geophysics, 2005, 20(2), 365-369.]

[2] ALMEIDA J A. Stochastic simulation methods for characterization of lithoclasses in carbonate eservoirs[J]. Earth Science Reviews, 2010, 101： 250-270.

[3] KUZNETSOVA A, ALMEIDA J A, LEGOINHA P. Stochastic simulation of the morphology of fluvial sand channel reservoirs[C]. Mathematics of Planet Earth. Berlin： Springer-Verlag, 2014： 689-693.

[4] AKI K. Analysis of the seismic coda of local earthquakes as scattered waves[J]. Journal of Geophysical Research, 1969, 74(2)： 615-631.

[5] AKI K, CHOUET B. Origin of coda waves： source, attenuation, and scattering effects[J]. Journal of Geophysical Research, 1975, 80：3322-3342.

[6] MICHAEL A H, SHEARER P M, EARLE P S. Seismic evidence for small-scale heterogeneity throughout the Earth’s mantle[J]. Nature, 1997, 387： 145-150.

[7] MARGERIN L G, NOLET G. Multiple scattering of high-frequency seismic waves in the deep Earth： PKP precursor analysis and inversion for mantle granularity[J]. Journal of Geophysical Research, 2003, 108(B11)： 2514-2529.

[8] VIDALE J E, EARLE P S. Fine-scale heterogeneity in the Earth’s inner core[J]. Nature, 2000, 404： 273-275.

[9] 王勤彩, 陳章立, ASANO Y等. 利用尾波包絡(luò)線反演方法研究伽師強(qiáng)震群區(qū)地殼的非均勻結(jié)構(gòu)[J]. 地球物理學(xué)報(bào), 2009, 52(1)：90-98. [WANG C Q, CHEN L Z, ASANO Y, et al. Imaging crustal heterogeneity in Jiashi strong earthquake swarm region by coda envelope inversion analysis[J]. Chinese Journal of Geophysics, 2009, 52(1)： 90-98.]

[10] FEHLER M, SATO H, HUANG L J. Envelope broadening of outgoing waves in 2-D random media： A comparison between the Markov approximation and numerical simulation[J]. Bulletin of the Seismological Society of America, 2000, 90(4)： 914-928.

[11] MARGERIN L, CAMPILLO M, TIGGELEN B A. Radiative transfer and diffusion of waves in a layered medium： New insight into coda Q[J]. Geophysical Journal International, 1998, 134(2)： 596-612.

[12] WEGLER U. Diffusion of seismic waves in a thick layer： Theory and application to Vesuvius volcano[J]. Journal of Geophysical Research, 2004, 109(B7)： B07303.

[13] SATO H. Broadening of seismogram envelopes in the randomly inhomogeneous lithosphere based on the parabolic approximation：Southeastern Honshu, Japan[J]. Journal of Geophysical Research, 1989, 94(B12)： 17735-17747.

[14] HOSHIBA M. Simulation of multiple scattered coda wave excitation based on the energy conservation law[J]. Physics of the Earth and Planetary Interiors, 1991, 67(1-2)： 123-136.

[15] ZENG Y, SU F, AKI K. Scattering wave energy propagation in a random isotropic scattering medium-I theory[J]. Journal of Geophysical Research, 1991, 96(B1)： 607-619.

[16] GUSEV A A, ABUBAKIROV I R. Simulated envelopes of non-isotropically scattered body waves as compared to observed ones：Another manifestation of fractal heterogeneity[J]. Geophysical Journal International, 1996, 127(1)： 49-60.

[17] HOSHIBA M. Estimation of non-isotropic scattering in western Japan using coda wave envelopes’ application of a multiple non-isotropic scattering model[J]. Journal of Geophysical Research, 1995, 100(B1)： 645-657.

[18] YOSHIMOTO K. Monte Carlo simulation of seismogram envelopes in scattering media[J]. Journal of Geophysical Research, 2000, 105(B3)： 6453-6161.

[19] WEGLER U, KORN M, PRZYBILLA J. Modeling full seismogram envelopes using radiative transfer theory with born scattering coeff i cients[J]. Pure and Applied Geophysics, 2006, 163(2)： 503-531.

[20] SATO H, FEHLER M, Maeda T. Seismic wave propagation and scattering in the heterogeneous earth[M]. Berlin： Springer-Verlag, 2012. [21] BURRIDGE R, PAPANICOLAOU G C. Transport equations for Stokes’ parameters from Maxwell’s equations in a random medium[J]. Journal of Mathematical Physics, 1975, 16： 2074-2085.

[22] RYTOV S M, KRAVTSOV Y A, TATARSKII V I. Principles of statistical radiophysics 4-Wave propagation through random media[M]. Berlin： Springer-Verlag, 1987.

[23] WEAVER R L. Diffusivity of ultrasound in polycrystals[J]. Journal of the Mechanics and Physics of Solids, 1990, 38(1)： 55-86.

[24] RYZHIK L V, PAPANICOLAOU G C, KELLER J B. Transport equation for elastic and other waves in random medium[J]. Wave Motion, 1996, 24(4)： 327-362.

[25] 奚先, 姚姚. 二維隨機(jī)介質(zhì)及波動(dòng)方程正演模擬[J]. 石油地球物理勘探, 2001, 36(5)： 546-552. [XI X, YAO Y. 2-D random media and wave equation forward modeling[J]. Oil Geophysical Prospecting, 2001, 36(5)： 546-552.]

[26] 奚先, 姚姚. 二維橫各向同性彈性隨機(jī)介質(zhì)中的波場(chǎng)特征[J]. 地球物理學(xué)進(jìn), 2004, 19(4)： 924-932. [XI X, YAO Y. The wave fi eld characteristics of 2-D transversely isotropic elastic random medium[J]. Progress in Geophysics, 2004, 19(4)： 924-932.]

[27] 劉永霞, 徐春明, 寧俊瑞. 不同模式自組織介質(zhì)中聲波傳播特性的比較研究[J]. 地球物理學(xué)報(bào), 2007, 50(3)： 830-836. [LIU Y X, XU C M, NING J R. Comparison of acoustic propagation in several different types self-organized media[J]. Chinese Journal of Geophysics, 2007, 50(3)： 830-836.]

[28] 郭乃川, 王尚旭, 董春暉, 等. 地震勘探中小尺度非均勻性的描述及長(zhǎng)波長(zhǎng)理論[J]. 地球物理學(xué)報(bào), 2012, 55(7)： 2385-2401. [GUO N C, WANG S X, DONG C H, et al. Description of small scale inhomogeneities in seismic prospecting and long-wavelength theory[J]. Chinese Journal of Geophysics, 2012, 55(7)： 2385-2401.]

[29] 范小平, 楊從杰, 李清河. 介質(zhì)非均勻性參數(shù)對(duì)散射包絡(luò)的影響[J]. 地球物理學(xué)進(jìn)展, 2013, 28(2)： 0687-0694. [FAN X P, YANG C J, LI Q H. The effect of inhomogeneity parameters on scattering wave envelopes[J]. Progress in Geophysics, 2013, 28(2)： 0687-0694.]

Monte-Carlo radiative transfer simulation of acoustic waves in two-dimensional random media and it’s application

MENG Xiangcui1,2, WANG Shangxu1,2, TANG Genyang1,2, LIU Zhiyuan3
1 CNPC Key Lab of Geophysical exploration, China University of Petroleum-Beijing, Beijing 102249, China 2 State Key Laboratory of Petroleum Resource and Prospecting, China University of Petroleum-Beijing, Beijing 102249, China 3 Sinopec Petroleum Exploration and Production Research Institute, Beijing 100083, China

In order to describe seismic energy transport in exploration fi elds full of small-scale heterogeneities, the acoustic radiative transfer theory in seismology is introduced into frequency zones of seismic exploration. In this paper, the acoustic radiative transfer theory based on Born scattering coeff i cients in 2-D random media are given, and the Monte-Carlo simulation is brief l y presented to solve the radiative transfer equation. In order to verify the radiative transfer method, Monte-Carlo solutions of the radiative transfer equation are compared with finite difference solutions of the acoustic wave equation, and we find a good coincidence of the two theories. Meanwhile, the Monte-Carlo radiative transfer simulation costs less time than the fi nite difference simulation of acoustic wave equations. Thus, it is a good option to apply Monte-Carlo radiative transfer simulation to model the energy transport of seismic scattering waves in complex heterogeneous reservoirs. Furthermore, we use a grid search procedure based on Monte-Carlo radiative transfer simulation to infer the stochastic parameters - the correlation length and the fl uctuation strength from the observed seismic scattering data acquired in the physical experiment. The estimated results agree with the true values of the physical sample, which indicate that the Monte-Carlo radiative transfer simulation can be used to invert the stochastic parameters of heterogeneous reservoirs.

Born scattering coeff i cient; random medium; radiative transfer theory; Monte-Carlo simulation

10.3969/j.issn.2096-1693.2017.02.019

(編輯 付娟娟)

孟祥翠, 王尚旭, 唐跟陽, 劉志遠(yuǎn). 二維隨機(jī)介質(zhì)中的聲波輻射傳輸蒙特卡洛數(shù)值模擬及應(yīng)用. 石油科學(xué)通報(bào), 2017, 02： 199-209

MENG Xiangcui, WANG Shangxu, TANG Genyang, LIU Zhiyuan. Monte-Carlo radiative transfer simulation of acoustic waves in two-dimensional random media and it’s application. Petroleum Science Bulletin, 2017, 02： 199-209.doi：10.3969/ j.issn.2096-1693.2017.02.019

*通信作者, mengxiangcui35105@163.com

2016-11-23

國(guó)家重點(diǎn)基礎(chǔ)研究發(fā)展計(jì)劃項(xiàng)目(2013CB228600)以及國(guó)家自然科學(xué)基金(41304042)共同資助