韓復(fù)興,宋炳宣,陳雨貝,吳林駿,黃夢婷,潘延杰

吉林大學(xué)地球探測科學(xué)與技術(shù)學(xué)院,長春 130026

0 引言

微動H/V譜比方法(又稱HVSR(horizontal to vertical spectral ratio)方法)是一種快速、經(jīng)濟(jì)的非侵入探測方法,它在城市地球物理勘探中發(fā)揮重要作用。1989年日本學(xué)者Nakamura[1]率先提出用同一地表測點(diǎn)地脈動水平分量與豎向分量的傅里葉幅值譜比值,即H/V譜比,來估計(jì)場地的共振頻率和放大因子。H/V譜比法認(rèn)為H/V譜比曲線的峰值頻率對應(yīng)于場地剪切波的共振頻率,通過該頻率可以獲得場地的放大系數(shù)和沉積層厚度,進(jìn)而推測地下地層構(gòu)造。

微動H/V譜比方法對于劃分土石分界面、檢索剪切波阻抗變化等信息實(shí)際應(yīng)用效果好,且具有便捷高效、勘探周期短的優(yōu)點(diǎn)。與野外地球物理勘探相比,城市地球物理探測一方面精度要求高,另一方面需要考慮城市復(fù)雜的噪聲和人文環(huán)境[2-3]。微動H/V譜比方法作為以非地震引起的微弱震動為觀測對象的非侵入探測方法,契合以上兩方面特點(diǎn),故在城市地球物理勘探中應(yīng)用較多。但城市中震動干擾因素眾多,如車輛震動噪聲、施工噪聲等,采集的一手譜比數(shù)據(jù)用于譜比計(jì)算效果較差;因此需要對原始的三分量數(shù)據(jù)進(jìn)行信噪分析,獲得精細(xì)化處理后的優(yōu)質(zhì)譜比曲線。

H/V譜比法的數(shù)據(jù)精細(xì)化處理環(huán)節(jié)可在時(shí)域和頻域上分別或結(jié)合進(jìn)行。文獻(xiàn)資料表明,目前國內(nèi)外基于時(shí)域識別瞬態(tài)干擾信號多采用長短時(shí)窗平均振幅比(short time average/long time average, STA/LTA)方法[4],如:2018年Setiawan等[5]在南澳大利亞州阿德萊德市進(jìn)行了環(huán)境噪聲測量,在此基礎(chǔ)上通過STA/LTA方法排除瞬態(tài)干擾,選擇環(huán)境噪聲中最穩(wěn)定的部分,并運(yùn)用H/V譜比法研究地震場地劃分;張若晗等[6]選用微動H/V譜比法對濟(jì)南中心城區(qū)的土石分界面展開研究,在計(jì)算每個(gè)時(shí)窗的H/V譜比之前先采用反觸發(fā)STA/LTA方法剔除瞬態(tài)干擾;在SESAME(site effects assessment using ambient excitations)團(tuán)隊(duì)的H/V譜比技術(shù)實(shí)施指南中[7],該算法也被推薦用于時(shí)窗選擇。

同樣是基于時(shí)間域,2016年D’Alessandro等[8]提出通過層次聚類算法實(shí)現(xiàn)時(shí)窗的自動選擇。該算法通過提取HVSR曲線的自洽聚類來對數(shù)據(jù)進(jìn)行處理,消除了STA/LTA方法中閾值和時(shí)窗長度選取的主觀因素影響,提高了數(shù)據(jù)分析的客觀性。

近年,基于頻域的降噪算法也逐漸得到研究者的關(guān)注。Cox等[9]提出了一種基于頻域的拒絕算法,在進(jìn)行計(jì)算時(shí),該方法可有效剔除峰值頻率與峰值頻率均值相差過大的時(shí)間窗口,對數(shù)據(jù)進(jìn)行降噪。Dal等[10]同樣通過譜比曲線峰值提出了對多峰譜比曲線進(jìn)行降噪的頻域方法。

此外,2015年Liu等[11]結(jié)合希爾伯特-黃變換適合處理非平穩(wěn)信號的特點(diǎn),提出采用希爾伯特-黃變換來對微動數(shù)據(jù)進(jìn)行分析。該方法通過經(jīng)驗(yàn)?zāi)B(tài)分解,將原始信號分解為固有模態(tài)函數(shù),并將頻率最高的固有模態(tài)函數(shù)去除,剩余部分作為微動信號的有效部分進(jìn)行希爾伯特-黃變換以實(shí)現(xiàn)對原始信號的去噪處理。

除了自動降噪算法,還有一些人工方法被用于降噪,例如:根據(jù)數(shù)據(jù)使用帶通濾波器,使用0.5～20.0 Hz的帶通濾波器對福州地區(qū)進(jìn)行微動探測[12];或?qū)D譜進(jìn)行目視檢查,手動刪除受局部干擾影響的部分[13]。

上述方法在實(shí)際運(yùn)用過程中均存在面對高噪聲數(shù)據(jù)信號提取困難、受主觀因素影響大、參數(shù)調(diào)整困難等問題。對此,本文提出一種基于XGBoost(extreme gradient boosting)的多重加權(quán)譜比方法對譜比曲線進(jìn)行預(yù)處理操作,從而實(shí)現(xiàn)較高噪聲下的微動信號提取。本方法在訓(xùn)練集建立后,只需將原始數(shù)據(jù)讀入軟件即可自動完成信號提取,在提取過程中無需調(diào)參,降低了微動數(shù)據(jù)信號提取的工作量,實(shí)現(xiàn)了微動信號噪聲壓制的智能化。

1 算法技術(shù)路線

1.1 多重加權(quán)譜比算法

微動數(shù)據(jù)可用時(shí)間序列X(H1,H2,V)表示(H1、H2為水平分量,V為垂直分量)。在數(shù)據(jù)處理過程中,X應(yīng)由X′(f,I,H/V)代替(f為樣本點(diǎn)頻率,I為樣本點(diǎn)傅里葉頻譜幅值,H/V為水平分量與垂直分量的傅里葉幅值比)。顯然,僅從X′的3個(gè)維度(f,I,H/V)難以直接判斷該采樣點(diǎn)是否為噪聲,故應(yīng)對X′進(jìn)行升維,在更高維空間中尋找噪點(diǎn)與信號的特征。

1.1.1 幅值加權(quán)譜比

幅值加權(quán)譜比可表示為

(1)

由原始譜比數(shù)據(jù)(圖1a)與幅值加權(quán)的譜比數(shù)據(jù)(圖1b)對比可見,具備穩(wěn)定信號且譜比值較大的數(shù)據(jù)位于幅值加權(quán)譜比曲線的頂端。

a. 原始譜比數(shù)據(jù);b. 幅值加權(quán)譜比;c. 頻率加權(quán)譜比;d. 多重加權(quán)譜比。

由地脈波穩(wěn)定、長周期的特點(diǎn)可知,加權(quán)譜比應(yīng)當(dāng)篩選出幅值穩(wěn)、譜比接近峰值的采樣點(diǎn)。故引入經(jīng)驗(yàn)常數(shù)λ,一方面改造了極點(diǎn),使程序簡單,另一方面符合加權(quán)譜比的意義。加權(quán)譜比是幅值賦權(quán)下的H/V,采樣點(diǎn)的加權(quán)譜比值應(yīng)由幅值與H/V值共同決定,加權(quán)后二者對加權(quán)譜比的數(shù)值影響相當(dāng)。

1.1.2 頻率加權(quán)譜比

常見H/V譜比曲線有個(gè)明顯的峰值頻率,單臺數(shù)據(jù)無法區(qū)分微動是由體波還是面波組成,很難判斷H/V峰值是由體波中的剪切波在松散沉積層的共振引起,還是由瑞利面波的極化作用或勒夫波的震相引起。通過Sylvette等[14]的模擬實(shí)驗(yàn)發(fā)現(xiàn),微動H/V譜比曲線峰值頻率與波場組成關(guān)系不大,而與松散沉積層的共振頻率相吻合,峰值頻率接近或等于松散沉積層的共振頻率。從場地效應(yīng)的傳遞函數(shù)出發(fā),可推導(dǎo)出H/V譜比曲線峰值頻率與松散覆蓋層的平均剪切波速和厚度的關(guān)系如下:

(2)

式中:f0為H/V譜比曲線峰值頻率;vS為覆蓋層平均剪切波速;Dov為松散覆蓋層厚度,可由測井等工程資料直接獲得。H/V譜比曲線關(guān)注峰值頻率,故給出頻率加權(quán)譜比:

(3)

如圖1c所示,在頻率加權(quán)譜比曲線中,兼具信號位于峰值信號附近且譜比值較大的數(shù)據(jù)處于加權(quán)譜比曲線的頂端。

1.1.3 多重加權(quán)譜比

由式(1)(3),可定義多重加權(quán)譜比:

(4)

如圖1d所示,在多重加權(quán)譜比曲線中,兼具信號穩(wěn)定、位于峰值信號附近且譜比值較大的數(shù)據(jù)處于加權(quán)譜比曲線的頂端。

由式(1)(3)(4)可對X′進(jìn)行升維,得到

由地脈波穩(wěn)定的特點(diǎn)與HVSR對峰值頻率的需要,可得信號點(diǎn)均位于X″中多重加權(quán)譜比曲線的頂端,故經(jīng)升維后的X″可在此高維空間中將信號與噪聲較好地分離。

1.2 XGBoost算法

(5)

XGBoost旨在通過迭代方式集成弱分類器以形成預(yù)測精度和魯棒性更高的模型。每一次迭代都是在前一步的基礎(chǔ)上增加一棵樹,而新增的樹視為擬合上次預(yù)測的殘差。迭代過程為:

(6)

目標(biāo)函數(shù)是研究一切機(jī)器學(xué)習(xí)問題的出發(fā)點(diǎn),求解目標(biāo)函數(shù)最低數(shù)值即為求解機(jī)器學(xué)習(xí)問題解。XGBoost的目標(biāo)函數(shù)由損失函數(shù)和正則項(xiàng)兩部分組成。正則項(xiàng)用于控制模型的復(fù)雜度,模型的預(yù)測精度由偏差和方差共同決定。由于本研究面向成分復(fù)雜、易出現(xiàn)極端數(shù)值的微動數(shù)據(jù),故需要利用正則項(xiàng)防止模型過擬合:

(7)

式中:T為節(jié)點(diǎn)數(shù)量;ωj為葉子節(jié)點(diǎn)j的分?jǐn)?shù);γ和ε均為系數(shù)。

綜上,XGBoost的目標(biāo)函數(shù)為

(8)

XGBoost采用增量訓(xùn)練的方式進(jìn)行迭代,以最大化減小目標(biāo)函數(shù)值。在t輪迭代之后,目標(biāo)函數(shù)更新為

(9)

為了使損失函數(shù)梯度收斂更快、更準(zhǔn)確,XGBoost利用泰勒展開式對所有二階可導(dǎo)的損失函數(shù)做近似替換,目標(biāo)函數(shù)近似為

(10)

將每個(gè)數(shù)據(jù)的損失函數(shù)疊加,去掉無關(guān)項(xiàng),得到

(11)

(12)

上述ωj相互獨(dú)立,當(dāng)新增樹的結(jié)構(gòu)已知時(shí),可以求得葉節(jié)點(diǎn)對應(yīng)的最優(yōu)權(quán)重和最優(yōu)目標(biāo)函數(shù):

(13)

(14)

在實(shí)際運(yùn)用中,只需保存訓(xùn)練好的XGBoost模型,并將待預(yù)測數(shù)據(jù)輸入使用即可。

1.3 訓(xùn)練集建立

幅值在幅值均值附近的信號符合地脈波穩(wěn)定的特征,實(shí)際應(yīng)用效果顯示,對干擾較大的數(shù)據(jù),幅值中位數(shù)對幅值均值的衡量效果優(yōu)于幅值算術(shù)平均值的衡量效果;因此,將傅里葉頻譜幅值位于幅值中位數(shù)鄰域的信號標(biāo)注為關(guān)注信號。頻率方面,用地層已知的HVSR曲線進(jìn)行訓(xùn)練集建立,將f0鄰域內(nèi)的數(shù)據(jù)標(biāo)注為關(guān)注信號,對基頻進(jìn)行重點(diǎn)提取。對地層未知的曲線可選取關(guān)注的頻率窗口,如本文采取頻率位于0.05～20.00 Hz之間的信號為有效信號,實(shí)現(xiàn)關(guān)注窗口的重點(diǎn)提取。綜合上述,幅值與頻率的選取可建立訓(xùn)練集。

1.4 算法適用范圍分析

顯然,本算法流程中多重加權(quán)譜比算法的建立效果決定了本算法的適用范圍。對于受低頻高振幅噪聲污染嚴(yán)重的微動數(shù)據(jù),如圖2所示,有效信號幅值遠(yuǎn)小于干擾數(shù)據(jù)幅值,加權(quán)譜比算法建立效果較差,難以提取有效信息。

a. 南北方向;b. 東西方向;c. 垂直方向。

2 實(shí)測微動數(shù)據(jù)信號提取分析

為了檢驗(yàn)基于XGBoost的多重加權(quán)譜比降噪方法的有效性,本文應(yīng)用該方法與STA/LTA方法分別對吉林大學(xué)朝陽校區(qū)的微動譜比數(shù)據(jù)進(jìn)行處理。朝陽校區(qū)臨近文化廣場,有多種文娛活動在此開展,噪聲復(fù)雜,干擾較大。圖3為三分量微動儀記錄的原始數(shù)據(jù),可見前72 s內(nèi)信號幅值巨大,72 s后信號幅值難以辨認(rèn)。圖4為較大干擾事件的局部數(shù)據(jù)放大圖,可見前72 s干擾嚴(yán)重。

a. 南北方向;b. 東西方向;c. 垂直方向。

a. 南北方向;b. 東西方向;c. 垂直方向。

圖5為STA/LTA方法篩選后的時(shí)間序列數(shù)據(jù),可見前72 s存在的大振幅干擾事件并沒有被去除,這是STA/LTA方法固有缺陷所導(dǎo)致的。在STA/LTA方法中,閾值的設(shè)定和時(shí)窗的選取具有主觀性且對降噪效果的影響較大。閾值和時(shí)窗主要基于統(tǒng)計(jì)規(guī)律來確定,當(dāng)高振幅噪聲持續(xù)時(shí)間較長、超過長時(shí)窗所截時(shí)段,且閾值選擇較高時(shí),該噪聲難以被去除。而時(shí)窗選擇過長,閾值選擇較低,可能導(dǎo)致有效信息不足。在該組數(shù)據(jù)的處理中,前72 s內(nèi)的噪聲持續(xù)時(shí)間長、振幅大,故STA/LTA方法難以去除,且過濾了較多有效信息,導(dǎo)致譜比曲線峰值不清晰。

a. 南北方向;b. 東西方向;c. 垂直方向。

分別對該信號利用STA/LTA方法與本文方法進(jìn)行信號提取,圖6為兩種方法提取結(jié)果的譜比曲線?？梢姳疚姆椒ɑl峰值與高頻的譜比峰值均得到有效提取,SLA/LTA方法難以識別曲線峰值。在面對較復(fù)雜環(huán)境時(shí),STA/LTA方法難以識別微動信號中的有效數(shù)據(jù),導(dǎo)致譜比曲線峰值無法識別,本文方法可清晰識別譜比曲線峰值。

圖6 STA/LTA方法(a)與多重加權(quán)降噪方法(b)譜比曲線對比

在智能化程度方面,相較于STA/LTA方法需要依靠經(jīng)驗(yàn)判斷時(shí)窗與閾值,本文方法在模型訓(xùn)練結(jié)束后,無需調(diào)整任何參數(shù)即可實(shí)現(xiàn)信號提取,實(shí)現(xiàn)信號提取智能化,增加提取結(jié)果客觀程度,減少因調(diào)參產(chǎn)生的主觀性和不確定性。

在算法效率方面,使用平臺為AMD RyzenTM7 5800H 3.2 GHz。由于本算法的計(jì)算時(shí)長主要取決于使用模型大小,以16 GB模型為例,本文方法用時(shí)728 s,STA/LTA方法用時(shí)28 s。雖然用時(shí)大幅上升,但STA/LTA方法無法處理上文中的高噪聲微動數(shù)據(jù),且本文方法使用時(shí)無需調(diào)參,更為智能化。

3 結(jié)論

1)本文針對復(fù)雜環(huán)境下高噪聲微動數(shù)據(jù)的降噪需求,提出了基于XGBoost的多重加權(quán)譜比降噪方法。

2)對高噪聲的實(shí)際數(shù)據(jù)應(yīng)用顯示,本文方法信號提取的準(zhǔn)度優(yōu)于STA/LTA方法。同時(shí),本文方法在模型訓(xùn)練完畢后無需調(diào)整參數(shù),在智能化程度上也優(yōu)于傳統(tǒng)STA/LTA方法。但本文方法訓(xùn)練模型較大,耗時(shí)較長,在干擾信號振幅大于微動數(shù)據(jù)振幅時(shí)方法失效。

針對本文方法的缺陷,未來將重點(diǎn)研究干擾信號振幅大于微動數(shù)據(jù)振幅的強(qiáng)干擾數(shù)據(jù)信號提取,并加入并行算法,使現(xiàn)有方法運(yùn)行效率更高。