基于黏彈性介質(zhì)波動(dòng)理論的頁(yè)巖超聲波數(shù)值模擬

陳 喬 徐烽淋 程 亮 劉 洪 簡(jiǎn) 旭 朱洪林 陳吉龍

1.中國(guó)科學(xué)院重慶綠色智能技術(shù)研究院 2.重慶市涪陵頁(yè)巖氣環(huán)保研發(fā)與技術(shù)服務(wù)中心3.中國(guó)石油川慶鉆探工程公司地質(zhì)勘探開(kāi)發(fā)研究院 4.重慶地質(zhì)礦產(chǎn)研究院5.“油氣藏地質(zhì)及開(kāi)發(fā)工程”國(guó)家重點(diǎn)實(shí)驗(yàn)室·西南石油大學(xué) 6. 重慶工商大學(xué)

0 引言

頁(yè)巖巖石的非均質(zhì)性及層理布局的隨機(jī)性，給頁(yè)巖中聲波傳播特性研究增加了較大的難度[1-7]。因此，學(xué)者們針對(duì)超聲波在層狀巖石中的傳播特征展開(kāi)了研究。在物理實(shí)驗(yàn)方面，Johns等[8]用X射線衍射和掃描電鏡技術(shù)得到層理是引起頁(yè)巖彈性波復(fù)雜和彈性波分裂特征的重要原因；Sondergeld等[9]利用超聲波透射實(shí)驗(yàn)的手段，分析了圍壓對(duì)不同層理角度頁(yè)巖的縱、橫波波速影響規(guī)律；Horne等[10]在建造斜井階段利用偶極子聲波測(cè)井?dāng)?shù)據(jù)來(lái)判別薄互層頁(yè)巖的彈性各向異性；Domnesteanu等[11]開(kāi)展了頁(yè)巖的聲波速度和衰減分析，結(jié)果表明頁(yè)巖各向異性的大小取決于層理方向與聲波傳播方向的夾角以及頁(yè)巖滲透率的方向；Szewczyk等[12]在一系列地震超聲室內(nèi)實(shí)驗(yàn)的基礎(chǔ)上對(duì)曼柯斯頁(yè)巖和皮埃爾頁(yè)巖的彈性應(yīng)力敏感度進(jìn)行了研究，認(rèn)為頁(yè)巖分別在地震波和超聲波頻率下具有不同的應(yīng)力敏感度。在國(guó)內(nèi)，鄧?yán)^新等[13]利用多頻率超聲波對(duì)層理發(fā)育的頁(yè)巖各向異性進(jìn)行了研究，給出了不同狀態(tài)下，樣品不同方向上聲波速度、各向異性參數(shù)隨壓力的變化規(guī)律；王倩等[14]通過(guò)頁(yè)巖的超聲波實(shí)驗(yàn)測(cè)量探索了不同層理傾角頁(yè)巖動(dòng)、靜態(tài)彈性模量和泊松比的關(guān)系。陳喬等[15]、熊健等[16]采用超聲波透射法系統(tǒng)地分析了下志留統(tǒng)龍馬溪組頁(yè)巖不同層理角度的巖石縱橫波速度、衰減以及頻率響應(yīng)，獲得了層理角度與聲波速度、衰減系數(shù)統(tǒng)計(jì)函數(shù)關(guān)系。

在數(shù)值模擬方面，Saenger等[17]利用交錯(cuò)有限差分法來(lái)研究孔隙度和孔洞密度與聲波的響應(yīng)關(guān)系，幾年后，又利用旋轉(zhuǎn)的交錯(cuò)差分網(wǎng)格對(duì)含黏彈性液體孔洞、無(wú)黏彈性天然氣條件下的波速進(jìn)行了模擬，并把模擬結(jié)果與真實(shí)測(cè)量結(jié)果進(jìn)行了對(duì)比分析[18]；Zhu、Baechle、Weger等[19-21]利用不同分辨率電鏡X射線分辨率為3.1 μs、EPMA分辨率為1 μs、SEM分辨率為0.1 μs和鑄體薄片進(jìn)行數(shù)字圖像分析，提取出孔隙大小、球形度、孔隙縱橫波、孔隙網(wǎng)絡(luò)的復(fù)雜程度等參數(shù)定量描述孔隙空間，在此基礎(chǔ)上進(jìn)行巖石聲波數(shù)值模擬；Elhusseiny等[22]和Xu等[23]利用人工模擬方法探討了巖石孔隙結(jié)構(gòu)的超聲波傳播特性；Tseng等[24]利用超聲物理模型測(cè)試了P—SV轉(zhuǎn)換波在TI介質(zhì)中的轉(zhuǎn)換點(diǎn)及其時(shí)程的數(shù)值計(jì)算精確度，通過(guò)費(fèi)曼原理和各向同性非雙曲時(shí)差方程兩種方法計(jì)算發(fā)現(xiàn)，P—SV轉(zhuǎn)換波在TI介質(zhì)各向同性平面?zhèn)鞑ブ芯哂邢嗤霓D(zhuǎn)換點(diǎn)和時(shí)程；國(guó)內(nèi)Yao等[25]利用“疊加型”多尺度隨機(jī)介質(zhì)模型、“區(qū)域多尺度隨機(jī)介質(zhì)模型”來(lái)模擬隨機(jī)溶洞介質(zhì)；魏建新等[26]以棍形、片形、球形、柱形、裂隙孔隙形狀形態(tài)的物理模型設(shè)計(jì)了填充物相同、橫向、縱向尺度不同的地震模型，在此基礎(chǔ)上，通過(guò)數(shù)值模擬計(jì)算，分析了由孔隙結(jié)構(gòu)變化引起的地震響應(yīng)特征；劉向君等[27-28]基于二維聲波數(shù)值模擬計(jì)算，開(kāi)展裂縫模型聲波衰減系數(shù)研究，同時(shí)還開(kāi)展了孔洞型儲(chǔ)層孔隙度預(yù)測(cè)模型研究。

綜上所述，國(guó)內(nèi)外學(xué)者主要通過(guò)超聲波透射實(shí)驗(yàn)獲得頁(yè)巖聲波傳播規(guī)律的定性認(rèn)識(shí)。同時(shí)，考慮到在物理模型實(shí)驗(yàn)過(guò)程中頁(yè)巖巖樣制備較為困難，國(guó)內(nèi)學(xué)者對(duì)作為油氣儲(chǔ)集層蓋層的頁(yè)巖地層聲波速度各向異性的巖石物理數(shù)據(jù)研究較少。僅使用實(shí)體物理模型研究復(fù)雜結(jié)構(gòu)的地層聲學(xué)特性是片面、非客觀的，因此，認(rèn)識(shí)和分析聲波在復(fù)雜地層介質(zhì)中的傳播規(guī)律及其響應(yīng)特征須結(jié)合數(shù)值模擬技術(shù)。然而，目前的數(shù)值模擬研究主要針對(duì)常規(guī)地層結(jié)構(gòu)的地震識(shí)別，以頁(yè)巖地層聲波測(cè)井為目的的聲學(xué)響應(yīng)鮮有報(bào)道，相關(guān)研究將頁(yè)巖當(dāng)作彈性介質(zhì)處理，未考慮頁(yè)巖的黏彈性特征。因此，筆者基于四川盆地渝東南下志留統(tǒng)龍馬溪組頁(yè)巖巖心，建立了不同角度的層理結(jié)構(gòu)模型，采用黏彈性介質(zhì)波動(dòng)理論和高階交錯(cuò)網(wǎng)格有限差分方法，開(kāi)展頁(yè)巖超聲波透射實(shí)驗(yàn)數(shù)值模擬研究，計(jì)算不同層理角度頁(yè)巖的聲學(xué)參數(shù)特征，總結(jié)頁(yè)巖聲波參數(shù)隨層理角度的變化規(guī)律。

1 數(shù)值模擬理論基礎(chǔ)

數(shù)值模擬中常用的波動(dòng)方程可分為二階彈性波動(dòng)方程和一階位移—應(yīng)力彈性波方程組2種。由于一階位移—應(yīng)力彈性波動(dòng)方程不需對(duì)彈性常數(shù)進(jìn)行空間微分，利用交錯(cuò)網(wǎng)格半程計(jì)算可以避免求解時(shí)間的高階微分，同時(shí)，網(wǎng)格引入的頻散較小，筆者采用一階位移—應(yīng)力彈性波動(dòng)方程來(lái)開(kāi)展數(shù)值模擬研究。

1.1 聲波波動(dòng)理論

聲波波動(dòng)方程所包括的本構(gòu)方程、運(yùn)動(dòng)微分方程和幾何方程揭示了物體內(nèi)部應(yīng)變、應(yīng)力和位移之間的相互關(guān)系，其中本構(gòu)方程主要反映的是物體的各向異性、黏滯性和彈性等性質(zhì)。Kelvin—Voigt黏彈性介質(zhì)模型能較好地描述聲波在介質(zhì)中的衰減[29]?；趹?yīng)力—應(yīng)變關(guān)系、本構(gòu)方程以及Kelvin—Voigt黏彈性模型，得到了二維黏彈性介質(zhì)的波動(dòng)方程，為了便于對(duì)其做有限差分，將二維黏彈性波動(dòng)方程轉(zhuǎn)化為一階速度—應(yīng)力方程，即

式中ρ表示密度，g/cm3；vx、vz分別表示速度分量，m/s；σxx、σzz、τxz、τzx分 別 表 示 應(yīng) 力 分 量，MPa；λ、μ分別表示拉梅系數(shù)；λ'、μ'分別表示黏滯系數(shù)，Pa·s；Qp、Qs分別表示縱波、橫波的品質(zhì)因子；vp、vs分別表示縱波和橫波的速度，m/s；ω表示圓頻率。

1.2 數(shù)值離散

根據(jù)泰勒展開(kāi)式的基本原理，速度分量vz在2階時(shí)間精度和2n階空間精度下的離散量為：

式中W表示速度分量（vz）的離散量；S、T分別表示應(yīng)力分量σzz、τxz的離散量；Cn表示差分階數(shù)；n表示空間域的階數(shù)；Δt表示時(shí)間步長(zhǎng)；Δx、Δz分別表示空間步長(zhǎng)；（i，j）表示網(wǎng)格節(jié)點(diǎn)標(biāo)號(hào)；k表示時(shí)間節(jié)點(diǎn)標(biāo)號(hào)。

1.3 穩(wěn)定性分析

為了避免高階有限差分的不穩(wěn)定性，需要選擇合理的數(shù)值模擬參數(shù)。對(duì)于二維黏彈性介質(zhì)，時(shí)間2階和空間4階差分的穩(wěn)定性條件為[30]：

1.4 振源函數(shù)

為了更好地模擬頁(yè)巖的超聲波特性，振源函數(shù)是模擬頁(yè)巖聲波響應(yīng)的關(guān)鍵，數(shù)值模擬中選取的聲源為超聲波探頭原始信號(hào)為：

式中R表示震源函數(shù)；x0、z0分別表示坐標(biāo)值。

1.5 邊界條件

在數(shù)值模擬的過(guò)程中，引入完全匹配層吸收模型，在給定的模型外面構(gòu)造有限厚度的吸收層，當(dāng)聲波傳入到吸收層時(shí)能量會(huì)隨傳播距離呈指數(shù)衰減，從而消除邊界反射對(duì)數(shù)值模擬的影響[31]。即

式中Gb表示邊界矩形區(qū)域的表達(dá)式；Gc表示邊界是個(gè)角區(qū)的函數(shù)表達(dá)式；α表示系數(shù)；N表示邊界網(wǎng)格的層數(shù)；（i0，j0）表示4個(gè)角的內(nèi)側(cè)頂點(diǎn)標(biāo)號(hào)。

2 數(shù)值模擬驗(yàn)證

2.1 模擬波形與實(shí)驗(yàn)波形驗(yàn)證

為了與物理實(shí)驗(yàn)相統(tǒng)一，針對(duì)室內(nèi)標(biāo)準(zhǔn)巖心的超聲波透射實(shí)驗(yàn)，筆者采用250 kHz的“探頭頻率”表示超聲波頻率大小。為了驗(yàn)證二維時(shí)域有限差分程序的正確性，計(jì)算了均勻地層條件下的聲波波形，并將有限差分結(jié)果與實(shí)軸積分結(jié)果進(jìn)行對(duì)比，對(duì)比結(jié)果如圖1所示。圖1中藍(lán)線為有限差分計(jì)算的結(jié)果，紅線表示實(shí)軸積分法的計(jì)算結(jié)果。有限差分法和實(shí)軸積分法的模擬計(jì)算結(jié)果基本一致，驗(yàn)證了有限差分?jǐn)?shù)值模擬方法的正確性。

圖1 二維有限差分與實(shí)軸積分波形對(duì)比圖

2.2 真實(shí)巖心數(shù)值計(jì)算及驗(yàn)證

考慮到層理角度為頁(yè)巖最顯著的結(jié)構(gòu)特征，因此，基于龍馬溪組真實(shí)頁(yè)巖巖心，建立不同層理角度的結(jié)構(gòu)模型開(kāi)展數(shù)值模擬（圖2）。從圖2可以看見(jiàn)，物理實(shí)驗(yàn)所用的真實(shí)巖心與數(shù)值模擬所用的數(shù)字巖心層理特征較為吻合，能真實(shí)反應(yīng)頁(yè)巖的層理特征。

巖石超聲波衰減系數(shù)既可用于分析巖石的層理構(gòu)造及分布，還可用于辨別巖石在地下所處的地質(zhì)環(huán)境，同時(shí)，針對(duì)巖石物理狀態(tài)的變化，衰減特性相較于波速特性更為敏感。衰減系數(shù)是表征衰減特性的關(guān)鍵參數(shù)。因此，筆者利用數(shù)值模擬不同頻率下巖心的衰減系數(shù)（圖3），隨著測(cè)試頻率的增加，超聲波衰減系數(shù)呈增大趨勢(shì)，數(shù)值模擬計(jì)算結(jié)果與陳喬物理實(shí)驗(yàn)[15]所得到總體規(guī)律較為吻合。

圖2 物理實(shí)驗(yàn)與數(shù)值計(jì)算模型對(duì)比圖

圖3 超聲波衰減特性結(jié)果對(duì)比圖

針對(duì)巖心重復(fù)性較差的緣故，上述數(shù)值模擬計(jì)算是基于理想的巖心層理角度模型。接下來(lái)，為了進(jìn)一步驗(yàn)證數(shù)值模擬結(jié)果的有效性，采用數(shù)值圖像處理方法提取出渝東南龍馬溪組真實(shí)巖心的層理角度，開(kāi)展不同層理角度的衰減系數(shù)數(shù)值計(jì)算。

圖4 實(shí)物巖樣照片

圖5 實(shí)物巖樣模型速度云圖

圖4是巖樣的實(shí)物照片（記為A1～A4），從照片可以看出，頁(yè)巖層理發(fā)育，層理角度在0°、45°、60°和90°變化。以圖5為研究對(duì)象，利用邊緣提取和Hough直線檢測(cè)方法將巖心的基質(zhì)骨架和層理區(qū)域進(jìn)行分類識(shí)別， 基于此，建立對(duì)應(yīng)的數(shù)字模型。將層理結(jié)構(gòu)識(shí)別后的數(shù)字圖像利用方形單元映射技術(shù)可以轉(zhuǎn)化成所需的有限差分網(wǎng)格數(shù)據(jù)，因?yàn)閿?shù)字圖像是由一個(gè)個(gè)像素點(diǎn)排列組合而成，每個(gè)像素點(diǎn)都可作為有限差分計(jì)算中的四邊形單元。其中，每個(gè)四邊形單元所對(duì)應(yīng)的節(jié)點(diǎn)坐標(biāo)又可通過(guò)圖像實(shí)際尺寸大小與像素之間的比例值，轉(zhuǎn)換成所對(duì)應(yīng)矢量空間的物理坐標(biāo)，由此，數(shù)值圖像便成功轉(zhuǎn)換為正方形的差分網(wǎng)格。最后，依據(jù)每個(gè)像素點(diǎn)所對(duì)應(yīng)的灰度值，為每個(gè)單元賦予相應(yīng)的特征參數(shù)，從而實(shí)現(xiàn)將數(shù)字圖像的層理結(jié)構(gòu)導(dǎo)入數(shù)字巖心模型（圖5）。從速度云圖來(lái)看，B1～B4的層理傾角分別分布在0°、45°、60°和90°，與實(shí)驗(yàn)巖樣吻合較好。最后，利用Matlab編寫代碼進(jìn)行數(shù)值模擬計(jì)算，從圖6所示的結(jié)果來(lái)看，衰減系數(shù)與層理角度都呈正相關(guān)規(guī)律變化，相同層理角度的衰減系數(shù)要小于模擬值，這是由于實(shí)物巖樣（A1～A4）所對(duì)應(yīng)的層理厚度略低于B1～B4，對(duì)應(yīng)的衰減系數(shù)的值及上升的幅度也低于理論模型，但是從兩者之間的平均值和標(biāo)準(zhǔn)差發(fā)現(xiàn)，數(shù)據(jù)偏離平均值程度較小，數(shù)據(jù)離散程度低，兩者的衰減系數(shù)變化趨勢(shì)基本一致。同時(shí)，本次實(shí)驗(yàn)結(jié)果與陳喬等[15]、范翔宇等[32]的室內(nèi)實(shí)驗(yàn)結(jié)果相比較，衰減系數(shù)隨著層理角度的編號(hào)趨勢(shì)也基本一致，進(jìn)一步說(shuō)明本數(shù)值模擬方法的有效性。

綜上所述，針對(duì)常規(guī)室內(nèi)超聲波透射法實(shí)驗(yàn)成本高、誤差較大的局限性，筆者提出了一種模擬頁(yè)巖超聲波透射實(shí)驗(yàn)的數(shù)值計(jì)算方法，該方法與實(shí)際物理模型實(shí)驗(yàn)結(jié)果吻合度較高，具有進(jìn)一步通過(guò)大規(guī)模數(shù)值模擬計(jì)算，開(kāi)展頁(yè)巖不同層理結(jié)構(gòu)聲學(xué)特征響應(yīng)研究的價(jià)值。

圖6 室內(nèi)巖樣實(shí)驗(yàn)與數(shù)值模擬計(jì)算所得衰減系數(shù)結(jié)果對(duì)比圖

3 數(shù)值計(jì)算分析

利用數(shù)值計(jì)算建立頁(yè)巖不同層理角度的數(shù)字模型如圖7所示。層理厚度均為0.2 mm，層理延伸方向與聲波激發(fā)方向夾角分別為0°、10°、20°、30°、40°、50°、60°、70°、80°、90°（記為 C1～ C10）。

以不同層理角度的頁(yè)巖數(shù)值模型為基礎(chǔ)開(kāi)展數(shù)值模擬計(jì)算，得到對(duì)應(yīng)的波形，從中提取巖石聲波的衰減系數(shù)和速度參數(shù)，分析層理角度對(duì)聲學(xué)特性的影響。

圖7 數(shù)字模型速度云圖

3.1 層理角度對(duì)聲波波速的影響

由圖8可知，隨著層理角度的增加，超聲波速度呈冪函數(shù)遞減規(guī)律。由數(shù)字模型速度云圖（圖7）可知，在層理密度相同的前提下，隨著層理角度的增大，在首波傳播路徑上的層理界面數(shù)目越來(lái)越多，導(dǎo)致聲波波速逐漸減小。

圖8 層理角度對(duì)波速的影響圖

3.2 層理角度對(duì)衰減系數(shù)的影響

由圖9可知：隨著層理角度增大，聲波衰減系數(shù)呈線性遞增。從數(shù)字模型速度云圖（圖7）可以看出：當(dāng)層理角度逐漸增大時(shí)，層理在水平方向阻礙超聲波在頁(yè)巖中的傳播路徑越來(lái)越多，并且在相同尺寸巖心中，層理數(shù)目隨著角度增大而增多，擴(kuò)展了超聲波的傳播路徑，導(dǎo)致巖石聲波能量消耗增大。同時(shí)，骨架與層理面之間的界面數(shù)量也隨之增多，進(jìn)而產(chǎn)生更多的反射、散射，導(dǎo)致超聲波穿透能量減弱，衰減系數(shù)增大。

4 結(jié)論

1）針對(duì)常規(guī)室內(nèi)超聲波透射法實(shí)驗(yàn)成本高、誤差較大的現(xiàn)狀，筆者基于黏彈性介質(zhì)聲波波動(dòng)理論，結(jié)合高階交錯(cuò)網(wǎng)格差分技術(shù)，建立了一種模擬頁(yè)巖超聲波透射實(shí)驗(yàn)的數(shù)值計(jì)算方法。

2）數(shù)值計(jì)算結(jié)果與物理實(shí)驗(yàn)得到的頁(yè)巖聲波特性變化趨勢(shì)吻合，同時(shí)，基于理想和真實(shí)巖心的數(shù)值模擬計(jì)算的衰減特性與物理實(shí)驗(yàn)結(jié)果一致，表明該方法具有較強(qiáng)的適應(yīng)性。

3）在層理尺度和密度恒定的條件下，隨著層理角度的增大，波速呈冪函數(shù)遞減，衰減系數(shù)呈線性增加。

4）該數(shù)值模擬方法可從微觀角度分析超聲波傳播規(guī)律，為進(jìn)一步研究分析頁(yè)巖層理結(jié)構(gòu)對(duì)聲學(xué)特性的影響提供了思路。

圖9 層理角度與衰減系數(shù)關(guān)系圖