楊長衛(wèi)，張凱文，吳東升，張志方，張 良，瞿立明

(1.西南交通大學(xué) 土木工程學(xué)院，四川 成都 610031；2.中國鐵路武漢局集團(tuán)有限公司，湖北 武漢 430071；3.西南交通大學(xué) 地球科學(xué)與環(huán)境工程學(xué)院，四川 成都 611756)

截至2022年底,中國高鐵運營里程突破4.2萬km,占據(jù)世界高鐵運營里程2/3以上[1]。高速鐵路的運行很大程度地緩解了我國交通運營壓力,為我國經(jīng)濟(jì)發(fā)展做出了巨大貢獻(xiàn)。地震作為一種隨機(jī)性強(qiáng)、破環(huán)性大的自然災(zāi)害,對高速鐵路的運營造成了極大的安全隱患[2]。運行中的高速鐵路列車在地震作用下可能有脫軌、碰撞等危害性極強(qiáng)的次生災(zāi)害發(fā)生,造成嚴(yán)重的人員傷亡及經(jīng)濟(jì)損失,許多研究學(xué)者研究了地震作用下的鐵路橋梁、軌道的運行安全[3-6]。我國在2017年建成了高鐵地震預(yù)警系統(tǒng)[7],該系統(tǒng)主要利用P波和S波、S波和電磁波的速度差來實現(xiàn)的[8-9]。地震P波自動識別作為高速鐵路地震預(yù)警系統(tǒng)的首要工作,其識別地震P波的精度與速度直接影響了后續(xù)地震預(yù)警系統(tǒng)工作震中定位[10]及震級估算[11-12]的準(zhǔn)確性。研究適應(yīng)我國高速鐵路地震預(yù)警的快速、精準(zhǔn)地震P波震相識別算法十分必要,能夠有效保障我國高速鐵路的安全運營。

數(shù)十年來,大量研究學(xué)者研究了提高地震P波自動識別的速度與精度的方法,主要形成了基于地震波長時窗均值與短時窗均值之比(short time average/ long time average,STA/LTA)[13-15]與基于赤池信息法則(Akaike information criterion,AIC)[16-17]的地震波識別方法。此外,也有國內(nèi)外研究學(xué)者提出了不同的地震P波震相初至識別方法,包括相鄰道互相關(guān)法[18]、分形維數(shù)方法[19]、小波變換方法[20-21]、數(shù)字圖像分割法[22]、CM及MCM算法[23]、人工智能算法[24-25]等。上述方法各有特點,但因存在局限性沒有作為地震P波的震相識別主流算法。已有研究結(jié)果顯示,STA/LTA方法在地震波振幅較大時的工作效果良好,AIC方法的識別的精準(zhǔn)度較高。STA/LTA方法由于工作原理方法,在地震P波振幅較小時的識別精度較差;而AIC方法通過建立自回歸模型,識別時間窗口長度內(nèi)的函數(shù)最小值作為地震波的精準(zhǔn)到達(dá)時刻,所需要的時間較長。據(jù)統(tǒng)計,僅在四川境內(nèi),2009年1月1日至2022年2月12日共發(fā)生3級以上地震 1 084次,其中5級以下的地震1 040次,占比95.94%。在小震頻發(fā)的背景下,為確保高速鐵路地震預(yù)警系統(tǒng)的良好工作及高速鐵路的安全運營,研究快速精準(zhǔn)的小震P波自動識別算法十分必要。分形理論應(yīng)用在地震P波震相識別最早由Boshetti在1996提出的,后續(xù)不斷有學(xué)者對分形理論應(yīng)用在地震波識別的方法進(jìn)行研究[26-28],提出STAFD/LTAFD、灰度邊緣識別等方法,并得出微震震相識別分形理論具有更高精確度的結(jié)論。該理論由于分形維數(shù)計算較慢的缺點,許多研究學(xué)者將計算一次分形維數(shù)的時間間隔設(shè)置為5個采樣點及以上以提高識別速度,識別精度仍有進(jìn)一步提高的空間。

本文在已有研究學(xué)者的基礎(chǔ)上,對分形維數(shù)應(yīng)用在地震P波震相識別進(jìn)行了進(jìn)一步的研究,優(yōu)化了分形維數(shù)的計算方法,提高了分形維數(shù)的計算速度,提升了識別精度。此外,本文在時效性和準(zhǔn)確性上與STA/LTA方法及AIC方法進(jìn)行對比。結(jié)果表明,在準(zhǔn)確性方面,本文提出方法的平均識別誤差達(dá)到了0.006 3 s,優(yōu)于STA/LTA方法且與AIC方法基本持平(<0.01 s)。在時效性方面,本文提出方法的平均識別時長為0.16 s,優(yōu)于AIC方法。這些結(jié)果滿足我國高速鐵路地震預(yù)警系統(tǒng)的要求,以期為我國高速鐵路地震預(yù)警系統(tǒng)提供參考并作為補(bǔ)充算法。

1 基于分形斜率的地震P波識別算法

1.1 分形理論

分形理論最早由美國數(shù)學(xué)家Mandelbrot在1975年提出。自此,分形理論被廣泛應(yīng)用解決不同的實際問題。分形理論改變了傳統(tǒng)的維數(shù)觀念,如點是一維的,面是二維的,體是三維的。在分形理論中,曲線的維數(shù)被認(rèn)為是1～2維的,而不規(guī)則平面的維數(shù)被認(rèn)為是2～3維的,即引入線、面的“粗糙度”的概念。分形理論在地震波震相識別的核心理念是,通過識別地震波到達(dá)前后曲線粗糙度變化來識別地震波的精準(zhǔn)到達(dá)時間。曲線分形維數(shù)的計算方法通常包括尺碼法和網(wǎng)格法,已有研究成果表明尺碼法更適合應(yīng)用在地震P波自動震相識別[29],本文同樣采用尺碼法進(jìn)行分形維數(shù)計算。尺碼法計算流程見圖1(a),其核心計算原理為對一段曲線,移動不同的尺碼ri覆蓋整條曲線,共移動了Ni次,則可得到曲線的近似長度Li,對不同的尺碼ri和近似長度Li取對數(shù)進(jìn)行一次擬合,見圖1(c),得到斜率F,則分形維數(shù)D計算式為

D=1-F

(1)

式中:D為分形維數(shù);F為擬合一次斜率。

在計算分形維數(shù)時,應(yīng)該保證在移動尺碼時,尺碼的兩端都在曲線上,見圖1(b)。同時,為得到曲線的近似長度Li,需要多次改變尺碼ri的長度,通常需要改變十余次尺碼長度才可以完成一次分形維數(shù)計算,這一行為嚴(yán)重降低了分形維數(shù)的計算速度。為提高分形維數(shù)的計算速度,改變尺碼ri長度的次數(shù)被降低至5次,減少了分形維數(shù)計算工作量,同時并沒有影響分形維數(shù)曲線的性質(zhì)。結(jié)果表明,在降低尺碼變化次數(shù)后,一次分形維數(shù)的計算速度由0.27 s提高至0.05 s。

1.2 計算原理

地震預(yù)警系統(tǒng)對自動震相識別算法的最大要求在于快速及時準(zhǔn)確。很多研究學(xué)者為實現(xiàn)分形理論在地震預(yù)警系統(tǒng)實現(xiàn)的可能,選擇間隔數(shù)個采樣點進(jìn)行一次分形維數(shù)計算,然而這一做法損失了識別的精度。基于改進(jìn)的分形維數(shù)算法,本文提出了連續(xù)采樣點計算分形維數(shù)的方法。連續(xù)地震波分形維數(shù)具體計算方法是選擇一定長度的時間窗口覆蓋地震波,計算一次分形維數(shù),并將分形維數(shù)標(biāo)記在時間窗口的右側(cè)時刻。隨著時間窗口的逐步移動,可以在地震波的每個時刻標(biāo)記上對應(yīng)分形維數(shù),即曲線的粗糙度,進(jìn)而形成連續(xù)的分形維數(shù)曲線。該方法在不損失計算速度的情況下,得到連續(xù)分形維數(shù)曲線,充分反映了地震波的細(xì)節(jié)特征,從而得到更加精準(zhǔn)的地震P波到達(dá)時間。分形維數(shù)連續(xù)與間斷計算曲線見圖2。圖2中,顯示了間隔5個采樣點計算一次分形維數(shù)得到的分形維數(shù)曲線與連續(xù)分形維數(shù)曲線的區(qū)別。

圖2 分形維數(shù)連續(xù)與間斷計算曲線

計算分形維數(shù)的時間窗口沿著地震波曲線移動,地震波可以被劃分為3個階段,見圖3(a),同樣地,在分形維數(shù)曲線中3個階段也有明顯的特征,見圖3(b)。第一階段為地震波到達(dá)之前,時間窗口內(nèi)的地震波信號僅有震前白噪聲組成,因此分形維數(shù)曲線十分平坦,記該段為CFD1;第二階段為地震波初至,時間窗口內(nèi)的地震波信號由部分震前白噪聲及地震波與白噪聲組成,曲線的粗糙度快速上升,因此分形維數(shù)曲線也快速上升,記該段為CFD2;第三階段為地震波完全到達(dá),時間窗口內(nèi)的地震波信號僅有地震波與白噪聲組成,分形維數(shù)曲線反映了地震波變化的細(xì)節(jié)特征,開始無規(guī)則波動,記該段為CFD3。其中第一階段與第二階段的臨界點記為點1,第二階段與第三階段的臨界點記為點2。

圖3 分形斜率識別原理

本文進(jìn)一步引入了分形斜率Ki去分析分形維數(shù)曲線的變化趨勢,其中斜率的計算式為

(2)

式中:Di為i時刻的分形維數(shù);t為時間窗口長度。

分形斜率即一個時間窗口的左右兩點的分形維數(shù)形成的斜率,反應(yīng)了分形維數(shù)曲線隨著時間窗口逐步移動的變化特征,進(jìn)而形成分形斜率曲線,進(jìn)一步反映了地震波的波動特點,分形斜率曲線可以劃分為四個階段,見圖3(c)。第一階段,分形斜率十分平坦,此階段時間窗口的左點與右點均在CFD1上,兩點的分形維數(shù)基本持平,從而分形斜率也基本持平;第二階段,分形斜率快速上升,此階段時間窗口的左點在CFD1上、右點在CFD2上,隨著時間窗口的移動,左點保持不變,右點快速上升,分形斜率也快速上升并上升;第三階段,分形斜率快速下降,此階段時間窗口的左點在CFD2上、右點在CFD3上,隨著時間窗口的移動,左點快速上升,而右點進(jìn)入波動段,從而分形斜率快速下降;第四階段,分形斜率不斷波動,此階段時間窗口的左點和右點均在CFD3上,反映了分形維數(shù)與地震波信號的變化特征。同樣地,臨界點1、臨界點2在分形斜率曲線中也是第一階段與第二階段、第二階段與第三階段的臨界點2。在分形維數(shù)斜率第二階段的末尾時刻,時間窗口的左點剛好在臨界點1上,右點剛好在臨界點2上。若時間窗口繼續(xù)推進(jìn),由于時間窗口左點落入到CFD2上,分形維數(shù)開始上升,則分形斜率開始下降,將該臨界點2時刻定義為極值時刻,而分形斜率也有極值的表現(xiàn)。本文通過識別分形斜率極值判斷地震波的精準(zhǔn)到達(dá)時刻。若分形斜率極值時刻為i,則地震波的精準(zhǔn)到達(dá)時刻為i-t。

本文通過設(shè)置極值時刻判斷標(biāo)準(zhǔn)和引入標(biāo)準(zhǔn)差去提高識別精度。首先,分形斜率極值時刻應(yīng)該大于相鄰4個采樣的分形斜率。其次,計算極值點時刻前15個采樣點至215個采樣點的分形斜率標(biāo)準(zhǔn)差,即CFD1階段的標(biāo)準(zhǔn)差,評估該階段的分形斜率曲線的波動程度。最后,要求極值時刻的分形斜率應(yīng)該大于分形斜率第一階段最大值數(shù)倍及以上。標(biāo)準(zhǔn)差的計算式為

(3)

根據(jù)震前白噪聲的分形斜率波動情況,可以將其分為3個波動程度水平,分別為:標(biāo)準(zhǔn)差s小于0.002、大于0.02和0.002～0.02之間。為進(jìn)一步減少識別誤差,當(dāng)標(biāo)準(zhǔn)差s小于0.002時,分形斜率的極值時刻應(yīng)該大于分形斜率第一階段最大值的6倍及以上;當(dāng)標(biāo)準(zhǔn)差s在0.002～0.02之間時,分形斜率的極值時刻應(yīng)該大于分形斜率第一階段最大值的4倍及以上;當(dāng)標(biāo)準(zhǔn)差s大于0.02時,分形斜率的極值時刻應(yīng)該大于分形斜率第一階段最大值的2倍及以上。本文提出的方法的具體識別流程見圖4。

圖4 識別流程

1.3 影響因素分析

時間窗長度的選擇直接影響了分形維數(shù)地震P波震相識別工作的速度與精度。時間窗長度過長導(dǎo)致分形維數(shù)計算速度降低,時間窗長度過短則導(dǎo)致地震波加速度曲線波動過于敏感容易誤觸發(fā)。為保證高速鐵路地震預(yù)警系統(tǒng)工作的時效性與準(zhǔn)確度,有必要研究確定合適的時間窗長度。

本文選取了100組日本K-Net &Kik-net地震臺網(wǎng)的地震數(shù)據(jù)進(jìn)行分析,其地震事件的震級均小于M 5,震中距小于100 km,采樣頻率為100 Hz,地震波方向均為UD方向。在進(jìn)行對比時,認(rèn)為人工檢測地震P波到達(dá)時間是精準(zhǔn)的。

本文提出的地震P波自動識別算法的精準(zhǔn)到達(dá)時刻為極值時刻減去時間窗長度,為了滿足高速鐵路地震預(yù)警系統(tǒng)的時效性,時間窗口長度的選擇不亦過大,保證得到地震波精準(zhǔn)到達(dá)時刻的速度滿足Q/CR 634—2018《高速鐵路地震預(yù)警監(jiān)測系統(tǒng)技術(shù)條件》[30]的技術(shù)要求。總共4個時間窗口長度被選擇進(jìn)行分形維數(shù)、分形斜率計算并進(jìn)行識別,包括5、10、15、20個采樣點,即0.05、0.10、0.15、0.20 s。按照本文提出的識別方法,可以得到4種不同時間窗口自動識別得到的地震波震相到達(dá)時刻,其識別結(jié)果與人工檢測的誤差統(tǒng)計見圖5。結(jié)果顯示,不同時間窗口長度下的識別誤差基本保持在±0.1 s以內(nèi),但也存在一定差異。不同時間窗口長度的識別平均誤差均在±0.05 s以內(nèi),以不同時間窗口長度識別標(biāo)準(zhǔn)差作為標(biāo)準(zhǔn)能夠更好的評價識別效果。顯然,15個采樣點的識別結(jié)果標(biāo)準(zhǔn)差更加優(yōu)異。

圖5 不同時間窗口長度識別誤差直方圖

2 案例分析

STA/LTA與AIC方法作為地震預(yù)警系統(tǒng)的主流地震P波震相識別方法,將本文提出的方法與這兩種方法進(jìn)行比較。其中,STA/LTA的計算參數(shù)按照馬強(qiáng)提出的進(jìn)行設(shè)置,短時間窗口的長度為50個采樣點,長時間窗口的長度為3 000個采樣點,閾值為10[31]。AIC的時間窗范圍為STA/LTA方法識別得到的到達(dá)時刻前300個采樣點及到達(dá)時刻后30個采樣點,計算共計330個采樣點的AIC函數(shù),識別最小值作為地震波的精準(zhǔn)到達(dá)時刻,其中,k時刻AIC函數(shù)fALC(k)的計算式為

fAIC(k)=k×lg{var[x(1,k)]}+(N-k-1)×lg{var[x(k+1,N)]}

(4)

式中:N為選定窗口中所有數(shù)據(jù)的數(shù)量;var[x(1,k)]、var[x(k+1,N)]均為兩個窗口數(shù)據(jù)采樣點的方差;x(i)為i時刻的地震波信號加速度值,i=1,2,…,N。

本文共選取368組地震波進(jìn)行了魯棒性測試,并與STA/LTA及AIC方法進(jìn)行對比,其選擇標(biāo)準(zhǔn)與之前相同。同樣,人工檢測的地震波到達(dá)時刻被認(rèn)為是準(zhǔn)確的。本文選擇的對比標(biāo)準(zhǔn)是平均誤差及標(biāo)準(zhǔn)差。三種不同方法的對比結(jié)果見表1,統(tǒng)計直方圖如圖6所示。根據(jù)結(jié)果顯示,AIC算法的平均識別誤差及標(biāo)準(zhǔn)差最低達(dá)到了0.003 4 s及0.025 6 s。STA/LTA算法的平均識別誤差為0.090 3 s,標(biāo)準(zhǔn)差為0.063 0 s。而本文提出的分形斜率算法的平均識別誤差為0.006 3 s、標(biāo)準(zhǔn)差為0.043 8 s,基本與AIC方法持平,優(yōu)于STA/LTA算法。

表1 不同識別算法識別結(jié)果 s

圖6 不同識別方法的識別誤差統(tǒng)計

本文提出的地震P波識別算法,是通過識別地震波分形斜率極值時刻來判斷地震波的精準(zhǔn)到達(dá)時刻,即在地震波到達(dá)一個時間窗口后加上分形維數(shù)運算時間可以得到地震波的精準(zhǔn)到達(dá)時刻。本文選取的時間窗口長度為0.15 s,而優(yōu)化分形維數(shù)計算后的一次計算時間為0.05 s,平均誤差為0.006 3 s,即平均在地震波到達(dá)0.16 s內(nèi)可以得到地震波的精準(zhǔn)到達(dá)時間。STA/LTA方法是識別均值觸發(fā)閾值的時刻作為地震波的精準(zhǔn)到達(dá)時刻,即平均在地震波到達(dá)后0.090 3 s可以得到地震波的精準(zhǔn)到達(dá)時間。而AIC算法需要在利用STA/LTA方法識別震相時刻后,后推0.3 s后進(jìn)行AIC函數(shù)計算,即地震波到達(dá)時刻后0.39 s后獲得地震波的到達(dá)時刻,同樣滯后于本文提出的分形斜率算法。本文提出的算法,在僅比STA/LTA算法延遲0.07 s的情況下,提高了0.1 s的識別精度,在與AIC算法相比提前了0.23 s獲得了與AIC算法基本相同的精準(zhǔn)度,為后續(xù)高速鐵路地震預(yù)警系統(tǒng)的震中定位及震級估算工作的進(jìn)行提供了參考,并且在識別速度上達(dá)到《高速鐵路地震預(yù)警監(jiān)測系統(tǒng)技術(shù)條件》[30]的技術(shù)標(biāo)準(zhǔn)。

3 結(jié)論

1)提出了一種基于分形維數(shù)的小震地震P波震相識別方法。與現(xiàn)有分形理論地震P波自動震相識別算法比較,本文改進(jìn)了分形維數(shù)的計算方法,在不改變分形維數(shù)曲線的特征的情況下,提高了分形維數(shù)的計算速度。

2)在改進(jìn)的分形維數(shù)計算方法的基礎(chǔ)上,改變了以往間隔計算分形維數(shù)的算法,采用連續(xù)滑移時間窗口進(jìn)行分形維數(shù)計算,得到連續(xù)分形維數(shù)曲線充分反映地震波的細(xì)節(jié)特征。此外,分形斜率被引入,并劃分為4個階段,分析臨界點的分形斜率特征,通過引入標(biāo)準(zhǔn)差和判斷標(biāo)準(zhǔn),識別分形斜率極值點識別判斷地震P波的精準(zhǔn)到達(dá)時刻。

3)此外,研究了在震級小于5的地震事件中,分形維數(shù)時間窗口長度的選擇的最優(yōu)值,為0.15 s。本文提出的算法與STA/LTA、AIC算法相比,識別精度的平均誤差達(dá)到了0.006 3 s,標(biāo)準(zhǔn)差達(dá)到了0.043 8 s,識別誤差優(yōu)于STA/LTA算法。在識別速度上,本文提出的算法平均可以在地震P波到達(dá)0.16 s后得到精準(zhǔn)到達(dá)時刻,優(yōu)于AIC算法,且可以滿足高速鐵路地震預(yù)警系統(tǒng)時效性與準(zhǔn)確度的要求。