場地土層速度結構的背景噪聲反演方法及其應用

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

1 北京工業(yè)大學城市與工程安全減災教育部重點實驗室，北京 100124 2 防災科技學院河北省地震災害防御與風險評價重點實驗室，河北三河 065201

0 引言

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

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

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

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

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

1 NHV全局反演方法

1.1 土層速度結構的NHV正演方法

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

散射場理論解是隨著背景噪聲成像的地震干涉理論發(fā)展起來的(秦彤威等, 2021)，地震干涉理論認為，完全散射波場中任意兩個觀測點獲得的位移記錄的互相關正比于兩觀測點之間的格林函數(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點j方向作用的一個單元簡諧點力δijδ(|x-xB|)exp(iωt)在xA點i方向引起的位移.在單元簡諧點力δi jδ(|x-xB|)exp(iωt)中，i為虛數(shù)單位，ω為圓頻率，t為時間.k=ω/β為S波波數(shù)，β為S波波速，Es=ρω2S2為平均S波能流密度，S2為S波平均譜密度，星號表示復共軛，Im表示給定的函數(shù)的虛部，尖括號表示平均.假定xA=xB，上式可以寫為

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

(2)

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

(3)

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

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

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

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

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

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

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

(5)

(6)

(7)

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

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

(8)

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

圖1 反演流程圖

(1)生成初始種群；

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

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

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

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

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

2 反演方法的驗證

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

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

(9)

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

表1 簡化的土層模型

2.2 目標NHV曲線隨機擾動后的反演實例

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

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

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

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

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

圖4 目標函數(shù)變化曲線(a)及反演過程中速度結構的變化(b)

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

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

2.3 對三類典型模型反演效果的驗證及不同方法的比較

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

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

圖6 不同方法對三類典型模型反演結果的比較

3 反演方法的應用實例

3.1 唐山響嘡豎向觀測臺陣場地及其噪聲觀測

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

圖7 唐山響嘡豎向臺陣位置示意圖

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

圖8 試驗場地3#測點的噪聲觀測記錄

3.2 反演獲得的土層速度結構模型與鉆孔的比較

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

圖9 試驗場地3#測點觀測記錄功率譜

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

(10)

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

圖10 唐山響嘡臺3#場地實測NHV

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

圖11 唐山響嘡臺3#場地基于傅里葉譜的NHV隨水平面方位角的變化

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

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

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

表6 初始模型和反演模型

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

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

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

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

4 結論

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

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