基于分形維數的高速鐵路地震預警系統地震P波震相識別方法

楊長衛，張凱文，吳東升，張志方，張 良，瞿立明

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

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

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

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

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

1.1 分形理論

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

D=1-F

(1)

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

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

1.2 計算原理

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

圖2 分形維數連續與間斷計算曲線

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

圖3 分形斜率識別原理

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

(2)

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

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

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

(3)

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

圖4 識別流程

1.3 影響因素分析

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

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

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

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

2 案例分析

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

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

(4)

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

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

表1 不同識別算法識別結果 s

圖6 不同識別方法的識別誤差統計

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

3 結論

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

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

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