雷宇航, 尹扶, 洪鶴庭, 李娛蘭, 王寶善
中國科學(xué)技術(shù)大學(xué)地球和空間科學(xué)學(xué)院, 合肥 230026
隨著全球城市化進程的不斷加快,城市地質(zhì)環(huán)境越來越受到關(guān)注.大多數(shù)基礎(chǔ)建設(shè)和市政公共設(shè)施都位于市區(qū)地表及淺層地下空間,了解不同地區(qū)地層的巖土特性對于場地評價、監(jiān)測城市地表下沉及維護地下管網(wǎng)運行等至關(guān)重要(Mei et al., 2015).同時地下空間的探測和開發(fā)對于推動城市由外延擴張式向內(nèi)涵提升式轉(zhuǎn)變,改善城市環(huán)境,建設(shè)宜居城市具有重要意義(楊文采等,2019).然而城市的發(fā)展給城市地區(qū)地球物理數(shù)據(jù)采集帶來了諸多困難,而傳統(tǒng)密集接收器的部署及維護又需要短時大量人力物力的支撐.
近年來,分布式聲波傳感(DAS,Distributed Acoustic Sensing)技術(shù)得到迅猛發(fā)展,為實現(xiàn)便捷低成本城市地下空間探測提供了可能(王寶善等,2021).DAS通過測量光纖中激光脈沖產(chǎn)生的Rayleigh后向散射信號的相位變化來獲得地震波在光纖上引起的軸向應(yīng)變(Lindsey et al., 2020).整條光纖可以等效成等間距排列的若干個水平分量應(yīng)變計,相比于傳統(tǒng)地震儀,DAS可以實現(xiàn)低成本、多尺度、可重復(fù)的超密集地震觀測.
與傳統(tǒng)速度或加速度地震儀不同,DAS信號是地震波在光纖軸向上引起的應(yīng)變或應(yīng)變率響應(yīng)(Blum et al., 2008),應(yīng)變與地表速度或位移存在數(shù)值關(guān)系(Paitz et al., 2021; Lior et al., 2021),同樣可以反映徑向分量的波場信息.已經(jīng)在主、被動源速度結(jié)構(gòu)成像(Luo et al., 2020; 宋政宏等, 2020; 林融冰等, 2020),垂直測深剖面(Mateeva et al., 2014; Lindsey et al., 2019),接收函數(shù)(Yu et al., 2019),地震監(jiān)測及預(yù)警(Lior et al., 2021; Nayak and Ajo-Franklin, 2021)等方面得到了廣泛有效的應(yīng)用.
利用DAS進行淺地表面波結(jié)構(gòu)成像中,提取頻散曲線是其中最為關(guān)鍵的環(huán)節(jié)之一,獲得高分辨率、高信噪的頻譜圖像對于提取頻散曲線及后續(xù)壓制反演結(jié)果的多解性至關(guān)重要.近年來,大量研究表明高模式頻散曲線的加入可以顯著提高反演地層的置信度及深度(Xia et al., 2003; 羅銀河等, 2008; Pan et al., 2019),而且高階頻散對于異常夾層的響應(yīng)更為顯著(Arai and Tokimatsu, 2004; Pan et al., 2019).因為面波信號具有顯著的時頻特性,常規(guī)頻率-波數(shù)變換法(F-K)、傾斜疊加法(τ-p)及相移法等經(jīng)典Radon變換在某一頻率或時刻存在多模式相速度的混疊嚴重影響了高階面波的分辨率,而Serdyukov等(2019)對經(jīng)過S變換的擬炮記錄進行的傾斜疊加說明傳統(tǒng)的波場變換也可以獲得高分辨率的多模式相速度曲線.
Wang等(2019)首先利用頻率-貝塞爾變換(F-J)從USAarry記錄的背景噪聲中提取到了高分辨率的多模式瑞雷波頻散曲線;基于該方法提取的多模式面波反演得到了北美大陸更高分辨率的三維速度模型,展示了其在大尺度地殼、巖石圈結(jié)構(gòu)探測方面卓有成效的結(jié)果(Wu et al., 2020);該方法不僅適用于天然地震記錄中高階面波相速度的提取(Li and Chen, 2020),在淺地表結(jié)構(gòu)成像中也進行了多次實例性研究(楊振濤等, 2019; 吳華禮等, 2019; 李雪燕等, 2020).Hu等(2020)進一步推導(dǎo)出多分量瑞雷波和勒夫波F-J頻譜計算公式并進行了實例驗證,Xi等(2021)提出改進的頻率-貝塞爾變換(MF-J)以消除垂直分量F-J變換造成的交叉假頻(cross artifact)現(xiàn)象.
與傳統(tǒng)地震儀不同,DAS記錄主要反映了沿光纜方向的應(yīng)變信號.為探索如何在光纖觀測數(shù)據(jù)中提取高階面波頻散,本文推導(dǎo)了適用于DAS系統(tǒng)地震記錄中提取多模式面波頻散曲線的MF-J變換公式并應(yīng)用于實際數(shù)據(jù).希望通過本文的研究為利用城市通信光纜進行低成本高分辨率地下結(jié)構(gòu)探測奠定基礎(chǔ).
本節(jié)簡單回顧了F-J方法,并針對DAS記錄水平分量的面波信號,推導(dǎo)出改進的F-J(MF-J)頻散能量譜計算公式.在假設(shè)的水平層狀各向同性均勻介質(zhì)地層中,地表距離震源r處的地震波垂直分量記錄uz(r,t)在頻率域可以表示為:
UZ(r,ω)=F(ω)G(r,ω),
(1)
式中,ω是頻率,F(xiàn)(ω)是地震子波,G(r,ω)為震源到檢波器之間的格林函數(shù).
對(1)式波場記錄進行矢量波數(shù)變換,利用貝塞爾函數(shù)的正交性(沈永歡等, 1992),最終可以得到Z分量F-J譜計算公式(Wang et al., 2019):
(2)
DAS的地震信號反映了沿光纖軸向上的地震動,所以理論上記錄的是徑向分量(RR),因為與速度存在數(shù)值關(guān)系(Paitz et al., 2021; Lior et al., 2021),可以用來表示經(jīng)驗格林函數(shù).根據(jù)Hu等(2020)的研究,我們將DAS地震記錄的F-J譜表示為:
×J2(kr)rdr=F(ω)g(ω,k),
(3)
(4)
故式(3)可以重寫成:
(5)
(6)
式中,*表示復(fù)共軛,I(ω,k)為MF-J頻散能量譜.
(7)
I(ω,c)近似離散積分參照Wang等(2019)中的式(A7)進行計算.
根據(jù)式(7)計算得到DAS面波MF-J頻譜圖后,對能量峰值進行識別獲得瑞雷波相速度,然后進行反演以構(gòu)建地層S波速度VS模型.
由于瑞雷波頻散曲線主要受地層VS及層厚度(h)的影響(見Xia et al., 1999),故一般將縱波速度(VP)、密度(ρ)通過經(jīng)驗公式表示成VS的函數(shù)(例如Gardner et al., 1974; Broche, 2005, 2008),反演過程只優(yōu)化VS、h對頻散曲線進行擬合.本文通過遺傳算法與阻尼最小二乘法聯(lián)合反演提取的多模式瑞雷波頻散曲線(Lei et al., 2019),該算法在數(shù)值模型及實際資料處理中都表現(xiàn)出良好的穩(wěn)定性及準確性.這里我們將目標函數(shù)定義為:
(8)
式中,Φ是評價函數(shù),N是參與反演的階數(shù),M是第i階模式下實測相速度的個數(shù),cc是理論頻散曲線,co是觀測頻散曲線.
為了衡量反演結(jié)果的穩(wěn)定性,定義了群體的不確定度:
(9)
我們設(shè)計了一個含低速夾層三層水平層狀模型進行數(shù)值測試.該模型地層參數(shù)如表1所示,采用爆炸源激發(fā)信號,地震信號通過布設(shè)在地表的90個檢波器接收,道間距為5 m,最小偏移距為10 m,采樣率為500 Hz,模擬炮記錄為水平分量地震信號,記錄時長3 s.
表1 三層模型地層參數(shù)Table 1 Parameter of the horizontally layered stratum model
圖1是通過F-J變換及MF-J變換得到的水平分量面波頻譜.根據(jù)F-J變換計算的頻譜在30 Hz以上可以看到三條明顯的交叉假頻(圖1a)(Wang et al., 2019),嚴重污染到高頻高階模式相速度的識別和提取.而根據(jù)式(7)計算的MFJ頻譜消除了上述交叉假頻,真實頻散能量并沒有受此影響且與理論頻散曲線重合,證明式(7)推導(dǎo)的水平分量MF-J譜的計算是正確的.
圖1 模擬地層地震記錄計算的頻散能量譜(a) F-J頻散能量圖; (b) MF-J頻散能量圖,黑色點線是正演的理論頻散曲線.Fig.1 Dispersion energy spectrograms of simulated model seismic record(a) The F-J spectrogram; (b) The MF-J spectrogram, and the black dotted line represents the forward theoretical dispersion curve.
本次實驗場地位于北京市海淀區(qū)西北部白家疃臺地上(北京國家地球觀象臺),觀測實驗于2020年9月30日到10月8日間進行.場地地形、光纖布置及成像區(qū)域如圖2所示.鉆井資料表明實驗區(qū)下方50 m以淺主要為洪積礫巖及風(fēng)化土層,更深處為奧陶紀灰?guī)r.
圖2 白家疃實驗場地地形圖及DAS觀測系統(tǒng)藍色實線為埋置的光纖,紅色三角形為主動源炮點位置(純色三角形為本文研究成像區(qū)域),光纖拐角及示例炮點位置通過金色序號標記.Fig.2 Baijia Tuan Experiment Field topographic map and DAS arrayThe blue line is the embedded fiber, all active source positions are indicated by red triangles (the pure color triangles are located in the imaging area in this study), the golden numbers are corners and examples of source location.
本次實驗采用了與宋政宏等(2020)和林融冰等(2020)研究相同的光纖,該光纖在實驗前增加了一段大約400 m長的內(nèi)圈路徑.光纖實際有效記錄長度為1040 m,道間距設(shè)置為4 m,通道數(shù)為260,數(shù)據(jù)采樣率為500 Hz.為了避免平均應(yīng)變的測量結(jié)果出現(xiàn)較大偏差,我們實驗中所考慮的最小視波長大約為6 m,明顯大于DAS系統(tǒng)2 m的標距,故平均效應(yīng)基本可以忽略(Paitz et al., 2021).
主動源地震通過70 kg的落錘激發(fā),炮間距10 m,最小偏移距在1 m左右,本文選取了光纖南面近直線上19炮的地震數(shù)據(jù)(圖2純色三角形區(qū)域),利用提取的多模式頻散曲線反演構(gòu)建淺地表地層VS模型.落錘震源可以近似為垂直的單力源,光纖沿震源徑向分布,所以DAS記錄數(shù)據(jù)以瑞雷波為主.
與其他面波成像方法(比如多道瞬態(tài)面波分析,空間自相關(guān)及密集臺陣法等)類似,利用DAS觀測面波淺地表成像時,同樣需要經(jīng)過數(shù)據(jù)采集、頻譜成像及頻散反演三個主要環(huán)節(jié).雖然DAS光纖沿線空間間隔提前設(shè)定,真實的地理位置需要通過錘擊實驗確定.但錘擊確定的炮點位置會存在一個道間距之內(nèi)的誤差.本次實驗DAS數(shù)據(jù)是連續(xù)采集的,沒有采用觸發(fā)計時與主動源聯(lián)動.這為直接進行F-J頻譜的計算造成了困難.考慮到時間和空間的不精確性,本文設(shè)計了一套適用于DAS觀測面波淺地表成像的處理流程(圖3).
圖3 DAS觀測面波淺地表成像處理流程Fig.3 Shallow structure imaging flowchart using DAS records
主動源數(shù)據(jù)采集于10月4日進行,震源初步激發(fā)時刻通過手機秒表計時得到,所以圖3中的輸入是以上估測時間基礎(chǔ)上向前、向后共截取20 s的連續(xù)波形記錄(圖4),主動源信號通過類似STA/LTA方法(Trnkoczy, 2009)識別并截取,設(shè)置短窗時長(STA)為0.1 s,長窗時長(LTA)根據(jù)單炮最大偏移距及信噪比人為調(diào)節(jié)在1 s左右,滑動步長為0.1 s,并依此設(shè)定圖3中win為0.1 s.
圖4 DAS實驗部分連續(xù)地震記錄及截取的單炮主動源信號右下角子圖是截取窗長1s紅色框線內(nèi)的地震記錄.Fig.4 DAS experiment continuous seismic record and the single shot signal intercepted from itThe subgraph in the right is the active seismic record within the red frame of the intercept window length of 1 s.
但STA/LTA獲取的時刻t與震源激發(fā)時刻并不對應(yīng).而在頻譜分析環(huán)節(jié),如果缺少準確的發(fā)震時刻及最小偏移距信息,F(xiàn)-J變換會對不同頻率波數(shù)進行累積,導(dǎo)致計算的頻譜能量發(fā)散嚴重甚至錯誤,這與傳統(tǒng)波場變換方法僅利用不同道之間的相位及走時差的疊加不同.針對這一問題,我們在計算F-J頻譜能量圖時,假定最小偏移距為1 m,在時間上以一定步長(s·dt)整體滑動(等同于激發(fā)時刻提前或延后)炮記錄,dt是采樣間隔.當假設(shè)的最小偏移距xo與滑動炮記錄相匹配時,F(xiàn)-J譜中的相速度能量峰值會通過疊加收斂達到極值(如圖5b和圖5e).
圖5 S25炮滑動不同步長地震記錄及對應(yīng)的F-J頻譜能量圖黑色點線是預(yù)測的理論頻散曲線, (d)、(e)、(f)分別是向上滑動12、向下滑動4、向下滑動12個采樣點的地震記錄,(a)、(b)、(c)分別是對應(yīng)的F-J譜.Fig.5 Seismic records and corresponding F-J energy spectrums by sliding different step length of the shot S25The black dotted line represents the predicted theoretical dispersion curve, and (d), (e), (f) are the seismic records of sliding up 12, sliding down 4, and sliding down 12 sampling points, respectively, (a), (b), (c) are the corresponding F-J energy spectrums, respectively.
圖6是選取的其中兩個單炮地震記錄及計算的F-J頻譜能量圖(S25,S35).圖6b和圖6e(依據(jù)式(2)J0項積分)可以看出在高頻部分(20~60 Hz)至少存在3個明顯的交叉假頻,這種假頻尤其使得高階頻散受損嚴重,基于此提取的頻散曲線會影響到后續(xù)反演結(jié)果的穩(wěn)定性.而通過MF-J計算的頻譜(式(7)),交叉假頻得到顯著壓制(圖6c,6f),增加了頻譜的信噪比及分辨率.在S25記錄中,可以明顯地分辨出基階、第2、第3及第4高階頻散曲線.盡管S35炮測線更短,分辨率稍差,但仍可以明顯地分辨出基階、第2及第3高階頻散曲線.理論上MF-J變換只是去除了F-J計算頻譜中存在的“交叉”假頻而不會對信號頻散形態(tài)產(chǎn)生影響(Wang et al., 2019; Xi et al., 2021),但本文依此計算的F-J(圖6b,6c)和MF-J的頻散(圖6e,6f)存在明顯差異,關(guān)于該問題會在討論部分進行詳細分析.
圖6 S25(左側(cè))及S35(右側(cè))原始炮記錄及計算的F-J(MF-J)頻譜能量圖黑色點線是預(yù)測的理論頻散曲線,(a)和(d)是DAS主動源地震記錄,(b)和(e)、(c)和(f)分別是J0項積分的F-J譜、H0項積分的MF-J譜.Fig.6 S25 (left) and S35 (right) are seismic records and calculated F-J (MF-J) energy spectrumThe black dotted line represents the predicted theoretical dispersion curve, (a) and (d) are the DAS active seismic records, (b) and (e), (c) and (f) are the F-J spectrum of the J0 integrated and the MF-J spectrum of the H0 integrated, respectively.
與常規(guī)使用垂直分量的瑞雷波不同,水平分量的各高階模態(tài)在某些條件下會更發(fā)育(Ikeda et al.,2015; Hu et al., 2020),反映了瑞雷波垂直、水平分量能量衰減上的差異.故DAS的頻譜記錄有助于瑞雷波各高階模態(tài)的分離及提取,也為反演深部地層參數(shù)時提供了更精確約束.可以看出,不同高階頻散明顯地在不同頻段相繼占據(jù)主導(dǎo),整體上呈現(xiàn)出“波浪”式頻譜(圖6).Mi等(2018)通過數(shù)值模擬證明頻譜能量上的不連續(xù)或“跳躍”是由于低速層導(dǎo)波無法穿透地表引起的.反演該區(qū)域速度剖面(圖7,圖9)存在明顯的異常速度夾層也證實了這一觀點,而鉆孔結(jié)果也表明地下的確存在洪積礫巖及風(fēng)化土層的堆積.
圖7 S25炮初始模型及反演結(jié)果(a) 初始群體、反演群體及平均模型的橫波速度分布; (b) 最終反演地層結(jié)構(gòu)及群體基于平均模型不確定度分析.Fig.7 Shot S25 initial model and inversion results(a) The shear wave velocity distribution of the initial group, the inversion group and the average model, (b) The final inversion of the stratigraphic structure and the uncertainty analysis of the group based on the average model.
本文使用GADLS方法(Lei et al., 2019)反演頻散曲線,設(shè)置的反演初始模型包含15個薄層(前4層厚度h為3 m,后11層厚度為6 m),群體規(guī)模60,迭代次數(shù)20,每次迭代過程中均選取了殘差最小的前20個優(yōu)勢個體進行局部線性優(yōu)化以引導(dǎo)群體進化方向.
圖7是S25炮記錄的反演結(jié)果,隨機產(chǎn)生的初始模型較好地均勻分布在可能存在的速度空間.平均模型(紅色三角形)總體上表現(xiàn)出速度遞增的特征,前4層出現(xiàn)的負速度梯度除了因為本次激發(fā)震源位于硬化道路上.由于缺少高頻相速度的約束,在各炮反演結(jié)果中出現(xiàn)了一定范圍的數(shù)值波動.圖7b是依據(jù)式(9)計算的關(guān)于平均模型的不確定度分析.因為增加了5個高階相速度的約束,群體收斂結(jié)果較好,群體距平均模型偏差在4%以內(nèi).第6層出現(xiàn)了一個明顯的高速異常夾層,收斂效果也較其他層更差,雖然該層群體速度收斂效果較差,可能會大于真實地層速度,但從最后的測線剖面來看也是準確反映了高速夾層的存在位置(圖9).
基于上文同樣的反演方法及參數(shù)設(shè)定,僅對基階頻散曲線進行反演結(jié)果如圖8所示.由于缺少高模態(tài)相速度約束,反演各層速度均與多模態(tài)反演結(jié)果(圖7)出現(xiàn)較大差異,多數(shù)層不確定度分析達到6%左右,反映了基階反演群體的穩(wěn)定性較差.同時基于基階頻散反演的結(jié)果能很好地擬合觀測的基階頻散但是無法擬合高階面波的頻散特征(圖7b).所以高階模態(tài)對于壓制反演多解性,提高地下結(jié)構(gòu)成像的準確性及分辨率至關(guān)重要.
圖8 S25基階模態(tài)頻散曲線反演結(jié)果(a) 只用基階(黑色)和使用所有階(紅色)獲得的地層模型結(jié)構(gòu)及基于基階平均模型的不確定度分析(綠色); (b) 觀測的頻散能量圖,黑色及紅色點劃線分別是由圖(a)中的黑色及紅色地層模型正演的頻散曲線.Fig.8 Shot S25 inversion result of fundamental mode dispersion curve(a) The shear wave velocity model obtained by using only the fundamental (black) and using all modes (red), and the uncertainty analysis based on the fundamental average model (green); (b) F-J spectrum observed and the dispersion curves forwarded according to the fundamental (black) and multimode (red) inversion model in (a).
圖9是光纖測線南面(圖2純色三角形區(qū)域)近直線上19炮記錄反演得到的VS速度剖面,與宋政宏等(2020)和林融冰等(2020)結(jié)果大約150~260號通道之間的區(qū)域重合.由于采集時間、解調(diào)設(shè)備存在差異,以及本文增加了高模態(tài)頻散曲線約束等原因,反演結(jié)果與前人研究尤其是深部地層VS存在些許差異.因為光纖呈方形分布及炮間距較大(10 m),不同于常規(guī)多道面波分析方法,本文將單炮反演結(jié)果近似為炮點處地層結(jié)構(gòu).因而,得到的圖9地層剖面可能會沿著東西方向有壓縮效果,且坐標較真實位置偏西一些.增加了高階面波相速度約束之后,本次實驗主動源面波成像的深度可以達到80 m.結(jié)合鉆孔信息做出如下解釋:0~3 m是硬化或夯實的路基,速度接近400 m·s-1,但橫向波動較大;3~22 m為低速潛水面,巖性為淤泥質(zhì)或粉質(zhì)黏土層,表現(xiàn)出負速度梯度特征;20~30 m之間是一個水平分布的等厚高速夾層,速度大于760 m·s-1,可能是快速堆積的洪積礫巖層;30~50 m之間的速度變化較為劇烈,在300~900 m·s-1均有分布,而且左半部存在一個寬度70 m左右的碗狀凹陷,可能經(jīng)歷了不同的沉積過程.
圖9 研究區(qū)南面沿光纖測線反演地層剖面Fig.9 Inversion of stratum profile along the optical fiber in the south of the study area
圖10 炸藥激發(fā)炮記錄F-J (MF-J)頻譜(a) F-J頻譜; (b) MF-J頻譜.Fig.10 F-J (MF-J) energy spectrum of seismic record generated by explosive sources(a) F-J energy spectrum; (b) MF-J energy spectrum.
圖11 S25(左側(cè))及S35(右側(cè))MF-J頻譜能量圖黑色點線是預(yù)測的理論頻散曲線,(a)和(b)、(c)和(d)分別是H0項積分的MF-J譜及H0+H2項積分的MF-J譜.Fig.11 S25 (left) and S35 (right) MF-J energy spectrumThe black dotted line represents the predicted theoretical dispersion curve, (a) and (b), (c) and (d) are the MF-J spectrum of the H0 integrated and the MF-J spectrum of the H0+H2 integrated, respectively.
在之前Sladen等(2019)、Lior等(2021)等人利用DAS技術(shù)對水下光纜進行的研究表明不均勻的光纜-地面耦合及水土相互作用極大地影響了對地震動的敏感性,從而限制了其進行地震監(jiān)測的可靠性.與之類似,陸上通信光纜通常穿插在地下PVC (Polyvinyl chloride, 聚氯乙烯材料)管道中,對地耦合條件可能更差,而且通信光纜中間分布的光交機(通常盤留不等長度的光纜)也增大了地震監(jiān)測的不確定性(王寶善等,2021),這可能是利用通信光纜進行地震監(jiān)測研究面臨最大的非技術(shù)性挑戰(zhàn)之一,而本文實驗所用的是專門布設(shè)在地下約30 cm的光纖,對地耦合條件較好,所以不會受此困擾.
由于DAS只能探測到沿光纖軸向上的地震動信號,而當城市地質(zhì)條件或街道規(guī)劃復(fù)雜時(光纖一般沿街道鋪設(shè)),單分量的坐標旋轉(zhuǎn)、儀器響應(yīng)的去除等都需要就近布設(shè)的地震儀的配合,但也會影響到最終結(jié)果的準確性(Spica et al., 2020; Paitz et al., 2021; Lior et al., 2021).
本文基于推導(dǎo)的水平分量MF-J變換從DAS系統(tǒng)記錄的主動源地震信號中提取了高分辨率的多模式瑞雷波相速度,并用于實際資料淺地表結(jié)構(gòu)成像.基于上述研究,得出以下結(jié)論:
F-J(MF-J)變換可以從DAS系統(tǒng)的地震記錄中提取到多模式面波頻散曲線,加入的高階面波信息約束很大程度上提高了僅依靠基階頻散時的探測深度,并進一步提高了反演地層模型的穩(wěn)定性及分辨率.
在利用DAS系統(tǒng)進行主動源面波探測時,需要確定準確的震源激發(fā)時刻及最小偏移距,否則會影響F-J(MF-J)頻散能量譜的分辨率,甚至導(dǎo)致頻散曲線出現(xiàn)較大偏差.
致謝謹此祝賀陳颙先生從事地球物理教學(xué)科研工作60周年.感謝上海樸牛科技有限公司提供DAS解調(diào)設(shè)備,感謝楊微、徐善輝、胡久鵬、冀戰(zhàn)波、喬寶平等同志在實驗中給予的幫助.