解滔 盧軍 杜學(xué)彬

1)中國地震臺網(wǎng)中心，北京 100045 2)甘肅省地震局，蘭州 730000

0 引言

直流視電阻率(以下簡稱視電阻率)是實(shí)施地震中短期預(yù)測的有效觀測手段，在50多年的觀測實(shí)踐中，記錄到了數(shù)十次在觀測站網(wǎng)內(nèi)及附近發(fā)生5～8級地震前突出的中期及短臨異常(錢家棟等，1985； 錢復(fù)業(yè)等，1998； 汪志亮等，2002； 杜學(xué)彬，2010)。與地震有關(guān)的異常變化形態(tài)表現(xiàn)為偏離之前的多年背景值變化范圍、持續(xù)時間為數(shù)月至兩年左右的下降或上升變化以及地震發(fā)生后的恢復(fù)過程，部分地震前1個月內(nèi)出現(xiàn)過加速下降變化，部分臺站觀測到“準(zhǔn)同震”階躍變化(錢復(fù)業(yè)等，1982； 錢家棟等，1985； 趙玉林等，2001；Lu et al，2016)。地震前的視電阻率異常以下降變化為主，對于7級以上地震，異常主要集中在距離震中約200km的范圍內(nèi)(錢家棟等，1995； 汪志亮等，2002； 杜學(xué)彬，2010)，許多大地震前出現(xiàn)與主壓應(yīng)力方位有關(guān)的各向異性變化(錢復(fù)業(yè)等，1996； 杜學(xué)彬等，2007； 解滔等，2018)，因而視電阻率異常主要反映近源區(qū)域介質(zhì)受孕震應(yīng)力的影響(趙玉林等，1996； 解滔等，2020b)。

從原始觀測數(shù)據(jù)中識別異常變化，是實(shí)施地震預(yù)測的重要環(huán)節(jié)。視電阻率觀測具有大尺度體積平均效應(yīng)，在測區(qū)環(huán)境未受顯著干擾的情況下，觀測值背景變化較為平穩(wěn)。因此，原始曲線分析是最為直接和有效的分析方法，可以直接識別出趨勢上升/下降、趨勢轉(zhuǎn)折、年變化畸變等較為顯著的變化。在分析變化幅度時，通常采用傅立葉滑動、距平/動態(tài)距平等方法進(jìn)行去年變化處理(杜學(xué)彬等，2017)。但是，在實(shí)際的預(yù)測實(shí)踐過程中，異常的分析和確認(rèn)十分困難，一方面是由于觀測數(shù)據(jù)受到不同程度的干擾，但更為核心的是異常形態(tài)、幅度和空間分布非常復(fù)雜。諸如1976年唐山MS7.8 地震和2008年汶川MS8.0 地震前幅度大且形態(tài)清晰的異常極為少見，更多地震前的異常十分微弱。從去年變后曲線上識別異常起始時間也存在一定的不確定性，不利于提取弱異常。此外，不同臺站地下地層結(jié)構(gòu)的放大系數(shù)存在較大差異，對相同裂隙體積變化的響應(yīng)本身也存在差異(趙玉林等，1983； 解滔等，2020a；Xie et al，2020)，采用統(tǒng)一的變化幅度作為異常閾值容易遺漏弱異常。杜學(xué)彬等(2001)提出歸一化變化速率方法，基于分析時段內(nèi)變化速率的均方差設(shè)定異常閾值，克服了不同臺站異常幅度之間存在差異的問題。但是，歸一化變化速率方法采用固定長度的窗口進(jìn)行滑動，變化速率也無法反映異常時段內(nèi)的變化幅度。因此，在開展異常分析時，變化速率不如變化幅度直觀。

本文提出一種自適應(yīng)計算相對變化幅度的方法，用于提取視電阻率中短期異常。該方法有效緩解了異常起始時間不確定、固定窗長無法獲取異常時段內(nèi)累計變化幅度的問題； 采用基于分析時段內(nèi)常態(tài)變幅均方差的異常閾值，也克服了不同臺站響應(yīng)能力存在差異的問題。

1 算法流程

自適應(yīng)變化幅度方法的計算過程包括年變化消除、趨勢提取、變化幅度計算、閾值線計算4個環(huán)節(jié)。為同步分析年變化形態(tài)類異常，算法中增加了年變化幅度計算。采用該方法對視電阻率數(shù)據(jù)進(jìn)行分析的流程如圖1所示。

圖1 自適應(yīng)變化幅度方法分析流程

1.1 年變化消除

采用傅立葉滑動方法提取觀測數(shù)據(jù)中的年變化信息。設(shè)X={x1，x2，…，xN}為視電阻率觀測值序列，數(shù)據(jù)長度為N，年變化周期為T(日均值：T=365天； 月均值：T=12個月)，記Y={y1，y2，…，yN}為觀測數(shù)據(jù)中的年變化成分。

對于年變化成分中的第k個值yk(k=T，T+1，…，N)，采用觀測數(shù)據(jù)X中的第k個值和之前T-1個值{xk-T+1，xk-T+2，…，xk}進(jìn)行擬合(趙躍辰等，1984)

(1)

其中，j=k-T+1，k-T+2，…，k。

取式(1)中零相位基波作為年變化成分，即j=k，此時有yk=ak，即

(2)

對于年變化成分中的前T-1個值，通常采用觀測值{x1，x2，…，x2T-2}進(jìn)行逆向擬合獲得。從觀測數(shù)據(jù)X中減去年變化Y，可得去年變化數(shù)據(jù)G{gi=xi-yi}。

這里以定西臺NS測道2016—2020年的月均值數(shù)據(jù)為例(圖2(a))，圖2(b)為年變化成分，圖2(c)為去年變后數(shù)據(jù)。

圖2 自適應(yīng)變化幅度方法提取視電阻率中短期異常分析過程示例(a)定西臺NS測道月均值；(b)年變化成分；(c)去年變數(shù)據(jù)；(d)去年變數(shù)據(jù)和趨勢變化；(e)相對變化幅度；(f)相對變化幅度(黑色實(shí)線)、2.5倍"平均均方差"閾值線(紅色虛線)和年變化幅度曲線(藍(lán)色實(shí)線)

1.2 趨勢提取

采用離散小波變換提取去年變數(shù)據(jù)G中的趨勢變化成分Q={q1，q2，…，qN}，選擇10階Daubichies小波基Db10對數(shù)據(jù)G進(jìn)行分解。在離散小波分析中，為了減少數(shù)據(jù)兩端的邊界效應(yīng)，分解時在數(shù)據(jù)兩端進(jìn)行鏡像延拓，DbN小波基的延拓長度為2N-1。對于存在趨勢上升或趨勢下降的月均值數(shù)據(jù)，采用Db10進(jìn)行分解時，單側(cè)延拓長度達(dá)19個月，會引起數(shù)據(jù)起始段和末尾段趨勢線擬合的嚴(yán)重失真。因此，在分析月均值數(shù)據(jù)時，先將月均值數(shù)據(jù)線性插值為日均值，待趨勢線擬合之后按對應(yīng)時間點(diǎn)抽取月均值的趨勢變化。對于日均值數(shù)據(jù)(含月均值插值為日均值的情況)，小波分解階數(shù)設(shè)定為10，取10階尺度部分作為趨勢變化成分。圖2(d)中藍(lán)色實(shí)線為定西臺NS測道的趨勢變化成分。

1.3 相對變化幅度

以去年變數(shù)據(jù)G與趨勢變化Q的每一次相交點(diǎn)作為計算變化幅度的起點(diǎn)，計算其之后數(shù)據(jù)相對于該相交點(diǎn)的變化幅度，記相對變化幅度數(shù)據(jù)為F={f1，f2，…，fN}。如圖2(d)所示，數(shù)據(jù)G與趨勢變化Q在gm位置存在相交點(diǎn)，下一個相鄰的相交點(diǎn)位于gn處。則gm至gn-1區(qū)間的相對變化幅度為

(3)

對于數(shù)據(jù)起始段，G與Q通常不存在相交點(diǎn)，記第一個相交點(diǎn)位于gk處(圖(2)d)，則g1至gk-1區(qū)間的相對變化幅度為

(4)

1.4 異常閾值線

在震情跟蹤工作中進(jìn)行異常分析時，通常關(guān)注分析數(shù)據(jù)末尾段是否存在異常。因此，若將基于變化幅度的均方差作為異常閾值，則分析位置的異常閾值應(yīng)只與之前的數(shù)據(jù)有關(guān)。此外，視電阻率中期異常持續(xù)時間通常在數(shù)月至兩年左右，且參與分析的數(shù)據(jù)長度通常為震前數(shù)年。隨著異常測值點(diǎn)的增加，相對變化幅度數(shù)據(jù)集F的均方差也會快速增加，造成異常閾值線的不穩(wěn)定。我們采用計算前向累計均方差和后向累計均方差的方式，可使得異常閾值線基本保持穩(wěn)定。

前向累計均方差按如下方式計算

(5)

后向累計均方差按如下方式計算

(6)

1.5 年變化幅度

年變化形態(tài)類異常也是視電阻率異常中非常重要的一種類型，通常表現(xiàn)為年變化幅度減小、增大或形態(tài)消失。為便于異常分析，在相對變化幅度分析的基礎(chǔ)上增加年變幅度分析。

對年變化序列Y={y1，y2，…，yN}進(jìn)行Hilbert變換，構(gòu)造年變化的解析信號序列S={s1，s2，…，sN}。年變化時間序列Y的Hilbert變換為

(7)

式(7)為時間域的褶積運(yùn)算，可先經(jīng)傅里葉變換至頻率域進(jìn)行乘積運(yùn)算后，再進(jìn)行傅里葉逆變換至?xí)r間域。Hilbert變換函數(shù)h(t)=1/(πt)的頻率響應(yīng)為

(8)

(9)

2 異常提取示例

分別以通渭臺在2013年甘肅岷縣-漳縣MS6.6 地震、柯坪臺在2020年伽師MS6.4 地震、石嘴山臺在2015年阿左旗MS5.8 地震前的觀測數(shù)據(jù)為例，提取地震前視電阻率的中短期弱幅度異常。

通渭臺距2013年甘肅岷縣-漳縣MS6.6 地震約125km，從原始月均值觀測曲線上可以看出，N20°W測道(圖3(a))在2013年存在較為突出的下降變化； EW方向長極距EW測道(圖3(b))和短極距EW′測道(圖3(c))的趨勢變化和年變化較為平穩(wěn)，均難以識別較為顯著的異常變化。采用月均值去年變處理后，可以分辨出N20°W測道的下降異常(圖3(d))，但EW和EW′測道仍然難以識別出異常(圖3(e)～(f))。其中，N20°W測道下降幅度為1.04%，略微超過傳統(tǒng)的異常閾值1%(杜學(xué)彬，2010)，而EW′測道的下降幅度僅0.37%，遠(yuǎn)小于1%。采用自適應(yīng)變化幅度方法進(jìn)行分析，N20°W測道(圖3(g))和EW′測道(圖3(i))變化幅度均超過異常閾值線，EW測道則位于閾值線之內(nèi)(圖3(h))。N20°W和EW′測道在震前出現(xiàn)異常變化，而EW測道則未出現(xiàn)異常，與已有的分析結(jié)果一致(杜學(xué)彬等，2013)。

圖3 2013年岷縣-漳縣 MS6.6 地震前通渭臺視電阻率異常變化(a)NW測道月均值； (b)EW測道月均值； (c)EW′測道月均值； (d)NW測道月均值去年變(藍(lán)色實(shí)線為趨勢變化)；(e)EW測道月均值去年變； (f)EW′測道月均值去年變； (g)NW測道相對變化幅度(紅色虛線為變化幅度的2.5倍均方差)； (h)EW測道相對變化幅度； (i)EW′測道相對變化幅度

圖4 南天山西段8次MS5.0以上地震前柯坪臺視電阻率異常變化(a)柯坪臺NS測道月均值；(b)月均值去年變；(c)相對變化幅度

表1 柯坪臺異常期間200km范圍內(nèi)5級以上地震

柯坪臺距離2020年伽師MS6.4 地震約172km，從NS測道月均值原始曲線上可以看出，在2018年和2019年，年變化幅度略微有所減小(圖4(a))，但是難以將其判定為異常。從去年變后曲線上來看(圖4(b))，由于整體存在趨勢下降，在此背景上識別出中短期異常仍然比較困難。采用自適應(yīng)變化幅度方法分析后，在2018年和2019年2次出現(xiàn)下降幅度超過閾值線的異常(圖4(c))，分別對應(yīng)于圖4(b)中的2次加速變化。2018年異常出現(xiàn)后，距離臺站124km發(fā)生新疆阿圖什MS5.1 地震，2019年異常出現(xiàn)后，臺站200km范圍內(nèi)發(fā)生7次5級以上地震，最大為伽師MS6.4 地震(表1)。2次異常的下降幅度均低于1%，最大幅度為0.886%。

圖5 2015年阿左旗 MS5.8 地震前石嘴山臺視電阻率異常變化(a)石嘴山臺NW測道月均值； (b)月均值去年變； (c)相對變化幅度

3 討論

異常分析方法的目的，是從原始觀測數(shù)據(jù)中將有別于背景常態(tài)變化的信息凸顯出來，其提取的異常信息是觀測數(shù)據(jù)本身的數(shù)據(jù)異常。而數(shù)據(jù)異常產(chǎn)生的原因很多，需要回歸到原始觀測曲線，進(jìn)一步分析該部分?jǐn)?shù)據(jù)異常的原因，只有在排除了環(huán)境干擾和觀測系統(tǒng)故障等原因之后，才能作為地震異常開展后續(xù)的工作。此外，觀測數(shù)據(jù)的變化形態(tài)是十分復(fù)雜的，即使排除了干擾因素的影響，從異常信度方面考慮，還需要該異常變化形態(tài)與已有地震異常相似，在此基礎(chǔ)上開展后續(xù)的地震預(yù)測工作，才能減少虛報。因此，在開展基于異常的地震預(yù)測之前，對提取的數(shù)據(jù)異常進(jìn)行甄別和解釋是十分重要的。

與地震晚期孕育過程有關(guān)的視電阻率異常類型，包括趨勢轉(zhuǎn)折、趨勢背景變化基礎(chǔ)上的下降/上升、年變形態(tài)畸變以及短臨階段的加速變化等。不同類型的異常通常不是單獨(dú)出現(xiàn)，而是經(jīng)常疊加在一起。任何類型的異常，都應(yīng)歸結(jié)為觀測值下降或上升變化，不同類型異常只是形態(tài)表象上的描述，在分析異常機(jī)理以及異常與地震之間的關(guān)系時，需要將異常還原為下降/上升變化以及與之對應(yīng)的變化幅度、起始時間和變化速率等。

自適應(yīng)變化幅度方法采用了基于分析時段內(nèi)背景常態(tài)變化幅度的異常閾值，有利于提取弱幅度異常，但與之相對應(yīng)的是，該方法對觀測數(shù)據(jù)穩(wěn)定性也提出了較高的要求。在進(jìn)行分析之前，需要通過預(yù)處理消除突跳、階躍變化及其他影響數(shù)據(jù)穩(wěn)定性的干擾變化。對于趨勢轉(zhuǎn)折變化，在轉(zhuǎn)折時段會偏離原有趨勢線，在采用該方法進(jìn)行分析時會超過異常閾值線。全國多數(shù)臺站存在不同程度的趨勢變化，而震例總結(jié)的結(jié)果顯示，趨勢轉(zhuǎn)折變化與臺站周圍地震之間的對應(yīng)關(guān)系并不理想，可能反映大尺度區(qū)域應(yīng)力場的調(diào)整(杜學(xué)彬，2010)，而與地震有關(guān)的異常更多表現(xiàn)為在趨勢背景變化基礎(chǔ)上的中短期變。此外，趨勢轉(zhuǎn)折可以從原始觀測數(shù)據(jù)上直觀地予以識別，無需采用其他分析方法。自適應(yīng)變化幅度方法并不適用于趨勢轉(zhuǎn)折類異常的分析，因此，對視電阻率進(jìn)行分析時，需要根據(jù)趨勢變化，選取趨勢變化較為一致的時段對觀測數(shù)據(jù)進(jìn)行分段處理。

由于該方法在利用傅立葉滑動方法去年變時，采用了正向擬合和逆向擬合相結(jié)合的方式，用于分析的數(shù)據(jù)長度至少需要2年。分析異常時，需要有正常的背景變化，視電阻率異常通常持續(xù)數(shù)月至2年尺度。因此，用于分析的數(shù)據(jù)總長度宜大于3年。

4 結(jié)論

本文提出一種自適應(yīng)變化幅度方法，用于分析視電阻率中短期異常。以去年變數(shù)據(jù)與趨勢背景變化的交點(diǎn)作為計算之后數(shù)據(jù)變化幅度的起點(diǎn)，克服了采用固定長度窗口滑動類方法無法獲取單次異常期間累計變化幅度的困難； 采用基于分析時段內(nèi)背景常態(tài)變化幅度的異常閾值標(biāo)準(zhǔn)，有利于識別弱幅度異常。該方法對變化幅度較為敏感，對觀測數(shù)據(jù)的穩(wěn)定性要求較高，分析之前需要進(jìn)行預(yù)處理，剔除突跳、階躍等干擾變化，并根據(jù)趨勢變化進(jìn)行分段處理。