海底不均勻沉積層斜向聲學原位縱波測量的仿真研究

王 瑩,陶春輝*,,張國堙,周建平,3,沈洪壘

(1.上海交通大學 海洋學院,上海 200030; 2.自然資源部海底科學重點實驗室,自然資源部第二海洋研究所,浙江 杭州 310012; 3.自然資源部海洋智能觀測技術創新中心,浙江 杭州 310012; 4.中國礦業大學 深地工程智能建造與健康運維全國重點實驗室,江蘇 徐州 221116)

0 引言

海底沉積物聲學特性通常用聲速、聲衰減系數等參數表示,可以用于反演沉積物的物理力學性質,在海洋工程地質評價、海洋聲場預報、潛在工程地質災害評估等領域有著廣泛的應用前景[1]。受沉積環境以及沉積物被改造、埋藏后壓實與固結等作用影響,從陸架、陸坡到深海盆地,不同沉積單元地層分層,導致沉積層垂向分布不均勻,其聲學特性存在著明顯的垂向變化[2]。另外,以花崗巖為基底的沉積層中常見孤石分布[3],會造成沉積層的橫向分布不均勻。針對這些不均勻沉積層開展聲學特性測量方法研究具有重要意義。

海底沉積物的聲學特性測量方法主要包括實驗室測量與原位測量兩種。實驗室測量即是在室內測量沉積物樣品,方法簡單高效,但是樣品脫離了海底原始的溫度與壓力環境,且在搬運過程中可能造成結構的變化,測量精度較低[4]。原位測量是將測量儀器放置在海底沉積物中進行測量,該方法保持了海底原始的溫度與壓力環境,且最大程度降低了對沉積物的擾動,可以獲得高精度的沉積物聲學特性[1,5-6]。根據聲源與接收換能器相對位置的不同,原位聲學測量方式可分為垂向測量、橫向測量、孔中測量以及斜向測量四大類[4](圖1)。橫向測量的發射與接收換能器在同一水平線上,由于受尺寸限制,現有的橫向聲學原位測量系統只能測所在層位一定水平范圍內的沉積物聲學特性,很難識別宏觀的沉積層空間不均勻性,代表性設備有美國的ISSAMS 系統[7]、ISSAP系統[8]、ACS系統[9]、AA沉積物聲衰減測量陣列[10],英國的SPADE探針[11]、SAPPA系統[12],我國自然資源部第一海洋研究所研發的HISAMS系統[13]和BISAMS系統[14]。垂向測量的發射與接收換能器在同一垂線上,其識別垂向不均勻性效果較好,但無法識別橫向不均勻性,代表性設備有美國的AL聲學長矛測量系統[15],我國自然資源部第二海洋研究所研制的多頻海底聲學原位測試系統MFI GeoA系統[6]和第二代MFI GeoA系統[16]?？字袦y量通常依靠鉆機進行裝置下放,可實現大深度測量,對水平方向的沉積空間不均勻性具有較好的識別效果,但是難以避免滑行波對測量結果的影響,由于探測波頻率通常較高,水平探測距離受到限制,典型代表有傳統聲波測井。斜向測量的發射與接收換能器呈一定角度傾斜,可獲取垂向和橫向信息,其識別垂向不均勻性的效果不如垂向測量好,代表性設備有美國的SAMS系統[17],該系統可以獲取橫向、縱向的沉積層信息,由于探測距離較短,且聲速計算方法過于簡單,對大深度沉積層空間不均勻性識別效果較差?？傮w而言,現有的原位聲學測量系統主要針對淺表層沉積物聲學特性測量,探測深度較淺,且受激發聲源頻帶和接收偏移距等客觀因素制約,橫向探測范圍有限。斜向測量能夠同時獲取沉積層垂向、橫向信息,是獲取不均勻沉積層聲學特性的有效方法,通過擴大探測距離并優化聲速計算方法,將有效推動大深度沉積層原位測量技術的發展。

(a)橫向測量

(b)垂向測量

(c)孔中測量

(d)斜向測量圖1 四種原位測量方式的示意圖Fig.1 Schematic diagram of four in situ measurement methods(圖片根據文獻[4]改繪。)(Figure was modified from reference[4].)

隨著海洋工程建設不斷發展,擴大探測范圍的趨勢明顯,由此帶來的潛在沉積層不均勻性問題對原位測量方式提出了新的挑戰。對此,本研究針對垂向百米和橫向十幾米甚至更遠的雙重探測范圍內極有可能出現淺層氣、孤石等地質異常體的情況,提出一種基于海底沉積層不均勻性的斜向聲學原位縱波測量方法,并利用COMSOL軟件構建沉積分層和含孤石地質異常體兩種不均勻沉積層模型,對分布式斜向聲學原位測量方法進行仿真研究,旨在解決斜向測量在獲取大范圍聲學特征時受識別地層不均勻性影響的問題。同時針對斜向測量模式對垂向不均勻性識別效果較差的特點,提出了基于等效偏移距校準的聲速計算方法,以期有效識別沉積層不均勻性。

1 斜向聲學原位縱波測量方法

斜向聲學原位縱波測量方法可同時獲取百米深度內沉積層垂向和橫向的不均勻性信息,其測量結果可通過聲速、聲衰減系數等表示。斜向聲學原位縱波測量系統包括甲板顯控單元與水下測量單元兩部分(圖2)。其中,甲板顯控單元負責數據采集與顯示,與水下測量單元通過電源線和信號傳輸線連接,可實時配置水下測量單元的工作參數,并顯示存儲采集的聲波信號。水下測量單元包括放置于海底表面的聲學發射換能器TX與搭載于探桿上的聲學接收換能器RX陣列。測量時,利用貫入裝置以預設速度v0貫入探桿,將接收換能器陣列全部剛好插入沉積物的時間標記為零時刻。繼續貫入深度Δh,發射一次信號并采集數據。每貫入Δh,重復上述操作。Δh取值越小,探測精度越高。

圖2 斜向聲學原位縱波測量系統及工作原理示意圖Fig.2 Oblique acoustic in situ longitudinal wave measurement system and schematic diagram of working principle

該系統尚在研發階段,缺乏實測數據,因此基于海底原位測量相關文獻及工程勘探資料,構建不同沉積層模型,對海底聲學原位測量時聲在不同沉積層結構中的傳播過程進行數值仿真,并根據仿真結果進行聲速反演,從而分析該方法對沉積層不均勻性的識別效果。

2 斜向聲學原位縱波測量的數值仿真

2.1 聲傳播特性仿真基礎

基于COMSOL Multiphysics壓力聲學模塊,使用有限元方法開展沉積物中聲傳播的仿真研究。利用微擾理論將聲學問題轉化為線性問題,在脫離背景屬性的情況下進行簡化分析?；谀芰?、質量和動量守恒,假設整個系統中的熱力學過程絕熱可逆,黏度和導熱系數忽略不計,材料屬性數據為常數且無熱源,將聲場描述為由波動方程控制的聲壓p,波動方程[18]如下:

(1)

2.2 沉積層模型構建

根據東海工程勘探資料構建均勻型沉積層、分層型沉積層、含孤石型沉積層三種模型。為模擬聲波在海底半無限連續空間的傳播過程,模型主體采用 50 m×60 m的矩形,兩側和底部設置厚度L0為10 m的完美匹配層(perfectly matched layer,PML),以避免模型邊界處的反射對仿真結果產生影響。沉積層中的聲波瞬態傳播過程與沉積物的聲阻抗(密度×速度)相關。

東海大陸架寬且沉積層厚[19],為使仿真更符合實際情況,搜集了東海某區域工程勘探中的140個土樣資料,其主要特征為:土質類型以淤泥、淤泥質黏土、淤泥質粉質黏土為主,少量為粉砂;干密度ρd為830～1 630 kg/m3;含水率ω為20%～80%。采用濕密度ρ作為仿真密度參數,根據公式ρ=ρd(1+ω) 換算,取值1 600～1 800 kg/m3。參考以往東海沉積物的取樣測量及原位測試測量結果,聲速范圍設定為 1 500～1 700 m/s[1,6,20-21]。模型基本參數見表1。

表1 模型基本參數表Tab.1 Model basic parameter table

設計裝置測量環境為淺海,測量時間一般控制在幾小時內,假設工作區域變化不大,水深變化不大。為更好地刻化海底原位環境對聲傳播的影響,考慮壓強對于聲速的影響,引入背景壓強p1(單位:atm,1 atm=101 325 Pa)。假設測量點位處水深為20 m,自海底面(y=0)起,p1可用如下公式表達:

(2)

式中:ρsediment為沉積物的密度,ρwater為上覆水層的密度,y對應圖2中的縱坐標。

為充分保證聲源發射的連續性,采用域聲源的方式來近似點聲源,其函數關系式定義為描述空間部分Gx,y與描述時間部分Gt的乘積形式。

Q=Gx,y·Gt

(3)

(4)

(5)

式中:Gx,y為二維平面上的高斯脈沖,Gt為雷克子波, dS表示聲源的范圍,x、y表示空間中任意點的橫、縱坐標,x0、y0表示聲源位置的橫、縱坐標,f0為發射信號主頻,t表示時間。

模型采用笛卡爾坐標系,起點坐標為(0, 0),聲源位置坐標為(35, 0);經多次試驗,dS取2e-4仿真效果較好;根據系統實際工作情況,頻率f0設為1 kHz。

2.3 網格剖分與求解

對均勻型與分層型沉積層模型形狀規則的物體,采用映射網格剖分比自由網格剖分質量更好;含孤石型沉積層模型,由于幾何形狀不規則,選用自由網格剖分效果更為理想。為保證精度,控制沉積物區域最大單元格為波長的1/10。PML區域形狀規則,對求解精度并無特殊要求,只需保證信號進入該區域后迅速衰減,因此,采用映射網格剖分,沿PML拉伸方向設置固定單元數為25。

采用向后差分法對0.04 s時間內聲波傳播過程進行仿真,時間步長設置為1×10-5s。

2.4 建?？煽啃则炞C

沿右側PML拉伸方向提取橫截線聲壓圖像(圖3),可以看出進入PML后,信號逐漸衰減,在到達壁面之前聲壓值已衰減為零,未產生反射信號。表明所采用的仿真模型能夠表征真實海底沉積層環境的半無限連續介質聲波傳播特征。

圖3 右側完美匹配層橫截線聲壓Fig.3 Right perfect match layer transversal sound pressure

3 結果與討論

3.1 三種沉積層模型仿真結果

3.1.1 均勻型沉積層模型

從均勻型沉積層模型二維平面內聲壓隨時間變化的圖像可以看出,自聲源發出脈沖信號后,以半球形不斷向外擴展(圖4a),進入PML后信號衰減劇烈(圖4b),傳播至邊界處無反射(圖4c),證明半無限連續的沉積層仿真可靠。

圖4 均勻沉積模型中的聲傳播Fig.4 Acoustic propagation in homogeneous sediments model

在對應均勻型沉積層模型中,聲源下方10、20、30、40 m的垂向距離處提取仿真數據,并繪制聲壓隨時間變化的圖像,作為接收換能器接收聲波信號的記錄(圖5)。

圖5 均勻沉積模型中聲源正下方接收器聲壓隨時間變化Fig.5 Receiver pressure directly below the acoustic source change with time in the homogeneous sediment model

3.1.2 分層型沉積層模型

在均勻型沉積層模型的基礎上構建分層型沉積層模型,上層厚20 m,下層厚40 m(圖6a)。根據勘探資料,通常淺層沉積物隨深度增加,密度和聲速呈增高趨勢。故模型上層設置為低速層,詳細材料參數取值參考沉積物Ⅰ(表1);下層為高速層,詳細材料參數取值參考沉積物Ⅱ(表1)。根據2.3節所述的網格剖分方法,采用映射網格對全部區域進行剖分,得到如圖6b所示的由規則矩形構成的疏密不一的網格。

圖6 分層型沉積層模型結構(a)及映射網格剖分(b)Fig.6 Model structure (a) and mapping grid subdivision (b) of stratified sedimentary layer model

完成瞬態求解后,可以得到分層型沉積層模型二維平面內聲壓時空分布規律,圖7分別表示在0.021、0.027、0.033 s三個時刻的瞬態聲壓分布。

圖7 分層型沉積層模型中的聲傳播Fig.7 Acoustic propagation in stratified sedimentary layer model

從圖7可以看出,區別于均勻介質,聲波傳播到低速層與高速層分界面時,一部分能量在界面處被反射(圖7a),開始反向傳播(圖7b),并在海底界面處形成二次反射(圖7c)。這表明在測量過程中,如果接收器位置恰好處于分界面附近,收到的波形信號很可能是因存在反射波干擾而失真的。

3.1.3 含孤石型沉積層模型

沉積層中出現與沉積物本身物理性質存在較大差異的地質體,會對聲傳播過程及其速度與衰減產生較大影響。為探明這種影響,構建含孤石型沉積層模型,并將其與均勻型沉積層模型仿真結果進行對比。根據東海某海域海上風電場項目工程地質勘察資料,在風電場區島嶼、暗礁等附近的殘積土、全風化和散體狀強風化巖體內多見球狀風化體(孤石)。孤石直徑一般在0.8～4.0 m。模型中以直徑4.0 m球狀中風化花崗巖為例,取密度為2 600 kg/m3,聲速為3 200 m/s[22]。孤石中心位置坐標設為(30, 20),即孤石設置在相距聲源水平距離5 m,垂向距離20 m處(圖8a)。孤石與沉積物接觸邊界的阻抗為8.32×106Pa·s·m-1。采用自由三角形網格對模型進行剖分,沉積物與孤石區域最大單元格為波長的1/10(圖8b)。

圖8 含孤石型沉積層模型結構(a)及自由網格剖分(b)Fig.8 Model structure (a) and free grid division (b) of a sedimentary layer model containing the boulder

對含孤石的沉積層模型求解后,在沿聲源發射水平距離10 m處的截線上,30、40 m深度處各取測點,分別繪制測點聲壓隨時間變化的圖像(圖9)作為接收換能器的波形信號結果?？梢钥闯?相比均勻型沉積層模型,接收到的聲壓幅值有明顯衰減,且波形也發生了明顯畸變。

圖9 含孤石型沉積層模型接收器聲壓Fig.9 Receiver sound pressure in a model of a sedimentary layer containing the boulder

圖10表示在0.015、0.017、0.033 s三個時刻的含孤石型沉積層模型的聲壓分布特征。從圖中可以看到由于沉積物中存在孤石,聲信號在傳播過程中一部分會發生透射穿過孤石(圖10a),而另一部分在孤石與沉積物接觸面發生反射(圖10b),并且回波至海底再次反射(圖10c),這也是導致信號波形畸變的主要原因。

圖10 含孤石型沉積層模型中的聲傳播Fig.10 Acoustic propagation in a sedimentary layer modle containing the boulder

3.2 聲速計算結果與討論

由于孤石等地質體與沉積層本身存在較大的物性差異,從聲學角度可將其視為沉積層中存在的物性不均勻結構體。如果測點附近存在這種物性不均勻結構體,聲學原位測量所得的聲速往往會與周圍沉積層存在顯著異常。另外,沉積層也存在垂向的分層物性分界面,可能也會使測量聲速結果異常。為了在貫入過程中實時觀測到這種異常,利用聲波的初至時刻直接計算得到地層聲速剖面。

3.2.1 傳統井下地震試驗聲速測量法

在斜向聲學原位測量方法中,聲速計算算法是決定測量精度的關鍵?，F有聲學原位測量系統貫入深度均較淺,忽略了沉積層不均勻性對聲速測量的影響。以美國華盛頓大學研制的沉積物聲學原位測量系統SAMS[17, 23]為例,其貫入深度為3 m,采用的聲速計算方法就是直接使用聲信號的傳播距離除以聲波的旅行時。對于大深度的聲學原位測量,由于未考慮沉積層的不均勻性,用該計算方法測量的結果準確度較低。地震波靜力觸探(seismic cone penetration test, SCPT)是一種大深度的地質勘探方法,它根據同一換能器兩次測量結果來進行聲速計算。假定到達h1深度發射換能器TX發射聲學信號,經過t1時間后信號到達接收換能器RX1。繼續貫入Δh后,接收換能器RX1到達RX1′位置處,再次發射聲學信號,經過t′1時間后信號到達接收換能器RX1′位置處(圖2)。美國材料與試驗學會(American Society of Testing Materials, ASTM)在井下地震試驗SCPT的標準試驗方法中使用的聲速計算公式[24]如下:

(6)

式中:cⅠ表示根據井下地震試驗SCPT的標準試驗方法計算得到的聲速,ΔL1是RX1下降前后兩次的測量過程中聲信號傳播直線距離之差,Δt1是前后兩次的測量聲波旅行時之差,L′1是后一次測量時接收換能器與發射換能器之間的直線距離,L1為前一次測量時接收換能器與發射換能器之間的直線距離,t′1為后一次信號發出到接收器接收所用時間,t1為前一次信號發出到接收器接收所用時間,x為聲源水平偏移距。

SCPT的聲速算法可應用于海底不均勻性沉積層的大深度測量中:根據沉積層模型的仿真結果,在不同聲源水平偏移距下,按測點間距Δh=1 m提取數據;參考李倩宇 等[25]提出的聲學原位信號自動拾取方法,采用長短時窗均值比(short term average/long term average, STA/LTA)結合峰值校正方法拾取聲波初至時刻;用公式(6)計算得到測點處聲速,并繪制聲速剖面圖像(圖11)。對于分層型沉積層模型:在x=0 m 的垂向測量時,公式(6)的計算結果與假定值基本一致;但在x=20 m的斜向測量時,在聲速假定值為1 700 m/s的沉積層中,聲速計算結果最高值為 1 934.4 m/s。由公式(6)計算得到的聲速結果在遇到沉積層分界面時異常偏高。對于含孤石型沉積層模型,在距離聲源0、10、20 m測量時均能明顯觀察到由于孤石存在導致的聲速曲線異常情況(圖11d～11f),分別在27、22、33 m處聲速開始出現異常。斜向測量聲速結果偏離假定值更大,更易識別。上述結果表明,斜向聲學原位縱波測量方法能識別到位于發射和接收范圍內的孤石,相比垂向測量在識別水平不均勻性方面具有明顯優勢。

圖11 不同模型在不同水平偏移距時的聲速剖面Fig.11 Sound velocity profiles of different models at different horizontal offsets

對比測點間距Δh=1 m時,不同聲源水平偏移距下根據分層型沉積層模型的仿真數據得到的聲速計算結果,隨著聲源水平偏移距增大,按照公式(6)計算得到的聲速結果在遇到沉積層分界面時的異常偏高愈明顯。原因是該公式忽略了復雜地質情況導致的聲線彎曲,將從聲源發射到接收點處的聲傳播路徑近似為直線。以本文提出的分層模型為例,聲源距接收器存在一定水平偏移距,根據斯涅爾定律,聲波從聲速小的介質進入聲速大的介質時,折射角大于入射角。此時,真實的聲傳播路徑與公式(6)表示的路徑長度不符,故依據有限元仿真結果代入公式(6)得到的計算聲速與假定聲速存在較大偏差。

3.2.2 基于水平偏移距校準的原位聲速測量法

在沉積層情況未知的條件下,為獲取更貼近真實情況的實時聲速剖面,在公式(6)的基礎上進行改進,引入權重系數α,以降低聲源水平偏移距增大加劇聲線彎曲帶來的測量結果異常突變效應。優化后的聲速計算公式如下:

(7)

式中:cⅡ表示優化后的聲速計算公式計算得到的聲速;α是權重系數,取值范圍為0～1,與聲源水平偏移距大小有關。通過多次仿真發現,聲源水平偏移距離為20 m,α取值為0.5時,聲速計算效果最優。

為進一步明確原位聲速測量的影響因素,通過調整每次下降的深度,得到不同測點間距Δh下的聲速剖面結構,并驗證公式(7)的計算效果(圖12)。對于分層型沉積層模型,在聲源水平偏移距x=20 m,Δh為2、1、0.5和0.2 m時,根據公式(7)計算的聲速范圍分別為1 494.8～1 722.9、1 480.8～1 720.8、1 430.6～1 767.3 和1 332.9～1 925.0 m/s,公式(7)的計算結果比公式(6)更貼近于假定沉積層聲速結構(圖12a～12d)。結果表明,不同Δh下,分層界面位置也隨之變化,Δh越大,沉積層分界面識別滯后的效應越明顯。

圖12 不同模型在不同測點間距時的聲速剖面Fig.12 Acoustic velocity profiles of different models at different measuring point spacing

上述規律不僅適用于分層型沉積層模型,對含孤石型沉積層模型同樣適用。對于在距聲源水平距離5 m、垂向距離20 m處存在孤石異常體的沉積層模型,當聲源水平偏移距x=20 m,Δh為2、1、0.5和0.2 m時,根據公式(7)計算的聲速范圍分別為 1 491.2～1 665.9、1 465.9～1 690.1、1 455.7～1 696.5 和 1 279.0～1 843.7 m/s(圖12e～12h)。聲速在約32 m深度處開始出現異常升高,這與孤石所在20 m深度相比是滯后的。說明斜向原位縱波測量方法能識別沉積層水平方向上的不均勻性,但當異常地質情況與測點水平距離較大時,識別結果呈現一定的滯后性,后期可以通過數據處理反演孤石真實層位。

兩種模型在不同Δh下的聲速計算結果均表明:隨著Δh的減小,計算聲速數據波動加劇,計算結果誤差明顯增大。推測這是由于Δh越小,聲信號到達接收點處的時間差越小,精度要求過高導致聲速計算結果誤差增大。這說明在實際工程測量中,若時間精度不能保證,加密測量反而可能會放大誤差。為優化上述結果,保證空間采樣率,同時降低計算誤差,自上而下每隔Δh依次取間隔1 m的兩個測點數據,根據公式(7)計算得到相應聲速結果(圖13):分層型沉積層模型測點間距為0.5 和0.2 m,計算聲速范圍分別為 1 462.5～1 734.5 和1 456.6～1 770.6 m/s;含孤石型沉積層模型測點間距為0.5 和0.2 m,計算聲速范圍分別為1 469.0～1 679.5 和1 461.5～1 693.8 m/s。與優化前相比,計算誤差明顯降低,曲線更平滑。將圖 13 與圖 12b和圖 12f對比,發現測點間距小的聲速結構更精細,這是由于在相同長度上的測點間距越小,測點數量越多,降低了數據隨機誤差。

圖13 優化后的不同模型在不同測點間距時的聲速剖面Fig.13 Optimized acoustic velocity profiles of different models at different measuring point spacing

4 結論與展望

本文提出一種針對海底沉積不均勻性的斜向聲學原位縱波測量方法,通過仿真海底沉積層聲波傳播的瞬態過程,對分層和含孤石兩類典型的沉積層不均勻性結構進行仿真,探究發射換能器與接收換能器水平距離、測點間距對于聲學原位測量結果的影響,結論如下。

1)沉積層空間結構的不均勻性導致信號異常衰減甚至波形變化,進而對聲學特性測量結果產生影響。測點間距Δh越小,獲取的沉積層聲速剖面結構越精細,但時間精度不夠會導致誤差放大,增大計算所用測點間隔并對每個測點的結果進行重新計算,可有效減少由于時間測量精度不夠帶來的誤差。

2)斜向測量聲源水平偏移距越大,沉積層垂向不均勻性造成的聲線彎曲越嚴重,增大了準確獲取沉積層聲速剖面的難度,改進后的聲速計算公式在一定程度上彌補了斜向測量在準確識別垂向不均勻性方面的不足,后續可考慮進行迭代反演獲取沉積層真實聲速結構。

3)斜向原位聲學探測方法對于沉積空間不均勻性具有良好的識別效果,可廣泛應用于海洋工程勘察,服務于海洋工程基礎設計,有望在天然氣水合物探測及大洋礦產資源調查等領域實現拓展應用。