張雁濱 蔣 駿 郭 斌 顧升宇 徐文杰 王 朝
1 華中科技大學(xué)物理學(xué)院,武漢市珞喻路1037號,430074
寬頻帶地震觀測技術(shù)提供了前所未有的更豐富的地球物理信號,例如在俯沖帶和板塊交界處觀測到的非火山震顫(nonvolcanic tremor)、間歇式震顫滑動(ETS)、低頻地震滑動(LFEs)、甚低頻地震(VLFE)、慢滑動事件(SSE)等[1-6]。在中國大陸地區(qū),寬頻帶地震計(jì)以及重力和傾斜儀普遍都能觀測到類似的信號,這里統(tǒng)稱為震顫波信號[7-10],與國際上在特殊地區(qū)觀測到的Tremor、ETS、LFEs、VLFE、SSE 等信號相比,既有類似之處,又有所不同。顯然,在分析研究大陸地震觀測中的震顫波與大規(guī)模天氣系統(tǒng)活動的關(guān)系、震顫波與構(gòu)造運(yùn)動背景及一些強(qiáng)震活動之間的關(guān)系、產(chǎn)生機(jī)理等問題時,對震顫波信號的定位是關(guān)鍵。
對能夠拾取初動的信號,只要滿足觀測條件,便可通過提取波初動的方法進(jìn)行較精確定位。而對小振幅、低信躁比、無明顯初動到時的弱信號,例如非火山間歇震顫信號,目前國際上采用的方法有顫動信號提取、信號振幅平方、采用互相關(guān)方法確定震顫信號包絡(luò)線的相對到時、以網(wǎng)格搜索法來進(jìn)行信號源的定位[1,11-12]。此外,Kao等[13]提出震源掃描算法(the source-scanning algorithm,SSA)進(jìn)行震顫源的定位。這些方法都主要針對間歇式的震顫信號(每個震顫的持續(xù)時間為s~min),研究的空間尺度范圍一般都在幾十km 或上百km,分布有非常密集的觀測臺網(wǎng)。
連續(xù)震顫波信號持續(xù)時間在幾h以上,觀測范圍可達(dá)幾千km,信噪比低、傳播路程長、不同臺站的記錄變異大,因此信號源的定位難度非常大。本文參考國際上對間歇式震顫的定位方法,在對連續(xù)震顫波信號的特征進(jìn)行分析的基礎(chǔ)上,根據(jù)實(shí)際觀測特點(diǎn),在信號的提取和識別及處理方面進(jìn)行整合和改進(jìn),得到可行的定位技術(shù)和方法,并通過對已知臺風(fēng)引發(fā)的連續(xù)震顫波的定位進(jìn)行方法檢驗(yàn)。
早先,我們參考文獻(xiàn)[11-12],對地震觀測數(shù)據(jù)中幾個已知的來自西太平洋臺風(fēng)產(chǎn)生的震顫波信號進(jìn)行定位方法的實(shí)驗(yàn)。在實(shí)驗(yàn)中發(fā)現(xiàn),由于連續(xù)震顫波持續(xù)震動幾十h不斷,觀測臺站的分布范圍大,信號傳遞時間長,分析所涉及信號的時間尺度范圍大,觀測中不同信號的互相關(guān)相差不明顯,給“同一信號包絡(luò)線”的相對到時判斷工作帶來困難。對此,采取了一系列處理手段:從連續(xù)的長時段中選擇有相似突變特征信號的時段討論分析①;在信號提取過程中,除在對信號特征分析的基礎(chǔ)上進(jìn)行濾波外,考慮到震顫波中含有多種復(fù)雜的信號成分,以及定位對信號時間變化的精準(zhǔn)要求,對信號識別提取以及分析方法進(jìn)行了改進(jìn)。利用小波分析在信號時頻分析上的優(yōu)勢,進(jìn)行震顫波中各信號的識別和分析②,比較不同臺站、不同支波的信號,進(jìn)行多個信號時間段的大致確定(進(jìn)一步縮小分析的時間尺度);在對包絡(luò)信號的互相關(guān)分析上,不僅對單個包絡(luò)信號,而且將一個確定時間窗中連續(xù)的幾個包絡(luò)信號作為一體來進(jìn)行相關(guān)分析等??傊?,在“更準(zhǔn)確地提取、識別相同的震顫波信號,提高信噪比”上做工作,最后以所確定時段內(nèi)的同一組信號包絡(luò)線的相對到時,利用網(wǎng)格搜索法進(jìn)行信號源的定位。以1 Hz采樣數(shù)據(jù),對幾個已知臺風(fēng)信號進(jìn)行定位方法實(shí)驗(yàn),定位結(jié)果與同時段臺風(fēng)中心位置的誤差范圍約40~100km①易鵬程.震顫信號的起始時間判斷與定位.華中科技大學(xué)物理學(xué)院畢業(yè)論文,2009年6月。指導(dǎo)教師:張雁濱,蔣駿。②劉偉.基于小波分析方法的震顫信號起始時間判斷與定位.華中科技大學(xué)物理學(xué)院畢業(yè)論文,2010年6月。指導(dǎo)教師:張雁濱,蔣駿。。根據(jù)臺風(fēng)的結(jié)構(gòu)和運(yùn)動狀況,對其產(chǎn)生的信號源的定位來說,這樣的結(jié)果是合理的。以上采用的方法最終需要獲得信號的相對到時。雖然實(shí)驗(yàn)結(jié)果都較理想,但在多臺站觀測信號的互相關(guān)分析以及“同一組包絡(luò)信號線的相對到時”的確定環(huán)節(jié),困難較大。例如實(shí)際工作中所遇到的最大問題是不相關(guān)的波形間相關(guān)系數(shù)也很高,從而需加入人為判斷,不利于方法的程式化。因此,采用震源掃描算法(SSA 法)來取代前面的到時判斷以及網(wǎng)格搜索定位環(huán)節(jié)。
SSA 法[13]的核心是通過觀測波形信號的相對振幅和到時來判定震源是否位于特定的時間和地點(diǎn),通過系統(tǒng)地掃描一個試探性的位置和發(fā)震時刻,得到整個震源序列分布,而不需要拾取到時,也不需要計(jì)算理論地震圖(圖1)。
圖1 震源掃描算法原理示意圖Fig.1 Schematic diagram to illustrate the concept of SSA
假設(shè)一個地震事件被N個臺站的地震儀記錄到。首先對每個數(shù)字地震信號的振幅歸一化,然后計(jì)算點(diǎn)η在時刻τ的亮度函數(shù)br(η,τ),定義如下[13]:
式中,un是在臺站n記錄的歸一化地震記錄,tηn為從η點(diǎn)到臺站n計(jì)算的某個最大震相的走時。如果所有的最大振幅都是由點(diǎn)η在時間τ產(chǎn)生,那么br(η,τ)=1(圖1)。同樣,br(η,τ)=0.1意味著位于點(diǎn)η、時間τ的震源產(chǎn)生的最大振幅在每個臺站只有平均10%被觀測到。通過系統(tǒng)地搜索η點(diǎn)和時間τ來尋找亮度函數(shù)的最大值,便可有效地重構(gòu)震源的時空分布。在通常情況下,br(η,τ)=1極少出現(xiàn),可以認(rèn)為,整個時空中“亮度”函數(shù)取最大時是震源發(fā)生的位置和事件發(fā)生的時刻。有關(guān)該方法的更多信息及應(yīng)用細(xì)節(jié)可參考文獻(xiàn)[3,13-14]。
在觀測范圍較大、臺網(wǎng)較稀疏的條件下,對大陸觀測中那些持續(xù)時間長、到時難以確定、不同臺站觀測信號的相關(guān)不明顯的連續(xù)震顫信號,在對信號進(jìn)行識別分析以及一系列處理的基礎(chǔ)上,采用SSA 法來取代到時確定、搜索定位環(huán)節(jié),以震源掃描的“最大亮度”值進(jìn)行連續(xù)震顫波信號的定位,可解決前述實(shí)際困難。
經(jīng)過反復(fù)應(yīng)用和檢驗(yàn),得到適合于連續(xù)震顫波信號的定位技術(shù)方案:
1)應(yīng)用頻譜特征、時頻分析、小波變換等多種方法,對觀測信號中的連續(xù)震顫信號進(jìn)行識別和提??;
2)采用“振幅平方法”、“差分法”等對信號進(jìn)行處理,提高其信噪比;
3)選擇包含有較強(qiáng)特征信號的時段(可以是多個)作為分析樣本信號,盡可能縮小連續(xù)震顫波分析討論的時間尺度;
4)對樣本信號取包絡(luò)線(也可采用對數(shù)據(jù)平滑的方法),即對波群信號進(jìn)行有關(guān)時間的分析,并根據(jù)具體觀測情況(臺站分布、掃描空間范圍),來選擇掃描所用信號的時間段(時間窗長);
5)分別對所選擇時段內(nèi)的信號進(jìn)行歸一化后,用SSA 法進(jìn)行空間掃描定位,給出“亮度”較大的結(jié)果。
基于上述基本原則和技術(shù)方法,對實(shí)際觀測中不同的信號進(jìn)行測試和檢驗(yàn)。這里給出其中幾個實(shí)例。
SSA已被引用于微震的定位[14-15],這里把上述技術(shù)方法用于稀疏臺站觀測、低采樣信號,對已知小地震進(jìn)行定位檢驗(yàn)。采用中國大陸CD2、KMI、GOM、GYA、LAS、LZH 和XAN 等7個臺站觀測的2008-05-27四川青川4.4級地震連續(xù)波形數(shù)據(jù)(圖2(a)),采樣率1Hz(對原始地震信號作重采樣處理)。需要說明的是,與連續(xù)震顫波信號不同,地震信號的信噪比高、持續(xù)時間短、信號清晰明確,但因觀測點(diǎn)不同,各臺觀測的地震波形有差別,所以這里只是對這個地震信號進(jìn)行了平方、取包絡(luò)及歸一化等處理(圖2(b)),然后用SSA 法,在95°~115°E、20°~40°N 空間范圍掃描定位。在一次掃描的基礎(chǔ)上縮小空間進(jìn)行二次掃描計(jì)算,當(dāng)取3.1km/s波速模型時(圖3),“亮度”最大值位于104.92°E、32.24°N(中國地震臺網(wǎng)中心公布的地點(diǎn)為104.74°E、32.24°N)。由于這里所用信號的采樣率相對較低,又對數(shù)據(jù)進(jìn)行了一系列處理,所以這個偏差是合理的。關(guān)于波速模型,考慮到掃描空間較大,而且對信號作了取包絡(luò)處理,波速模型的準(zhǔn)確性可能帶來偏差,所以在實(shí)驗(yàn)中,還采用了不同的面波波速模型(2.8~3.5km/s)進(jìn)行計(jì)算比較。結(jié)果顯示,采用不同波速會帶來幾十km的偏差,波速在3.0km/s左右時結(jié)果偏差小。
圖2 2008-05-27各臺站觀測的青川4.4地震信號(1Hz采樣)Fig.2 The observation of Qingchuan 4.4earthquake on May 27th,2008(1Hz sampling)
圖3 2008-05-27青川4.4地震SSA 掃描圖Fig.3 Application of SSA to locate Qingchuan 4.4 earthquake on May 27th,2008
將上述技術(shù)用于已知臺風(fēng)引發(fā)的連續(xù)震顫波信號來進(jìn)行定位實(shí)驗(yàn)。選擇兩個發(fā)生在西太平洋、運(yùn)動路徑比較有代表性的臺風(fēng)——2006-05強(qiáng)臺風(fēng)珍珠(Chanchu)和2007-07碧利斯(Bilis)為研究對象,采用中國大陸GZH、CD2、WHN、GYA、GUL、CNS、HEF、QZH、WZH 和SSE 等臺站的地震觀測記錄。根據(jù)對各個臺站觀測信號的小波分析結(jié)果發(fā)現(xiàn),臺風(fēng)引發(fā)的連續(xù)震顫波信號主要分布在d1、d2和d3層,而且與d1和d2層信號比,d3層信號(1Hz采樣數(shù)據(jù),8~16s周期段)雖相對弱,但在較大的觀測范圍內(nèi),不同臺站的d3層信號特征表現(xiàn)一致性強(qiáng),信號在傳播空間上衰減相對緩慢,所以采用d3層的信號進(jìn)行震顫波源的定位。
圖4 臺站分布示意圖Fig.4 Schematic diagram of station distribution
對2006-05-14~19“珍珠”引發(fā)的連續(xù)震顫波的d3信號,選其中信號較強(qiáng)時段的數(shù)據(jù)取振幅平方(以05-16 20:00、05-07 00:00各1h的信號為例,見圖5),根據(jù)需要確定掃描時間窗長(這里選1 000s),再對所選時間窗內(nèi)的信號(例如圖5(a)中陰影部分,可任意前后挪動這個窗來選不同時段的信號)經(jīng)取包絡(luò)線、歸一化處理后以SSA方法進(jìn)行震源掃描。在一次掃描的基礎(chǔ)上,縮小掃描的空間為110°~120°E、15°~25°N,進(jìn)行二次定位。用圖5(a)、5(b)所示數(shù)據(jù)中多個不同時段的信號,以不同的波速(采用面波波速模型2.8~3.5km/s)進(jìn)行掃描計(jì)算實(shí)驗(yàn),所得結(jié)果穩(wěn)定、地點(diǎn)集中,與臺風(fēng)路徑相比較,所得地點(diǎn)均分布在距臺風(fēng)中心幾十到100余km 的范圍內(nèi),并能正確反映臺風(fēng)中心隨時間的運(yùn)動變化趨勢,結(jié)果可信。因篇幅有限,這里僅給出對應(yīng)圖5(a)、5(b)所示數(shù)據(jù)中多個時段掃描結(jié)果中的一個圖(圖6(a)、6(b),波速采用3.0km/s)。
圖5 各臺站觀測的震顫波取平方后的信號Fig.5 The signals after square processing
圖6 臺風(fēng)“珍珠”的連續(xù)震顫波信號震源掃描結(jié)果圖例Fig.6 The map by SSA of Chanchu’s tremor wave
將同樣的方法用于2006-07臺風(fēng)“碧利斯”引發(fā)的連續(xù)震顫波信號。同樣采用db小波分層進(jìn)行信號的識別,采用d3層的信號,選取其中07-13 12:00、07-14 04:00的數(shù)據(jù)為例(圖7(a)、7(b)),取掃描時間窗長1 000s,采用與前面相同的方法,對這2h的數(shù)據(jù)分時段進(jìn)行信號處理后,用不同的波速(采用面波波速模型2.8~3.5km/s)進(jìn)行震源掃描實(shí)驗(yàn)計(jì)算,掃描后的定位結(jié)果也同樣能反映臺風(fēng)的運(yùn)動趨勢,結(jié)果也穩(wěn)定分布在距臺風(fēng)中心20~100 多km 的范圍內(nèi)(圖8(a)、8(b))。并且,對“碧利斯”的觀測而言,其臺站的空間分布相對要比“珍珠”的合理,掃描的結(jié)果也優(yōu)于“珍珠”。
以上實(shí)驗(yàn)中,因?yàn)閷?shù)據(jù)進(jìn)行了小波變換、平方、取包絡(luò)或平滑等一系列處理,可考慮這些信號是以波群的“視速度”傳播。若視其在地表傳播,以面波波速2.5~3.5km/s進(jìn)行測試,得到的結(jié)果基本分布在距該時段臺風(fēng)運(yùn)動路徑中心的20~180km 范圍內(nèi),說明這個波速模型是可行的。臺風(fēng)的本身結(jié)構(gòu)半徑在500~1 000km,即使是臺風(fēng)的環(huán)狀區(qū)域(最大風(fēng)速和最強(qiáng)暴雨集中出現(xiàn)的區(qū)域),也是位于距臺風(fēng)中心約幾十到100 多km 半徑范圍內(nèi),加上大陸的觀測臺站均分布在臺風(fēng)的西北部,位置并不十分理想,綜上因素,定位的結(jié)果是合理的。此外,對連續(xù)震顫波中不同時段的信號以相同方法進(jìn)行定位,最后所得結(jié)果也都是穩(wěn)定的,地點(diǎn)相對集中,定位結(jié)果也能反映臺風(fēng)的運(yùn)動趨勢,達(dá)到對此類大規(guī)模天氣運(yùn)動與大陸地震觀測中出現(xiàn)的震顫波相關(guān)分析中對信號源定位的要求。這也表明,對于連續(xù)的震顫波信號用該方法得到的定位結(jié)果是可信的。
圖7 各臺站觀測的震顫波取平方后的信號Fig.7 The signals after square processing
本文參考國際上對間歇式震顫的定位方法,根據(jù)連續(xù)震顫波的特點(diǎn),在準(zhǔn)確識別提取信號、提高信噪比、縮小分析討論的時間尺度等原則基礎(chǔ)上,對信號的識別及處理技術(shù)進(jìn)行整合和改進(jìn),得到對這類信號所采用的定位途徑和方法流程:識別分析-信號處理-掃描定位。對已知臺風(fēng)引發(fā)的連續(xù)震顫波定位檢驗(yàn)結(jié)果表明,這些技術(shù)方法是可行的,結(jié)果合理可信。
對連續(xù)震顫波信號的識別提取是其中關(guān)鍵的一環(huán),信號若識別準(zhǔn)確,再經(jīng)過有效的處理,即使對信號以SSA 方法進(jìn)行無特定信號的“盲目掃描”,也能獲得可信的穩(wěn)定結(jié)果。該定位方法和技術(shù)已實(shí)現(xiàn)計(jì)算機(jī)程式化。
以上對連續(xù)震顫波定位所采用的技術(shù)方法還有可提升的空間,例如在信號的識別環(huán)節(jié),可以采用時頻分析技術(shù)來提高對“相同信號”的確定;數(shù)據(jù)處理中可用取包絡(luò)線的方法,也可嘗試采用數(shù)據(jù)平滑等多種方式;根據(jù)實(shí)際觀測的信號情況來選擇合適的掃描時間窗的長度等。由于連續(xù)震顫波的持續(xù)時間長,所用數(shù)據(jù)的采樣率一般為1 Hz,嘗試用采樣率為20Hz的連續(xù)信號來進(jìn)行定位時,準(zhǔn)確度雖有提高,但在“信號識別”及“相關(guān)分析”環(huán)節(jié)會帶來更多的困難。
[1]Obara K.Nonvolcanic Deep Tremor Associated with Subduction in Southwest Japan[J].Science,2002,296(5 565):1 679-1 681
[2]Rogers G,Dragert H.Episodic Tremor and Slip on the Cascadia Subduction Zone:The Chatter of Silent Slip[J].Science,2003,300(5 627):1 942-1 943
[3]Kao H,Shah S J,Herb D,et a1.A Wide Depth Distribution of Seismic Tremors along the Northern Cascadia Margin[J].Nature,2005,436(7 052):841-844
[4]Shelly D R,Beroza G C,Ide S,et al.Low-Frequency Earthquakes in Shikoku,Japan,and Their Relationship to Episodic Tremor and Slip[J].Nature,2006,442:188-191
[5]Ito Y,Obara K,Shiomi K,et al.Slow Earthquakes Coincident with Episodic Tremors and Slow Slip Events[J].Science,2007,315(5 811):503-506
[6]Shelly D R,Beroza G C,Ide S.Non-Volcanic Tremor and Low-Frequency Earthquake Swarms[J].Nature,2007,446(7 133):305-307
[7]張雁濱,蔣駿,廖盈春,等.寬頻地震計(jì)及傾斜、重力儀對長周期波動信號的綜合觀測[J].地震學(xué)報(bào),2008,30(6):626-633(Zhang Yanbin,Jiang Jun,Liao Yingchun,et al.Joint Observation of Long Period Tremor Signals with Broadband Seismometer,Tiltmeter and Gravimeter[J].Acta Seismologica Sinica,2008,30(6):629-635)
[8]張雁濱,蔣駿,李勝樂,等.熱帶氣旋引起的震顫波[J].地球物理學(xué)報(bào),2010,53(2):335-341(Zhang Yanbin,Jiang Jun,Li Shengle,et al.The Tropic Cyclone Aroused the Tremor[J].Chinese J Geophys,2010,53(2):335-341)
[9]蔣駿,張雁濱,林鋼,等,固體潮觀測中的震顫異常波[J].地球物理學(xué) 報(bào),2012,55(2):462-471(Jiang Jun,Zhang Yanbin,Lin Gang,et al.The Tidal Instruments Recorded Abnormal Tremor Wave[J].Chinese J Geophys,2012,55(2):462-471)
[10]張雁濱,蔣駿,李才媛,等,昆侖山強(qiáng)震前的震顫波并非源自慢地震[J].地球物理學(xué)報(bào),2013,56(3):869-877(Zhang Yanbin,Jiang Jun,Li Caiyuan,et al.The Tremor Wave before the Kunlun Strong Earthquake is not the Slow Earthquake Event[J].Chinese J Geophys,2013,56(3):869-877)
[11]Rubinstein J L,Rocca M L,Vidale J E,et al.Tidal Modulation of Nonvolcanic Tremor[J].Science,2008,319(5 860):186-189
[12]Gomberg J,Rubinstein J L,Peng Z G,et al.Widespread Triggering of Nonvolcanic Tremor in California[J].Science,2008,319(5 860):173
[13]Kao H,Shan S J,The Source-Scanning Algorithm:Mapping the Distribution of Seismic Sources in Time and Space[J].Geophys J Int,2004,157(2):589-594
[14]李文軍,陳棋福.用震源掃描算法(SSA)進(jìn)行微震的定位[J].地震,2006,26(3):107-115(Li Wenjun,Chen Qifu.Micro-Earthquake Location by Means of Source-Scanning Algorithm[J].Earthquake,2006,26(3):107-115)
[15]李文軍,李麗,陳棋福.用震源掃描算法(SSA)研究列車源的運(yùn)動[J].地球物理學(xué)報(bào),2008,51(4):1 146-1 151(Li Wenjun,Li Li,Chen Qifu.Research on the Movement of Vibration Source of Train by Means of SSA[J].Chinese J Geophys,2008,51(4):1 146-1 151)