榮棉水, 王繼鑫, 李小軍，*, 劉奧懿, 孔小山, 李航

1 北京工業(yè)大學(xué)城市與工程安全減災(zāi)教育部重點(diǎn)實(shí)驗(yàn)室，北京 100124 2 防災(zāi)科技學(xué)院河北省地震災(zāi)害防御與風(fēng)險(xiǎn)評(píng)價(jià)重點(diǎn)實(shí)驗(yàn)室，河北三河 065201

0 引言

大量震害和理論研究表明，區(qū)域近地表土層速度結(jié)構(gòu)對(duì)地震動(dòng)特性有顯著的影響，主要表現(xiàn)為對(duì)地震動(dòng)峰值和地震波中不同頻率成分的放大或縮小，從而直接影響到地震災(zāi)害的嚴(yán)重程度(師黎靜, 2007; 李小軍等, 2020)，準(zhǔn)確確定近地表土層速度結(jié)構(gòu)對(duì)防震減災(zāi)極為重要(馮少孔, 2003).

由于鉆孔、P-S波測(cè)井和地震勘探方法存在高成本、破壞環(huán)境、有地形時(shí)難開展等缺點(diǎn)，近幾十年來，利用地表背景噪聲觀測(cè)反演土層速度結(jié)構(gòu)的間接方法得到迅速發(fā)展.如歐洲自2001年起就設(shè)立了Site Effects Assessment using Ambient Excitations(簡(jiǎn)化為SESAME)計(jì)劃(Bard et al., 2004)，致力于發(fā)展基于背景噪聲的場(chǎng)地效應(yīng)估計(jì)和沉積層速度結(jié)構(gòu)反演方法.美、日、墨西哥等許多國家學(xué)者都相繼研發(fā)操作簡(jiǎn)便、成本相對(duì)較低的背景噪聲反演方法(Cho, 2019; Pia-Flores et al., 2017; Thomas et al., 2020).我國也有許多地震學(xué)家將噪聲密集臺(tái)陣方法應(yīng)用于淺層速度結(jié)構(gòu)的確定，取得了很好的效果(Li et al., 2020; 雷宇航等, 2021; 王寶善等, 2021).綜合來看，具體到土層的速度結(jié)構(gòu)的背景噪聲觀測(cè)記錄反演，主要發(fā)展有三類方法：①用單點(diǎn)背景噪聲的水平與豎向譜比法(均簡(jiǎn)稱NHV)測(cè)定的卓越周期簡(jiǎn)單推算土層模型參數(shù)(王未來等, 2011)；②在觀測(cè)區(qū)架設(shè)小孔徑地震臺(tái)陣，通過測(cè)定背景噪聲面波相速度頻散曲線反演場(chǎng)地下方隨深度變化的速度結(jié)構(gòu)(Tian et al., 2020; 王繼鑫等, 2020)；③基于面波理論(van Ginkel et al., 2020)或散射場(chǎng)理論(Sánchez-Sesma et al., 2011)建立NHV與地下土層結(jié)構(gòu)的聯(lián)系，以單點(diǎn)或多點(diǎn)臺(tái)陣方式反演土層速度結(jié)構(gòu)以及在此基礎(chǔ)上發(fā)展的NHV與面波的聯(lián)合反演.第一類方法是根據(jù)卓越周期Tp和場(chǎng)地覆蓋土層厚度h來估算地基上層的平均剪切波速度(VS=4h/Tp)，僅能反映場(chǎng)地土層的總體平均特性(Arai and Tokimatsu, 2004)，第二類方法在實(shí)際應(yīng)用時(shí)需要將背景噪聲探測(cè)儀器布設(shè)成一定尺度的臺(tái)陣.相比較而言，第三類方法既沒有臺(tái)陣布設(shè)的限制，又能方便地與面波頻散形成聯(lián)合反演，應(yīng)用更為靈活.在第三類方法中，單點(diǎn)NHV反演方法是多點(diǎn)反演的基礎(chǔ)，因而成為當(dāng)前研究的重點(diǎn).

根據(jù)對(duì)引起地表NHV的噪聲波場(chǎng)認(rèn)識(shí)的差異，單點(diǎn)NHV反演方法大體分為三類：①認(rèn)為NHV主要受背景噪聲中體波成分影響，如Herak(2008)、Bignardi等(2016)用垂直入射的S波與P波傳遞函數(shù)的比值來計(jì)算NHV并據(jù)此發(fā)展了利用背景噪聲反演淺層速度結(jié)構(gòu)的ModelHVSR和OpenHVSR方法；②認(rèn)為NHV主要受背景噪聲中面波成分影響，如Arai和Tokimatsu(2004)在面波理論基礎(chǔ)上發(fā)展的背景噪聲反演方法，張立(2009)發(fā)展的基階瑞利面波NHV反演剪切波速結(jié)構(gòu)的遺傳算法等；③認(rèn)為NHV是包括體、面波兩種不同類型波場(chǎng)共同作用的結(jié)果.Sánchez-Sesma等(2011)的研究表明，水平層狀場(chǎng)地背景噪聲波場(chǎng)可用散射方程來描述，該方程借助波動(dòng)傳播的平均能量密度建立了NHV與格林函數(shù)的關(guān)系式，首次將背景噪聲波場(chǎng)考慮為含有體波、面波的共同作用，從理論上建立了利用基于背景噪聲的地球物理特征(此處為NHV)反演場(chǎng)地土層模型參數(shù)的基礎(chǔ).隨后García-Jerez等(2013)利用邊界線積分對(duì)NHV的正演計(jì)算進(jìn)行了改進(jìn)，并結(jié)合蒙特卡洛或模擬退火建立了NHV反演方法.

從以上研究現(xiàn)狀可知，當(dāng)用不同的正演理論來解釋噪聲波場(chǎng)時(shí)，所發(fā)展的NHV反演方法也相應(yīng)地存在本質(zhì)上的差異.大量觀測(cè)和理論研究均表明背景噪聲波場(chǎng)是面波和體波共同作用的結(jié)果，基于散射場(chǎng)理論的NHV更全面考慮了不同波動(dòng)類型，顯然更符合觀測(cè)實(shí)際，因此再回顧前述的三類背景噪聲反演方法，現(xiàn)有基于散射場(chǎng)理論的NHV反演方法及實(shí)用的反演工具的開發(fā)就顯得非常匱乏.已有的基于散射場(chǎng)理論NHV的反演方法如García-Jerez等(2016)發(fā)展的HV-Inv，其反演算法采用蒙特卡洛或模擬退火全局反演，仍存在計(jì)算效率不高、不利于快速求解等突出的問題，因此在考慮不同類型波場(chǎng)共同作用的前提下，發(fā)展快速有效的全局NHV反演方法是當(dāng)前亟需開展的一項(xiàng)重要工作.

針對(duì)這一研究現(xiàn)狀，本文基于散射場(chǎng)理論NHV的正演計(jì)算方法，引入遺傳與模擬退火相結(jié)合的反演框架，提出了一種場(chǎng)地土層速度結(jié)構(gòu)的背景噪聲全局反演方法，利用多個(gè)算例驗(yàn)證了方法的合理性和適用性，并在唐山響嘡豎向觀測(cè)臺(tái)站場(chǎng)地開展了本文方法的實(shí)例應(yīng)用.

1 NHV全局反演方法

1.1 土層速度結(jié)構(gòu)的NHV正演方法

NHV是場(chǎng)地地表觀測(cè)到的背景噪聲的水平和豎向分量的比值，前文已述，基于對(duì)背景噪聲波場(chǎng)的不同認(rèn)識(shí)，NHV的理論計(jì)算可主要分為體波垂直入射解(Herak, 2008; Bignardi et al., 2016)、隨機(jī)點(diǎn)源面波解(Arai and Tokimatsu, 2000)、散射場(chǎng)理論解(Sánchez-Sesma et al., 2011)三類.如何解釋不同類型波場(chǎng)作用下的NHV一直是許多研究者尋求解決的難點(diǎn)問題，僅考慮單獨(dú)某一種類型的波場(chǎng)(如體波或者面波)來解釋NHV仍存在不足之處，散射場(chǎng)理論解全面考慮了背景噪聲波場(chǎng)中體波和面波的共同作用，理論和觀測(cè)結(jié)果均證實(shí)了用其解釋NHV的合理性和適用性(Sánchez-Sesma et al., 2011)，本文采用該理論解作為NHV正演方法，對(duì)其簡(jiǎn)要介紹如下：

散射場(chǎng)理論解是隨著背景噪聲成像的地震干涉理論發(fā)展起來的(秦彤威等, 2021)，地震干涉理論認(rèn)為，完全散射波場(chǎng)中任意兩個(gè)觀測(cè)點(diǎn)獲得的位移記錄的互相關(guān)正比于兩觀測(cè)點(diǎn)之間的格林函數(shù)(Shapiro et al., 2005)，在頻率域中可以表示為以下公式(Sánchez-Sesma and Campillo, 2006; Perton et al., 2009; Sánchez-Sesma et al., 2011)：

(1)

其中，格林函數(shù)Gij(xA,xB,ω)表示在xB點(diǎn)j方向作用的一個(gè)單元簡(jiǎn)諧點(diǎn)力δijδ(|x-xB|)exp(iωt)在xA點(diǎn)i方向引起的位移.在單元簡(jiǎn)諧點(diǎn)力δi jδ(|x-xB|)exp(iωt)中，i為虛數(shù)單位，ω為圓頻率，t為時(shí)間.k=ω/β為S波波數(shù)，β為S波波速，Es=ρω2S2為平均S波能流密度，S2為S波平均譜密度，星號(hào)表示復(fù)共軛，Im表示給定的函數(shù)的虛部，尖括號(hào)表示平均.假定xA=xB，上式可以寫為

=-2πμEsk-1×Im[Gm m(xA,xA,ω)],

(2)

其中，μ=ρβ2，E(xA)為能量密度，下標(biāo)m表示不同方向.該公式表明當(dāng)激發(fā)力源和觀測(cè)接收點(diǎn)重合時(shí)，能量密度與格林函數(shù)虛部成正比.當(dāng)取x=0表示地表觀測(cè)點(diǎn)，下標(biāo)用1和2表示兩個(gè)水平方向，3表示豎直方向，則水平和豎向譜比可寫為

(3)

上式中格林函數(shù)虛部是與各土層的VS、VP、h、ρ相關(guān)的函數(shù)，由于其函數(shù)形式既涉及各層VS、VP、h、ρ組成的矩陣求解，又涉及同一土層的波數(shù)積分，函數(shù)形式較為復(fù)雜，對(duì)此，文獻(xiàn)Sánchez-Sesma等(2011)附錄中給出了分層場(chǎng)地格林函數(shù)虛部計(jì)算的詳細(xì)推導(dǎo)過程和計(jì)算表達(dá)式，篇幅所限本文不再贅述.但無論格林函數(shù)虛部的計(jì)算表達(dá)式如何復(fù)雜，只要各層土體特性參數(shù)確定，則格林函數(shù)虛部也唯一地確定，這樣上式就將地表平均測(cè)量值與地下介質(zhì)的固有特性聯(lián)系起來，本文即采用該式正演計(jì)算已知場(chǎng)地土層模型的理論NHV.

1.2 目標(biāo)函數(shù)的確定

反演的過程就是以地表觀測(cè)NHV曲線為目標(biāo)提取地下一維速度結(jié)構(gòu)的過程，為此需要確定一目標(biāo)函數(shù)來描述反演模型的理論NHV與實(shí)際觀測(cè)的吻合程度.作為對(duì)已有目標(biāo)函數(shù)的改進(jìn)，我們引入了斜率相關(guān)項(xiàng)并構(gòu)建了一種新形式的目標(biāo)函數(shù)(Rong et al., 2020)，其定義如下：

該式第一項(xiàng)代表模擬曲線和觀測(cè)曲線幅值和譜形的差異，第二項(xiàng)表示模擬曲線和觀測(cè)曲線之間斜率比的差異.采用兩項(xiàng)乘積形式的目標(biāo)函數(shù)將使得反演曲線需要同時(shí)擬合NHV曲線及其一階導(dǎo)數(shù)曲線，從而實(shí)現(xiàn)多目標(biāo)反演，這種多目標(biāo)反演有利于結(jié)果的快速收斂，同時(shí)可以減少反演的多解性.在等式(4)中，NHVt(f)為理論NHV，NHVo(f)為觀測(cè)NHV，f為頻率，dNHVt(f)/df和dNHVo(f)/df分別代表NHVt(f)、NHVo(f)對(duì)頻率f的一階導(dǎo)數(shù).P=(p1,p2,p3,…,pN)T是模型的參數(shù)向量，pi={VSi,VPi,hi,ρi}是第i層的參數(shù).VSi、VPi、hi、ρi分別代表橫波速度、縱波速度、厚度和密度，N表示土層總數(shù).當(dāng)N較大時(shí)，Φ(P)為多參數(shù)、多極值的非線性函數(shù).

1.3 遺傳模擬退火混合反演算法

遺傳算法(GA)和模擬退火算法(SA)混合的全局優(yōu)化算法已被用于土層模型參數(shù)的地震動(dòng)觀測(cè)記錄反演(榮棉水等, 2018)，這種混合算法利用了模擬退火算法局部搜索能力強(qiáng)的優(yōu)勢(shì)避免遺傳算法陷入局部極值而過早收斂，同時(shí)也利用遺傳算法全局搜索能力強(qiáng)的優(yōu)勢(shì)有效克服了模擬退火算法全局搜索能力差的缺點(diǎn)，二者取長補(bǔ)短提高了反演的效率與準(zhǔn)確性.為適應(yīng)背景噪聲數(shù)據(jù)的反演，本文對(duì)該混合方法進(jìn)行了兩方面的調(diào)整：一是不再以目標(biāo)函數(shù)的精度作為反演是否結(jié)束的判斷條件，而是改用預(yù)先設(shè)定的遺傳迭代數(shù)來控制反演是否結(jié)束；二是不再根據(jù)目標(biāo)函數(shù)的精度來進(jìn)行遺傳迭代過程中的早熟判斷.這些調(diào)整進(jìn)一步提高了反演的效率.算法簡(jiǎn)述如下：

首先用pi=(VSi,VPi,hi,ρi)表示第i個(gè)土層所需考慮的參數(shù)，i=1，2，3，…，N.設(shè)置全局反演參數(shù)，這些參數(shù)包括遺傳種群大小(M)、變量維數(shù)(反演中考慮的最大變量數(shù))、用于表示單個(gè)個(gè)體的二進(jìn)制位數(shù)、SA初始溫度T0和在特定溫度下的迭代代數(shù)(l).然后初始化遺傳種群模型，計(jì)算初始種群的目標(biāo)函數(shù)值和其中每個(gè)個(gè)體的相對(duì)適應(yīng)度值.模型種群中的每個(gè)個(gè)體用一串0和1的二進(jìn)制碼表示，其相對(duì)適應(yīng)度值表達(dá)式為

(5)

(6)

(7)

適應(yīng)度值是從上一遺傳代進(jìn)化到下一遺傳代時(shí)對(duì)個(gè)體進(jìn)行優(yōu)勝劣汰的依據(jù)，根據(jù)適應(yīng)度值進(jìn)行遺傳的選擇、重組或變異操作.種群中個(gè)體的適應(yīng)度值越大，其被復(fù)制并傳遞到下一代的概率就越高.完成遺傳操作后父代種群就可以進(jìn)化生成子代種群.

完成遺傳操作后，對(duì)于新生成的子代種群，再引入SA操作以減少GA的過早收斂.在這個(gè)操作中，我們先比較子代種群中每個(gè)個(gè)體(Φchild)與其父代種群中相應(yīng)個(gè)體(Φparent)的目標(biāo)函數(shù)，然后計(jì)算它們的差值(ΔΦ=Φchild-Φparent).當(dāng)差值ΔΦ<0時(shí)該個(gè)體可繼續(xù)遺傳到下一代，當(dāng)差值大于零時(shí)，則給定一個(gè)隨機(jī)擾動(dòng)值，按照當(dāng)前模擬退火的溫度值計(jì)算接受新個(gè)體的概率值exp(-ΔΦ/Tk)，如果該個(gè)體的接受概率大于隨機(jī)擾動(dòng)值，則接受該子代個(gè)體加入當(dāng)前的種群，參與后續(xù)的遺傳進(jìn)化，否則不接受該子代個(gè)體.接受概率計(jì)算式中Tk是SA第k次迭代的溫度，Tk隨著迭代代數(shù)(l)的增加而減少，其按下式計(jì)算獲得：

(8)

式中C0為溫度衰減系數(shù)，通常是一個(gè)小于1的常數(shù).通過遺傳和模擬退火操作的組合，反演循環(huán)計(jì)算設(shè)定的最大遺傳代數(shù)，并給出所反演的參數(shù)模型空間的最優(yōu)解.為了更清晰地展示本文混合GA和SA全局優(yōu)化算法的完整框架，用圖1給出了反演計(jì)算流程圖，從該圖也能看出反演涉及的主要步驟如下：

圖1 反演流程圖

(1)生成初始種群；

(2)設(shè)置全局反演參數(shù)(模擬退火的初始溫度、最大遺傳代數(shù)、交叉率、變異率等)；

(3)根據(jù)遺傳代數(shù)進(jìn)行循環(huán)計(jì)算.計(jì)算當(dāng)前遺傳代下所有個(gè)體在當(dāng)前溫度下的NHVs和目標(biāo)函數(shù)，之后根據(jù)目標(biāo)函數(shù)計(jì)算種群的適應(yīng)度值；

(4)執(zhí)行每一代循環(huán)的GA操作，進(jìn)行選擇、重組和變異；

(5)執(zhí)行SA操作，比較子代個(gè)體(Φchild)和父代個(gè)體(Φparent)的目標(biāo)函數(shù)，判斷在當(dāng)前溫度下子代個(gè)體是否能依據(jù)給定的隨機(jī)擾動(dòng)以可接受的概率被替換，依據(jù)最大遺傳代數(shù)不斷循環(huán)；

(6)輸出最優(yōu)解，結(jié)束反演過程.

2 反演方法的驗(yàn)證

2.1 反演參數(shù)的敏感性及反演策略的確定

前文已述對(duì)于每一層土體，都有VP、VS、h和ρ四個(gè)參數(shù)，當(dāng)需考慮的土層數(shù)目較大時(shí)，模型空間將變得很大，會(huì)顯著降低反演計(jì)算的效率.進(jìn)行參數(shù)的敏感性分析可以明確對(duì)NHV起主要影響作用的模型參數(shù)，忽略影響極小的參數(shù)，從而降低模型空間的維度，節(jié)約計(jì)算資源.為了量化參數(shù)的敏感性，本文采用Arai和Tokimatsu(2004)參數(shù)敏感性的定義，用一無量綱偏導(dǎo)數(shù)的絕對(duì)值來表征反演模型參數(shù)的敏感度，其公式表示如下：

(9)

圖2 簡(jiǎn)化模型的參數(shù)敏感性

表1 簡(jiǎn)化的土層模型

2.2 目標(biāo)NHV曲線隨機(jī)擾動(dòng)后的反演實(shí)例

為了驗(yàn)證本文方法，從多個(gè)KiK-net臺(tái)站中選擇了符合一維場(chǎng)地假設(shè)(Thompson et al., 2012)的FKSH14臺(tái)站場(chǎng)地作為一實(shí)例.KiK-net網(wǎng)站提供了該臺(tái)站場(chǎng)地地面以下106 m的PS波測(cè)井記錄，考慮到使算例模型中理論計(jì)算的NHV曲線能完整地體現(xiàn)該場(chǎng)地的深部和淺部特征，引用日本地震危險(xiǎn)性信息站(J-SHIS)給出的FKSH14臺(tái)站深部模型將其速度剖面延伸至地震基巖，但由于J-SHIS模型淺部信息較粗略，該模型第一層土厚度為170 m，我們將106～170 m范圍的土層分為兩層，其中上層的波速與PS測(cè)井揭示的工程基巖相同，下層的波速與J-SHIS模型的最上層相同，最終給出的FKSH14臺(tái)站模型見表2所示.雖然本文重點(diǎn)關(guān)注土層的模型參數(shù)，但對(duì)于理論算例而言，給出兼具深、淺部特征的速度剖面可以為反演提供更合理的初始模型，也有利于獲得深部特征對(duì)NHV影響的初步認(rèn)識(shí).

表2 FKSH14的PS測(cè)井、J-SHIS模型和下至地震基巖的綜合速度剖面

在反演算法的檢驗(yàn)研究中，一般都需對(duì)反演的目標(biāo)曲線增加一定程度的隨機(jī)擾動(dòng)以考慮反演目標(biāo)曲線的不確定性，使得反演算法的檢驗(yàn)更符合實(shí)際情況.在FKSH14場(chǎng)地速度剖面的理論NHV曲線上添加30 dB的隨機(jī)擾動(dòng)以模擬目標(biāo)NHV曲線的不確定性，并將添加隨機(jī)擾動(dòng)后的NHV曲線作為反演目標(biāo)譜，如圖3a所示.不同頻率下添加隨機(jī)擾動(dòng)后的目標(biāo)譜與理論NHV譜的差異如圖3b所示，可見目標(biāo)譜與理論NHV曲線差異在20%范圍內(nèi).確定場(chǎng)地反演目標(biāo)NHV譜之后，在開始反演過程之前，應(yīng)預(yù)先確定所有模型參數(shù)的搜索范圍(即為模型空間).由于本文僅關(guān)注場(chǎng)地的土層速度結(jié)構(gòu)，將J-SHIS模型揭示的深部速度結(jié)構(gòu)取為固定值，將FKSH14場(chǎng)地上部四層土體取為可變值.本文對(duì)上部四層土體所有模型參數(shù)的搜索范圍都是根據(jù)初始模型速度剖面確定.例如，F(xiàn)KSH14的壓縮、剪切波速度、厚度、密度分別設(shè)置為VP0、VS0、ρ0和h0，每層參數(shù)搜索范圍為0.65VP0～1.35VP0、0.65VS0～1.35VS0、0.65h0～1.35h0，每層土體密度取為固定值.

圖3 反演目標(biāo)曲線、FKSH14模型理論曲線及二者的差異

反演中將遺傳代數(shù)取為100，每一遺傳代的個(gè)體數(shù)為200.對(duì)于GA操作，代溝為0.9，交叉率為0.7，變異率為0.01.SA的初始溫度設(shè)置為10 ℃，溫度衰減系數(shù)C0為0.99.當(dāng)遺傳代數(shù)達(dá)到100代時(shí)，反演過程自動(dòng)結(jié)束，得到與最小目標(biāo)函數(shù)相對(duì)應(yīng)的反演模型作為最佳反演模型.圖4給出了反演各代的平均目標(biāo)函數(shù)與遺傳代數(shù)的關(guān)系，該圖表明隨著遺傳代數(shù)的增加，目標(biāo)函數(shù)快速降低，遺傳40代之后目標(biāo)函數(shù)曲線趨緩，此后即使遺傳代數(shù)繼續(xù)增加，目標(biāo)函數(shù)也沒有明顯的降低，這也說明經(jīng)過100代的遺傳，獲得的反演結(jié)果穩(wěn)定在最佳反演模型附近.圖4b給出了反演過程中所有的速度結(jié)構(gòu)隨模型出現(xiàn)先后的變化，該圖表明隨著遺傳代數(shù)的增加，反演逐步收斂到目標(biāo)模型，反演最佳模型與目標(biāo)模型吻合良好.圖5給出了各土層VS、h隨遺傳迭代次數(shù)增大的變化關(guān)系，該圖清晰地顯示隨著迭代次數(shù)增大，各層VS、h很快收斂至某一固定值附近，表明了反演算法的有效性.初始模型、反演模型各參數(shù)的均值和方差均示于表3，由該表可知反演模型土層的參數(shù)與表2目標(biāo)模型極為接近，尤其是各土層VS、h值幾乎重現(xiàn)了目標(biāo)模型.

圖4 目標(biāo)函數(shù)變化曲線(a)及反演過程中速度結(jié)構(gòu)的變化(b)

圖5 各層VS與h隨迭代次數(shù)增大的變化關(guān)系

表3 FKSH14的初始模型和反演模型

2.3 對(duì)三類典型模型反演效果的驗(yàn)證及不同方法的比較

為了進(jìn)一步驗(yàn)證本文方法對(duì)幾類不同模型的反演效果，根據(jù)場(chǎng)地剪切波速變化的一般規(guī)律構(gòu)建了三類典型模型，分別是剪切波速隨深度逐漸變大的遞增層模型、軟夾層和硬夾層模型，模型參數(shù)列于表4中.反演目標(biāo)譜分別取為模型理論NHV曲線添加30 dB隨機(jī)擾動(dòng)后的譜曲線，如圖6a中黑色粗實(shí)線所示.反演中將表4所示模型剪切波速2000 m·s-1及以上的基巖層各參數(shù)取為固定值，上部三層土體的縱波、橫波速度及厚度取為可變值，每層參數(shù)搜索范圍為0.65VP0～1.35VP0、0.65VS0～1.35VS0、0.65h0～1.35h0，每層土體密度取為固定值.利用本文提出的反演方法和HV-Inv法分別反演100代.對(duì)于本文方法，與上一節(jié)中反演參數(shù)設(shè)置相似，將遺傳代數(shù)取為100，每一遺傳代的個(gè)體數(shù)為200，GA操作中代溝為0.9，交叉率為0.7，變異率為0.01，SA的初始溫度設(shè)置為10 ℃，SA溫度衰減系數(shù)C0為0.99.HV-Inv反演設(shè)置中采用全局的Monte Carlo法，迭代次數(shù)設(shè)置為100，擾動(dòng)范圍設(shè)置為5%(此值控制當(dāng)前模型周圍的最大擾動(dòng)范圍，在此范圍內(nèi)產(chǎn)生下一個(gè)模型).當(dāng)代數(shù)達(dá)到100代時(shí)，反演過程自動(dòng)結(jié)束，得到與最小目標(biāo)函數(shù)相對(duì)應(yīng)的反演模型作為最佳反演模型.不同方法反演目標(biāo)函數(shù)隨遺傳代數(shù)的變化關(guān)系如圖6b所示，由該圖可知本文反演方法目標(biāo)函數(shù)下降速度更快，更易趨近于穩(wěn)定的最佳模型，說明本文反演方法較HV-Inv方法效率更高.圖6c給出了經(jīng)過相同代數(shù)的迭代之后，不同方法反演獲得的最佳模型，可知本文方法與HV-Inv方法具有幾乎相當(dāng)?shù)姆囱菪Ч容^而言，本文方法對(duì)巖土分界面的反演更為準(zhǔn)確，針對(duì)軟夾層模型，本文方法獲得的反演模型更接近目標(biāo)模型.從圖6a不同方法獲得的最佳模型理論NHV曲線對(duì)目標(biāo)NHV曲線的擬合情況來看，本文方法與HV-Inv方法獲得的最佳模型的NHV曲線均能非常好地?cái)M合目標(biāo)NHV，兩種方法對(duì)目標(biāo)NHV曲線的擬合效果基本相當(dāng).這一結(jié)果進(jìn)一步表明本文方法適用于不同類型模型的反演.

表4 遞增層模型、軟夾層和硬夾層模型參數(shù)

圖6 不同方法對(duì)三類典型模型反演結(jié)果的比較

3 反演方法的應(yīng)用實(shí)例

3.1 唐山響嘡豎向觀測(cè)臺(tái)陣場(chǎng)地及其噪聲觀測(cè)

本文選取唐山響嘡豎向觀測(cè)臺(tái)陣場(chǎng)地的3#測(cè)井作為觀測(cè)場(chǎng)地，該測(cè)井穿透強(qiáng)風(fēng)化層并到達(dá)弱風(fēng)化層，井深47 m，其鉆孔數(shù)據(jù)如表5所示(周雍年等, 2005).測(cè)井位于河北省唐山市灤縣響嘡鎮(zhèn)杜峪村(39.70°N，118.74°E)，地理位置示于圖7.

圖7 唐山響嘡豎向臺(tái)陣位置示意圖

筆者于2019年7月2—3日在3#測(cè)井場(chǎng)地布設(shè)了背景噪聲觀測(cè)臺(tái)陣以1000 Hz的采樣率進(jìn)行了連續(xù)24 h的背景噪聲探測(cè)實(shí)驗(yàn).本文截取了最接近3#測(cè)井的觀測(cè)點(diǎn)凌晨2∶00—3∶00的三分量記錄，并進(jìn)行了去均值、去趨勢(shì)項(xiàng)等預(yù)處理.圖8給出了預(yù)處理之后的三分量噪聲波形記錄及其振幅譜.

圖8 試驗(yàn)場(chǎng)地3#測(cè)點(diǎn)的噪聲觀測(cè)記錄

3.2 反演獲得的土層速度結(jié)構(gòu)模型與鉆孔的比較

本文利用開源軟件GEOPSY(http:∥www.geopsy.org-SESAME Project)處理噪聲觀測(cè)數(shù)據(jù)，該軟件可在若干可選擇的時(shí)間窗內(nèi)計(jì)算三個(gè)分量的功率譜密度，在頻域中對(duì)記錄的波場(chǎng)進(jìn)行統(tǒng)計(jì)分析.這些時(shí)間窗口的寬度取決于目標(biāo)頻帶和記錄長度.其中時(shí)間窗口的長度定義了可以進(jìn)行合理分析的最小頻率(Bard et al., 2004)，由于SESAME標(biāo)準(zhǔn)建議每個(gè)窗口至少有10個(gè)周期，因此最小頻率(以Hz為單位)定義為10/lw(lw為窗口長度，以s為單位)；對(duì)于每個(gè)窗口，應(yīng)用5%的錐度，以避免可能的截?cái)嘈?yīng)；將反觸發(fā)算法(Withers et al., 1998)應(yīng)用于噪聲記錄序列，以便在計(jì)算中僅包括具有準(zhǔn)平穩(wěn)振幅的信號(hào)，選擇標(biāo)準(zhǔn)是基于短時(shí)窗(STA)和長時(shí)窗(LTA)平均振幅的比較，STA、LTA、Min STA/LTA和Max STA/LTA的典型值分別取1、30、0.2和2.5 s，這一比率的上限旨在排除例如靠近傳感器的操作人員產(chǎn)生的尖峰；對(duì)于所考慮的所有情況，時(shí)間窗口重疊5%；頻率分析中采用了25%平滑常數(shù)的平滑函數(shù)(Picotti et al., 2017; Konno and Ohmachi, 1998)以確保低頻和高頻點(diǎn)的數(shù)量恒定.計(jì)算所得三分量功率譜密度如圖9所示.

圖9 試驗(yàn)場(chǎng)地3#測(cè)點(diǎn)觀測(cè)記錄功率譜

傳統(tǒng)的NHV曲線是根據(jù)傅里葉振幅譜(FAS; 例如Bard et al., 2004)或響應(yīng)譜(Zhu et al., 2020)計(jì)算.此處，文中根據(jù)功率譜密度(PSD)估算NHV曲線，其中通過矢量求和對(duì)兩個(gè)水平分量求平均值，然后從PSD中獲得噪聲的NHV(Arai and Tokimatsu, 2004)：

(10)

式中，下標(biāo)EW，NS和UD分別表示東西，南北和豎向分量.由于噪聲觀測(cè)結(jié)果易受環(huán)境、場(chǎng)地、儀器多因素影響，因此對(duì)觀測(cè)NHV值還需檢驗(yàn)其可靠性.歐盟SESAME項(xiàng)目(Bard et al., 2004)研究給出了檢驗(yàn)NHV可靠性和峰值清晰度的標(biāo)準(zhǔn)(詳見SESAME, 2004 “2.1結(jié)果可靠性標(biāo)準(zhǔn)”).本文依據(jù)該標(biāo)準(zhǔn)對(duì)唐山響嘡臺(tái)3#場(chǎng)地實(shí)測(cè)NHV(見圖10)進(jìn)行檢驗(yàn)，結(jié)果證實(shí)實(shí)際觀測(cè)獲得的NHV曲線滿足其可靠性和峰值清晰度準(zhǔn)則，可以認(rèn)為唐山響嘡臺(tái)3#場(chǎng)地實(shí)測(cè)NHV結(jié)果用于場(chǎng)地反演是可靠的.

圖10 唐山響嘡臺(tái)3#場(chǎng)地實(shí)測(cè)NHV

許多研究證實(shí)NHV的空間變化與地下介質(zhì)的橫向不均勻性有密切聯(lián)系(Macau et al., 2015; Bonnefoy-Claudet et al., 2009; Matsushima et al., 2014; Uebayashi et al., 2012).本文通過旋轉(zhuǎn)實(shí)測(cè)噪聲記錄的東西、南北水平分量計(jì)算得到了NHV隨水平面方位角的變化(見圖11)，結(jié)果顯示不同的水平面方位角情形下NHV頻譜曲線較為一致，NHV頻譜曲線隨水平面方位角的變化是微小的，表明該場(chǎng)地在東西、南北兩個(gè)水平方向僅表現(xiàn)出弱各向異性，將其近似為一維場(chǎng)地是合適的.圖12給出了地表觀測(cè)NHV與鉆孔模型正演NHV的比較，雖然觀測(cè)NHV與鉆孔模型的理論NHV譜形接近，但圖中也能看到二者在1.5 Hz以上頻段仍有明顯差異，據(jù)此可以推測(cè)反演模型的自振周期等特征可能接近鉆孔模型，而分層波速等細(xì)部特征可能與鉆孔存在明顯差別.

圖11 唐山響嘡臺(tái)3#場(chǎng)地基于傅里葉譜的NHV隨水平面方位角的變化

圖12 地表觀測(cè)NHV與鉆孔模型正演NHV的比較及不同方法對(duì)地表觀測(cè)NHV曲線的擬合

對(duì)土層速度結(jié)構(gòu)的反演需要事先設(shè)定一個(gè)假設(shè)的初始速度結(jié)構(gòu)模型.本文采用王繼鑫等(2020)提取的3.5～20 Hz實(shí)測(cè)頻散曲線，以0.1 Hz為步長采用改進(jìn)半波長法估算出試驗(yàn)場(chǎng)地剪切波速剖面.由于本文主要關(guān)注場(chǎng)地的土層，巖土工程一般將剪切波速>500 m·s-1的層位定義為軟基巖，借助上述改進(jìn)半波長法估算的剪切波速剖面，得出剪切波速<500 m·s-1的覆蓋土層厚度為29.2 m.再結(jié)合實(shí)測(cè)NHV的最高頻率(20 Hz)以及頻散曲線該頻率下的相速度(193 m·s-1)，可以簡(jiǎn)單估算NHV曲線最高頻率所對(duì)應(yīng)的最小土層厚度為9.65 m，即用NHV反演所能分辨的土層厚度將不小于9.65 m.鑒于此，從工程應(yīng)用角度對(duì)初始模型進(jìn)行了一定簡(jiǎn)化，將覆蓋土層總厚度取為30 m，并按10 m(高于最小分辨率)一層分為3層.根據(jù)此分層給出初始速度結(jié)構(gòu)模型如表6所示.

隨后分別利用本文提出的反演方法和HV-Inv法反演場(chǎng)地土層速度結(jié)構(gòu).對(duì)于本文方法，采用第2.3節(jié)理論算例中相同的反演參數(shù)，將遺傳代數(shù)取為200，每一遺傳代的個(gè)體數(shù)為200，GA操作中代溝值為0.9，交叉率為0.7，變異率為0.01，SA的初始溫度設(shè)置為10℃，SA溫度衰減系數(shù)C0為0.99.當(dāng)遺傳代數(shù)達(dá)到200代時(shí)，反演過程自動(dòng)結(jié)束，得到與最小目標(biāo)函數(shù)相對(duì)應(yīng)的反演模型即為最佳反演模型.對(duì)于HV-Inv法反演，首先采用全局的Monte Carlo抽樣，迭代次數(shù)設(shè)置為100，擾動(dòng)范圍設(shè)置為10%(此值控制當(dāng)前模型周圍的最大擾動(dòng)范圍，在此范圍內(nèi)產(chǎn)生下一個(gè)模型)，而后采用修正的模擬退火(Modified SA)，其初始溫度設(shè)置為6.0℃，溫度衰減系數(shù)C0為0.9，當(dāng)遺傳代數(shù)達(dá)到500代時(shí)，反演過程自動(dòng)結(jié)束，得到與最小目標(biāo)函數(shù)相對(duì)應(yīng)的反演模型即為最佳反演模型.反演過程中不同方法的目標(biāo)函數(shù)隨迭代次數(shù)的變化曲線如圖13所示，由該圖可知本文方法收斂速度遠(yuǎn)快于HV-Inv方法，反演經(jīng)100代之后本文方法的反演結(jié)果即趨于穩(wěn)定.圖14給出了不同方法反演的模型與鉆孔數(shù)據(jù)的比較.VP、VS、h隨迭代次數(shù)的變化情況示于圖15，限于篇幅，僅給出了第一層土體參數(shù)的結(jié)果，該結(jié)果能清晰地顯示反演參數(shù)隨迭代次數(shù)增大快速收斂至最佳值的過程.我們也很容易注意到圖15所示的土層厚度并不隨迭代次數(shù)增加而變化，這主要是因?yàn)槿狈s束土層厚度的信息，將土層厚度取為了固定值.兩種反演方法獲得的最佳擬合模型示于圖14中，將表5所示的鉆孔模型采用等效波速法折算到相應(yīng)層厚作為反演結(jié)果的參照.由該圖可知，兩種反演方法獲得的四層土體速度結(jié)構(gòu)模型總體上與鉆孔揭示的速度結(jié)構(gòu)較為接近，本文方法相比HV-Inv方法能更準(zhǔn)確地獲得基巖層剪切波速.對(duì)第一和第三層土體的波速不同方法的反演結(jié)果接近，且與鉆孔數(shù)據(jù)較為接近，對(duì)第二層土體本文方法與鉆孔數(shù)據(jù)的差別要遠(yuǎn)小于HV-Inv方法且反演結(jié)果與鉆孔數(shù)據(jù)十分接近.圖12給出了反演模型對(duì)應(yīng)的NHV曲線對(duì)目標(biāo)曲線的擬合結(jié)果，該結(jié)果表明本文反演模型能很好地?cái)M合目標(biāo)曲線，反演模型可視為更符合地表觀測(cè)結(jié)果的土體等效速度結(jié)構(gòu)模型，其模型參數(shù)的不確定性可用表6中的方差來描述.

表6 初始模型和反演模型

圖13 歸一化目標(biāo)函數(shù)隨迭代次數(shù)的變化曲線

圖14 不同方法反演的最佳模型及與鉆孔數(shù)據(jù)的比較

圖15 參數(shù)隨迭代次數(shù)變化

綜合圖12—15的結(jié)果我們可以看到，無論是采用本文方法還是HV-Inv方法，反演過程本質(zhì)上都是使得反演模型的NHV曲線不斷趨近目標(biāo)NHV曲線的過程，而目標(biāo)NHV與鉆孔模型的理論NHV在幅值和頻譜上本身就有一定差別(如圖12)，從這個(gè)角度來看，反演得到的最佳模型與鉆孔有差別屬于正?，F(xiàn)象，鉆孔數(shù)據(jù)僅能作為反演結(jié)果的參照，但從圖11揭示的結(jié)果，將該場(chǎng)地視為一維場(chǎng)地是合理的，鉆孔在很大程度上能較好地代表背景噪聲臺(tái)站布設(shè)的場(chǎng)地特征，因此反演模型與鉆孔越接近表明反演效果越好，圖14的結(jié)果恰好證實(shí)本文反演方法的效果要優(yōu)于已有的HV-Inv方法，尤其對(duì)于巖土分界面的揭示及基巖剪切波速的確定，本文方法更有明顯的優(yōu)勢(shì).

4 結(jié)論

場(chǎng)地土層是各類工程建設(shè)極為關(guān)注的近地表介質(zhì)層，利用地表觀測(cè)的背景噪聲反演場(chǎng)地土層速度結(jié)構(gòu)的研究仍較少.本文將已有的遺傳模擬退火混合算法與基于散射場(chǎng)假定的NHV正演計(jì)算方法相結(jié)合，提出了一種基于背景噪聲波場(chǎng)NHV反演場(chǎng)地土層速度結(jié)構(gòu)的全局反演方法，通過簡(jiǎn)化模型的參數(shù)敏感性分析，揭示了場(chǎng)地土層的VP、VS及厚度是影響場(chǎng)地地表NHV的主要因素，而密度參數(shù)的影響幾乎可以忽略.理論算例驗(yàn)證了本文全局反演方法的可靠性.隨后本文利用唐山響嘡豎向臺(tái)陣場(chǎng)地的背景噪聲觀測(cè)獲得了面波頻散曲線，為土層速度結(jié)構(gòu)反演提供了初始模型信息.基于該初始模型，采用本文反演方法獲得了較為接近鉆孔測(cè)井?dāng)?shù)據(jù)的地下速度結(jié)構(gòu).

理論和實(shí)際應(yīng)用算例以及與同類方法的比較均表明本文方法是基于地表背景噪聲觀測(cè)反演場(chǎng)地土層速度結(jié)構(gòu)的合理、可靠的方法，可為近地表土層模型參數(shù)的確定提供一種低成本且有效的途徑，但由于土層速度結(jié)構(gòu)的背景噪聲反演研究尚處于起步階段，目前反演獲得的土層速度結(jié)構(gòu)的細(xì)致程度仍取決于與面波頻散曲線相關(guān)的初始模型的精細(xì)程度以及地表NHV曲線的頻段范圍，這方面的研究仍有大量基本問題如模型參數(shù)對(duì)背景噪聲波場(chǎng)的影響規(guī)律、反演模型的分辨率及其不確定性評(píng)價(jià)等都亟待解決，對(duì)于這些問題尚需進(jìn)一步開展更深入的研究工作.