2011年日本東北大地震特低頻地磁信號的分形標度特征研究

榮揚名，王 橋，丁 霞，黃清華

北京大學地球與空間科學學院地球物理學系，北京 100871

1 引 言

大量與地震相關的電磁觀測資料和研究成果在一定程度上證實了地震電磁現象的客觀存在性，地震電磁觀測作為一種前兆監測方法已在地震預測研究中發揮了一定的積極作用[1-12］.地震電磁現象中的特低頻（Ultra-Low-Frequency，ULF，根據國際電聯和國內通信界稱謂的特低頻為小于3Hz的頻帶范圍）電磁現象更是被認為是一種有潛力的前兆現象而受到了廣泛的關注[4-12］，也帶動了關于特低頻電磁現象物理機制的理論模型、實驗和數值模擬等研究[1-2，13-15］，只是由于實際地下結構和孕震過程的復雜性，當前對包括特低頻在內的地震電磁現象的認識還比較有限，尚有待進一步的研究.

另一方面，研究表明地震和其它一些自然災害系統呈現出典型的非線性自組織臨界狀態（selforganized criticality，SOC）[16-17］.自組織臨界性的基本圖像是，大自然在某處永久的偏離平衡，卻又被組織在一種穩定狀態中——一種臨界狀態：各種事件都能按照完全確定的統計規律——冪律發生.自組織臨界性理論是一種新型的統計理論，不依賴于系統的初始條件和任何細節部分.系統演化行為遵從冪律分布，而服從冪律分布的動力學系統表明系統內部存在自相似性.因此，自組織臨界性被猜測是相互作用的多體系統所具有的典型行為，它無論是在時間上還是空間上都具有豐富的分形結構.國內外已經在地磁和地電時間序列的分形分析方法和實際數據分析上做了大量的工作[9，18-19］，相關的數學處理手段和科學分析得到了很大程度的發展.本文采用的去傾擾動分析方法（Detrended Fluctuation Analysis，DFA）即是一種有效的數學分析方法，最先在20世紀90年代由Peng用于醫學研究，用來分析心跳間隔時間序列、DNA序列和其他生理信號記錄的分形標度特征[20-21］.近幾年，意大利的Telesca等研究了DFA方法在處理地學信號記錄方面的有效性，研究對象包括與地震和火山活動相關的地電、地磁和自電位信號記錄，得到了信號分形標度異常與相應地震火山活動比較好的關聯性[22-24］.然而，以往的研究結果大多限于單臺數據，且缺乏統計顯著性方面的檢驗，這在一定程度上制約了相關結論的可信度.

2011年3月11日在北緯38.103°，東經142.861°發生的日本東北部Mw9.0級大地震是有現代地震記錄以來最大的地震之一，并伴隨著一系列的余震活動.位于該地區附近的電磁臺站采集到了震前一段時期和余震持續期間的ULF地磁信號，這無疑為進一步探討ULF地磁信號異常與地震活動之間的關聯提供了一次良好的機會.

本文嘗試運用DFA方法處理2010年1月1日到2011年4月30日之間485天ULF地磁記錄的秒值（1Hz）和分鐘值（1／60Hz）地磁記錄，數據來源于大地震發生區域附近的 Kakioka（KAK，36.232°N，140.186°E）、Haramachi（HAR，37.615°N，140.953°E）和 Yokohama（YOK，40.993°N，141.240°E）三個地磁臺站.我們基于以上方法研究了24h單位長度的地磁時間序列分形標度特征，結合前人的研究，定量分析了非均勻標度特征的時間變化，提出了一種能反映地磁三分量非均勻標度特征同步變化的新指標，并據此探討了特低頻地磁信號分形標度特征變化與日本東北大地震之間的可能關聯性，并對得到的信號異常進行統計顯著性方面的檢驗，以提高異常的可信度.

2 觀測資料和分析方法

本文研究了來自日本東北大地震區域附近三個地磁觀測臺站的ULF地磁信號記錄，分別是作為國際地磁基準臺的KAK的1Hz數據，以及來源于另外兩個地磁連續觀測臺站HAR和YOK的分鐘值數據（1／60Hz），數據記錄時間為2010年1月1日至2011年4月30日共485天.三個臺站的分布見圖1，所選取的研究區域和時間的地震事件也在圖1標出.

我們采用的DFA方法是處理時間序列分形特征的一種統計學方法，它能去除人為造成的標度關聯性的一些信息.

第一步：對于時間序列x（i），i=1，2，…，N，求得序列的平均值，用下式求得一個新的序列：

第二步：將總長為N的積分時間序列分成等長，取每份為n的窗，n的取值范圍從4到N／4.

圖1 臺站及震中分布圖.圓圈標示的是研究區域中2010／01／01—2011／04／30期間大于4級的地震，五角星表示2011／03／11 M9.0.東北大地震的震中（38.103°N，142.861°E），三角形為三個臺站 KAK、HAR和 YOK的位置Fig.1 Distribution of the geomagnetic stations and the epicenter of earthquakes in the investigated region.The circles indicate the earthquakes with magnitude≥4.0which occurred from 2010／01／01to 2011／04／30；the triangles show three stations；the star shows the epicenter of Tohoku earthquake（38.103°N，142.861°E）

第三步：對每個長為n的數據序列進行最小二乘法線性擬合，代表此數據窗的傾向.線性擬合后得到的y軸坐標記作yn（k）.進行去傾處理，即在每個窗里用y（k）減去相應的局部傾向yn（k），然后得到用下式計算的均方根

這里N是序列的總長度.

第四步：對每個窗口標度重復第三步的運算，我們可以得到F（n）-n的一一對應關系，代表了每個窗口標度n的平均的擾動情況；

第五步：如果F（n）對n的關系符合冪律，即

將F（n）-n畫在雙對數坐標系中應呈一直線，此直線的擬合斜率便是目標時間序列對應的標度特征指數α.

對于白噪聲隨機序列，α值理論上為0.5，但如果序列過短，初始擬合斜率可能跟0.5有所偏差，當序列長度足夠大時，會無限逼近0.5；當α＞0.5時表明此序列具有強的長程相關性，即相對均值較大的值后邊出現較大的值幾率比較大，反之亦然；當α＜0.5時表明一種弱的長程相關性，即對均值較大的值后邊出現較小的值幾率比較大，反之亦然；α=1表明此序列是典型的閃爍噪聲或者說1／f信號過程；α=1.5表明是典型的布朗隨機過程[21］.作為一個例子，圖2給出了MATLAB產生的布朗隨機過程序列數值模擬的結果，序列長為217，擬合結果與理論值α=1.5吻合的很好.

但對于實際數據的時間序列，往往得不到單一

(1)確定整體教學目標：學生完成該課程的學習應達到的具體知識、能力、素質目標，并分別指出中職、高職階段應分別達到的目標，教學目標應該量化并具有可執行性。

由本文的具體數據，即每日的地磁東西分量（EW）、南北分量（NS）和垂向分量（Z）的三分量序列，可以得到三個β序列，這三個序列各包含地磁分量的分形信息，如果有大的擾動應該也包含一定的異常.為了既凸現三分量同步的異常信息，又能抑制各分量中可能出現的零星干擾，我們提出了一種新的標度不穩定性指標βXYZ，具體將每日三分量的標度不穩定性指數β取乘積，以此來表征總的標度變化效應，即βXYZ=βEW＊βNS＊βZ.

3 結 果

我們用DFA方法處理了KAK臺站485天1Hz的地磁數據，數據包括地磁場東西分量（EW）、南北分量（NS）和垂向分量（Z）的三分量時間序列.圖3給出了運用DFA方法獲得的雙對數坐標系中標度與擾動狀況冪律擬合的例子.從圖3a中可以看到對不同標度的擬合斜率比較一致，說明對所取的時間標度存在相近的擾動狀況，即分形特征比較一致，這種只用單一斜率就能描述整體擬合的情形對大多數天數都適用；而從3b和3c兩幅圖中看出，部分擬合曲線分為明顯的兩段或三段，中間存在轉折點（crossover），反應了不同標度擾動情況的非均勻性，而且這些轉折點對應的標度以及各段的擬合斜率也不盡一致；更有極少情況，如3d圖所示，雖然可以看出大致的線性擬合趨勢，但是擬合點分布過于分散，擬合相關系數相對較小.

結合以上分析，從每日的DFA冪律擬合結果并不是都能得到唯一對應的標度指數α，不同的標度范圍擬合得到的α可能不同.正因為存在這種不一致性，我們如方法中所述選擇分析擾動變化曲線的標度不穩定指數β，這樣能更有效地反映非線性系統的非均勻標度特征變化.

為了更好地分析β值的變化，可取β的滑動平均值（本文取11天滑動平均）作為參考指標.圖4（a-c）分別給出了KAK、HAR和YOK三個臺站EW地磁分量標度不穩定指數β滑動平均值的變化，為便于比較，取其均值加2倍均方根作為異常的判定閾值（圖中用虛線表示）.結果表明三個臺站的地磁EW分量對應的滑動平均β值都在2011年3月11日大規模地震事件前25～50天左右出現了異常增大現象，其中，KAK異常出現時間為2011／2／15—2011／2／18，HAR 和 YOK 為 2011／1／20—2011／1／30.其他分量所對應的β值在震前這段時期內也基本維持在相對較高的水平.

需要指出的是，上述三個臺站中HAR臺EW分量的β值在其他時間也出現了兩次相對較弱異常（圖4b），分別為2010／4／10—2010／4／11的異常1以及2010／10／12的異常2，且這兩次異常的持續時間較短（均不超過兩天），此外，同時期的全球地磁活動Ap值（圖4d）顯示2010／4／5-6兩天的Ap值分別高達55和44，異常1可能與此地磁活動有一定的關聯，因為圖4b反映的是11天滑動平均的結果，至于異常2，考慮到其只在一個臺站的一個分量圖中出現，持續時間也較短，可能是由于臺站附近局部范圍的擾動所致.

圖3 DFA方法處理KAK三分量地磁數據得到的冪律擬合示例（a）用一個標度特征指數就可以描述數據的標度特征.（b）-（c）需要兩或三個分段特征指數才能描述標度特征.（d）圖雖也能看出較好的擬合趨勢，但數據點的分布過于離散.Fig.3 Examples of power law fitting curves for the investigated data of KAK（a）Single scaling exponent；（b）Two scaling exponents；（c）Three scaling exponents；（d）Single power law fitting of scattered data.

考慮到單一地磁分量的異常可能并不足以反映非線性系統整體的異常情況，本文提出了綜合考慮三分量的標度不穩定指數β并取其乘積的一種新指標，即βXYZ=βEW＊βNS＊βZ，該指標能有效地反映地磁三分量內在的非均勻標度特征的同步變化.圖5給出了三個臺站11天滑動平均結果.從圖中可以看出，單一分量中出現的個別異常基本得到了較好的抑制，而圖4中顯示的三個臺在日本東北大地震前25～50天左右出現的異常在圖5中也有明顯體現.KAK臺站滑動平均值（圖5a）在2011／3／11—2011／3／21期間出現了非常顯著的異常變化，這可能與日本東北大地震及其余震等地震活動有關.事實上，YOK臺也出現了類似的可能與該地震余震活動相關的異常變化，只是幅度不如KAK顯著.考慮到KAK臺在主震及余震期間出現的異常過于顯著，可能會在一定程度上掩蓋震前的可能的變化，為此，我們重新分析了KAK臺截至日本東北大地震前的變化，即不分析地震及余震持續期間的βXYZ而關注震前的βXYZ并作相應的滑動平均處理（圖5d），結果表明在2011／2／13—2011／2／23出現了超出閾值的情況，即震前25天左右開始βXYZ值相對背景值明顯增大，同時期Ap值（圖4d）相對較小，表明該異常并非全球地磁擾動所致.HAR臺滑動平均值在2011／1／24—2011／1／30出現了明顯超出閾值的異常現象（圖5b），該異常大約始于東北大地震46天前.YOK 臺站滑動平均值在2011／1／20—2011／1／30期間出現了較明顯的異常（圖5c），該異常在東北大地震前50天左右出現，與HAR臺同期的異常吻合較好.總之，無論是從地磁單分量的變化特征（圖4），還是從同時考慮三分量信息的綜合指標（圖5）來看，三個臺在日本東北大地震前25～50天左右都出現了非均勻標度特征的顯著異常變化，這種基本一致的異常現象可能反映了該地震對地磁時間序列內在的非線性系統的影響.

除了在日本東北大地震前25～50天左右以及余震期間出現的異常之外，HAR臺在2010／11／11（圖5b），KAK臺在2010／3／6（圖5d）也分別出現了非均勻標度特征的異常，但這兩次異常僅在單臺出現，而且相應的地磁單分量在對應的時間段并未出現顯著異常（圖4），因此，這兩次異常的可信度并不太高.

4 討 論

研究表明，地震等構造活動可能導致ULF地磁信號的 分 形 特 征 異 常 變 化[9，19，25］.在 地 震 活 動 水平相對較低的階段，構造系統的自組織臨界性質能令ULF地磁信號的分形標度特征變化維持在比較穩定均勻的程度，而在強震等來臨前該非線性系統會逐漸失穩，進而導致ULF信號出現了較大的分形標度異常變化.運用統計分形分析方法可從實際觀測到的地磁信號中提取出這種分形異常特征，有助于加深對地震及其孕育過程的理解.

鑒于以往的研究大多限于單臺數據，且缺乏統計顯著性方面的檢驗，相關結論的可信度有待進一步的研究，本文提出了一種能反映地磁三分量非均勻標度特征同步變化的新指標，并據此探討了ULF地磁信號分形標度特征變化與日本東北大地震之間的可能關聯性.結果表明，震區附近三個臺的ULF地磁信號在日本東北大地震前25～50天左右都出現了非均勻標度特征的顯著異常變化，其中，HAR和YOK兩個臺站的異常起始時間比較一致（圖4，圖5），而KAK略晚，這可能與地震孕育過程本身的復雜性有關，也可能受到實際條件下地震電磁信號選擇性特征的影響[27-29］.

圖6 磁靜日隨機合成數據統計檢驗結果（a）合成EW分量β值11天滑動平均結果，虛線表示與實際數據得到的閾值（0.785）比較；（b）合成三分量數據的βXYZ11天滑動平均結果與實際閾值（0.315）比較.Fig.6 Stochastic test using the synthetic geomagnetic data with randomized disturbances（a）The 11-day-running-mean ofβobtained from the synthetic data of EW component；（b）The 11-day-running-mean of obtained from the synthetic geomagnetic data of three components.

本文通過地磁單分量、三分量以及多臺檢測到的所有非均勻標度特征異常的綜合分析，發現在日本東北大地震前25～50天左右以及余震期間出現的異常在多方面有基本一致的變化，相比之下，圖4和圖5中其他所有的異常由于僅在單臺或單分量序列中出現而缺乏可信度，因此，在整個研究期間內，相對可信的異常只出現在震前25～50天左右以及余震期間，該異常現象可能反映了該地震對地磁時間序列內在的非線性系統的影響.

在目前尚難以用確定的物理模型定量解釋地震前兆異常的背景下，對得到的信號異常進行統計顯著性方面的檢驗是提高異常可信度的一種有效的方法.為了檢驗本文得到的震前非均勻標度特征異常是否是一種隨機的頻發異常，本文借鑒了統計地震學中常用的隨機統計檢驗的思路[30］，來考察所得異常的顯著性.我們通過將本文方法應用于隨機合成的地磁數據的方式來進行隨機統計檢驗.具體選取某一磁靜日（例如，2011／01／09）YOK臺站的地磁記錄數據序列，每次生成相當于其2倍標準差的隨機擾動噪聲序列并與原始數據相加得到一天的隨機合成數據，重復生成485次隨機噪聲與原始信號相加即得到485天（本文研究的時間范圍）的隨機合成數據，然后運用DFA方法得到隨機合成數據的標度不穩定指數β和βXYZ滑動平均結果，并與YOK實際485天數據得到的β和βXYZ異常閾值標準（圖4c和5c）進行比較，判斷分形標度特征的異常情況，看是否會有超出閾值的情況.圖6（a、b）分別為485天隨機合成數據得到的β及βXYZ滑動平均結果.該模擬顯示隨機序列的結果并未產生如實際數據（圖4，圖5）顯示的異常變化.我們重復進行了100次類似的隨機檢驗，證實了類似圖6的結論，即隨機模擬結果均未超出實際數據的異常閾值.因此，本文得到的震前非均勻標度特征的異常變化并非由于隨機擾動所致的頻發異常，而是一種具有良好統計顯著性的可信異常，該異常可能反映了日本東北大地震對ULF地磁信號內在的非線性系統的影響.

5 結 論

本文提出了一種能反映地磁三分量非均勻標度特征同步變化的新指標，對2011年3月11日M9.0日本東北大地震震中附近三個地磁臺站16個月的ULF地磁三分量記錄進行DFA處理，得到了每日地磁時間序列的分形標度特征，并提取了各分量標度特征及三分量標度特征同步變化的異常.與地震活動的對比分析表明，三個臺站地磁記錄的標度不穩定指數滑動平均值在強震前25～50天左右都出現了比較一致的異常增大，顯示ULF地磁信號分形標度特征異常與地震活動具有較好的關聯性.進一步對所得異常進行了隨機統計檢驗，結果表明本文得到的分形標度異常并非隨機擾動所致，因而具有良好的統計顯著性.

致 謝 感謝兩位審稿人對本文提出的寶貴意見和建議.本文所使用的HAR和YOK兩個臺站的地磁觀測資料來自http：／／www.gsi.go.jp／ENGLISH／index.html KAK臺站的觀測資料來自http：／／www.kakioka-jma.go.jp.

（References）

[1]Varotsos P.The Physics of Seismic Electric Signals.Terra Pub，Tokyo，2005.

[2]Park S K，Johnston M J S， Madden T R ，et al.Electromagnetic precursors to earthquakes in the ULF band—a review of observations and mechanisms.Rev.Geophys.，1993，31：117-132.

[3]Huang Q H.Retrospective investigation of geophysical data possibly associated with the Ms8.0Wenchuan earthquake in Sichuan，China.J.Asian Earth Sci.，2011，41：421-427.

[4]Huang Q H.Rethinking earthquake-related DC-ULF electromagnetic phenomena：towards a physics-based approach.Natural Hazards and Earth System Sciences，2011，11：2941-2949.

[5]Han P，Hattori K，Huang Q，et al.Evaluation of ULF electromagnetic phenomena associated with the 2000Izu Islands earthquake swarm by wavelet transform analysis.Natural Hazards and Earth System Sciences，2011，11：965-970.

[6]Molchanov O A，Kopytenko Y A，Voronov P M，et al.Results of ULF magnetic field measurements near the epicenters of the Spitak（Ms=6.9）and Loma Prieta（Ms=7.1）earthquakes：comparative analysis.Geophys.Res.Lett.1992，19：1495-1498.

[7]Fraser-Smith A C，Bernardi A，McGill P R，et al.Lowfrequency magnetic field measurements near the epicenter of the Ms7.1Loma Prieta earthquake.Geophys.Res.Lett.，1990，17：1465-1468.

[8]Hayakawa M，Kawate R，Molchanov O A，et al.Results of ultra-low-frequency magnetic field measurements during the Guam earthquake of 8August 1993.Geophys.Res.Lett.，1996，23：241-244.

[9]Hayakawa M，Itoh T，Smirnova N.Fractal analysis of geomagnetic ULF data associated with the Guam earthquake on August 8，1993.Geophys.Res.Lett.1999，26：2797-2800.

[10]Kawate R， Molchanov O A， Hayakawa M.Ultra-low-frequency magnetic fields during the Guam earthquake of 8 August 1993and their interpretation.Phys.Earth Planet.Inter.1998，105：229-238.

[11]Du A，Huang Q，Yang S.Epicenter location by abnormal ULF electromagnetic emissions.Geophys.Res.Lett.，2002，29（10）：1455，641-643.

[12]杜愛民，周志堅，徐文耀等.新疆和田ML7.1地震前ULF電磁輻射的激發機理.地球物理學報，2004，47（5）：832-837.Du A M，Zhou Z J，Xu W Y，et al.Generation mechanisms of ULF electromagnetic emissions before the ML=7.1 earthquake at Hotan of Xinjiang.Chinese J.Geophys.（in Chinese），2004，47（5）：832-837.

[13]Zhao G Z，Zhan Y，Wang L F，et al.Electromagnetic anomaly before earthquakes measured by EM experiment.Earthquake Science，2009，22：395-402.

[14]Huang Q.One possible generation mechanism of co-seismic electric signals.Proceedings of the Japan Academy，2002，78：173-178.

[15]Ren H X，Chen X F，Huang Q H.Numerical simulation of co-seismic electromagnetic fields associated with seismic waves due to finite faulting in porous media.Geophysical Journal International.2012，188：925-944.

[16]Bak P，Tang C，Wiesenfeld K.Self-organized criticality：An explanation of 1／fnoise.Phys.Rev.Lett.，1987，59：381-384.

[17]Bak P.How Nature Works（The Science of Self-organized Critiality），Oxford University Press，1997：210.

[18]Higuchi T.Approach to an irregular time series on the basis of the fractal theory.Physica D，1988，31：277-283.

[19]Gotoh K，Hayakawa M，Smirnova N A，et al.Fractal Analysis of Seismogenic ULF Emissions.Physics and Chemistry of the Earth，2004，29：419-424.

[20]Peng C-K，Havlin S，Stanley H E，et al.Quantification of Scaling exponents and crossover phenomena in nonstationary heartbeat time series.Chaos，1995，5（1）：82-87.

[21]Peng C K，Buldyrev S V， Havlin S，et al.Mosaic organization of DNA nucleotides.Phys.Rev.，E，1994，49：1685-1689.

[22]Telesca L，Hattori K.Non-uniform scaling behavior in ultralow-frequency（ULF）earthquake-related geomagnetic signals.Physica A，2007，384：522-528.

[23]Telesca L，Balasco M，Lapenna V.Investigating the timecorrelation properties in self-potential signals recorded in a seismic area of Irpinia，southern Italy.Chaos，Solitions and Fractals，2007，32：199-211.

[24]Balasco M，Lapenna V，Romano G，et al.Extracting quantitative dynamics in Earth′s apparent resistivity time series by using the detrended fluctuation analysis.Physica A，2007，374：380-388.

[25]Telesca L，Lapenna V，Macchiato M，et al.Investigating non-uniform behavior in Ultra Low Frequency （ULF）earthquake-related geomagnetic signals.Earth and Planetary Sciences Letters，2008，268：219-224.

[26]Viswanathan G M，Peng C-K，Stanley H E，et al.Deviations from uniform power law scaling in nonstationary time scales.Phys.Rev.，E.，1997，55：845-849.

[27]黃清華，劉濤.新島臺地電場的潮汐響應與地震.地球物理學報，2006，49：1745-1754.Huang Q H，Liu T.Earthquakes and tide response of geoelectric potential field at the Niijima station.Chinese J.Geophys.，2006，49：1745-1754.

[28]黃清華，林玉峰.地震電信號選擇性數值模擬及可能影響因素.地球物理學報，2010，53：535-543.Huang Q H，Lin Y F.Numerical simulation of selectivity of seismic electric signal and its possible influences.Chinese J.Geophys.（in Chinese），2010，53：535-543.

[29]Huang Q H，Lin Y F.Selectivity of seismic electric signal（SES）of the 2000Izu earthquake swarm：a 3DFEM numerical simulation model.Proceedings of the Japan Academy，2010，86：257-264.

[30]Huang Q H.Search for reliable precursors：A case study of the seismic quiescence of the 2000western Tottori prefecture earthquake.J.Geophys.Res.，2006，111：B04301.