基于黏彈性介質波動理論的頁巖超聲波數值模擬

陳 喬 徐烽淋 程 亮 劉 洪 簡 旭 朱洪林 陳吉龍

1.中國科學院重慶綠色智能技術研究院 2.重慶市涪陵頁巖氣環保研發與技術服務中心3.中國石油川慶鉆探工程公司地質勘探開發研究院 4.重慶地質礦產研究院5.“油氣藏地質及開發工程”國家重點實驗室·西南石油大學 6. 重慶工商大學

0 引言

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

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

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

1 數值模擬理論基礎

數值模擬中常用的波動方程可分為二階彈性波動方程和一階位移—應力彈性波方程組2種。由于一階位移—應力彈性波動方程不需對彈性常數進行空間微分，利用交錯網格半程計算可以避免求解時間的高階微分，同時，網格引入的頻散較小，筆者采用一階位移—應力彈性波動方程來開展數值模擬研究。

1.1 聲波波動理論

聲波波動方程所包括的本構方程、運動微分方程和幾何方程揭示了物體內部應變、應力和位移之間的相互關系，其中本構方程主要反映的是物體的各向異性、黏滯性和彈性等性質。Kelvin—Voigt黏彈性介質模型能較好地描述聲波在介質中的衰減[29]?；趹Α獞冴P系、本構方程以及Kelvin—Voigt黏彈性模型，得到了二維黏彈性介質的波動方程，為了便于對其做有限差分，將二維黏彈性波動方程轉化為一階速度—應力方程，即

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

1.2 數值離散

根據泰勒展開式的基本原理，速度分量vz在2階時間精度和2n階空間精度下的離散量為：

式中W表示速度分量（vz）的離散量；S、T分別表示應力分量σzz、τxz的離散量；Cn表示差分階數；n表示空間域的階數；Δt表示時間步長；Δx、Δz分別表示空間步長；（i，j）表示網格節點標號；k表示時間節點標號。

1.3 穩定性分析

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

1.4 振源函數

為了更好地模擬頁巖的超聲波特性，振源函數是模擬頁巖聲波響應的關鍵，數值模擬中選取的聲源為超聲波探頭原始信號為：

式中R表示震源函數；x0、z0分別表示坐標值。

1.5 邊界條件

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

式中Gb表示邊界矩形區域的表達式；Gc表示邊界是個角區的函數表達式；α表示系數；N表示邊界網格的層數；（i0，j0）表示4個角的內側頂點標號。

2 數值模擬驗證

2.1 模擬波形與實驗波形驗證

為了與物理實驗相統一，針對室內標準巖心的超聲波透射實驗，筆者采用250 kHz的“探頭頻率”表示超聲波頻率大小。為了驗證二維時域有限差分程序的正確性，計算了均勻地層條件下的聲波波形，并將有限差分結果與實軸積分結果進行對比，對比結果如圖1所示。圖1中藍線為有限差分計算的結果，紅線表示實軸積分法的計算結果。有限差分法和實軸積分法的模擬計算結果基本一致，驗證了有限差分數值模擬方法的正確性。

圖1 二維有限差分與實軸積分波形對比圖

2.2 真實巖心數值計算及驗證

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

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

圖2 物理實驗與數值計算模型對比圖

圖3 超聲波衰減特性結果對比圖

針對巖心重復性較差的緣故，上述數值模擬計算是基于理想的巖心層理角度模型。接下來，為了進一步驗證數值模擬結果的有效性，采用數值圖像處理方法提取出渝東南龍馬溪組真實巖心的層理角度，開展不同層理角度的衰減系數數值計算。

圖4 實物巖樣照片

圖5 實物巖樣模型速度云圖

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

綜上所述，針對常規室內超聲波透射法實驗成本高、誤差較大的局限性，筆者提出了一種模擬頁巖超聲波透射實驗的數值計算方法，該方法與實際物理模型實驗結果吻合度較高，具有進一步通過大規模數值模擬計算，開展頁巖不同層理結構聲學特征響應研究的價值。

圖6 室內巖樣實驗與數值模擬計算所得衰減系數結果對比圖

3 數值計算分析

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

以不同層理角度的頁巖數值模型為基礎開展數值模擬計算，得到對應的波形，從中提取巖石聲波的衰減系數和速度參數，分析層理角度對聲學特性的影響。

圖7 數字模型速度云圖

3.1 層理角度對聲波波速的影響

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

圖8 層理角度對波速的影響圖

3.2 層理角度對衰減系數的影響

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

4 結論

1）針對常規室內超聲波透射法實驗成本高、誤差較大的現狀，筆者基于黏彈性介質聲波波動理論，結合高階交錯網格差分技術，建立了一種模擬頁巖超聲波透射實驗的數值計算方法。

2）數值計算結果與物理實驗得到的頁巖聲波特性變化趨勢吻合，同時，基于理想和真實巖心的數值模擬計算的衰減特性與物理實驗結果一致，表明該方法具有較強的適應性。

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

4）該數值模擬方法可從微觀角度分析超聲波傳播規律，為進一步研究分析頁巖層理結構對聲學特性的影響提供了思路。

圖9 層理角度與衰減系數關系圖