基于地震動(dòng)信號(hào)分析的地質(zhì)災(zāi)害過(guò)程重構(gòu)方法研究與應(yīng)用*

嚴(yán) 炎 崔一飛 周開來(lái) 尹述遙 田 鑫 田繼楓

(①西南交通大學(xué)土木工程學(xué)院，成都 610031，中國(guó))(②清華大學(xué)水沙科學(xué)與水利水電工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 100084，中國(guó))(③西南交通大學(xué)高速鐵路線路工程教育部重點(diǎn)實(shí)驗(yàn)室，成都 610031，中國(guó))(④西南交通大學(xué)陸地交通地質(zhì)災(zāi)害防治技術(shù)國(guó)家工程實(shí)驗(yàn)室，成都 610031，中國(guó))(⑤中國(guó)地震局地球物理研究所，北京 100081，中國(guó))

0 引 言

中國(guó)西南山區(qū)地震多發(fā)、地勢(shì)陡峻、地形地質(zhì)復(fù)雜、降雨豐沛，從而導(dǎo)致地質(zhì)災(zāi)害頻發(fā)，嚴(yán)重?fù)p害了人民的生命財(cái)產(chǎn)安全(崔鵬等，2008；殷躍平，2008；黃潤(rùn)秋，2009)。如2010年四川綿竹“8·13”群發(fā)性泥石流(許強(qiáng)，2010；蘇鵬程等，2011)、2017年四川茂縣“6·24”新磨特大滑坡(許強(qiáng)等，2017；曾慶利等，2018；Yan et al.，2020a)。然而由于地質(zhì)災(zāi)害的強(qiáng)破壞性和突發(fā)性，導(dǎo)致現(xiàn)有的監(jiān)測(cè)方法難以獲取完整的災(zāi)害過(guò)程數(shù)據(jù)，這就使得對(duì)地質(zhì)災(zāi)害演進(jìn)過(guò)程的系統(tǒng)性研究工作難以開展。因此，亟需探尋一種有效的方法用于研究災(zāi)害的演進(jìn)過(guò)程。

在地質(zhì)災(zāi)害研究的早期，大部分研究采用的是現(xiàn)場(chǎng)調(diào)查的方法，但其浪費(fèi)人力、物力和時(shí)間，且應(yīng)用范圍較小，具有一定的局限性。隨著科技的進(jìn)步，一種新的技術(shù)受到廣大研究者們的青睞——遙感技術(shù)，其獲取的遙感影像用于災(zāi)害的早期識(shí)別可對(duì)災(zāi)害應(yīng)急指導(dǎo)提供有效支持(許強(qiáng)等，2019)，通過(guò)定期監(jiān)測(cè)的遙感數(shù)據(jù)可描述災(zāi)害體的變形過(guò)程(Ouyang et al.，2019)，再配合詳細(xì)的災(zāi)后現(xiàn)場(chǎng)調(diào)查資料能夠?qū)崿F(xiàn)對(duì)災(zāi)害過(guò)程的判斷，基于這些優(yōu)點(diǎn)，該技術(shù)已被廣泛運(yùn)用于地質(zhì)災(zāi)害研究中(鄧輝，2007；何秋玲，2019；許強(qiáng)，2020)。但是，對(duì)于很多新生的地質(zhì)災(zāi)害，往往缺乏定期的遙感或現(xiàn)場(chǎng)觀測(cè)調(diào)查資料，而僅通過(guò)災(zāi)害發(fā)生前后的資料難以對(duì)演進(jìn)過(guò)程形成清晰準(zhǔn)確的認(rèn)知。近年來(lái)，計(jì)算機(jī)技術(shù)的迅速發(fā)展使得數(shù)值模擬逐漸成為一種較成熟的災(zāi)害運(yùn)動(dòng)過(guò)程還原方法，通過(guò)調(diào)整相關(guān)參數(shù)，根據(jù)數(shù)值模擬的結(jié)果便可以清晰準(zhǔn)確地認(rèn)識(shí)真實(shí)的災(zāi)害演進(jìn)過(guò)程。數(shù)值模擬方法包括離散元法(DEM)、有限元法(FEM)、光滑粒子流動(dòng)力學(xué)法(SPH)等，國(guó)內(nèi)外學(xué)者運(yùn)用這些方法已開展了大量研究對(duì)災(zāi)害過(guò)程進(jìn)行模擬(Huang et al.，2012；Moretti et al.，2012；周禮等，2019)。數(shù)值模擬結(jié)果的合理性可通過(guò)與災(zāi)害現(xiàn)場(chǎng)的最終狀態(tài)的吻合度來(lái)確定，但模擬的中間過(guò)程可能與實(shí)際仍存在一定偏差。鑒于現(xiàn)有方法的局限性，需要提出一種能夠用于完整、準(zhǔn)確地分析災(zāi)害過(guò)程的方法。

隨著對(duì)地質(zhì)災(zāi)害研究的不斷深入，發(fā)現(xiàn)大部分的地質(zhì)災(zāi)害(落石、滑坡、泥石流等)在其發(fā)生過(guò)程中都伴隨著強(qiáng)烈的地表或地表以下震動(dòng)，災(zāi)害的類型、規(guī)模、位置等不同則會(huì)導(dǎo)致地震動(dòng)信號(hào)的特征有所差異，因此通過(guò)對(duì)災(zāi)害產(chǎn)生的地震動(dòng)信號(hào)分析，便可以較為完整且準(zhǔn)確地分析災(zāi)害過(guò)程相關(guān)特征信息。一些研究者基于地震動(dòng)信號(hào)對(duì)落石災(zāi)害進(jìn)行研究(Vilajosana et al.，2008；Hibert et al.，2011；Farin et al.，2015)發(fā)現(xiàn)落石產(chǎn)生的地震動(dòng)信號(hào)可以提供與其有關(guān)的動(dòng)力學(xué)和位置信息，判斷出落石產(chǎn)生的原因，通過(guò)地震動(dòng)信號(hào)分析還能估算出落石的體積。也有一些研究對(duì)滑坡產(chǎn)生的長(zhǎng)周期和短周期的地震動(dòng)信號(hào)進(jìn)行解釋分析，揭示了邊坡破壞的動(dòng)力特性，并對(duì)滑坡動(dòng)態(tài)過(guò)程的不同階段進(jìn)行了定性描述(Chen et al.，2014；Hibert et al.，2014)。在20世紀(jì)末，便有研究者注意到了泥石流災(zāi)害與地震動(dòng)的相關(guān)性，并通過(guò)地震動(dòng)數(shù)據(jù)對(duì)流體流速進(jìn)行估算(Arattano，1999)，在此基礎(chǔ)上對(duì)泥石流災(zāi)害的研究日漸增多(Burtin et al.，2016；Schimmel et al.，2016)，這些研究表明基于地震動(dòng)信號(hào)分析可以對(duì)泥石流災(zāi)害的地貌事件進(jìn)行探測(cè)、定位和量化，同時(shí)對(duì)地震動(dòng)信號(hào)進(jìn)行時(shí)頻分析以識(shí)別地貌變化過(guò)程，通過(guò)對(duì)泥石流流域地震動(dòng)信號(hào)的分析還能得出其相關(guān)運(yùn)動(dòng)特征信息。此外，Dammeier et al.(2011)還探究了災(zāi)害特征與地震動(dòng)信號(hào)參數(shù)之間的關(guān)系，如地震動(dòng)信號(hào)持續(xù)時(shí)間、信號(hào)峰值包絡(luò)線、能量包絡(luò)線等，為災(zāi)害過(guò)程的反演提供了思路。這些研究均驗(yàn)證了基于地震動(dòng)信號(hào)進(jìn)行地質(zhì)災(zāi)害過(guò)程研究的可行性，且為后期的深入研究奠定了基礎(chǔ)。

本研究提出了一套基于地震動(dòng)信號(hào)的災(zāi)害重構(gòu)思路和方法體系：首先從地質(zhì)災(zāi)害發(fā)生地的周圍地震臺(tái)站提取出地震動(dòng)信號(hào)，同時(shí)對(duì)災(zāi)害后現(xiàn)場(chǎng)進(jìn)行調(diào)查，獲取一些基本資料；接著運(yùn)用BP-filter、EMD、FFT、STFT、PSD對(duì)信號(hào)進(jìn)行處理與分析，從而獲取時(shí)頻特征、能量分布和動(dòng)力特征等信號(hào)特征信息；然后結(jié)合現(xiàn)場(chǎng)調(diào)查資料反演災(zāi)害的基本特征；基于此最后運(yùn)用數(shù)值模擬方法實(shí)現(xiàn)災(zāi)害過(guò)程的重構(gòu)。該方法已運(yùn)用于堰塞湖、滑坡、泥石流3種不同災(zāi)害的研究中，且成功實(shí)現(xiàn)了災(zāi)害過(guò)程的重構(gòu)，基于此可完整且準(zhǔn)確地研究地質(zhì)災(zāi)害過(guò)程，并有效地指導(dǎo)地質(zhì)災(zāi)害的監(jiān)測(cè)預(yù)警和快速評(píng)估工作。

1 研究方法

基于地震動(dòng)信號(hào)的災(zāi)害重構(gòu)思路和方法體系如圖1所示。從地震臺(tái)網(wǎng)獲取的地震動(dòng)信號(hào)包含了很多與地質(zhì)災(zāi)害無(wú)關(guān)的地震動(dòng)信號(hào)噪聲，無(wú)法直接用于災(zāi)害過(guò)程分析，故通過(guò)BP-filter對(duì)其進(jìn)行去噪處理；然后利用EMD分離主特征分量；接著進(jìn)行地震動(dòng)信號(hào)時(shí)頻分析：運(yùn)用快速傅里葉變換FFT、短時(shí)傅里葉變換STFT、功率譜密度計(jì)算PSD等方法分析信號(hào)得出地震動(dòng)信號(hào)特征，再結(jié)合現(xiàn)場(chǎng)調(diào)查資料反演出災(zāi)害基本特征，最后通過(guò)數(shù)值模擬實(shí)現(xiàn)災(zāi)害過(guò)程的重構(gòu)。

圖1 基于地震動(dòng)信號(hào)的災(zāi)害過(guò)程重構(gòu)流程圖Fig.1 Flow chart of hazards process reconstruction based on seismic signal process

1.1 地震動(dòng)信號(hào)處理與分析方法

在實(shí)際地震動(dòng)信號(hào)監(jiān)測(cè)過(guò)程中，地震動(dòng)信號(hào)往往存在許多背景噪聲和其他事件引起的噪聲，這些噪聲會(huì)影響分析結(jié)果的準(zhǔn)確性，因此本研究以代替帶通濾波器和經(jīng)驗(yàn)?zāi)B(tài)分解為例，對(duì)原始地震動(dòng)信號(hào)進(jìn)行去噪處理；然后對(duì)獲得的去噪地震動(dòng)信號(hào)進(jìn)行快速傅里葉變換、短時(shí)傅里葉變換、希爾伯特-黃變換、小波變換和功率譜計(jì)算等方法處理，獲取其時(shí)域圖、時(shí)頻譜和能量譜后進(jìn)行分析。本研究以下列方法為例，旨在形成一套完整的地震動(dòng)信號(hào)處理分析方法，對(duì)災(zāi)害過(guò)程進(jìn)行重構(gòu)。

1.1.1 帶通濾波器(BP-filter，Band-Pass filter)

當(dāng)干擾信號(hào)和有效信號(hào)在不同的頻率范圍內(nèi)時(shí)，可以通過(guò)帶通濾波器抑制干擾信號(hào)。通過(guò)識(shí)別有效信號(hào)的頻率范圍，該頻帶內(nèi)的信號(hào)的幅度保持不變，并且該頻帶外的幅度被設(shè)置為0，排除了干擾信號(hào)，得到特定的頻段的災(zāi)害相關(guān)地震動(dòng)信號(hào)(Stein et al.，1983)。

該信號(hào)降噪方法適用于地質(zhì)災(zāi)害發(fā)生過(guò)程中伴隨有其他較強(qiáng)地震動(dòng)信號(hào)影響(其產(chǎn)生的頻率與目標(biāo)災(zāi)害頻率的差異較大)的情況。

1.1.2 經(jīng)驗(yàn)?zāi)B(tài)分解(EMD，Empirical Model Decomposition)

由于地質(zhì)災(zāi)害產(chǎn)生的地震動(dòng)信號(hào)屬于典型的非線性非平穩(wěn)信號(hào)(噪聲影響較大)，而EMD可以自適應(yīng)地將非平穩(wěn)、非線性信號(hào)分解為一組正交的組件，將主特征分量與其余分量分離，得到清晰的時(shí)域曲線，從而實(shí)現(xiàn)地震動(dòng)信號(hào)的降噪(Han et al.，2013)。EMD分解信號(hào)過(guò)程如式(1)所示。

(1)

任何信號(hào)x(t)都能在時(shí)域內(nèi)被EMD分解得到各IMF分量；IMF(Intrinsic Mode Function)代表本征模函數(shù)，每個(gè)IMF分量都包含了各自的頻率和振幅；n為模態(tài)數(shù)；rn(t)為去掉高頻成分的新信號(hào)。

該方法適用于地質(zhì)災(zāi)害發(fā)生過(guò)程中伴隨有其他頻率相近的地震動(dòng)信號(hào)影響的情況。

1.1.3 快速傅里葉變換(FFT，F(xiàn)ast Fourier Transform)

FFT可以獲取信號(hào)的頻譜，分析可得頻率隨時(shí)間的變化、各成分出現(xiàn)的時(shí)間和各時(shí)刻的瞬時(shí)頻率及其幅值等信息。但其缺陷在于對(duì)時(shí)域相差很大的兩個(gè)非平穩(wěn)信號(hào)，其頻譜卻可能一樣(Nussbaumer，1981)。

1.1.4 短時(shí)傅里葉變換(STFT，Short-Time Fourier Transform)

STFT能夠在保留信號(hào)的時(shí)域分布信息的同時(shí)又獲得信號(hào)的頻域信息(Liu et al.，2016)。對(duì)收集到的地震動(dòng)信號(hào)進(jìn)行STFT變換，得到其時(shí)頻譜，然后進(jìn)行時(shí)頻分析，獲取相應(yīng)的時(shí)頻特征。SFTF將信號(hào)進(jìn)行時(shí)頻域轉(zhuǎn)換的變換公式如式(2)所示。

STFT{x(t)}(m，ω)≡X(m，ω)

(2)

式中：x[n]為離散時(shí)間序列；ω[n]為時(shí)頻局部窗口函數(shù)；m為窗口起始時(shí)間；ω為角頻率(ω=2πf)。

功率譜密度(PSD，Power Spectral Density)反映了各時(shí)刻地震動(dòng)信號(hào)一定頻率范圍內(nèi)的能量分布情況(Martin，2001)。將地震動(dòng)信號(hào)分為不同頻帶，分別計(jì)算其功率譜密度，由此可以了解災(zāi)害過(guò)程地震動(dòng)信號(hào)高、中、低頻各成分在災(zāi)害過(guò)程中的統(tǒng)計(jì)特性及不同時(shí)刻的信號(hào)功率分布。

Tsai et al.(2012)提出了一種基于地震動(dòng)信號(hào)的泥沙運(yùn)移PSD評(píng)估模型，計(jì)算公式如式(3)：

(3)

式中：vc為瑞利波的相速度；vu為1iHz組的速度；Vp為石塊的體積；n為石塊的數(shù)量；n/ti為泥沙的流量；χβ為無(wú)量綱參數(shù)，表達(dá)如下：

在此基礎(chǔ)上，Lai et al.(2018)提出了一種簡(jiǎn)化的地震動(dòng)信號(hào)PSD評(píng)估模型，如式(4)。

(4)

式中：ξ為地震動(dòng)速度相關(guān)參數(shù)，其值約為0.25～0.5；r0為地震臺(tái)站到災(zāi)害地之間的平均距離；Q為瑞利波質(zhì)量因素；P是作為f的函數(shù)的地震動(dòng)信號(hào)PSD值；D為顆粒直徑；u為平均流速；L和W分別為運(yùn)動(dòng)物質(zhì)的長(zhǎng)度和寬度。

地震動(dòng)信號(hào)的功率譜密度可以歸納表示為式(5)。

(5)

式中：fmin和fmax分別表示最低和最高頻率；t代表時(shí)間；S(t，f)表示時(shí)頻功率譜。

FFT、STFT、PSD 3種方法基于降噪后獲得的地震動(dòng)信號(hào)，可形成一套完整的信號(hào)處理分析方法，適用于常見地質(zhì)災(zāi)害(滑坡、泥石流、堰塞湖、落石、崩塌等)。

1.2 災(zāi)害過(guò)程重構(gòu)方法

在反演災(zāi)害基本特征的基礎(chǔ)上運(yùn)用數(shù)值模擬方法對(duì)災(zāi)害過(guò)程進(jìn)行模擬，將得到的模擬結(jié)果與基于地震動(dòng)信號(hào)處理與分析后獲取的信號(hào)特征進(jìn)行耦合對(duì)比分析，最終實(shí)現(xiàn)對(duì)整個(gè)災(zāi)害過(guò)程重構(gòu)。常用的數(shù)值模擬方法有：離散元法(DEM，Discrete Element Method)(劉春等，2019；Hu et al.，2020)、有限元法(FEM，F(xiàn)inite Element Method)(鄭穎人等，2005；Zienkiewicz et al.，2005)和光滑粒子流體動(dòng)力學(xué)法(SPH，Smoothed Particle Hydrodynamics)(Liu et al.，2010；劉謀斌等，2011)等。

2 研究方法應(yīng)用

2.1 堰塞湖應(yīng)用

通過(guò)對(duì)堰塞湖潰決過(guò)程中提取到的地震動(dòng)信號(hào)進(jìn)行處理與分析，系統(tǒng)地研究堰塞湖潰決過(guò)程，從而獲得的災(zāi)害過(guò)程特征服務(wù)于堰塞湖的監(jiān)測(cè)預(yù)警工作。

本研究以中興大學(xué)在中國(guó)臺(tái)灣南投縣惠蓀林場(chǎng)的蘭島溪原位實(shí)驗(yàn)為例(Feng，2012；Yan et al.，2017)。蘭島溪上游修筑有一條數(shù)百米長(zhǎng)的排砂道，實(shí)驗(yàn)過(guò)程中，通過(guò)調(diào)控上游來(lái)水，進(jìn)行堰塞壩潰決觀察實(shí)驗(yàn)及其相關(guān)研究，實(shí)驗(yàn)區(qū)域位置及傳感器布置情況如圖2所示。

圖2 堰塞湖潰決原位實(shí)驗(yàn)概況圖Fig.2 Overview of large-scale dammed lake bursting physical modellinga.實(shí)驗(yàn)區(qū)域示意圖；b.傳感器布設(shè)示意圖；c.上游堰塞壩體與下游堰塞壩體正視圖

為了去除背景噪聲和其他事件干擾噪聲的影響，本研究以BP-filter對(duì)采集到的地震動(dòng)信號(hào)進(jìn)行去噪處理，得到去噪的地震動(dòng)信號(hào)后進(jìn)行EMD處理。EMD可以自適應(yīng)地將堰塞湖潰決的地震動(dòng)信號(hào)以及砂石撞擊、水流沖擊等噪聲分離出來(lái)，便于我們單獨(dú)對(duì)各個(gè)分量進(jìn)行分析，故對(duì)上述去噪后信號(hào)進(jìn)行EMD處理，根據(jù)信號(hào)的殘差情況選擇IMF分量數(shù)量，得到結(jié)果如圖3所示。

圖3 EMD處理后得到的IMF分量圖Fig.3 IMF component diagram obtained after EMD processing

通過(guò)進(jìn)行FFT處理，可以生成地震動(dòng)信號(hào)的頻譜，以便我們獲取其中各個(gè)頻率成分及其能量，進(jìn)而解釋過(guò)程中發(fā)生的事件，以及這些事件對(duì)地震動(dòng)帶來(lái)的影響及其大小。本次FFT處理后得到的頻譜中最大頻率為50Hz，如圖4所示。

圖4 FFT處理后得到的頻譜Fig.4 Spectrum obtained after FFT processing

運(yùn)用STFT對(duì)去噪信號(hào)進(jìn)行時(shí)域-頻域聯(lián)合變換，得出本次實(shí)驗(yàn)采集的地震動(dòng)信號(hào)范圍為0～50Hz，如圖5a所示。將0～10Hz和10～45Hz范圍內(nèi)的時(shí)頻譜進(jìn)行放大，如圖5b、圖5c所示。經(jīng)分析后可按頻率由低到高將地震動(dòng)信號(hào)劃分為3個(gè)主要部分：即低頻(0～1.5Hz)、中頻(1.5～10Hz)、高頻(10～45Hz)。

圖5 地震動(dòng)信號(hào)的STFT時(shí)頻譜Fig.5 STFT time-frequency spectrum of seismic signal

基于堰塞湖潰決的3個(gè)階段(潰決之前、潰決時(shí)、潰決之后趨于穩(wěn)定)，對(duì)更大范圍內(nèi)的時(shí)間段信號(hào)進(jìn)行上述處理，最終劃分時(shí)間段后得到3個(gè)階段的時(shí)頻譜如圖6所示。

圖6 堰塞湖潰決過(guò)程STFT時(shí)頻譜(0～50Hz)Fig.6 STFT time-frequency spectrum in the process of dammed lake failurea.第1階段(潰決之前)；b.第2階段(潰決時(shí))；c.第3階段(潰決之后趨于穩(wěn)定)

由圖6a可得：第1階段能量集中在2350～2500s期間，2500s之后信號(hào)的能量成分大幅度衰減，開始出現(xiàn)時(shí)頻譜中的雜亂特征，頻率寬度集中于10～40Hz，但隨著時(shí)間的推移，頻寬移動(dòng)向低頻，由此可推斷該階段對(duì)應(yīng)堰塞湖潰決前水流匯聚于上游壩體且水位不斷升高至漫頂過(guò)程中產(chǎn)生的地震動(dòng)信號(hào)。由圖6b可得：隨著災(zāi)害過(guò)程的演進(jìn)，時(shí)頻譜整體已呈現(xiàn)出雜亂的特征，并且主要能量所在的最高頻率由30Hz升至40Hz，最后又降為30Hz(流速增大，頻寬變寬)，基于此變化，推斷出2650～3100s期間內(nèi)出現(xiàn)了一次促使高頻成分增加的事件，由此可確定該階段對(duì)應(yīng)于堰塞湖潰決時(shí)水流通過(guò)潰口向下游泄水并侵蝕潰口過(guò)程中產(chǎn)生的地震動(dòng)信號(hào)。由圖6c可得：該階段相比于前兩階段能量較低，雖然時(shí)頻譜同樣呈現(xiàn)雜亂的特征，但頻寬卻很窄，且能量隨著時(shí)間的推移還有繼續(xù)減弱的趨勢(shì)，最終趨于穩(wěn)定狀態(tài)，由此推斷最后階段應(yīng)為堰塞湖潰決后湖面趨于穩(wěn)定情況下水流動(dòng)所產(chǎn)生的地震動(dòng)信號(hào)。

在上述堰塞湖潰決3個(gè)過(guò)程中的時(shí)頻分析的基礎(chǔ)上，對(duì)低頻(0～1.5Hz)、中頻(1.5～10Hz)、高頻(10～45Hz)3個(gè)主要部分的信號(hào)源進(jìn)行分析，結(jié)果如下：

對(duì)于低頻成分(0～1.5Hz)，結(jié)合視頻記錄的資料可得：上游與下游潰壩過(guò)程時(shí)段如圖7a所示，上游壩體的大片土體從壩體脫離并滑落至河床，同時(shí)其潰口在2708～2722s和2804～2851s這兩個(gè)時(shí)間段內(nèi)發(fā)生了多次大規(guī)模的壩體崩塌事件，而下游壩體潰口發(fā)生的較大規(guī)模壩體崩塌事件集中在2851～2947s。從圖7a中可看出發(fā)生壩體崩塌事件在低頻信號(hào)下對(duì)應(yīng)的時(shí)間段呈現(xiàn)出高于穩(wěn)定狀態(tài)且各不相同的地震動(dòng)能量。由此證明壩體潰決過(guò)程中土體崩塌事件引起的地震動(dòng)造成了本實(shí)驗(yàn)地震動(dòng)信號(hào)中低頻成分(0～1.5Hz)的產(chǎn)生?；诖?，當(dāng)利用地震動(dòng)信號(hào)對(duì)其他堰塞湖災(zāi)害事件進(jìn)行分析時(shí)，若低頻成分出現(xiàn)上述情況，則可判斷對(duì)應(yīng)時(shí)段發(fā)生多次規(guī)模較大壩體崩塌事件。

對(duì)于中頻成分(1.5～10Hz)，根據(jù)本實(shí)驗(yàn)中記錄的視頻影響資料發(fā)現(xiàn)：一共出現(xiàn)3次石塊滑落事件(圖7b)，其中上游壩體石塊A、石塊B和下游壩體石塊C分別從壩頂?shù)袈洳_擊河床。由圖7b中可明顯看出，對(duì)應(yīng)出現(xiàn)了3次能量激增的現(xiàn)象，與視頻影像對(duì)比后發(fā)現(xiàn)時(shí)間吻合，由此可證明能量的短時(shí)激增是由于石塊滑落所引起。故本方法體系應(yīng)用于其他堰塞湖事件時(shí)，若中頻成分中某時(shí)刻出現(xiàn)能量激增現(xiàn)象，則可推斷此時(shí)發(fā)生壩體石塊滑落事件。

圖7 地震動(dòng)信號(hào)各頻帶能量與影像對(duì)比圖Fig.7 Energy of each frequency band of seismic signal and corresponding site picturesa.低頻成分(0～1.5Hz)能量與影像對(duì)比圖；b.中頻成分(1.5～10Hz)能量與影像對(duì)比圖；c～d.高頻成分(10～45Hz)功率譜密度與水位變化對(duì)比圖(顏色由藍(lán)色變?yōu)辄S色表示取樣點(diǎn)的時(shí)間先后順序)

對(duì)于高頻成分(10～45Hz)，我們分頻段(10～30Hz和30～45Hz)計(jì)算信號(hào)的功率譜密度(PSD)，并與水位進(jìn)行關(guān)聯(lián)性分析(圖7c～7d)。由圖7c～7d可得：前期PSD能量隨著水位的增加而增加，對(duì)應(yīng)上游來(lái)水進(jìn)入上游壩，水流以較快的速度在河床上流動(dòng)使砂礫之間產(chǎn)生較強(qiáng)的碰撞，推移質(zhì)運(yùn)動(dòng)加強(qiáng)，從而導(dǎo)致地震動(dòng)信號(hào)10～45Hz成分在該時(shí)間段的頻譜能量較強(qiáng)(圖6a)；后期PSD能量隨著水位的降低而降低，對(duì)應(yīng)水流速度隨著水位的降低相應(yīng)逐漸降低的現(xiàn)象，這導(dǎo)致了潰口推移質(zhì)運(yùn)動(dòng)逐漸減弱，10～45Hz成分能量因此持續(xù)下降。上述結(jié)果表明了堰塞湖潰決的地震動(dòng)信號(hào)高頻成分能量大小與其過(guò)程中水位變化的關(guān)系成正比，并進(jìn)一步推斷高頻成分的變化與推移質(zhì)的運(yùn)動(dòng)相關(guān)?；诖耍?dāng)對(duì)其他堰塞湖潰決災(zāi)害事件進(jìn)行分析時(shí)，通過(guò)其高頻成分的能量變化和對(duì)水位變化的分析，可推斷堰塞湖潰決過(guò)程的演進(jìn)情況。

綜上所述，將本研究思路和方法體系應(yīng)用于其他堰塞湖潰決災(zāi)害研究時(shí)，在根據(jù)上述堰塞湖地震動(dòng)信號(hào)的分析處理得到其信號(hào)特征后，便能反演出堰塞湖潰決全過(guò)程的基本特征，同時(shí)獲取各信號(hào)源成分等信息，基于這些特征信息，再結(jié)合數(shù)值模擬方法，就能實(shí)現(xiàn)對(duì)堰塞湖潰決過(guò)程的重構(gòu)。

2.2 滑坡應(yīng)用

滑坡屬于常見的地質(zhì)災(zāi)害，影響范圍廣且造成損失大，故本文選取兩次滑坡事件進(jìn)行研究，分析其發(fā)生過(guò)程中產(chǎn)生的地震動(dòng)信號(hào)，實(shí)現(xiàn)滑坡過(guò)程的重構(gòu)，從而能更好地對(duì)滑坡災(zāi)害進(jìn)行監(jiān)測(cè)、預(yù)警和評(píng)估。

本研究以兩次滑坡事件為例，分別為：“6·24”茂縣疊溪鎮(zhèn)新磨村高速遠(yuǎn)程滑坡(Fan et al.，2017；Meng et al.，2018；Yan et al.，2020a)和“7·23”水城縣山體滑坡事件(王立朝，2019；Yan et al.，2020b；鄭光等，2020)。兩起滑坡事件的位置及其附近地震臺(tái)網(wǎng)的分布如圖8a～圖8b所示。研究采用EMD對(duì)兩起滑坡事件附近的地震臺(tái)站的地震動(dòng)信號(hào)進(jìn)行降噪處理，過(guò)濾信號(hào)噪音并提取滑坡地震動(dòng)信號(hào)的IMF分量；采用FFT和STFT方法分析了滑坡地震動(dòng)信號(hào)的時(shí)域特征和時(shí)頻域特征；最后結(jié)合數(shù)值模擬方法重構(gòu)了滑坡的整個(gè)過(guò)程，確定了各階段滑坡體的運(yùn)動(dòng)特征等信息。

圖8 滑坡位置概況示意圖Fig.8 Schematic diagram of landslide location overviewa.新磨滑坡地理位置及地震臺(tái)站分布圖；b.水城滑坡地理位置及地震臺(tái)站分布圖

通過(guò)EMD處理后得新磨滑坡中MXI臺(tái)站的IMF4分量如圖9a所示，由圖中兩次信號(hào)振幅波動(dòng)(①與②)可推斷出現(xiàn)了兩次地震動(dòng)事件。水城滑坡中XUW臺(tái)站的IMF1分量如圖10a所示，基于振幅波動(dòng)推斷只出現(xiàn)了一次地震動(dòng)事件。分別進(jìn)行FFT變換，得到新磨滑坡和水城滑坡的頻譜分別如圖9b、圖10b所示(由于水城滑坡中信號(hào)在4.5～6Hz之間出現(xiàn)較大異常，其信號(hào)特征最為明顯，故選取該頻段的地震動(dòng)信號(hào)進(jìn)行分析)，根據(jù)圖中信息便可得出能量主要集中的頻段。通過(guò)STFT獲得兩次滑坡地震動(dòng)信號(hào)的時(shí)頻譜，分析可知新磨滑坡兩次主信號(hào)頻率范圍在0～2Hz與0～4Hz之間(圖9c)，而水城滑坡的地震動(dòng)信號(hào)時(shí)頻譜顯示能量峰值呈現(xiàn)“W”形變化(圖10c)。我們選取水城滑坡的地震動(dòng)信號(hào)進(jìn)行PSD計(jì)算，得到結(jié)果如圖10d所示，根據(jù)圖中PSD曲線的振幅大小和趨勢(shì)，便可估算出滑坡的速度和能量大小(成正比關(guān)系)。

圖9 新磨滑坡地震動(dòng)信號(hào)處理結(jié)果圖Fig.9 Seismic signal processing result figure of Xinmo landslidea.EMD分解后得到的MXI臺(tái)站IMF4成分結(jié)果圖；b.MXI臺(tái)站FFT頻譜；c.MXI臺(tái)站STFT時(shí)頻譜

圖10 水城滑坡地震動(dòng)信號(hào)處理結(jié)果圖(4.5～6Hz)Fig.10 Seismic signal processing result diagram of Shuicheng landslide(4.5～6Hz)a.EMD處理后得到的XUW臺(tái)站IMF1成分結(jié)果圖；b.XUW臺(tái)站FFT頻譜；c.XUW臺(tái)站STFT時(shí)頻譜；d.XUW臺(tái)站PSD功率譜密度圖

基于上述地震動(dòng)信號(hào)的處理與分析，將兩次滑坡事件各自獲得的地震動(dòng)信號(hào)特征與數(shù)值模擬分別進(jìn)行耦合分析，對(duì)災(zāi)害過(guò)程進(jìn)行重構(gòu)，得到主要結(jié)論如下：

(1)新磨滑坡過(guò)程大致可分為5個(gè)階段：靜止階段、滑移階段、過(guò)渡階段、夾帶運(yùn)移階段和沉積階段(圖9a Ⅰ～Ⅴ)，其中首次從地震動(dòng)特征中識(shí)別出過(guò)渡階段；由時(shí)頻譜(圖9c)中可見，過(guò)渡階段的低頻特性介于兩個(gè)高頻特征之間(分別代表滑移階段和夾帶運(yùn)移階段)。結(jié)合數(shù)值模擬結(jié)果，分析發(fā)現(xiàn)過(guò)渡階段的主要運(yùn)動(dòng)特征是滑坡整體速度的相對(duì)穩(wěn)定和坡積物開始被物源土體刮鏟，與之前研究得出的結(jié)論不同。

(2)水城滑坡過(guò)程大致可分為4個(gè)階段：斷裂過(guò)渡階段、加速啟動(dòng)階段、分叉刮鏟階段和沉積階段(圖10a Ⅰ～Ⅳ)。本研究首次確定出滑坡初始斷裂點(diǎn)并非位于滑源區(qū)，而是位于滑動(dòng)面的上部的位置。

2.3 泥石流應(yīng)用

泥石流演進(jìn)過(guò)程與地震動(dòng)信號(hào)具有很高的相關(guān)性，在可能發(fā)生泥石流的流通區(qū)安裝地震動(dòng)信號(hào)監(jiān)測(cè)設(shè)備，實(shí)時(shí)監(jiān)控，通過(guò)分析地震動(dòng)信號(hào)來(lái)進(jìn)行預(yù)警，可降低下游因泥石流造成的經(jīng)濟(jì)損失和人員傷亡。

本研究以2014年6月24日記錄的東川蔣家溝泥石流為例，研究其災(zāi)害演進(jìn)過(guò)程。通過(guò)安裝在蔣家溝岸邊的地震動(dòng)信號(hào)監(jiān)測(cè)設(shè)備記錄了災(zāi)害演進(jìn)期間觀測(cè)到的泥石流(圖11a)。監(jiān)測(cè)儀器包括MEMS三軸數(shù)字加速度計(jì)芯片及地震動(dòng)信號(hào)監(jiān)測(cè)探頭(圖11b)。以加速度計(jì)作為傳感元件，記錄泥石流經(jīng)過(guò)岸邊處的重力加速度并傳回觀測(cè)站。

圖11 東川泥石流地震動(dòng)信號(hào)監(jiān)測(cè)位置及時(shí)頻譜Fig.11 Monitoring location and time-frequency spectrum of Dongchuan debris flow seismic signala.地震動(dòng)信號(hào)監(jiān)測(cè)位置照片；b.監(jiān)測(cè)設(shè)備現(xiàn)場(chǎng)布置照片；c.地震動(dòng)信號(hào)時(shí)域曲線；d.地震動(dòng)信號(hào)FFT頻譜；e.地震動(dòng)信號(hào)STFT時(shí)頻譜

采用FFT、STFT對(duì)地震動(dòng)信號(hào)整段進(jìn)行分析，獲取信號(hào)的時(shí)域圖和時(shí)頻譜。由時(shí)域圖中信息可知，重力加速度異常僅出現(xiàn)在第30s左右且持續(xù)時(shí)間約1s(圖11c)，信號(hào)包括兩段呈紡錘狀的強(qiáng)信號(hào)段，即信號(hào)先呈現(xiàn)增強(qiáng)，達(dá)到最大值后減弱，然后逐漸恢復(fù)至正常狀態(tài)。研究發(fā)現(xiàn)該信號(hào)包含有兩次泥石流事件的疊加，兩次事件間隔時(shí)間為0.6s(圖11c)。FFT分析的結(jié)果表明信號(hào)在頻譜中呈現(xiàn)出能量的全頻段分布(0～200Hz)，優(yōu)勢(shì)段(70～200Hz)具有明顯的特征能量較強(qiáng)(圖11d)。由STFT獲得的時(shí)頻譜也顯示出相似的變化規(guī)律，在持續(xù)時(shí)間1s內(nèi)有較強(qiáng)的地震動(dòng)信號(hào)變化，而未記錄到信號(hào)異常的時(shí)間段內(nèi)信號(hào)頻譜呈白噪化(圖11e)。

通過(guò)FFT與STFT處理得到的時(shí)頻特征，再結(jié)合東川站觀測(cè)的泥石流資料進(jìn)行對(duì)比分析，得到下列主要結(jié)論，最終實(shí)現(xiàn)了對(duì)其災(zāi)害過(guò)程的重構(gòu)。

(1)根據(jù)圖11c中出現(xiàn)的較大重力加速度異常現(xiàn)象，可推斷流體對(duì)所流經(jīng)的監(jiān)測(cè)斷面基巖施加了一個(gè)較大的作用力，結(jié)合現(xiàn)場(chǎng)記錄資料可得本次泥石流為黏性泥石流且攜帶顆粒物含量大。

(2)通過(guò)圖11c中出現(xiàn)的重力加速度異常值快速變化，衰減迅速的現(xiàn)象，可推斷流體流過(guò)監(jiān)測(cè)斷面的速度很快，對(duì)監(jiān)測(cè)斷面施加的作用力時(shí)間極短，結(jié)合監(jiān)測(cè)資料得出本次泥石流的規(guī)模較小且為斷斷續(xù)續(xù)的陣性流。

3 結(jié)論與展望

(1)本研究將BP-filter、EMD、FFT、STFT、PSD等信號(hào)處理方法聯(lián)合運(yùn)用，并結(jié)合現(xiàn)場(chǎng)調(diào)查資料和數(shù)值模擬最終形成一套系統(tǒng)且完整的基于地震動(dòng)信號(hào)重構(gòu)災(zāi)害的研究方法和體系，將其應(yīng)用于上述堰塞湖、滑坡和泥石流3種典型地質(zhì)災(zāi)害得出的結(jié)果與真實(shí)情況具有較高的吻合度，證明了本方法的可靠性與可行性。

(2)本研究提出的基于地震動(dòng)信號(hào)重構(gòu)地質(zhì)災(zāi)害的方法和思路體系，應(yīng)用于上述3種常見地質(zhì)災(zāi)害分析其過(guò)程，得出主要結(jié)論如下：

1)堰塞湖：低、中、高頻成分的信號(hào)源分別來(lái)自：壩體崩塌滑坡、壩體石塊滑落、推移質(zhì)的運(yùn)動(dòng)；在此基礎(chǔ)上經(jīng)探究后得出堰塞湖壩體潰決先兆特征。

2)滑坡：將新磨滑坡劃分為5個(gè)階段，首次從地震動(dòng)信號(hào)分析中識(shí)別出過(guò)渡階段；水城滑坡劃分為4個(gè)階段，首次確定了滑坡初始斷裂點(diǎn)位于滑動(dòng)面的上部，而并非處于滑源區(qū)。

3)泥石流：獲取了強(qiáng)信號(hào)段呈紡錘體狀、時(shí)頻譜中能量呈現(xiàn)全頻段分布、優(yōu)勢(shì)頻段(70～200Hz)能量較強(qiáng)、重力加速度異常值大而持續(xù)時(shí)間短等特征信息，并在此基礎(chǔ)上推斷出此次泥石流呈現(xiàn)黏性、陣性流等特征。

(3)我們?cè)诤罄m(xù)的研究中會(huì)逐步減少對(duì)現(xiàn)場(chǎng)調(diào)查資料的依賴，甚至無(wú)需現(xiàn)場(chǎng)調(diào)查，僅通過(guò)結(jié)合地震動(dòng)信號(hào)的處理分析和數(shù)值模擬就能形成一套完整的方法體系來(lái)系統(tǒng)地重構(gòu)地質(zhì)災(zāi)害過(guò)程，提升災(zāi)害重構(gòu)效率的同時(shí)還能保證其結(jié)果的準(zhǔn)確性。此外，將本方法體系應(yīng)用于不同地質(zhì)災(zāi)害時(shí)，能準(zhǔn)確地獲得災(zāi)害發(fā)生全過(guò)程的一些特征信息，基于此能有助于對(duì)地質(zhì)災(zāi)害的監(jiān)測(cè)預(yù)警等工作的指導(dǎo)。