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

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

0 引言

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

在數(shù)值模擬方面，Saenger等[17]利用交錯(cuò)有限差分法來研究孔隙度和孔洞密度與聲波的響應(yīng)關(guān)系，幾年后，又利用旋轉(zhuǎn)的交錯(cuò)差分網(wǎng)格對(duì)含黏彈性液體孔洞、無黏彈性天然氣條件下的波速進(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ì)算精確度，通過費(fèi)曼原理和各向同性非雙曲時(shí)差方程兩種方法計(jì)算發(fā)現(xiàn)，P—SV轉(zhuǎn)換波在TI介質(zhì)各向同性平面?zhèn)鞑ブ芯哂邢嗤霓D(zhuǎn)換點(diǎn)和時(shí)程；國內(nèi)Yao等[25]利用“疊加型”多尺度隨機(jī)介質(zhì)模型、“區(qū)域多尺度隨機(jī)介質(zhì)模型”來模擬隨機(jī)溶洞介質(zhì)；魏建新等[26]以棍形、片形、球形、柱形、裂隙孔隙形狀形態(tài)的物理模型設(shè)計(jì)了填充物相同、橫向、縱向尺度不同的地震模型，在此基礎(chǔ)上，通過數(shù)值模擬計(jì)算，分析了由孔隙結(jié)構(gòu)變化引起的地震響應(yīng)特征；劉向君等[27-28]基于二維聲波數(shù)值模擬計(jì)算，開展裂縫模型聲波衰減系數(shù)研究，同時(shí)還開展了孔洞型儲(chǔ)層孔隙度預(yù)測(cè)模型研究。

綜上所述，國內(nèi)外學(xué)者主要通過超聲波透射實(shí)驗(yàn)獲得頁巖聲波傳播規(guī)律的定性認(rèn)識(shí)。同時(shí)，考慮到在物理模型實(shí)驗(yàn)過程中頁巖巖樣制備較為困難，國內(nèi)學(xué)者對(duì)作為油氣儲(chǔ)集層蓋層的頁巖地層聲波速度各向異性的巖石物理數(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í)別，以頁巖地層聲波測(cè)井為目的的聲學(xué)響應(yīng)鮮有報(bào)道，相關(guān)研究將頁巖當(dāng)作彈性介質(zhì)處理，未考慮頁巖的黏彈性特征。因此，筆者基于四川盆地渝東南下志留統(tǒng)龍馬溪組頁巖巖心，建立了不同角度的層理結(jié)構(gòu)模型，采用黏彈性介質(zhì)波動(dòng)理論和高階交錯(cuò)網(wǎng)格有限差分方法，開展頁巖超聲波透射實(shí)驗(yàn)數(shù)值模擬研究，計(jì)算不同層理角度頁巖的聲學(xué)參數(shù)特征，總結(jié)頁巖聲波參數(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)方程來開展數(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ù)泰勒展開式的基本原理，速度分量vz在2階時(shí)間精度和2n階空間精度下的離散量為：

式中W表示速度分量（vz）的離散量；S、T分別表示應(yīng)力分量σzz、τxz的離散量；Cn表示差分階數(shù)；n表示空間域的階數(shù)；Δt表示時(shí)間步長；Δx、Δz分別表示空間步長；（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ù)

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

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

1.5 邊界條件

在數(shù)值模擬的過程中，引入完全匹配層吸收模型，在給定的模型外面構(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)證

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

巖石超聲波衰減系數(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ì)算是基于理想的巖心層理角度模型。接下來，為了進(jìn)一步驗(yàn)證數(shù)值模擬結(jié)果的有效性，采用數(shù)值圖像處理方法提取出渝東南龍馬溪組真實(shí)巖心的層理角度，開展不同層理角度的衰減系數(shù)數(shù)值計(jì)算。

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

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

圖4是巖樣的實(shí)物照片（記為A1～A4），從照片可以看出，頁巖層理發(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)又可通過圖像實(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）。從速度云圖來看，B1～B4的層理傾角分別分布在0°、45°、60°和90°，與實(shí)驗(yàn)巖樣吻合較好。最后，利用Matlab編寫代碼進(jìn)行數(shù)值模擬計(jì)算，從圖6所示的結(jié)果來看，衰減系數(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)一步說明本數(shù)值模擬方法的有效性。

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

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

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

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

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

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

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

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

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

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

由圖9可知：隨著層理角度增大，聲波衰減系數(shù)呈線性遞增。從數(shù)字模型速度云圖（圖7）可以看出：當(dāng)層理角度逐漸增大時(shí)，層理在水平方向阻礙超聲波在頁巖中的傳播路徑越來越多，并且在相同尺寸巖心中，層理數(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ù)，建立了一種模擬頁巖超聲波透射實(shí)驗(yàn)的數(shù)值計(jì)算方法。

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

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

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

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