頁(yè)巖層理結(jié)構(gòu)對(duì)超聲波特性響應(yīng)分析及應(yīng)用

徐烽淋 ，陳喬 ，朱洪林，王丹，陳吉龍，劉璞，姚光華，張闊，霍振永

（1. 重慶市涪陵頁(yè)巖氣環(huán)保研發(fā)與技術(shù)服務(wù)中心，重慶 408000；2. 中國(guó)科學(xué)院重慶綠色智能技術(shù)研究院，重慶 400714；3. 油氣藏地質(zhì)及開(kāi)發(fā)工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，成都 610500；4. 重慶礦產(chǎn)資源資源開(kāi)發(fā)有限公司，重慶 401120；5. 中國(guó)石油西南油氣田公司蜀南氣礦，四川瀘州 646300）

0 引言

頁(yè)巖復(fù)雜層理結(jié)構(gòu)會(huì)導(dǎo)致巖石的超聲波特性變得復(fù)雜，無(wú)法準(zhǔn)確求取儲(chǔ)集層的力學(xué)參數(shù)。Josh等[1]通過(guò)分析頁(yè)巖氣區(qū)塊的相關(guān)資料發(fā)現(xiàn)，頁(yè)巖復(fù)雜層理結(jié)構(gòu)導(dǎo)致聲波測(cè)井資料存在各向異性，Vernik等[2]、鄧?yán)^新等[3]通過(guò)頁(yè)巖波速的各向異性也得到類(lèi)似的結(jié)論。Horne等[4]利用偶極子聲波測(cè)井?dāng)?shù)據(jù)來(lái)判斷層狀頁(yè)巖彈性各向異性，估算出各向異性參數(shù)與鄰井預(yù)測(cè)值一致；張白紅等[5]測(cè)量了2種巖塊垂直和平行層理兩個(gè)方向的超聲波速，數(shù)據(jù)表明平行于層理面方向及其法線(xiàn)方向具有明顯的波速各向異性。Yan等[6]、程禮軍等[7-8]等分別研究頁(yè)巖在單軸、三軸壓縮條件下的波速變化，發(fā)現(xiàn)隨著軸向載荷的增加，縱、橫波速度呈先增后降的趨勢(shì)。Chen等[9]、熊健等[10]和吳濤等[11]通過(guò)超聲波透射實(shí)驗(yàn)發(fā)現(xiàn)，波速隨層理角度的增加或頻率的減小而減?。凰p系數(shù)隨角度或頻率的增加而增大，平行于層理方向的波速與層理密度呈線(xiàn)性正相關(guān)，而垂直于層理方向的波速與層理密度呈線(xiàn)性負(fù)相關(guān)。Schn等[12]基于超聲波實(shí)驗(yàn)建立了頁(yè)巖彈性性質(zhì)的函數(shù)形式，并通過(guò)該函數(shù)計(jì)算速度和各向異性。

綜上所述，國(guó)內(nèi)外學(xué)者主要利用露頭或者井下巖心制作不同層理結(jié)構(gòu)的物理模型，通過(guò)超聲波透射實(shí)驗(yàn)總結(jié)頁(yè)巖聲波特征響應(yīng)的定性認(rèn)識(shí)[13-14]?？紤]到頁(yè)巖制樣困難，同時(shí)，由于實(shí)驗(yàn)研究對(duì)頁(yè)巖的大量耗損導(dǎo)致無(wú)法開(kāi)展大規(guī)模的物理模型實(shí)驗(yàn)，致使研究結(jié)果具有一定的局限性。針對(duì)上述問(wèn)題，本文以頁(yè)巖層理結(jié)構(gòu)為研究對(duì)象，基于波動(dòng)理論，對(duì)超聲波透射實(shí)驗(yàn)過(guò)程進(jìn)行數(shù)值模擬，計(jì)算頁(yè)巖不同層理結(jié)構(gòu)模型的超聲波響應(yīng)規(guī)律，利用灰色系統(tǒng)理論篩選對(duì)層理結(jié)構(gòu)敏感的聲學(xué)參數(shù)，進(jìn)一步建立頁(yè)巖的動(dòng)態(tài)力學(xué)參數(shù)模型，研究結(jié)果對(duì)提高利用聲波資料預(yù)測(cè)巖石力學(xué)參數(shù)的準(zhǔn)確性具有重要意義。

1 數(shù)值模擬方法

1.1 波動(dòng)理論

彈性理論主要研究物體受力與形變的關(guān)系，其波動(dòng)方程由本構(gòu)方程、運(yùn)動(dòng)平衡微分方程和幾何方程共同建立。在彈性各向同性介質(zhì)中，二維一階速度-應(yīng)力彈性波動(dòng)方程為：

在時(shí)間上，采用二階精度的時(shí)間差分方程解的近似值，在空間采用四階差分方程解的近似值，利用交錯(cuò)網(wǎng)格格式，取

可推導(dǎo)出波動(dòng)方程的差分格式：

1.2 波動(dòng)方程求解

1.2.1 邊界條件

在計(jì)算區(qū)域中左、右邊界與實(shí)際條件基本相符，按照常規(guī)反射邊界來(lái)處理，而上、下邊界和 4個(gè)邊角點(diǎn)設(shè)置吸收條件。

1.2.2 初始條件

當(dāng)t≤0時(shí)，速度和位移為零，即

1.2.3 振源條件

1.2.4 收斂條件

分析誤差傳播和積累的情況，是保證差分格式收斂性的基礎(chǔ)。董良國(guó)等[15]詳細(xì)推導(dǎo)了交錯(cuò)網(wǎng)格高階差分方程的穩(wěn)定性條件。當(dāng)介質(zhì)為均勻各向同性時(shí)，時(shí)間精度為二階，空間精度為四階的穩(wěn)定條件為：

實(shí)驗(yàn)室超聲波測(cè)試采用透射法進(jìn)行測(cè)量。但由于頁(yè)巖的非均質(zhì)性及層理分布的不確定性，給研究層理結(jié)構(gòu)對(duì)超聲波傳播特性的影響帶來(lái)了較大困難，而且透射實(shí)驗(yàn)中巖心與探頭接觸不平整以及人為讀取初始時(shí)間等也會(huì)造成誤差。數(shù)值模擬基于被測(cè)巖心的縱向剖面建立模型，以超聲波探頭的激發(fā)信號(hào)為振源，基于波動(dòng)理論推導(dǎo)出一階位移-應(yīng)力彈性波方程組，并開(kāi)展交錯(cuò)網(wǎng)格半程計(jì)算，用Matlab編程實(shí)現(xiàn)超聲波透射實(shí)驗(yàn)的模擬。

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

1.3.1 波形對(duì)比驗(yàn)證

為了使模型與超聲波透射實(shí)驗(yàn)相統(tǒng)一，將探頭頻率設(shè)為250 kHz。其中入射波信號(hào)顯示為初始波形，當(dāng)入射信號(hào)透過(guò)巖心到達(dá)接收端時(shí)，實(shí)驗(yàn)室可獲取透射信號(hào)，顯示為實(shí)驗(yàn)波形，而數(shù)值計(jì)算的結(jié)果顯示為模擬波形（見(jiàn)圖1）。波形對(duì)比結(jié)果顯示，模擬波形與實(shí)驗(yàn)波形到達(dá)接收端的時(shí)間均為0.9×10-5s左右，電壓均為-5～5 V，二者范圍相近（見(jiàn)圖1）。

將模擬波形和實(shí)驗(yàn)波形進(jìn)行相關(guān)性計(jì)算，將模擬波形信號(hào)設(shè)為X，實(shí)驗(yàn)波形信號(hào)設(shè)為Y：

X和Y的標(biāo)準(zhǔn)差分別為：

協(xié)方差：

圖1 標(biāo)準(zhǔn)頁(yè)巖超聲波透射實(shí)驗(yàn)與數(shù)值模擬結(jié)果對(duì)比圖

通過(guò)相關(guān)性計(jì)算得到模擬波形與實(shí)驗(yàn)波形之間的相關(guān)系數(shù)為 0.85，吻合度較高，所以該數(shù)值模擬方法可對(duì)超聲波透射實(shí)驗(yàn)進(jìn)行有效的模擬。

1.3.2 數(shù)值計(jì)算驗(yàn)證

考慮到層理角度為頁(yè)巖最顯著的結(jié)構(gòu)特征，因此，結(jié)合實(shí)際巖心建立類(lèi)似不同層理角度的理想巖心模型，在此基礎(chǔ)上開(kāi)展數(shù)值模擬計(jì)算。圖2是渝東南地區(qū)下志留統(tǒng)龍馬溪組露頭巖心，編號(hào)記為A1、A2、A3、A4。層理角度為層理方向與超聲波傳播方向的夾角，其中 A1的層理角度為 0°，A2為 30°，A3為 60°，A4為 90°（見(jiàn)圖2）。為了進(jìn)一步驗(yàn)證數(shù)值模擬方法的有效性，本文還以A1、A2、A3、A4巖心為模板，建立數(shù)值模型B1、B2、B3、B4（見(jiàn)圖3）。在模型中對(duì)速度參數(shù)進(jìn)行設(shè)置，將超聲波在層理中的速度值設(shè)置為1 900 m/s，將其在骨架中的速度值設(shè)置為4 200 m/s，該模型為理想的層理角度模型，其角度與巖心 A1、A2、A3、A4一一對(duì)應(yīng)。在此基礎(chǔ)上，對(duì)超聲波透射實(shí)驗(yàn)與數(shù)值模擬結(jié)果進(jìn)行對(duì)比驗(yàn)證。

圖2 渝東南地區(qū)下志留統(tǒng)龍馬溪組露頭巖心照片（據(jù)參考文獻(xiàn)[9]修改）

通過(guò)超聲波透射實(shí)驗(yàn)完成對(duì)巖心A1、A2、A3、A4在不同探頭頻率條件下衰減系數(shù)的計(jì)算，結(jié)果顯示，在同一層理角度下，隨著探頭頻率的增大，衰減系數(shù)整體呈現(xiàn)出增大趨勢(shì)（見(jiàn)圖4a）。同時(shí)，通過(guò)數(shù)值模擬也完成了對(duì)模型B1、B2、B3、B4在相同條件下的計(jì)算（見(jiàn)圖4b），數(shù)值模擬計(jì)算得到的規(guī)律與超聲波透射實(shí)驗(yàn)規(guī)律基本吻合[16]，即在同一層理角度下，衰減系數(shù)與探頭頻率均呈正相關(guān)。

圖3 渝東南地區(qū)龍馬溪組巖心數(shù)值模型（據(jù)參考文獻(xiàn)[9]修改）

圖4 超聲波衰減特性結(jié)果對(duì)比圖（樣品數(shù)4個(gè)）

圖5 渝東南地區(qū)下志留統(tǒng)龍馬溪組露頭巖心照片

圖6 渝東南地區(qū)下志留統(tǒng)龍馬溪組巖心層理提取結(jié)果圖

上述數(shù)值計(jì)算是基于理想的層理角度模型。為了進(jìn)一步驗(yàn)證上述數(shù)值模擬結(jié)果的可靠性，利用數(shù)值圖像處理技術(shù)提取出真實(shí)巖心的層理角度，然后再開(kāi)展不同層理角度的衰減系數(shù)數(shù)值計(jì)算。圖5是渝東南地區(qū)下志留統(tǒng)龍馬溪組露頭巖心照片，編號(hào)記為C1、C2、C3、C4，其中 C1的層理角度為 0°，C2為 30°，C3為 45°，C4為 90°（見(jiàn)圖5）。

以C1、C2、C3、C4巖心為研究對(duì)象，對(duì)照片中的層理和基質(zhì)骨架區(qū)域進(jìn)行邊緣提取、Hough直線(xiàn)檢測(cè)[17]，建立4個(gè)數(shù)值模擬模型，編號(hào)分別記為D1、D2、D3、D4（見(jiàn)圖6），其層理角度與 C1、C2、C3、C4一一對(duì)應(yīng)。由于數(shù)字圖像是由矩形排列的像素點(diǎn)組成，每個(gè)像素點(diǎn)是一個(gè)小正方形，這些小正方形均可作為有限差分計(jì)算中的四邊形單元，故可將層理結(jié)構(gòu)提取后的數(shù)字圖像通過(guò)方形單元映射的方法可以很方便地轉(zhuǎn)化成有限差分網(wǎng)格數(shù)據(jù)。首先，每個(gè)單元的 4個(gè)節(jié)點(diǎn)坐標(biāo)通過(guò)圖像實(shí)際尺寸與像素尺寸之間的比例關(guān)系，轉(zhuǎn)換為相應(yīng)矢量空間的物理坐標(biāo)，整個(gè)圖像就可轉(zhuǎn)換為正方形有限差分網(wǎng)格，此法即為方形單元映射法。然后根據(jù)每個(gè)像素點(diǎn)所屬灰度值，給每個(gè)單元賦予相應(yīng)的參數(shù)，從而把數(shù)字圖像表征的層理結(jié)構(gòu)引入到數(shù)值模型中。最后，利用Matlab算法開(kāi)展數(shù)值模擬計(jì)算，并將計(jì)算結(jié)果與超聲波透射實(shí)驗(yàn)結(jié)果對(duì)比。從衰減系數(shù)結(jié)果圖來(lái)看，數(shù)值模擬的衰減系數(shù)與層理角度呈正相關(guān)，這與超聲波透射實(shí)驗(yàn)中巖心的衰減系數(shù)變化規(guī)律基本一致（見(jiàn)圖7）。

圖7 實(shí)驗(yàn)結(jié)果與數(shù)值模擬結(jié)果圖（樣品數(shù)4個(gè)）

綜上所述，針對(duì)超聲波透射實(shí)驗(yàn)成本高，由于實(shí)驗(yàn)過(guò)程中巖心與探頭接觸不平整、人為數(shù)據(jù)讀取等造成誤差較大的現(xiàn)狀[18]，本文提出的頁(yè)巖超聲波透射實(shí)驗(yàn)的數(shù)值模擬方法與實(shí)際超聲波透射實(shí)驗(yàn)結(jié)果吻合度較高，因此，可依托此數(shù)值模擬計(jì)算手段，開(kāi)展頁(yè)巖層理結(jié)構(gòu)聲學(xué)特性研究。

2 頁(yè)巖層理結(jié)構(gòu)對(duì)超聲波特征的影響

2.1 層理結(jié)構(gòu)對(duì)超聲波特性的影響

圖8 不同層理厚度對(duì)超聲波特性的影響（模型個(gè)數(shù)為22個(gè)）

頁(yè)巖的顯著特征是層理發(fā)育，表征層理結(jié)構(gòu)的參數(shù)主要包括層理厚度、層理密度和層理角度。為了開(kāi)展與超聲波透射實(shí)驗(yàn)條件相吻合的數(shù)值模擬計(jì)算，本文將模型中的巖心柱塞高度設(shè)置為50 mm，直徑為25 mm，頁(yè)巖骨架中的縱波速度設(shè)為4 200 m/s，橫波速度設(shè)為2 500 m/s；頁(yè)巖層理中的縱波速度為1 900 m/s，橫波波速為1 100 m/s；頁(yè)巖骨架密度為2.6 g/cm3，層理密度為 2.3 條/cm2。提取實(shí)驗(yàn)室探頭入射波作為初始波，保持超聲波傳播方向不變，將柱塞縱向剖面區(qū)域劃分成長(zhǎng)度為250 mm、寬為125 mm的網(wǎng)格，空間網(wǎng)格間距為0.2 mm，時(shí)間步長(zhǎng)為20 ns。通過(guò)不斷調(diào)整柱塞的層理模型參數(shù)來(lái)實(shí)現(xiàn)頁(yè)巖結(jié)構(gòu)的變化，完成對(duì)超聲波透射實(shí)驗(yàn)的模擬，為探索層理結(jié)構(gòu)對(duì)超聲波特性的影響奠定基礎(chǔ)。

2.1.1 層理厚度對(duì)超聲波特性的影響

模型參數(shù)設(shè)置中，我們將平直且與層面平行的層理稱(chēng)為水平層理，定義平直且與層面垂直的層理為垂直層理。在層理角度和密度相同的情況下，將垂直層理厚度設(shè)為 0，0.2，0.4，0.6，0.8，1.0，1.2，1.4，1.6，1.8，2.0 mm，此組模型編號(hào)記為E1—E11。而另一組模型中水平層理厚度與垂直層理厚度設(shè)為相同，編號(hào)記為F1—F11。兩組模型中超聲波傳播方向均為由其設(shè)置的巖心柱塞底端向頂端傳播。本文從中選取了層理厚度結(jié)構(gòu)特征比較明顯的 E2、E10和 F2、F10模型，分析兩組模型不同層理厚度對(duì)超聲波特性的影響，并通過(guò)其波場(chǎng)快照?qǐng)D進(jìn)一步開(kāi)展層理厚度對(duì)超聲波特性影響的微觀(guān)分析。其中，E2的層理厚度為0.2 mm，E10的層理厚度為1.8 mm，F(xiàn)2層理厚度為0.2 mm，F(xiàn)10層理厚度為1.8 mm。由波速數(shù)值計(jì)算結(jié)果可知，隨著層理厚度的增加，無(wú)論水平層理還是垂直層理，波速均呈線(xiàn)性遞減，水平層理的遞減幅度更大（見(jiàn)圖8a）。衰減系數(shù)計(jì)算結(jié)果顯示，隨著層理厚度增加，無(wú)論水平層理還是垂直層理，衰減系數(shù)都呈增大趨勢(shì)，其中水平層理的衰減系數(shù)以?xún)绾瘮?shù)規(guī)律遞增，而垂直層理以指數(shù)函數(shù)規(guī)律遞增（見(jiàn)圖8b）。超聲波主頻和主振幅計(jì)算結(jié)果顯示，隨著層理厚度的變化，水平層理中的超聲波主頻、垂直層理中的超聲波主頻與垂直層理中的超聲波主振幅的變化無(wú)規(guī)律，而水平層理中的超聲波主振幅隨層理厚度的增大整體呈下降趨勢(shì)（見(jiàn)圖8c）。

前文預(yù)設(shè)模型中 E10垂直層理厚度比 E2的垂直層理厚度大，超聲波由巖心底端到頂端傳播。波場(chǎng)快照?qǐng)D顯示，相同時(shí)刻下的超聲波信號(hào)在E10中傳播所能到達(dá)的距離比在E2中的更短（見(jiàn)圖9a—圖9b），所以導(dǎo)致超聲波波速在E10中也較?。ㄒ?jiàn)圖8a）。超聲波在F2與F10中的水平層理傳播也呈同樣規(guī)律，即隨層理厚度增大，在相同時(shí)刻下超聲波信號(hào)在F10中傳播距離也更短，說(shuō)明超聲波波速在F10中也較小。圖9中使用不同的聲場(chǎng)幅值來(lái)表征超聲波透射能量大小，直觀(guān)可見(jiàn)E10的顏色總體較E2更淺，說(shuō)明隨著垂直層理厚度增大，E10中透射超聲波能量更?。ㄒ?jiàn)圖9b），所以其衰減系數(shù)也較大（見(jiàn)圖8b）。在模型F2和F10中的水平層理中也表現(xiàn)出與E2、E10同樣的變化規(guī)律，即隨層理厚度增大，透射超聲波能量減?。ㄒ?jiàn)圖9c—圖9d），衰減系數(shù)增大（見(jiàn)圖8b）。

圖9 層理厚度模型的波場(chǎng)快照

2.1.2 層理密度對(duì)超聲波特性的影響

在層理角度和厚度恒定的條件下，通過(guò)調(diào)整層理數(shù)目來(lái)模擬層理密度變化，即巖心剖面單位面積上層理的數(shù)目。將垂直層理厚度設(shè)為0.4 mm，層理數(shù)目設(shè)為 0、1、2、4、5、6、10、12、15、18、20 條，對(duì)應(yīng)的層理密度為 0、0.08×104、0.16×104、0.32×104、0.40×104、0.48×104、0.80×104、0.96×104、1.20×104、1.44×104和 1.60×104條/m2，此組模型編號(hào)記為 G1—G11。另一組模型編號(hào)記為 H1—H11，為水平層理，并將層理密度設(shè)置與上述一致。從中選取了能較好反映層理密度結(jié)構(gòu)特征的 G3、G7和 H3、H7模型分析兩組模型不同層理密度對(duì)超聲波特性的影響，并通過(guò)其波場(chǎng)快照?qǐng)D，開(kāi)展層理密度對(duì)超聲波特性影響的微觀(guān)分析。其中，G3的層理密度為0.16×104條/m2，G7的層理密度為 0.80×104條/m2，H3層理密度為 0.16×104條/m2，H7層理密度為 0.80×104條/m2。

由波速計(jì)算結(jié)果可知，隨著層理密度的增加，無(wú)論水平層理還是垂直層理，波速均呈線(xiàn)性遞減，垂直層理下降幅度相對(duì)較?。ㄒ?jiàn)圖10a）。衰減系數(shù)數(shù)值計(jì)算結(jié)果顯示，衰減系數(shù)隨著層理密度的增加呈線(xiàn)性遞增（見(jiàn)圖10b）。超聲波主頻和主振幅數(shù)值計(jì)算結(jié)果顯示，隨著層理密度的增加，垂直層理中的超聲波主頻和主振幅變化無(wú)規(guī)律；但是在水平層理中，隨著層理密度的增加，層理密度小于 0.3×104條/m2時(shí)，主頻和主振幅呈下降趨勢(shì)，層理密度大于 0.3×104條/m2后開(kāi)始緩慢上升直到基本保持不變（見(jiàn)圖10c）。

圖10 不同層理密度對(duì)超聲波特性的影響（模型個(gè)數(shù)為22個(gè)）

前文已對(duì)模型預(yù)設(shè)層理密度，由于模型中G7層理密度比G3層理密度大，H7層理密度比H3層理密度大，超聲波由巖心底端到頂端傳播，在G7中超聲波通過(guò)層理與骨架的界面數(shù)目要比G3多，波場(chǎng)快照?qǐng)D顯示，相同時(shí)間下的超聲波信號(hào)在 G7中傳播的距離比在 G3中的更短（見(jiàn)圖11），所以超聲波速度也較?。ㄒ?jiàn)圖10a）。超聲波在 H7與 H3中的水平層理傳播也呈同樣規(guī)律，即在相同時(shí)間下其超聲波信號(hào)傳播距離也更短，說(shuō)明超聲波波速也較小。在衰減特性方面，隨著層理密度的增加，G7中層理?xiàng)l數(shù)比G3更多，非均質(zhì)性更強(qiáng)。由于層理與骨架之間的模量不等，在骨架與層理之間產(chǎn)生一定的應(yīng)變相位差，從而在界面處產(chǎn)生內(nèi)摩擦，使機(jī)械能轉(zhuǎn)為熱能，極大消耗了能量[19]。圖11中使用不同的聲場(chǎng)幅值來(lái)表征超聲波透射能量大小，直觀(guān)可見(jiàn)G7顏色整體較G3更淺（見(jiàn)圖11）；因此，G7中的超聲波透射能量要明顯低于G3，G7對(duì)應(yīng)的衰減系數(shù)更大（見(jiàn)圖10b）。H7與H3類(lèi)似，也呈同樣規(guī)律，即隨著層理密度的增大，H7超聲波透射能量更小，說(shuō)明衰減系數(shù)更大。

圖11 層理密度模型的波場(chǎng)快照?qǐng)D

2.1.3 層理角度對(duì)超聲波特性的影響

在層理厚度和密度恒定的條件下，將層理角度為0°、10°、20°、30°、40°、50°、60°、70°、80°、90°的模型編號(hào)記為 I1—I10。文中選取了層理角度特征比較明顯的I3、I8，并通過(guò)其波場(chǎng)快照?qǐng)D進(jìn)一步開(kāi)展層理厚度對(duì)超聲波特性影響的微觀(guān)分析。其中，模型 I3的層理角度為30°，I8的層理角度為80°。

數(shù)值計(jì)算結(jié)果顯示，隨著層理角度的增加，波速呈冪函數(shù)遞減規(guī)律（見(jiàn)圖12a），而衰減系數(shù)呈線(xiàn)性遞增規(guī)律（見(jiàn)圖12b）。

由于模型預(yù)設(shè)I8的層理角度大于I3的層理角度，超聲波由底端向頂端傳播，在層理密度相同的前提下，I8的層理界面數(shù)目多于I3，波場(chǎng)快照?qǐng)D顯示，在相同時(shí)刻下I8中超聲波信號(hào)傳播所能到達(dá)的距離比在I3中的更短（見(jiàn)圖13a—圖13b），對(duì)應(yīng)的波速也就較?。ㄒ?jiàn)圖12a）。在衰減特性方面，由于I8中的骨架與層理面之間的界面?zhèn)€數(shù)多于 I3，導(dǎo)致阻礙超聲波傳播的路徑也多于I3，圖13中使用不同的聲場(chǎng)幅值來(lái)表征超聲波透射能量大小，直觀(guān)可見(jiàn)I8的顏色總體較 I3更淺，說(shuō)明透射超聲波能量更小，所以衰減系數(shù)也較大，對(duì)應(yīng)的衰減系數(shù)更大（見(jiàn)圖12b）。

圖12 不同層理角度對(duì)超聲波特性的影響圖（模型個(gè)數(shù)為11）

圖13 層理角度模型的波場(chǎng)快照?qǐng)D

2.2 超聲波特性統(tǒng)計(jì)分析

在頁(yè)巖數(shù)值模擬計(jì)算結(jié)果的基礎(chǔ)上，利用灰色系統(tǒng)關(guān)聯(lián)法[20]定量分析聲學(xué)參數(shù)對(duì)層理結(jié)構(gòu)的敏感性。首先，將超聲波速度、衰減系數(shù)等聲學(xué)參數(shù)作為參考序列，將層理角度、密度和厚度變化等層理結(jié)構(gòu)變化參數(shù)分別作為比較序列：

對(duì)序列中的數(shù)據(jù)進(jìn)行均值無(wú)量綱化處理，即：

求取序列的關(guān)聯(lián)系數(shù)：

獲取整個(gè)數(shù)列的關(guān)聯(lián)度：

聲學(xué)參數(shù)中不僅速度對(duì)層理角度和層理密度的變化較為敏感（關(guān)聯(lián)度值均為0.80以上），衰減系數(shù)對(duì)層理厚度的變化也較敏感（關(guān)聯(lián)度值可達(dá)0.89），而主振幅、主頻與層理結(jié)構(gòu)相關(guān)聯(lián)度則相對(duì)較?。P(guān)聯(lián)度均在 0.60以下）（見(jiàn)表1）。因此，在進(jìn)一步利用聲學(xué)參數(shù)表征頁(yè)巖層理結(jié)構(gòu)的過(guò)程中，除了選用常規(guī)敏感因子速度之外，還需配合衰減系數(shù)來(lái)綜合描述。

表1 層理結(jié)構(gòu)參數(shù)和各個(gè)影響因子的關(guān)聯(lián)度

3 頁(yè)巖動(dòng)態(tài)力學(xué)參數(shù)模型研究

3.1 頁(yè)巖動(dòng)態(tài)力學(xué)參數(shù)模型建立

利用巖石波速建立動(dòng)態(tài)力學(xué)參數(shù)模型是現(xiàn)場(chǎng)油氣井開(kāi)發(fā)中最常用的方法。其中各向同性模型和橫向各向同性模型表達(dá)式如下所示：

各向同性巖石力學(xué)參數(shù)預(yù)測(cè)模型：

橫向各向同性巖石力學(xué)參數(shù)預(yù)測(cè)模型[21]：

通過(guò)數(shù)值模擬研究發(fā)現(xiàn)，僅通過(guò)波速建立頁(yè)巖的動(dòng)態(tài)力學(xué)參數(shù)模型是不可行的，還需配合對(duì)巖石各向異性特征敏感的衰減系數(shù)，方能達(dá)到更好效果。利用實(shí)驗(yàn)測(cè)得的渝東南彭水地區(qū)井下龍馬溪組頁(yè)巖的波速和衰減系數(shù)，基于麥夸特與全局優(yōu)化算法[22]在考慮量綱的基礎(chǔ)上，建立了該地區(qū)頁(yè)巖的動(dòng)態(tài)力學(xué)參數(shù)模型：

3.2 模型室內(nèi)對(duì)比驗(yàn)證

采用渝東南ZY1井和YY1井的井下巖心開(kāi)展基礎(chǔ)物性實(shí)驗(yàn)獲取巖心密度，通過(guò)巖石力學(xué)實(shí)驗(yàn)獲取彈性模量和泊松比，利用超聲波透射實(shí)驗(yàn)獲取縱、橫波速度（見(jiàn)表2）。在此基礎(chǔ)上，對(duì)本文模型進(jìn)行驗(yàn)證。將表2中實(shí)驗(yàn)結(jié)果代入由公式（14）、（15）所建立的各向同性模型，以及由公式（16）、（17）所建立的橫向各向同性模型，開(kāi)展模型的對(duì)比分析。

表2 井下巖心實(shí)驗(yàn)結(jié)果

驗(yàn)證結(jié)果表明：利用本文模型計(jì)算得到的動(dòng)、靜態(tài)力學(xué)參數(shù)相關(guān)性整體較好，其中在ZY1井中靜態(tài)泊松比與動(dòng)態(tài)泊松比相關(guān)系數(shù)可達(dá) 0.83，靜態(tài)彈性模量與動(dòng)態(tài)彈性模量相關(guān)系數(shù)可達(dá)0.84。在YY1井中靜態(tài)泊松比與動(dòng)態(tài)泊松比相關(guān)系數(shù)可達(dá) 0.82，靜態(tài)彈性模量與動(dòng)態(tài)彈性模量相關(guān)系數(shù)可達(dá)0.81（見(jiàn)圖14）。本文的新模型模擬效果優(yōu)于通常使用的各向同性模型和橫向各向同性模型，其中各向同性模型是指不同方向上物理性質(zhì)相同的介質(zhì)，為均質(zhì)模型；橫向各向同性模型是指垂直于軸線(xiàn)的任何方向上彈性相同，為非均質(zhì)模型。究其原因，是因?yàn)楦飨蛲阅Ｐ椭会槍?duì)各向同性的巖石，未考慮頁(yè)巖的層理發(fā)育，而橫向各向同性模型將頁(yè)巖等效為橫向各向同性體處理，模型結(jié)構(gòu)過(guò)于理想化，影響準(zhǔn)確性，而本文模型是基于表征頁(yè)巖層理結(jié)構(gòu)最為敏感的波速和衰減參數(shù)，利用通用全局優(yōu)化算法求得，模型可較好地反映頁(yè)巖層理結(jié)構(gòu)的影響，而且本文計(jì)算得到的巖石動(dòng)態(tài)模量一般大于靜態(tài)模量，動(dòng)態(tài)泊松比與靜態(tài)泊松比基本為1︰1關(guān)系，與Kuhlman等得到的關(guān)系基本一致[23]。

3.3 模型現(xiàn)場(chǎng)驗(yàn)證

收集渝東南地區(qū)ZY2井的測(cè)井資料開(kāi)展巖石力學(xué)動(dòng)態(tài)參數(shù)剖面反演，獲取靜態(tài)力學(xué)參數(shù)剖面結(jié)果。從圖15中看出，在反演深度范圍內(nèi)，動(dòng)態(tài)彈性模量數(shù)值普遍大于靜態(tài)彈性模量數(shù)值，動(dòng)態(tài)泊松比與靜態(tài)泊松比數(shù)值很接近，基本為 1︰1，這與前文所述模型計(jì)算結(jié)果保持一致。同時(shí)，取ZY2井龍馬溪組不同深度的10塊巖心開(kāi)展巖石力學(xué)實(shí)驗(yàn)，將實(shí)驗(yàn)結(jié)果與現(xiàn)場(chǎng)測(cè)井資料的反演數(shù)據(jù)進(jìn)行驗(yàn)證。

圖14 動(dòng)態(tài)力學(xué)參數(shù)模型的對(duì)比驗(yàn)證圖（R—相關(guān)系數(shù)）

圖15 ZY2井下志留統(tǒng)龍馬溪組巖石力學(xué)參數(shù)預(yù)測(cè)結(jié)果圖

由表3可見(jiàn)，10個(gè)不同深度反演得到的力學(xué)參數(shù)值與對(duì)應(yīng)井下巖心的實(shí)驗(yàn)測(cè)試值基本一致。通過(guò)計(jì)算兩者的相對(duì)誤差可知，泊松比的相對(duì)誤差小于 10%，平均相對(duì)誤差為 6.94%；彈性模量最大相對(duì)誤差為21.25%，平均誤差為15.30%。可見(jiàn)，利用測(cè)井資料反演的力學(xué)剖面能準(zhǔn)確地反映地層力學(xué)參數(shù)真實(shí)情況，本文建立的動(dòng)態(tài)力學(xué)參數(shù)模型對(duì)層理發(fā)育的頁(yè)巖地層有較好的適應(yīng)性。

表3 現(xiàn)場(chǎng)靜態(tài)力學(xué)參數(shù)驗(yàn)證誤差分析表

4 結(jié)論

基于波動(dòng)理論，利用交錯(cuò)網(wǎng)格有限差分法，獲取了一種頁(yè)巖超聲波透射實(shí)驗(yàn)的數(shù)值模擬方法。該方法模擬波形與實(shí)驗(yàn)波形相關(guān)系數(shù)大于 80%，表明數(shù)值模擬方法可有效模擬超聲波透射實(shí)驗(yàn)。

超聲波參數(shù)中，聲波速度是用來(lái)表征頁(yè)巖層理結(jié)構(gòu)的常規(guī)敏感因子，但還需要利用衰減系數(shù)的歸一化結(jié)果來(lái)綜合描述頁(yè)巖層理，可使結(jié)果更加準(zhǔn)確。

利用頁(yè)巖超聲波透射實(shí)驗(yàn)的波速和衰減系數(shù)建立了渝東南下志留統(tǒng)龍馬溪組超聲波預(yù)測(cè)頁(yè)巖力學(xué)參數(shù)模型，此模型計(jì)算得到的動(dòng)、靜力學(xué)參數(shù)相關(guān)性?xún)?yōu)于傳統(tǒng)模型，利用模型和測(cè)井資料反演而獲取的巖石力學(xué)剖面預(yù)測(cè)值與實(shí)驗(yàn)值吻合較好。

符號(hào)注釋?zhuān)?/p>

A，B——速度u和ν的離散值，無(wú)因次；Cn(k)——待定差分系數(shù)，無(wú)因次；D——應(yīng)力τxx的離散值，無(wú)因次；D0——標(biāo)準(zhǔn)差，無(wú)因次；E——應(yīng)力τxz的離散值，無(wú)因次；E0——期望值，無(wú)因次；Et——彈性模量，MPa；Eβ——任意角度β方向的彈性模量，MPa；E1、E2——表示平行和垂直各向同性面的彈性模量，MPa；F——應(yīng)力τzz的離散值，無(wú)因次；G2——各向同性面內(nèi)的剪切模量，MPa；i，j，k——空間和時(shí)間網(wǎng)格點(diǎn)，無(wú)因次；p1、p2、p3，q1、 q2、q3——模型系數(shù)，無(wú)因次；ri——Xi(k)與Yi(k)的關(guān)聯(lián)度，無(wú)因次；R——相關(guān)系數(shù)，無(wú)因次；R1——聲波震動(dòng)位移，m；Δt——時(shí)間網(wǎng)格間距，無(wú)因次；u——x方向的速度分量，m/s；ν——z方向的速度分量，m/s；U——x方向的位移函數(shù)，無(wú)因次；V——z方向的位移函數(shù)，無(wú)因次；Vp、Vs——巖心縱、橫波速度，m/s；ν1、ν2——分別表示平行和垂直各向同性面的泊松比，無(wú)因次；Xi——表示比較序列，無(wú)因次；Yi——表示參考序列，無(wú)因次；Δx，Δz——空間網(wǎng)格間距，無(wú)因次；β——層理面與巖樣端面間的夾角，°；λ，μ——拉梅常數(shù)，無(wú)因次；ρj——巖心密度，g/cm3；ρij——表示質(zhì)點(diǎn)在x、z坐標(biāo)軸中密度的離散值，g/cm3；τxx、τxz、τzz——應(yīng)力分量，MPa；χ——?dú)w一化衰減系數(shù)，無(wú)因次。