基于電-力-聲多物理場耦合數(shù)值模型的含水合物多孔介質聲速和衰減特性研究*

徐 源，張歡歡，邢蘭昌?，魏 偉，韓維峰，魏周拓，鄧少貴

基于電-力-聲多物理場耦合數(shù)值模型的含水合物多孔介質聲速和衰減特性研究*

徐 源1，張歡歡1，邢蘭昌1?，魏 偉2，韓維峰2，魏周拓3，鄧少貴3

（1. 中國石油大學（華東）控制科學與工程學院，山東 青島 266580；2. 中國石油勘探開發(fā)研究院新能源研究所，河北 廊坊 065007；3. 中國石油大學（華東）地球科學與技術學院，山東 青島 266580）

在物理模擬實驗中對水合物微觀賦存模式和飽和度進行準確控制和評價尚存在技術困難，僅依賴實驗技術研究含水合物沉積物聲學特性、建立儲層參數(shù)解釋模型存在局限性。采用基于有限元的數(shù)字巖石物理技術，針對懸浮、接觸和膠結三種典型的水合物微觀賦存模式分別建立多孔介質的三維電?力?聲多物理場耦合模型，考察了微觀賦存模式和水合物飽和度對多孔介質聲速和衰減的影響規(guī)律，對比了聲速數(shù)值模擬與理論模型計算結果，建立了聲波衰減參數(shù)與水合物飽和度之間的關系式。研究結果表明：（1）對于三種水合物賦存模式，由于水合物相比孔隙水具有更高的彈性模量，多孔介質的聲速隨著水合物飽和度的增大而增大；水合物的存在導致聲波在傳播過程中遇到更多不連續(xù)的聲阻抗界面，聲衰減隨著水合物飽和度的增大而近似線性增大；（2）懸浮和接觸賦存模式條件下，水合物飽和度對多孔介質的聲速和衰減影響規(guī)律基本一致；對于相同的水合物飽和度，膠結模式條件下含水合物多孔介質具有更高的聲速和更小的聲衰減；（3）通過合理選擇參數(shù)值，利用權重方程與Lee改進的Biot-Gassmann Theory（BGTL）模型估算的含懸浮和接觸模式水合物多孔介質的聲速較為準確；通過等效介質理論模型C計算的含膠結模式水合物多孔介質的聲速更為準確。研究結果可為獲取復雜條件下含水合物沉積物的聲學特性提供數(shù)值建模方法，為基于聲波測井數(shù)據(jù)的水合物儲層精細評價提供理論支撐。

天然氣水合物；數(shù)字巖石物理；微觀賦存模式；水合物飽和度；聲速；聲衰減

0 引 言

天然氣水合物是由水和天然氣分子在低溫高壓環(huán)境下結合形成的與冰相似的化合物，其廣泛分布于海洋沉積物和凍土區(qū)域中[1-4]，是一種極具開發(fā)前景且儲量巨大的新型清潔能源[5-7]。地球物理測井是準確識別和評價水合物儲層的重要手段，其中聲波測井是評價水合物飽和度（h）的重要方法之一[8-9]。研究含水合物沉積物的聲學響應特性，深入探究水合物飽和度、微觀賦存模式等因素對沉積物聲波傳播速度及其能量衰減的影響規(guī)律，將為聲波測井數(shù)據(jù)解釋、水合物儲層識別與評價提供理論基礎。

天然氣水合物在儲層中具有多種微觀賦存模式，不同的微觀賦存模式會對沉積物的聲學、電學、力學、熱學等物性參數(shù)產(chǎn)生不同程度的影響[10-13]。為了準確評價儲層中水合物的飽和度，在研究水合物飽和度與沉積物中聲速和衰減之間的定量關系時，必須考慮水合物微觀賦存模式對沉積物聲速及衰減的影響規(guī)律。國內外有許多學者采用室內模擬實驗方法研究了含水合物（模擬）沉積物的聲學特性[14-20]。文獻[14-17]研究了水合物的存在對不同類型沉積物（固結、松散）聲速的影響。PRIEST等[18-19]研究了不同體系（過量氣、定量氣、溶解氣）中水合物的生成對沉積物聲速的影響規(guī)律。SAHOO等[21]和LIU等[22]分別研究了聲波頻率和水合物飽和度對沉積物聲波衰減特性的影響規(guī)律。雖然研究者提出了制備具有不同水合物飽和度、不同微觀分布模式的含水合物沉積物樣品的實驗方法，但是對于所測試的整個樣品而言，實驗中所得到的飽和度、微觀分布模式是否準確、是否具有重復性、是否方便控制，仍然是值得進一步探討的問題。而數(shù)值模擬方法能夠準確、方便地控制水合物飽和度、微觀賦存模式等影響因素，同時數(shù)值模擬方法具有重復性好、時間成本低的顯著優(yōu)勢。

數(shù)值模擬方法在巖石物理領域具有廣泛的應用，逐漸形成了數(shù)字巖石物理技術。利用數(shù)字巖石物理技術不僅能夠實現(xiàn)巖石孔隙尺度建模，也能夠研究微觀因素對巖石聲速、電阻率、滲透率等物理參數(shù)的影響[23]。針對巖石聲學特性的數(shù)值模擬方法主要包括旋轉交錯網(wǎng)格有限差分方法、聲格子方法和有限元方法。SAENGER等[24]采用旋轉交錯網(wǎng)格有限差分方法模擬出聲波在各向異性、黏彈性介質中的傳播過程。DEL VALLE-GARCíA等[25]利用聲格子方法研究了聲波在多孔介質中的傳播過程和規(guī)律。MAKARYNSKA等[26]利用有限元方法研究了不同含水飽和度巖石的聲學特性。ZHU等[27]利用有限差分方法實現(xiàn)了油氣儲層介質彈性波場的數(shù)值模擬。孫建孟等[28]、劉向君等[29]采用有限元方法計算了數(shù)字巖心的彈性參數(shù)。YANG等[30]采用旋轉交錯網(wǎng)格有限差分方法模擬了地震波在海底天然氣水合物沉積地層中的傳播。卜慶濤等[31]基于實驗條件對水合物模擬儲層開展了聲學特性數(shù)值模擬，發(fā)現(xiàn)水合物儲層聲速的實驗與數(shù)值模擬結果對應關系良好。目前，基于數(shù)字巖石物理技術研究微觀賦存模式、水合物飽和度等因素對含水合物沉積物聲學特性的影響規(guī)律鮮有報道。

本文采用基于有限元的數(shù)字巖石物理技術，針對懸浮、接觸和膠結三種水合物微觀賦存模式分別建立多孔介質的三維電?力?聲多物理場耦合模型，研究了微觀賦存模式和水合物飽和度對多孔介質聲速和衰減的影響規(guī)律。研究結果可為含水合物多孔介質聲學特性數(shù)值模擬提供建模方法，為聲波測井在水合物儲層精細評價中的應用提供理論依據(jù)和模型支撐。

1 含水合物多孔介質多物理場有限元建模

1.1 幾何結構與材料

采用多物理場耦合計算平臺（COMSOL Multiphysics）來建立含水合物多孔介質的有限元數(shù)值模型。數(shù)值模型的幾何結構包含壓電晶片和多孔介質兩部分。如圖1所示，兩塊壓電晶片位于多孔介質兩側，分別用于模擬發(fā)射和接收聲波。兩塊壓電晶片的幾何尺寸和材料相同，晶片厚度為1 mm，邊長為1.624 mm。水飽和多孔介質為正方體，邊長為1.624 mm。幾何結構包含64個大小相同的球形石英砂顆粒，半徑為0.21 mm，石英砂顆粒之間規(guī)則相切排布。多孔介質的孔隙度為相互連通的孔隙體積與總體積的比值，本文所構建的多孔介質幾何結構的孔隙體積約為1.8 mm3，孔隙度為0.42。

依據(jù)ECKER[32]提出的懸浮、接觸和膠結三種典型的水合物微觀賦存模式（如圖2所示），構建含水合物多孔介質模型的三維幾何結構。圖3所示為三種水合物微觀賦存模式條件下的含水合物多孔介質的三維立體結構。含懸浮模式水合物多孔介質截面中，代表水合物的9個球體均勻分布于孔隙空間中，水合物與石英砂顆粒不接觸。含接觸模式水合物多孔介質截面中，代表水合物的9個球體與石英砂顆粒相接觸（相交）。含膠結模式水合物多孔介質截面中，代表水合物的薄層包裹在石英砂顆粒外表面，相鄰石英砂顆粒表面的水合物層相互接觸（相交），將石英砂顆粒膠結在一起。

圖1 水飽和多孔介質數(shù)值模型幾何結構

圖2 沉積物中水合物的三種微觀賦存模式（改自文獻[32]）

圖3 含不同微觀賦存模式水合物的多孔介質模型三維立體結構：（a）懸浮模式；（b）接觸模式；（c）膠結模式

為得到不同的水合物飽和度，對于懸浮和接觸模式，保持代表水合物的球體個數(shù)不變，改變球體的直徑；對于膠結模式，改變石英砂顆粒表面水合物層的厚度。受到三維模型幾何結構的限制（保持水合物為單一球形顆粒），懸浮模式條件下模型所模擬的水合物飽和度存在上限，當水合物與石英砂相切時水合物飽和度為23% ～ 24%?？紤]到地層中水合物飽和度過高時，水合物的賦存模式傾向于其他模式（如接觸模式），因此選取含2%、4%、6%、8%、10%、12%、14%、16%、18%、20%、22%等11個水合物飽和度的多孔介質進行數(shù)值求解。懸浮與接觸模式條件下，以上11種水合物飽和度對應的水合物顆粒半徑的范圍為68.47 ～ 152.14 μm。膠結模式條件下對應的水合物層厚度范圍為1.16 ～ 13.18 μm。

在分別定義石英砂、孔隙水、水合物和壓電晶片等材料的基礎上，根據(jù)表1設定多孔介質中各組分的聲速和密度等參數(shù)。壓電晶片材料為鋯鈦酸鉛（PZT-5H），密度為7 500 kg/m3。

表1 多孔介質各組分聲學特性參數(shù)[33]

1.2 物理場控制方程與邊界條件

采用靜電、固體力學、壓力聲學多物理場耦合方法實現(xiàn)數(shù)值建模方案。采用聲?結構邊界耦合壓力聲學和固體力學兩個物理場，基于壓電效應實現(xiàn)靜電和固體力學兩個物理場的耦合。

靜電物理場中控制方程為高斯方程，如式（1）所示。

式中：為電位移，C/m2；為電場強度，V/m；為電勢，V；v為自由電荷密度，C/m3。

固體力學物理場控制方程為運動方程和壓電材料本構方程，計算得到位移、應力和應變值，方程如式（2）所示。

式中：為壓電晶片材料的密度，kg/m3；為位移場，m；V為體積力，N/m3；tt為壓電晶片振動加速度，m/s2；為應力，Pa；為應變，無量綱；E與分別為壓電晶片材料的彈性矩陣（Pa）和耦合矩陣（C/m2）；0為真空介電常數(shù)，8.854 × 10?12F/m；S為壓電晶片材料的相對介電常數(shù)，無量綱。

壓力聲學物理場控制方程為時域波動方程，如式（3）所示。

式中：0為聲波傳播介質的密度，kg/m3；為介質的聲速，m/s；t為總聲壓，Pa。

以上三個物理場的邊界條件設定如下：①對于靜電物理場，將用作發(fā)射聲波的壓電晶片內端面定義為電勢邊界，電勢為矩形波調制諧波后的電壓脈沖，脈沖幅值0為100 V，頻率為；將用作接收聲波的壓電晶片內端面定義為懸浮電位邊界（floating potential）；將兩塊壓電晶片的外端面定義為接地邊界。②對于固體力學物理場，將位移場初始值設為0，將兩塊壓電晶片外端面定義為輥支承邊界，防止其沿厚度方向移動。③對于壓力聲學物理場，為保證聲波在抵達多孔介質幾何外邊界時沿法向產(chǎn)生的反射最少，將多孔介質的外邊界定義為平面輻射波邊界。

聲?結構（本模型中的多孔介質?壓電晶片）邊界條件如式（4）所示。

式中：為邊界上的單位法向量，無量綱；A為邊界上單位面積所受的力，Pa。

1.3 網(wǎng)格劃分與求解器

使用自由四面體網(wǎng)格對壓電晶片幾何結構進行網(wǎng)格劃分。采用三種數(shù)量級別的網(wǎng)格、網(wǎng)格尺寸及數(shù)量如表2所示。利用穩(wěn)態(tài)（靜態(tài)）求解器對以上描述靜電、固體力學物理場的微分方程進行求解。

表2 壓電晶片部分網(wǎng)格尺寸及數(shù)量

使用自由四面體網(wǎng)格對多孔介質幾何結構進行網(wǎng)格劃分。采用多種數(shù)量級別的網(wǎng)格、網(wǎng)格尺寸及數(shù)量如表3所示。使用與瞬態(tài)求解器（廣義方法）求解上述靜電、固體力學、壓力聲學物理場中的微分方程，選擇并行稀疏直接求解方法（multi-frontal massively parallel sparse direct solver, MUMPS）。在獲取多孔介質聲速的計算模型中，將求解時長設定為2 μs，時間步長為1 ns；在獲取多孔介質聲衰減的計算模型中，將求解時長設定為20 μs，時間步長為10 ns。使用組合參數(shù)掃描方法對表3中不同數(shù)量級別的網(wǎng)格依次進行計算。

表3 多孔介質部分網(wǎng)格尺寸及數(shù)量

1.4 網(wǎng)格選取和首波到達時刻提取方法

使用表2所示三種數(shù)量級別的網(wǎng)格測試了壓電晶片中心頻率的模型計算值對網(wǎng)格的依賴性，根據(jù)測試結果選取合適的網(wǎng)格尺寸進行后續(xù)研究。根據(jù)計算結果得到三種網(wǎng)格下壓電晶片振動位移隨激勵電信號頻率的變化關系，最大振動位移對應的頻率即為壓電晶片中心頻率。結果表明，壓電晶片的中心頻率在三種數(shù)量級別網(wǎng)格下均為642.5 kHz，由于網(wǎng)格數(shù)量成倍增加，可認為基于這三種網(wǎng)格所得到的結果均為正確值。為減少計算時間，針對壓電晶片選擇單元數(shù)量為5724的網(wǎng)格級別，即最大網(wǎng)格尺寸為200 μm，最小網(wǎng)格尺寸為40 μm。

圖4 不同網(wǎng)格條件下接收的聲波

將激勵信號頻率設定為642.5 kHz，使用表3中不同級別網(wǎng)格，計算了水飽和多孔介質接收聲波的首個波谷到達時刻，測試了模型計算值對網(wǎng)格的依賴性。圖4所示為接收聲波信號，信號首個波谷到達時刻隨網(wǎng)格數(shù)量變化情況如圖5所示。結果表明，當網(wǎng)格數(shù)量達到588 993后，首個波谷到達時刻的模型計算值已經(jīng)不隨網(wǎng)格數(shù)量的增加而變化。為了平衡計算精度和計算時間，選擇單元數(shù)量為588 993的網(wǎng)格級別，對應的最大和最小網(wǎng)格尺寸分別為60 μm和20 μm。

圖5 首個聲波波谷到達時刻隨網(wǎng)格數(shù)量的變化

從接收聲波信號中準確提取首波到達時刻是獲取多孔介質聲速的前提。將整個多孔介質幾何結構設定為同種材料（即賦予相同的聲速值），結合多孔介質幾何尺寸能夠計算出首波到達時刻理論值。為了建立首波到達時刻提取方法，首先賦予多孔介質中各組分相同的聲速和密度值，然后保持密度不變，在1 500 ～ 4 500 m/s范圍內設定一系列材料聲速值進行數(shù)值求解。經(jīng)過數(shù)值求解得到樣品的接收聲波信號如圖6所示，首先從接收信號中提取波形的首個波谷到達時刻，然后計算其與首波到達時刻理論值0的差值Δ=–0。圖7所示為Δ隨的變化情況。

可以將圖7所示Δ隨的變化趨勢劃分為三個區(qū)間，兩個轉折點對應的首個波谷到達時刻分別為1.256 μs和1.457 μs。針對三個區(qū)間內Δ隨的變化規(guī)律，分別采用線性函數(shù)進行擬合，擬合結果如式（5）所示。

采用水飽和多孔介質時，接收聲波的首個波谷到達時刻t的模型計算值為1.515 μs，結合上述線性擬合式［式（5）］，可以得到水飽和多孔介質的Δt為0.604 μs，因此其首波到達時刻的真實值t0為0.911 μs。已知多孔介質樣品沿聲波傳播方向的長度L為1.624 mm，因此能夠計算出水飽和多孔介質的聲速Vp為1 782.7 m/s。在后續(xù)的模型求解數(shù)據(jù)分析中，采用相同的方法對所獲取的首波到達時刻t進行校正，從而得到真實的首波時刻t0和多孔介質聲速。

圖7 Δt隨首個波谷到達時刻變化

2 含水合物多孔介質聲速特性

圖8所示為懸浮模式、接觸模式和膠結模式條件下含水合物多孔介質接收的聲波信號。利用所接收的聲波信號計算出首個波谷到達時刻，從而得到三種微觀賦存模式下不同水合物飽和度的多孔介質樣品的聲速。含水合物多孔介質樣品的聲速隨水合物飽和度的變化情況如圖9所示。

分析圖8和圖9可知：

（1）三種賦存模式對應的含水合物多孔介質的聲速均隨水合物飽和度增大而增大。隨著水合物飽和度的增大，更多的孔隙水轉化為水合物，由于水合物比孔隙水具有更高的體積模量，因此多孔介質整體的體積模量和剪切模量增大，從而使多孔介質具有更高的聲速。

（2）含懸浮和接觸模式水合物多孔介質的聲速隨水合物飽和度的增大而增大的規(guī)律基本一致，聲速增長速率穩(wěn)定，與DUCHKOV等[17]采用模擬實驗方法得到的聲速與飽和度之間的關系一致。當水合物飽和度為22%時，兩種賦存模式下多孔介質的聲速分別達到2 090.1 m/s和2 117.3 m/s。參見圖2和圖3，接觸模式可通過將懸浮模式下的水合物顆粒運移到與顆粒接觸而得到，因此水合物飽和度變化所引起的多孔介質彈性模量的變化差異不大。

（3）含膠結模式水合物多孔介質中水合物飽和度由0增大到2%時，聲速出現(xiàn)大幅增長現(xiàn)象，由1 782.7 m/s增長為1 956.6 m/s；水合物飽和度由2%增大到22%時，聲速隨水合物飽和度的增長速率趨于穩(wěn)定，飽和度為22%時聲速達到2 343.4 m/s。在相同水合物飽和度下，含膠結模式水合物多孔介質的聲速顯著高于含懸浮模式和接觸模式水合物多孔介質[16]。膠結模式條件下的水合物使得多孔介質的體積模量和剪切模量大幅增長，因此多孔介質聲速大幅提升。隨著水合物飽和度進一步增大，水合物飽和度變化成為聲速變化的主導因素，因此聲速增長速率趨于穩(wěn)定（與懸浮和接觸模式類似）。

圖8 接收的聲波信號（時間長度為2 μs）：（a）懸浮模式；（b）接觸模式；（c）膠結模式

圖9 多孔介質聲速隨水合物飽和度變化

采用時間平均方程[34]、伍德方程[35]、權重方程[36]、Lee改進的Biot-Gassmann（Biot-Gassmann Theory modified by Lee, BGTL）模型[37]和等效介質理論模型[32,38]等對多孔介質的聲速進行計算，計算時采用的介質聲速、密度、剪切模量和體積模量等參數(shù)參考表1。將三種微觀賦存模式對應的含水合物多孔介質聲速數(shù)值模擬結果與理論模型計算結果進行了對比，如圖10所示。其中基于等效介質理論模型A、B和C的聲速計算結果分別與含懸浮、接觸、膠結模式水合物多孔介質數(shù)值模擬結果進行了比較。

圖10 聲速數(shù)值模擬與理論模型計算結果對比：（a）懸浮模式；（b）接觸模式；（c）膠結模式

分析圖10可知：

（1）對于懸浮模式（圖10a），聲速數(shù)值模擬結果低于時間平均方程計算結果，高于伍德方程計算結果；權重方程（取1.9、取0.8）計算出的聲速與數(shù)值模擬結果具有較高的一致性；BGTL模型（取0.9、取1.4）聲速計算結果與數(shù)值模擬結果基本一致；等效介質理論模型A計算出的聲速隨水合物飽和度增長速率略低于數(shù)值模擬結果。

（2）對于接觸模式（圖10b），聲速數(shù)值模擬結果低于時間平均方程計算結果，高于伍德方程計算結果；權重方程（取1.9、取1）計算出的聲速與數(shù)值模擬結果具有較高的一致性；BGTL模型（取1、取1.5）計算結果與數(shù)值模擬結果一致性較好；等效介質理論模型B計算出的聲速隨水合物飽和度增長速率略低于數(shù)值模擬結果。

（3）對于膠結模式（圖10c），聲速數(shù)值模擬結果低于時間平均方程計算結果，高于伍德方程計算結果；權重方程（取1.7、取2.9）計算出的聲速與數(shù)值模擬結果差異最??；BGTL模型（取1、取0.6）計算結果在水合物飽和度大于10%時與數(shù)值模擬結果幾乎一致，但在小于10%飽和度區(qū)間內，BGTL模型計算結果顯著高于數(shù)值模擬結果；等效介質理論模型C計算結果與數(shù)值模擬結果高度一致。

（4）對于三種賦存模式對應的含水合物多孔介質，在所研究的水合物飽和度區(qū)間內，基于時間平均方程的多孔介質聲速預測值偏高，而伍德方程預測值偏低[39]；權重方程和BGTL模型對預測含懸浮和接觸模式水合物多孔介質聲速具有優(yōu)勢；等效介質理論模型C對預測含膠結模式多孔介質聲速具有獨特優(yōu)勢。

3 含水合物多孔介質聲衰減特性

聲波在多孔介質中傳播時，在聲阻抗不連續(xù)的界面處會產(chǎn)生折射與反射現(xiàn)象，從而導致聲波能量衰減。本節(jié)討論多孔介質中水合物飽和度和賦存模式對聲衰減的影響。經(jīng)過數(shù)值求解，得到針對含懸浮模式、接觸模式和膠結模式水合物多孔介質所接收到的聲波信號，如圖11所示。研究聲波信號的衰減時，將信號的振幅（以最大電壓表示）作為考察指標。基于信號振幅的降低定義了用于量化聲波能量衰減的參數(shù)a，a由公式（6）計算得到。

從圖11所示接收聲波信號中，提取三種賦存模式對應的含水合物多孔介質在不同水合物飽和度條件下的信號振幅，計算得到多孔介質聲波衰減參數(shù)Ma。圖12展示了含水合物多孔介質聲波衰減參數(shù)Ma隨水合物飽和度變化的情況。

圖12 聲波衰減參數(shù)Ma隨水合物飽和度變化

分析圖12可知：

（1）三種水合物賦存模式對應的含水合物多孔介質聲衰減參數(shù)a隨著水合物飽和度的增大都呈現(xiàn)出近似線性增長的規(guī)律[22]。與水飽和多孔介質相比較，水合物的存在導致聲波在傳播過程中遇到更多不連續(xù)的聲阻抗界面，從而加劇聲波能量的衰減。采用線性函數(shù)對兩者關系進行擬合，得到式（7）。

（2）含懸浮和接觸模式水合物多孔介質的聲衰減參數(shù)a隨水合物飽和度的變化規(guī)律基本一致，這與以上聲速隨水合物飽和度變化規(guī)律類似，分析其原因為懸浮和接觸模式條件下的水合物對多孔介質內部聲阻抗結構影響的差異性不顯著。

（3）在相同的水合物飽和度條件下，含膠結模式水合物多孔介質的聲衰減參數(shù)a顯著低于懸浮和接觸模式。與懸浮和接觸模式相比較，膠結模式條件下多孔介質內部組分界面出現(xiàn)差異。主要差異在于：前兩種模式中，水直接與石英砂接觸，而在膠結模式中水合物與石英砂接觸，同時水與石英砂不接觸?？紤]到水合物與石英砂之間的聲阻抗差異小于水與石英砂，因此膠結模式條件下多孔介質聲衰減參數(shù)a較低，即衰減較小。

4 結 論

基于電?力?聲多物理場耦合建模方法，針對懸浮、接觸和膠結三種水合物微觀賦存模式，分別建立了含水合物多孔介質的有限元數(shù)值模型，分析了微觀賦存模式和水合物飽和度對多孔介質聲速及衰減的影響規(guī)律，對比了聲速數(shù)值模擬與理論模型計算結果，建立了聲波衰減參數(shù)與水合物飽和度之間的關系式。得到如下結論：

（1）對于三種水合物賦存模式，由于水合物相比孔隙水具有更高的彈性模量，多孔介質的聲速隨著水合物飽和度的增大而增大；水合物的存在導致聲波在傳播過程中遇到更多不連續(xù)的聲阻抗界面，聲衰減隨著水合物飽和度的增大而近似線性增大。

（2）懸浮和接觸賦存模式條件下，水合物飽和度對多孔介質的聲速和衰減影響規(guī)律基本一致；對于相同的水合物飽和度，膠結模式條件下含水合物多孔介質具有更高的聲速和更小的聲衰減。

（3）通過合理選擇參數(shù)值，利用權重方程與BGTL模型估算的含懸浮和接觸模式水合物多孔介質的聲速較為準確；通過等效介質理論模型C計算的含膠結模式水合物多孔介質的聲速更為準確。

[1] LI X S, XU C G, ZHANG Y, et al. Investigation into gas production from natural gas hydrate: a review[J]. Applied energy, 2016, 172: 286-322. DOI: 10.1016/j.apenergy. 2016.03.101.

[2] YE J L, QIN X W, XIE W W, et al. The second natural gas hydrate production test in the South China Sea[J]. China geology, 2020, 3(2): 197-209. DOI: 10.31035/ cg2020043.

[3] LIU J W, LI X S. Recent advances on natural gas hydrate exploration and development in the South China Sea[J]. Energy & fuels, 2021, 35(9): 7528-7552. DOI: 10.1021/ acs.energyfuels.1c00494.

[4] 蘇正, 劉麗華. 南海北部陸坡天然氣水合物成藏特征研究進展分析[J]. 新能源進展, 2020, 8(1): 35-41. DOI: 10.3969/j.issn.2095-560X.2020.01.006.

[5] 周守為, 李清平, 呂鑫, 等. 天然氣水合物開發(fā)研究方向的思考與建議[J]. 中國海上油氣, 2019, 31(4): 1-8. DOI: 10.11935/j.issn.1673-1506.2019.04.001.

[6] LI J F, YE J L, QIN X W, et al. The first offshore natural gas hydrate production test in South China Sea[J]. China geology, 2018, 1(1): 5-16. DOI: 10.31035/cg2018003.

[7] KOH D Y, KANG H, LEE J W, et al. Energy-efficient natural gas hydrate production using gas exchange[J]. Applied energy, 2016, 162: 114-130. DOI: 10.1016/j. apenergy.2015.10.082.

[8] 鐘廣法, 張迪, 趙巒嘯. 大洋鉆探天然氣水合物儲層測井評價研究進展[J]. 天然氣工業(yè), 2020, 40(8): 25-44. DOI: 10.3787/j.issn.1000-0976.2020.08.002.

[9] SHANKAR U, RIEDEL M. Gas hydrate saturation in the Krishna-Godavari basin from P-wave velocity and electrical resistivity logs[J]. Marine and petroleum geology, 2011, 28(10): 1768-1778. DOI: 10.1016/j. marpetgeo.2010.09.008.

[10] KOU X, LI X S, WANG Y, et al. Effects of gas occurrence pattern on distribution and morphology characteristics of gas hydrates in porous media[J]. Energy, 2021, 226: 120401. DOI: 10.1016/j.energy.2021.120401.

[11] DAI S, SANTAMARINA J C, WAITE W F, et al. Hydrate morphology: physical properties of sands with patchy hydrate saturation[J]. Journal of geophysical research: solid earth, 2012, 117(B11): B11205. DOI: 10.1029/2012JB009667.

[12] 寧伏龍, 梁金強, 吳能友, 等. 中國天然氣水合物賦存特征[J]. 天然氣工業(yè), 2020, 40(8): 1-24. DOI: 10.3787/ j.issn.1000-0976.2020.08.001.

[13] XING L C, QI S Y, XU Y, et al. Numerical study on complex conductivity characteristics of hydrate-bearing porous media[J]. Journal of natural gas science and engineering, 2021, 95: 104145. DOI: 10.1016/j.jngse. 2021.104145.

[14] WINTERS W J, WAITE W F, MASON D H, et al. Methane gas hydrate effect on sediment acoustic and strength properties[J]. Journal of petroleum science and engineering, 2007, 56(1/3): 127-135. DOI: 10.1016/j.petrol. 2006.02.003.

[15] YUN T S, FRANCISCA F M, SANTAMARINA J C, et al. Compressional and shear wave velocities in uncemented sediment containing gas hydrate[J]. Geophysical research letters, 2005, 32(10): L10609. DOI: 10.1029/2005GL022607.

[16] HU G W, YE Y G, ZHANG J, et al. Acoustic properties of hydrate-bearing unconsolidated sediments measured by the bender element technique[J]. Chinese journal of geophysics, 2012, 55(6): 635-647. DOI: 10.1002/cjg2. 1758.

[17] DUCHKOV A D, DUCHKOV A A, PERMYAKOV M E, et al. Acoustic properties of hydrate-bearing sand samples: laboratory measurements (setup, methods, and results)[J]. Russian geology geophysics, 2017, 58(6): 727-737. DOI: 10.1016/j.rgg.2016.09.029.

[18] PRIEST J A, BEST A I, CLAYTON C R I. A laboratory investigation into the seismic velocities of methane gas hydrate-bearing sand[J]. Journal of geophysical research: solid earth, 2005, 110(B4): B04102. DOI: 10.1029/ 2004JB003259.

[19] PRIEST J A, REES E V L, CLAYTON C R I. Influence of gas hydrate morphology on the seismic velocities of sands[J]. Journal of geophysical research: solid earth, 2009, 114(B11): B11205. DOI: 10.1029/2009JB006284.

[20] REN S R, LIU Y J, LIU Y X, et al. Acoustic velocity and electrical resistance of hydrate bearing sediments[J]. Journal of petroleum science and engineering, 2010, 70(1/2): 52-56. DOI: 10.1016/j.petrol.2009.09.001.

[21] SAHOO S K, NORTH L J, MARíN-MORENO H, et al. Laboratory observations of frequency-dependent ultrasonic P-wave velocity and attenuation during methane hydrate formation in Berea sandstone[J]. Geophysical journal international, 2019, 219(1): 713-723. DOI: 10.1093/gji/ggz311.

[22] LIU Z C, NING F L, HU G W, et al. Characterization of seismic wave velocity and attenuation and interpretation of tetrahydrofuran hydrate-bearing sand using resonant column testing[J]. Marine and petroleum geology, 2020, 122: 104620. DOI: 10.1016/j.marpetgeo.2020.104620.

[23] 趙建鵬, 陳惠, 李寧, 等. 三維數(shù)字巖心技術巖石物理應用研究進展[J]. 地球物理學進展, 2020, 35(3): 1099-1108. DOI: 10.6038/pg2020DD0486.

[24] SAENGER E H, BOHLEN T. Finite-difference modeling of viscoelastic and anisotropic wave propagation using the rotated staggered grid[J]. Geophysics, 2004, 69(2): 583-591. DOI: 10.1190/1.1707078.

[25] DEL VALLE-GARCíA R, SáNCHEZ-SESMA F J. Rayleigh waves modeling using an elastic lattice model[J]. Geophysical research letters, 2003, 30(16): 1866. DOI: 10.1029/2003GL017600.

[26] MAKARYNSKA D, GUREVICH B, CIZ R, et al. Finite element modelling of the effective elastic properties of partially saturated rocks[J]. Computers & geosciences, 2008, 34(6): 647-657. DOI: 10.1016/j.cageo.2007.06.009.

[27] ZHU X, MCMECHAN G A. Numerical simulation of seismic responses of poroelastic reservoirs using Biot theory[J]. Geophysics, 1991, 56(3): 328-339. DOI: 10.1190/1.1443047.

[28] 孫建孟, 閆國亮, 姜黎明, 等. 基于數(shù)字巖心研究流體性質對裂縫性低滲透儲層彈性參數(shù)的影響規(guī)律[J]. 中國石油大學學報(自然科學版), 2014, 38(3): 39-44. DOI: 10.3969/j.issn.1673-5005.2014.03.006.

[29] 劉向君, 朱洪林, 梁利喜. 基于微CT技術的砂巖數(shù)字巖石物理實驗[J]. 地球物理學報, 2014, 57(4): 1133- 1140. DOI: 10.6038/cjg20140411.

[30] YANG J J, HE B S, ZHANG J Z. Multicomponent seismic forward modeling of gas hydrates beneath the seafloor[J]. Applied geophysics, 2014, 11(4): 418-428. DOI: 10.1007/s11770-014-0465-x.

[31] 卜慶濤, 劉圣彪, 胡高偉, 等. 含水合物沉積物聲學特性——實驗模擬與數(shù)值模擬的對比分析[J]. 海洋地質前沿, 2020, 36(9): 56-67. DOI: 10.16028/j.1009-2722. 2020.079.

[32] ECKER C. Seismic characterization of methane hydrate structures[M]. Stanford University?ProQuest Dissertations Publishing, 1998, 9837195.

[33] LEE M W, COLLETT T S. Pore- and fracture-filling gas hydrate reservoirs in the gulf of Mexico gas hydrate joint industry project leg II green canyon 955 H well[J]. Marine and petroleum geology, 2012, 34(1): 62-71. DOI: 10.1016/j.marpetgeo.2011.08.002.

[34] WYLLIE M R J, GREGORY A R, GARDNER G H F. An experimental investigation of factors affecting elastic wave velocities in porous media[J]. Geophysics, 1958, 23(3): 459-493. DOI: 10.1190/1.1438493.

[35] PEARSON C F, HALLECK P M, MCGUIRE P L, et al. Natural gas hydrate deposits: a review of in situ properties[J]. The journal of physical chemistry, 1983, 87(21): 4180-4185. DOI: 10.1021/j100244a041.

[36] WOOD W T, STOFFA P L, SHIPLEY T H. Quantitative detection of methane hydrate through high-resolution seismic velocity analysis[J]. Journal of geophysical research: solid earth, 1994, 99(B5): 9681-9695. DOI: 10.1029/94jb00238.

[37] LEE M W. Biot–Gassmann theory for velocities of gas hydrate-bearing sediments[J]. Geophysics, 2002, 67(6): 1711-1719. DOI: 10.1190/1.1527072.

[38] HELGERUD M B, DVORKIN J, NUR A, et al. Elastic-wave velocity in marine sediments with gas hydrates: Effective medium modeling[J]. Geophysical research letters, 1999, 26(13): 2021-2024. DOI: 10.1029/1999GL900421.

[39] 雷亞妮, 吳時國, 孫金, 等. 海洋天然氣水合物綜合地球物理識別方法評述[J]. 中南大學學報(自然科學版), 2022, 53(3): 864-878. DOI: 10.11817/j.issn.1672-7207. 2022.03.009.

Sound Velocity and Attenuation Characteristics of Hydrate-Bearing Porous Media Based on Numerical Model of Electrical-Mechanical-Acoustic Multi-Physics-Field Coupling

XU Yuan1, ZHANG Huan-huan1, XING Lan-chang1, WEI Wei2, HAN Wei-feng2, WEI Zhou-tuo3, DENG Shao-gui3

(1. College of Information and Control Engineering, China University of Petroleum (East China), Qingdao 266580, Shandong, China; 2. Department of Alternative Energy, PetroChina Research Institute of Petroleum Exploration & Development, Langfang 065007, Hebei, China; 3. School of Geosciences, China University of Petroleum (East China), Qingdao 266580, Shandong, China)

There are still technical difficulties in accurately controlling and evaluating the micro-distribution mode and saturation of hydrate in physical simulation experiments. Limitations exist in the experimental technologies for investigating the acoustic characteristics of hydrate-bearing sediments and establishing interpretation models of reservoir parameters. Three-dimensional numerical models were established for hydrate-bearing porous media based on digital rock physics technology. For the three kinds of hydrate micro-distribution modes (suspension, contact, and cementation), finite-element models were established individually based on the method of electrical-mechanical-acoustic multi-physics-field coupling. The effects of micro-distribution mode and hydrate saturation on sound velocity and attenuation of porous media were examined. The results of sound velocity from the numerical and theoretical models were compared. The relationship between acoustic attenuation parameters and hydrate saturation was established. It was demonstrated that: (1) for the three micro-distribution modes, the sound velocity of porous media increased with the increase of hydrate saturation due to the higher elastic modulus of hydrate compared with pore water; the existence of hydrate caused the sound wave to encounter more discontinuous acoustic impedance interfaces in the process of propagation, thus the sound attenuation increased approximately linearly with the increase of hydrate saturation; (2) under the conditions of suspension and contact modes, the influence of hydrate saturation on sound velocity and attenuation in porous media was quite similar; under the condition of cementation mode, the sound velocity was higher and the sound attenuation was lower for the porous media with the same hydrate saturation; (3) based on a reasonable selection of parameters, the sound velocity of hydrate porous media with suspension and contact modes of hydrate estimated by the weighted equation and Biot-Gassmann Theory modified by Lee (BGTL) model was more accurate; the sound velocity of porous media containing cemented hydrate calculated by the effective medium theoretical model C was more accurate. This study may provide a numerical modeling method for obtaining the acoustic characteristics of hydrate-bearing sediments under complex conditions, and provide a theoretical support for the fine evaluation of hydrate reservoirs based on acoustic logging data.

natural gas hydrate; digital rock physics; micro-distribution mode; hydrate saturation; sound velocity; sound attenuation

2095-560X（2022）05-0400-10

TK01；P631；P315

A

10.3969/j.issn.2095-560X.2022.05.002

2022-04-29

2022-07-25

中石油重大科技項目（ZD2019-184-001）；中央高?；究蒲袠I(yè)務費專項資金項目（20CX05005A）；中國石油科技創(chuàng)新基金項目（2018D-5007-0214）；山東省自然科學基金項目（ZR2019MEE095）

邢蘭昌，E-mail：xinglc@upc.edu.cn

徐 源（1998-），男，碩士研究生，主要從事天然氣水合物儲層物性參數(shù)檢測技術研究。

張歡歡（1997-），女，碩士研究生，主要從事天然氣水合物儲層物性參數(shù)檢測技術研究。

邢蘭昌（1983-），男，博士，副教授，主要從事天然氣水合物、CO2等多相流相關檢測理論與方法、多物理場耦合數(shù)值模擬方法、智能感知與檢測技術、計算機測控系統(tǒng)研究工作。