方躍龍，胡高偉，劉昌嶺?，趙仕俊

（1. 中國石油大學(xué)（華東）信息與控制工程學(xué)院，山東 青島 266555；2. 國土資源部天然氣水合物重點(diǎn)實(shí)驗(yàn)室，山東 青島266071；3. 青島海洋地質(zhì)研究所，山東 青島 266071；4. 中國石油大學(xué)（華東）石油儀器儀表研究所，山東 東營 257061）

含水合物松散沉積物聲速剖面成像技術(shù)研究*

方躍龍1,2，胡高偉2,3，劉昌嶺2,3?，趙仕俊1,4

（1. 中國石油大學(xué)（華東）信息與控制工程學(xué)院，山東 青島 266555；2. 國土資源部天然氣水合物重點(diǎn)實(shí)驗(yàn)室，山東 青島266071；3. 青島海洋地質(zhì)研究所，山東 青島 266071；4. 中國石油大學(xué)（華東）石油儀器儀表研究所，山東 東營 257061）

本文基于射線理論的超聲層析成像技術(shù)，采用自適應(yīng)消噪技術(shù)除去天然氣水合物（NGH）實(shí)驗(yàn)?zāi)Ｐ凸逃性肼?，運(yùn)用直射線追蹤法和基于聯(lián)合迭代重建算法（Simultaneous Iterative Reconstruction Techniques, SIRT）的迭代重建法實(shí)現(xiàn)超聲層析成像的正演和反演，針對(duì)水合物實(shí)驗(yàn)獲取的波形，得到NGH在松散沉積物中形成時(shí)的地層聲速剖面結(jié)構(gòu)圖像。結(jié)果表明，使用超聲層析成像技術(shù)獲取的含NGH沉積物的二維聲速剖面結(jié)構(gòu)能夠準(zhǔn)確反應(yīng)不同時(shí)刻沉積物中NGH飽和度及分布情況，并且縱橫波速度剖面結(jié)構(gòu)的變化趨勢(shì)大體相同。同時(shí)本文也結(jié)合常用的NGH聲速預(yù)測(cè)模型，針對(duì)實(shí)驗(yàn)數(shù)據(jù)，對(duì)聲速與NGH飽和度之間的關(guān)系進(jìn)行了討論。

天然氣水合物；聲學(xué)速度剖面；自適應(yīng)消噪；層析成像；超聲波換能器

0 引 言

天然氣水合物（Natural gas hydrate, NGH）又稱可燃冰，具有燃燒值高、清潔無污染的特點(diǎn)，其物性分析與勘探開發(fā)技術(shù)受到人們的廣泛關(guān)注[1]。NGH與孔隙流體相比，具有較高的彈性模量，含NGH沉積層在地震和聲波測(cè)井剖面上表現(xiàn)為高速度異常[2]，因此研究NGH的聲學(xué)特性能夠反映巖性、NGH分布及含量等重要信息。在使用地震波技術(shù)進(jìn)行海底NGH的勘探中，常用到地震剖面上識(shí)別的似海底地震反射面（Bottom-Simulating Reflectors, BSR）判定NGH的分布情況。

由于受地層條件的限制，對(duì)NGH各項(xiàng)物性參數(shù)進(jìn)行原位測(cè)量十分困難，在實(shí)驗(yàn)室對(duì)人工合成NGH（含NGH沉積物）的基本物性參數(shù)進(jìn)行測(cè)量是比較經(jīng)濟(jì)有效的方法。目前，含NGH沉積物的聲學(xué)特性實(shí)驗(yàn)研究可分為靜態(tài)測(cè)試和實(shí)時(shí)測(cè)試兩種，如Winter等[3]利用GHASTLI裝置對(duì)水/冰、冰凍孔隙水、野外沉積物樣品、凍結(jié)沉積物樣品、以及不同孔隙中的沉積物合成的NGH樣品進(jìn)行了縱波的測(cè)量；Priest等[4,5]利用共振柱技術(shù)研究了“過量氣 +定量水”、“過量水 + 定量氣”體系下生成的NGH對(duì)沉積物縱橫波速和衰減的影響；業(yè)渝光等[6]和胡高偉等[7]使用縱橫波一體化探針測(cè)量了固結(jié)和松散沉積物中NGH生成過程和分解過程的聲波隨飽和度的變化情況，這些研究所用到的實(shí)驗(yàn)裝置多集中于一維模型。由于NGH在沉積物中的形成具有隨機(jī)分布特征，單一層位研究的疊加難以闡明NGH對(duì)沉積物速度剖面結(jié)構(gòu)的影響，因而需要在實(shí)驗(yàn)室內(nèi)通過模擬實(shí)驗(yàn)獲取NGH生成/分解過程中含NGH沉積層的縱橫波層速度特征和與之對(duì)應(yīng)的NGH飽和度數(shù)據(jù)。

計(jì)算機(jī)層析成像（Computed Tomography, CT）技術(shù)是指利用某種射線源，通過物體外部檢測(cè)出的信號(hào)，依照一定的物理和數(shù)學(xué)關(guān)系，生成圖像，重現(xiàn)物體內(nèi)部（剖面）特征[8]。由于超聲波具有頻率高、指向性好、能量強(qiáng)等特點(diǎn)，因此人們將超聲波作為發(fā)射源的檢測(cè)技術(shù)與射線CT技術(shù)相結(jié)合，在無損檢測(cè)領(lǐng)域具有良好的應(yīng)用。例如王振宇等[9]對(duì)已知缺陷的混凝土構(gòu)件的層析成像進(jìn)行了試驗(yàn)研究，可以識(shí)別出具有一定置信度的缺陷區(qū)域；Deidda等[10]在混凝土質(zhì)量無損檢測(cè)中采用基于SIRT法成像，在兩個(gè)相互垂直的方向上進(jìn)行透射法聲波檢測(cè)，獲取了二維聲學(xué)速度剖面，取得了良好的效果。

Netzeband等[11]在秘魯Yaquina盆地利用正演模擬得到BSR界面上方NGH穩(wěn)定帶聲速的二維圖像，并根據(jù)二維圖像對(duì)NGH穩(wěn)定帶中游離氣體的分布進(jìn)行了分析。然而在NGH模擬實(shí)驗(yàn)中尚未有這方面的研究。因此，本文將層析成像技術(shù)引入到實(shí)驗(yàn)室的NGH聲學(xué)測(cè)試中，發(fā)展了含NGH沉積物縱橫波速度剖面結(jié)構(gòu)的探測(cè)技術(shù)。通過實(shí)驗(yàn)獲得了NGH生成過程中含NGH沉積物不同層位的縱橫波速度，并利用自適應(yīng)消噪技術(shù)對(duì)噪聲信號(hào)進(jìn)行濾波。在此基礎(chǔ)上，采用超聲層析成像技術(shù)，繪制出不同時(shí)刻下的含NGH沉積物縱橫波速度二維剖面結(jié)構(gòu)。同時(shí)，利用時(shí)域反射技術(shù)（Time domain reflectometry, TDR）探測(cè)了含NGH沉積物不同層位的NGH飽和度。在此基礎(chǔ)上將兩者結(jié)合，可準(zhǔn)確獲取NGH生成/分解過程中含NGH沉積物的速度剖面結(jié)構(gòu)和其對(duì)應(yīng)層位的NGH飽和度。研究成果有望為NGH海上地球物理勘探提供一種新的方法，所得的數(shù)據(jù)處理結(jié)果能直觀地反映沉積物中含NGH飽和度及沉積物中NGH的分布情況。

1 含NGH沉積物的超聲探測(cè)實(shí)驗(yàn)

1.1 實(shí)驗(yàn)裝置描述

為獲取含NGH松散沉積物中的聲學(xué)速度剖面，本文使用含NGH沉積物速度剖面結(jié)構(gòu)特性研究的實(shí)驗(yàn)裝置（中國專利號(hào)：ZL201220232164.3）來模擬NGH的生成實(shí)驗(yàn)。實(shí)驗(yàn)裝置結(jié)構(gòu)如圖1所示，由NGH二維實(shí)驗(yàn)?zāi)Ｐ停ǜ邏悍磻?yīng)釜）、壓力控制系統(tǒng)、制冷系統(tǒng)和計(jì)算機(jī)采集系統(tǒng)組成。反應(yīng)釜內(nèi)徑為300 mm，內(nèi)部高度為400 mm，設(shè)計(jì)壓力為15 MPa，由增壓設(shè)備控制通入的甲烷氣壓力，并由一壓力傳感器測(cè)量反應(yīng)釜內(nèi)壓力的大小，壓力測(cè)量精度為±0.25%。實(shí)驗(yàn)溫度用循環(huán)水浴制冷系統(tǒng)調(diào)節(jié)。計(jì)算機(jī)采集系統(tǒng)包括聲學(xué)測(cè)試采集系統(tǒng)和TDR測(cè)試采集系統(tǒng)。

圖1 含NGH沉積物速度剖面結(jié)構(gòu)特性研究實(shí)驗(yàn)裝置簡圖Fig. 1 Experimental equipment for the study of velocity profile structure of hydrate in sediments

沿圓柱體反應(yīng)釜縱向徑面自上而下布設(shè)4對(duì)新型彎曲元換能器[13]，分別標(biāo)為AE、BF、CG、DH，具有純縱波片和彎曲元片，可單獨(dú)收發(fā)縱波信號(hào)和橫波信號(hào)。與4對(duì)聲波換能器對(duì)應(yīng)的每一層位設(shè)置雙棒型TDR探針，可分別測(cè)量4個(gè)對(duì)應(yīng)層位的含水量和NGH飽和度。8組溫度傳感器安放于釜底，利用探針的長度差分別測(cè)量沉積物4個(gè)不同層位的溫度，測(cè)量精度為 ±0.1%。反應(yīng)釜內(nèi)部結(jié)構(gòu)如圖2所示，可以看出超聲波探針?biāo)谄拭姹荛_了溫度探針，因而可認(rèn)為所布設(shè)的溫度探針不影響超聲波的測(cè)量。

圖2 反應(yīng)釜內(nèi)部結(jié)構(gòu)圖Fig. 2 Inner structure of the reactor vessel

1.2 實(shí)驗(yàn)材料及實(shí)驗(yàn)過程

本文所使用的松散沉積物為0.15～0.3 mm粒徑的砂（孔隙度為38.33%）和0.3～0.6 mm粒徑的砂（孔隙度為39.88%）。在實(shí)驗(yàn)過程中為加快NGH的生成，使用十二烷基硫酸鈉溶液（SDS溶液）作為實(shí)驗(yàn)溶液。具體實(shí)驗(yàn)步驟如下：

（1）實(shí)驗(yàn)?zāi)Ｐ椭袖佋O(shè)介質(zhì)。第1層加0.15～0.3 mm粒徑的砂并夯實(shí)，然后加入SDS水溶液直到飽和；第2層加0.3～0.6 mm粒徑的砂并夯實(shí)，然后加入SDS水溶液直到飽和；第3、4層鋪設(shè)同第1、2層。

（2）排除模型內(nèi)空氣。先向釜中通入一定壓力的甲烷氣，然后將其放出，從而將釜中的空氣排出。

（3）加壓。向釜中緩慢加入甲烷氣直到加到實(shí)驗(yàn)設(shè)定壓力為止。

（4）分步降溫來生成NGH。啟動(dòng)循環(huán)水浴制冷系統(tǒng)，先降到9.5℃左右，保持24 h后，再降到8.5℃左右，再保持24小時(shí)，以后每次降1℃，一直降到2.5℃左右為止。

整個(gè)實(shí)驗(yàn)過程中，計(jì)算機(jī)系統(tǒng)實(shí)時(shí)記錄溫度、壓力、超聲波形和TDR波形等數(shù)據(jù)。

1.3 超聲波信號(hào)采集

彎曲元換能器經(jīng)波形發(fā)生器觸發(fā)高斯脈沖信號(hào)后產(chǎn)生超聲波，聲波信號(hào)在經(jīng)功率放大器放大后由16位高速數(shù)據(jù)采集器采集，其采樣頻率為10 MHz。在實(shí)驗(yàn)過程中，縱波頻率設(shè)定為40 kHz，橫波頻率設(shè)定在20 kHz?？v波和橫波信號(hào)分兩個(gè)通道分別采集。每個(gè)換能器發(fā)射的信號(hào)可由4個(gè)換能器同時(shí)進(jìn)行接收，從而在每次測(cè)量中可獲取32個(gè)波形。為避免聲波換能器之間的干擾，以切換采集的方式，每隔一定時(shí)間分別記錄每個(gè)層位的縱橫波信號(hào)，如圖3所示，為其中一個(gè)發(fā)射換能器與4個(gè)接收換能器之間的傳播波形，上圖為采集的橫波波形，下圖為縱波波形。

圖3 采集到的縱橫波波形圖像Fig. 3 Collected images of P and S-wave

1.4 水合物的飽和度測(cè)量

利用TDR技術(shù)測(cè)量NGH的飽和度[13,14]。依據(jù)TDR波形獲取的介電常數(shù)和其與介質(zhì)含水量之間的關(guān)系式，可以計(jì)算出水合物的飽和度：

其中，Ф表示孔隙度，θv表示介質(zhì)含水量。

使用介電常數(shù)計(jì)算介質(zhì)含水量時(shí)可直接套用Wright等[14]的經(jīng)驗(yàn)公式。

利用TDR探針和上述公式，可獲取NGH生成過程中飽和度的變化情況，可以驗(yàn)證NGH飽和度與超聲波速度的關(guān)系。

2 聲波信號(hào)的自適應(yīng)消噪技術(shù)

由于聲波沿釜壁傳播的速度遠(yuǎn)大于在沉積物中的傳播速度，接收到的波形可能含有沿釜壁傳播的雜波，從而對(duì)提取首波信息帶來干擾，因此需要設(shè)計(jì)濾波器對(duì)噪聲信號(hào)進(jìn)行消除。自適應(yīng)濾波器能夠根據(jù)輸入信號(hào)統(tǒng)計(jì)特性的變化自動(dòng)調(diào)整其結(jié)構(gòu)參數(shù)，從而滿足某種最佳準(zhǔn)則（如最小均方誤差準(zhǔn)則、線性均方誤差準(zhǔn)則等）要求。因此，自適應(yīng)濾波器能夠很好地解決這一問題。

自適應(yīng)干擾抵消器是自適應(yīng)濾波器的一種應(yīng)用，其基本結(jié)構(gòu)如圖4所示。

圖4 自適應(yīng)抗干擾抵消的一般結(jié)構(gòu)Fig. 4 General structure of adaptive noise cancellation

期望信號(hào)d(n)是有用信號(hào)s(n)與干擾信號(hào)n1(n)之和，n2(n)是與n1(n)相關(guān)的一個(gè)干擾信號(hào)。由圖可看出噪聲n2(n)經(jīng)自適應(yīng)濾波器輸出y(n)，再從原始輸入中減去該輸出，產(chǎn)生了系統(tǒng)的輸出e(n)，即為對(duì)有用信號(hào)的最佳估計(jì)，干擾信號(hào)n1(n)也就得到一定程度的抵消。

最小均方（LMS）自適應(yīng)算法是一種常用的自適應(yīng)濾波算法，它是基于最小均方誤差準(zhǔn)則，在梯度法的基礎(chǔ)上，通過改進(jìn)均方差梯度的估計(jì)值的計(jì)算方法[16]。LMS算法可用下面一組遞推公式表示，即：

其中：W(n)為濾波器系數(shù)向量；X(n)是由輸入信號(hào)組成的一組輸入向量；y(n)是輸出信號(hào)；d(n)是期望信號(hào)；e(n)是誤差信號(hào)；μ為權(quán)矢量更新時(shí)的步長因子，μ越大，則算法收斂越快，但同時(shí)收斂后的誤差信號(hào)也越大；μ越小，則算法收斂速度越慢，但收斂后的誤差信號(hào)也相應(yīng)減小，穩(wěn)態(tài)性更好。為解決這一矛盾，需對(duì)步長因子加以改進(jìn)，從而提高收斂速度，對(duì)原算法加以改進(jìn)的歸一化LMS算法（NLMS算法）即可解決這一問題。

在NLMS算法中原算法步長因子μ可用下列步長因子代替：

即對(duì)公式（2）的算法做如下改進(jìn)：

其余不變。

為驗(yàn)證自適應(yīng)消噪技術(shù)的有效性，本文對(duì)實(shí)驗(yàn)裝置獲取反應(yīng)釜充滿水時(shí)聲波傳播的波形和空釜時(shí)聲波波形（即噪聲信號(hào)），運(yùn)用matlab編寫GUI圖形進(jìn)行分析驗(yàn)正。其GUI圖形主界面如圖5所示，左邊上下兩幅圖是期望信號(hào)水的波形和對(duì)應(yīng)的頻譜圖，中間兩幅圖為噪聲信號(hào)波形和對(duì)應(yīng)的頻譜圖，右邊兩幅圖是濾波后波形和對(duì)應(yīng)的頻譜圖。

利用傅里葉?小波變換法（FFT-WT）即可獲取波形的首波信息，其詳細(xì)原理和方法可以參見文獻(xiàn)[17]。根據(jù)濾波后的波形頻譜圖，可看出波形主頻在40 kHz附近；再使用以Gabor函數(shù)（中心頻率，7 MHz）作為母小波的小波分析軟件來確定聲波傳播時(shí)間：通過色譜圖的方式來顯示頻率隨時(shí)間變化的分布特征，同時(shí)分析某時(shí)刻波形頻率的小波系數(shù)，來確定此時(shí)刻的頻率分布特征。如圖6所示，可知時(shí)間點(diǎn)為200.1 μs時(shí)主頻在40 kHz處最強(qiáng)，與濾波后波形的頻譜圖結(jié)果相符，由此可推斷200.1 μs處為縱波到達(dá)時(shí)間點(diǎn)。

由于反應(yīng)釜直徑為300 mm，可知此時(shí)波速為1 499.3 m/s，基本符合聲波在水中傳播的速度（1 500 m/s）。證明本文采用的自適應(yīng)消噪技術(shù)的可行性。

胡高偉等[6]通過對(duì)南海沉積物的縱橫波測(cè)量結(jié)果的比較，證明了FFT-WT法在提取聲波的有效性。在NGH生成實(shí)驗(yàn)中，含NGH沉積層的縱橫波在經(jīng)過濾波處理后的首波信息可采用上述方法進(jìn)行提取。

圖5 自適應(yīng)消噪GUI界面Fig. 5 GUI interface of adaptive noise cancellation

圖6 水中超聲波形的小波分析圖Fig. 6 Wavelet analysis chart of ultrasonic waveform in water

3 超聲層析成像技術(shù)

超聲層析成像技術(shù)的數(shù)學(xué)基礎(chǔ)是Radon變換及逆變換。成像過程包括四個(gè)部分：數(shù)據(jù)的預(yù)處理，正演（即射線追蹤），反演，以及對(duì)反演結(jié)果的評(píng)價(jià)，其中正演和反演方程的求解是其核心內(nèi)容。本文利用新型彎曲元探針，采用自適應(yīng)消噪技術(shù)除去雜波，并結(jié)合FFT-WT法即可準(zhǔn)確獲得沉積物不同層位之間的縱橫波速度，實(shí)現(xiàn)了對(duì)數(shù)據(jù)的預(yù)處理。然后再利用層析成像技術(shù)即可得到沉積物中含NGH飽和度及沉積物中NGH的分布。

3.1 層析成像的射線追蹤法

射線追蹤是層析成像中正演和反演的關(guān)鍵步驟，它可以用矩陣形式來表示：

其中，T為投影向量，D為射線路徑矩陣，X為圖像向量。對(duì)于超聲層析成像，正演的目的在于計(jì)算超聲波從發(fā)射點(diǎn)到接受點(diǎn)的傳播路徑和走時(shí)，即構(gòu)建射線路徑矩陣D，反演的目的在于將實(shí)際的走時(shí)資料帶入反演算法中，得到速度模型，即已知投影向量T，求解公式（5）得到圖像向量X。

直射線法是構(gòu)建射線路徑矩陣中較為常用的方法，它以波在均勻介質(zhì)中傳播的路徑是直線傳播為理論基礎(chǔ)，具體算法如下。

由一個(gè)源點(diǎn)出發(fā)，計(jì)算從源點(diǎn)到其相鄰所有網(wǎng)格點(diǎn)的走時(shí)和射線長度。設(shè)源點(diǎn)坐標(biāo)A(x1,y1)，位于網(wǎng)格單元第p行與第q列，接受點(diǎn)坐標(biāo)為B(x2,y2)，位于網(wǎng)格單元的第k列與第l列。

設(shè)直射線在網(wǎng)格單元(m,n)的長度為l(m,n)

然后把除源點(diǎn)之外的所有網(wǎng)格點(diǎn)相繼當(dāng)作次級(jí)源，再計(jì)算從各個(gè)次級(jí)源點(diǎn)到其相鄰網(wǎng)格點(diǎn)的走時(shí)和射線長度，將每次計(jì)算出走時(shí)加上從波源到次級(jí)源走時(shí)，作為波源點(diǎn)到該網(wǎng)格點(diǎn)的走時(shí)[17]。對(duì)于其他發(fā)射源點(diǎn)，重復(fù)上述工作，直到選取完所有發(fā)射點(diǎn)，即可獲得每個(gè)發(fā)射點(diǎn)到其他接受點(diǎn)的射線路徑矩陣。

由于本文采用四對(duì)收發(fā)聲學(xué)探針（即圖7中的AE、BF、CG、DH），每對(duì)探針中發(fā)射探針與接受探針相距300 mm，相鄰每組探針之間的垂直距離為65 mm，所以可以構(gòu)建模型為260 mm×300 mm的矩形區(qū)域。在上面可劃分出4×4個(gè)單元格，射線條數(shù)為4×4條，每個(gè)單元格縱橫邊長分別為65 mm和75 mm（見圖7）。通過程序迭代得到所有的射線距離，再進(jìn)行格式變化可得到一個(gè)16×16的射線路徑矩陣。

圖7 層析成像模型的單元格與射線示意圖Fig. 7 Schematic diagram of the cells with the tomography model

3.2 SIRT算法

層析成像的反演過程是根據(jù)投影值重建物體內(nèi)部的分布圖像，也是實(shí)現(xiàn)Radon逆變換的過程。級(jí)數(shù)展開法是最常用的反演算法，包括反投影法和迭代類型的算法。反投影法（Back Projection technique, BPT）由于分辨率過低，一般都與迭代算法相結(jié)合，用來修改迭代算法所需要的初值。

BPT算法的基本原理可參見文獻(xiàn)[18]，利用其對(duì)迭代類算法初值的修正公式是：

其中，ai是第i條射線穿過第j個(gè)像元的距離，bi為投影向量。

較為常用的迭代類型算法有聯(lián)合迭代重建算法（SIRT），其基本原理是利用重建圖像中的某一像素所有通過它的射線的修正值來確定這個(gè)像素的修正值，以消除某些干擾因素[19]。本文所用到的SIRT算法如下。

利用BPT算法計(jì)算迭代算法所需的慢度初值xj(0)，（j=1, 2,…），再計(jì)算其估計(jì)值：

求得觀測(cè)值與估計(jì)值誤差：

再利用以下公式對(duì)第j個(gè)像素的慢度值sk進(jìn)行校正：

其中，Nj表示穿過第j個(gè)像素的射線數(shù)目，約束條件為smin＜s(k+1)＜smax。

設(shè)б為某個(gè)殘差值，當(dāng)|s(k+1)?sk|≤б時(shí)，可停止計(jì)算，否則重復(fù)反演過程。這時(shí)所得的圖像即為層析成像的結(jié)果。

4 實(shí)驗(yàn)結(jié)果分析

4.1 聲速與NGH飽和度關(guān)系的分析

根據(jù)TDR探針和彎曲元探針在反應(yīng)釜的位置，以每對(duì)探針為一層，將反應(yīng)釜從上到下分為1到4層。利用TDR技術(shù)即可得到松散沉積物的NGH飽和度。本文利用獲取的實(shí)驗(yàn)數(shù)據(jù)建立了各層松散沉積物中聲速與飽和度的對(duì)應(yīng)關(guān)系，如圖8所示。隨著NGH的生成，松散沉積物的聲速也隨之增強(qiáng)，縱橫波的增長趨勢(shì)呈對(duì)應(yīng)關(guān)系。對(duì)于不同尺寸的沉積物顆粒粒徑，第1、3層（粒徑0.3～0.6 mm）的初始縱橫波速度要大于第2、4層（粒徑0.15～0.30 mm）的速度，根據(jù)王開林等的研究，粒徑較細(xì)的沉積物縱波速度較低[20]。隨著NGH的生成，這種差異逐漸減小，聲速主要依賴于孔隙中生成的NGH。

圖8 實(shí)測(cè)各層沉積物的NGH飽和度與縱橫波聲速對(duì)應(yīng)關(guān)系圖Fig. 8 Relationships between the measured NGH saturation and velocity in various sediments

利用NGH飽和度進(jìn)行沉積物聲波速度的估算方法主要有經(jīng)驗(yàn)公式法（包括時(shí)間?平均方程、伍德方程、權(quán)重方程和BGTL理論等），等效介質(zhì)法（包含等效介質(zhì)A、B、C三種模式下對(duì)應(yīng)的公式），以及K-T方程法[21]。本文利用第四層沉積物的彈性模量和物理參數(shù)（詳見參考文獻(xiàn)[22]），繪制了利用各公式估算的沉積物聲速與實(shí)測(cè)聲速的隨飽和度的變化曲線，如圖9、圖10所示。

圖9 各模型計(jì)算的縱波速度值與實(shí)測(cè)值比較Fig. 9 Comparisons between P-wave velocity calculated by each model and the measured values

圖10 各模型計(jì)算的橫波速度值與實(shí)測(cè)值比較Fig. 10 Comparisons between S-wave velocity calculated by each model and the measured values

分析圖8、圖9、圖10可以看出：

（1）利用時(shí)間平均方程預(yù)測(cè)的縱波速度偏高，而利用伍德方程預(yù)測(cè)的縱波速度偏低；

（2）利用權(quán)重方程（W=1，n=2）估計(jì)的縱橫波速度與實(shí)測(cè)結(jié)果最為接近；

（3）利用BGTL理論（G=0.7，n=0.1）所測(cè)的縱波速度與實(shí)測(cè)值相近，而橫波速度在飽和度大于12%時(shí)與實(shí)測(cè)較為相符；

（4）對(duì)于介質(zhì)模型，模型B所預(yù)測(cè)的縱橫波速度最為接近，說明NGH在沉積物中的形成與顆粒接觸的微觀分布模式有關(guān)，但縱波速度略低于實(shí)測(cè)速度，而橫波速度在飽和度大于20%時(shí)才比較吻合；

（5）K-T方程能較準(zhǔn)確地預(yù)測(cè)縱波速度，而橫波速度預(yù)測(cè)并不準(zhǔn)確。

由于權(quán)重方程和BGTL理論具有調(diào)節(jié)因子，在含NGH松散沉積物的聲速預(yù)測(cè)上有良好的適應(yīng)性。在今后利用聲速預(yù)測(cè)飽和度的估算中應(yīng)較多考慮以上兩種模型。

4.2 聲學(xué)速度剖面結(jié)構(gòu)圖像的分析

本文選取了0 h、48 h、96 h和197 h時(shí)間點(diǎn)，利用計(jì)算機(jī)系統(tǒng)所記錄的上述時(shí)刻四對(duì)探針之間的縱橫波速度，依據(jù)前述超聲層析成像原理，分別繪制出縱橫波聲速的二維剖面圖，如圖11所示。從圖中能直觀辨別NGH未生成時(shí)、開始生成時(shí)、正在生成時(shí)和完全生成時(shí)的時(shí)間點(diǎn)和NGH在松散沉積物中的分布情況，所得的剖面圖像與Netzeband等[11]獲得的地層速度剖面圖像相似。同時(shí)利用TDR所測(cè)得的NGH飽和度數(shù)據(jù)加以說明（4個(gè)時(shí)刻不同層位的聲速、飽和度數(shù)據(jù)見表5）。

在二維速度剖面圖像中各單元格的灰度值0～255對(duì)應(yīng)的縱波波速為1.78～2.61 km/s，橫波為0.54～1.21 km/s。聲速剖面圖各灰度值表示的波速已在圖11中標(biāo)明。單元格的灰度值越高（即顏色越淺），表示波速度越快，NGH的含量也就越高。

圖11 0 h、48 h、96 h和197 h的二維縱橫波聲速剖面Fig. 11 Two dimensional P and S wave velocity profiles at 0 h, 48 h, 96 h and 197 h

結(jié)合上述時(shí)刻的聲速剖面結(jié)構(gòu)圖以及各層NGH飽和度數(shù)據(jù)進(jìn)行分析，可以看出：

（1）隨著NGH的逐漸生成，飽和度增加，剖面結(jié)構(gòu)中的各單元格的灰度值逐漸增加，表明松散沉積物中的縱橫波波速在增加；

（2）隨著NGH的增多，縱波剖面圖中灰度值較高的單元格逐漸多分布在1、2層，說明1、2層NGH含量較多，因?yàn)?、2層位于沉積層上方，更容易得到甲烷氣體，這與實(shí)際情況相吻合；

（3）實(shí)驗(yàn)結(jié)束時(shí)，第2層中部的灰度值最淺，說明此處水合物的含量最大，而結(jié)合飽和度數(shù)據(jù)，可知此時(shí)第2層的水合物飽和度最大，因此可認(rèn)為各單元格灰度值與飽和度的變化趨勢(shì)相符；

（4）由于同一時(shí)刻各個(gè)單元格灰度所代表的速度差異并不明顯，可以認(rèn)為NGH在反應(yīng)釜內(nèi)的分布較為均勻；

（5）由橫波速度剖面也可得到類似的結(jié)論，說明橫波速度剖面結(jié)構(gòu)特征的變化趨勢(shì)與縱波速度剖面基本相同；

（6）因模擬介質(zhì)的粒徑差異，其橫波速度與縱波速度在NGH生成初期差異較大，隨著NGH飽和度增加，其差異性呈減小趨勢(shì)。

（7）實(shí)驗(yàn)過程中，含NGH細(xì)粒徑松散沉積物（0.15～0.3 mm粒徑的砂，孔隙度為39.88%）的聲波速度低于粗粒徑（粒徑0.3～0.6 mm，孔隙度為38.33%），這是受沉積物孔隙度的影響，還是NGH在不同松散沉積物中的生長機(jī)制的影響，有待進(jìn)一步研究。

因此，利用本文提出超聲層析成像技術(shù)得到的縱橫波速剖面圖像與松散沉積物各層位的NGH飽和度狀況相吻合，能夠真實(shí)反映沉積物各層中的NGH含量。將超聲層析成像技術(shù)運(yùn)用于含NGH沉積物的聲學(xué)特性研究，利用該技術(shù)能使用較少的超聲波探頭得到較為準(zhǔn)確的二維聲速剖面圖像。該技術(shù)工程化后，在海底層NGH地球物理勘探中具有較強(qiáng)的適用性。

5 總結(jié)與展望

本文利用松散沉積物中NGH生成實(shí)驗(yàn)的聲波數(shù)據(jù)，通過自適應(yīng)消噪，結(jié)合FFT-WT法獲取了聲波數(shù)據(jù)的首波信息，基于射線追蹤法實(shí)現(xiàn)了層析成像的正演；再通過SIRT算法，完成層析成像的反演，從而獲取了含NGH松散沉積物的聲速剖面圖像，重建了松散沉積物內(nèi)部的NGH分布情況。結(jié)合不同層位飽和度的數(shù)值變化，該聲速剖面圖像能夠反應(yīng)沉積物中NGH的分布情況，效果良好。該技術(shù)在海底層NGH地球物理勘探中具有較強(qiáng)的適用性。

同時(shí)，本文利用實(shí)驗(yàn)所獲得的數(shù)據(jù)，結(jié)合常用的NGH聲速預(yù)測(cè)模型對(duì)NGH飽和度與縱橫波速之間的關(guān)系進(jìn)行了討論。通過對(duì)相關(guān)文獻(xiàn)資料分析，發(fā)現(xiàn)權(quán)重方程和BGTL理論在NGH飽和度預(yù)測(cè)上應(yīng)用較為廣泛。然而，受實(shí)驗(yàn)條件限制，本文僅利用現(xiàn)有文獻(xiàn)中的沉積層物理參數(shù)對(duì)NGH的飽和度進(jìn)行了討論，在今后的研究中應(yīng)對(duì)夯實(shí)后各層沉積物的物理力學(xué)參數(shù)進(jìn)行實(shí)際測(cè)量，以保證實(shí)驗(yàn)結(jié)果的準(zhǔn)確性。

由于實(shí)驗(yàn)裝置局限性，本文僅獲取了含NGH松散沉積物的二維聲速剖面，但可以在與本文所用超聲探頭所在剖面的垂直方向安裝新的超聲探頭，同時(shí)對(duì)算法加以改進(jìn)，以獲取含NGH沉積物的三維聲速剖面圖像，從而能立體且更加直觀地反映出NGH在松散沉積物中的分布情況，這也是今后的研究方向之一。

此外，相關(guān)實(shí)驗(yàn)研究發(fā)現(xiàn)在SDS溶液中形成的水合物不是很致密[23]，這對(duì)聲速測(cè)量可能會(huì)有一定影響，在今后的研究中也需考慮這方面的影響。

[1] Makogon Y F, Holditch S A, Makogon T Y. Natural gas-hydrates—A potential energy source for the 21st century[J]. Journal of petroleum science and engineering, 2007, 56: 14-31.

[2] 梁勁, 王明君, 陸敬安, 等. 南海北部神狐海域含天然氣水合物沉積層的速度特征[J]. 天然氣工業(yè), 2013, 33(7): 29-35.

[3] 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: 127-135.

[4] 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). DOI: 10.1029/2009JB006284.

[5] Best A I, Priest J A, Clayton C R I, et al. The effect of methane hydrate morphology and water saturation on seismic wave attenuation in sand under shallow sub-seafloor conditions[J]. Earth and Planetary Science Letters, 2013, 368: 78-87.

[6] 業(yè)渝光, 張劍, 胡高偉, 等. 天然氣水合物飽和度與聲學(xué)參數(shù)響應(yīng)關(guān)系的實(shí)驗(yàn)研究[J]. 地球物理學(xué)報(bào), 2008, 51(4): 1156-1164.

[7] 胡高偉, 業(yè)渝光, 張劍, 等. 南海沉積物中天然氣水合物飽和度與聲學(xué)特性的關(guān)系[J]. 石油學(xué)報(bào), 2013, 34(6): 1112-1118.

[8] 王浩全, 韓焱, 殷黎. 基于超聲CT的混凝土質(zhì)量陣列檢測(cè)方法研究[J]. 計(jì)算機(jī)工程與應(yīng)用, 2010, 46(15): 239-241.

[9] 王振宇, 劉國華. 混凝土構(gòu)件層析成像的試驗(yàn)研究[J].土木工程學(xué)報(bào), 2005, 38(6): 110-114.

[10] Deidda G P, Ranieri G. Seismic tomography imaging of an unstable embankment[J]. Engineering geology, 2005, 82(1): 32-42.

[11] Netzeband G L, Hübscher C P, Gajewski D, et al. Seismic velocities from the Yaquina forearc basin off Peru: evidence for free gas within the gas hydrate stability zone[J]. International Journal of Earth Sciences, 2005, 94: 420-432.

[12] 胡高偉, 業(yè)渝光, 張劍, 等. 基于彎曲元技術(shù)的含水合物松散沉積物聲學(xué)特性研究[J]. 地球物理學(xué)報(bào), 2012, 55(11): 3762-3773.

[13] Huisman J A, Hubbard S S, Redman J D, et al. Measuring soil water content with ground penetrating radar: a review[J]. Vadose Zone Journal, 2003, 2: 476-491.

[14] Wright J F, Nixon F M, Dallimore S R, et al. A method for direct measurement of gas hydrate amounts based on the bulk dielectric properties of laboratory test media[C]. Fourth International Conference on Gas Hydrate, Yokohama, 2002, 745-749.

[15] 陳云敏, 周燕國, 黃博. 利用彎曲元測(cè)試砂土剪切模量的國際平行試驗(yàn)[J]. 巖土工程學(xué)報(bào), 2006, 7(28): 874-880.

[16] 王海峰, 陳偉, 黃秋元, 等. 基于LMS算法自適應(yīng)噪聲抵消器的分析研究[J]. 計(jì)算機(jī)與數(shù)字工程, 2009, 37(3): 85-87.

[17] 毛偉偉, 于素萍. 走時(shí)層析成像正反演方法研究[J].西安郵電學(xué)院學(xué)報(bào), 2010, 15(1): 123-126.

[18] 許令周, 范宜仁, 夏衛(wèi)芳. ART算法中射線矩陣加權(quán)的BPT算法[J]. 計(jì)算物理, 2011, 28(5): 725-729.

[19] 劉良瓊, 劉江平, 張英德. 井間地震射線走時(shí)層析成像數(shù)值模擬[J]. 工程地球物理學(xué)報(bào), 2004, 1(5): 441-446.

[20] 王開林, 楊圣奇, 蘇承東. 不同粒徑大理巖樣聲學(xué)特性的研究[J]. 煤田地質(zhì)與勘探, 2004, 32(2): 30-32.

[21] 孫春巖, 章明昱, 牛濱華, 等. 天然氣水合物微觀模式及其速度參數(shù)估算方法研究[J]. 地學(xué)前緣, 2003, 10(1): 191-198.

[22] 胡高偉, 業(yè)渝光, 張劍, 等. 沉積物中天然氣水合物微觀分布模式及其聲學(xué)響應(yīng)特征[J]. 天然氣工業(yè), 2010, 30(3): 120-124, 140.

[23] 張琳, 王樹立, 周詩崠, 等. 表面活性劑用于促進(jìn)氣體水合物生成研究的進(jìn)展[J]. 應(yīng)用化學(xué), 2014, 31(5): 505-512.

Study of the Technology to Obtain the Sound Velocity Profile of Hydrate Bearing Unconsolidated Sediments

FANG Yue-long1,2, HU Gao-wei2,3, LIU Chang-ling2,3, ZHAO Shi-jun1,4
(1. College of Information and Control Engineering, China University of Petroleum, Shandong Qingdao 266555, China; 2. Key Laboratory of Gas Hydrate, Ministry of Land and Resources , Shandong Qingdao 266071, China; 3. Qingdao Institute of Marine Geology, Shandong Qingdao 266071, China; 4. Petroleum Instrument Research Institute, China University of Petroleum, Shandong Dongying 257061, China)

According to the ultrasonic tomography technology based on ray theory, the adaptive noise cancellation was adopted to remove the premise of autoclave under the inherent noise. The straight ray-tracing method and the iterative algorithm based on SIRT algorithm were applied to realize forward modeling and inversion of the ultrasonic tomography. Then the acoustic velocity profile distribution images of hydrate formation in loose sediments were obtained. The results showed that the two-dimensional sound velocity profile structures of hydrate bearing sediments obtained by using ultrasonic tomography technology could accurately reflect the distribution of hydrate in sediments at different time and saturation. The P-wave velocity profile structures fit in well with S-wave velocity profile structures. Besides, with the experimental data, several velocity models were used to validate the relationship between the sonic velocity and the saturation of hydrate.

natural gas hydrate; acoustic velocity profile; adaptive noise cancellation; computed tomography; ultrasonic transducer

TK01

A doi：10.3969/j.issn.2095-560X.2015.04.011

2095-560X（2015）04-0309-10

方躍龍（1989-），男，碩士，主要從事水合物探測(cè)技術(shù)實(shí)驗(yàn)研究。

2015-05-17

2015-06-29

國家自然科學(xué)基金（41104086）

? 通信作者：劉昌嶺，E-mail：qdliuchangling@163.com

劉昌嶺（1966-），男，博士，研究員，主要從事天然氣水合物研究。