成國(guó)輝，羅 勇，劉 斌，匡翠林

(1. 長(zhǎng)沙市規(guī)劃勘測(cè)設(shè)計(jì)研究院，湖南 長(zhǎng)沙 410007； 2. 長(zhǎng)沙理工大學(xué)交通運(yùn)輸工程學(xué)院測(cè)繪工程系，湖南 長(zhǎng)沙 410114； 3. 中南大學(xué)地球科學(xué)與信息物理學(xué)院，湖南 長(zhǎng)沙 410083)

瞬態(tài)形變是由地殼構(gòu)造運(yùn)動(dòng)(如震后地殼松弛形變、板塊邊界走滑斷層、消減帶慢地震等)導(dǎo)致的位移變化[1]。GPS觀測(cè)精度突破毫米級(jí)使得對(duì)坐標(biāo)序列中的瞬態(tài)形變信息進(jìn)行探測(cè)成為可能。GPS坐標(biāo)時(shí)間序列中瞬態(tài)形變探測(cè)在早期主要是基于物理模型來(lái)反演坐標(biāo)序列中的瞬態(tài)變形信號(hào)[2- 4]，該類方法實(shí)現(xiàn)起來(lái)稍復(fù)雜，且適用于監(jiān)測(cè)網(wǎng)較小的情況，在斷層幾何模型或格林函數(shù)相關(guān)的先驗(yàn)信息不準(zhǔn)確的情況下，給形變信息的提取帶來(lái)較大的偏差。南加州地震中心(southern California earthquake center,SCEC)于2010年提出基于SCIGN(sdouthern California integrated GPS network)監(jiān)測(cè)網(wǎng)開展瞬態(tài)形變信息探測(cè)的研究[5]，各學(xué)者陸續(xù)提出并發(fā)展了一些探測(cè)瞬態(tài)形變信息的算法：網(wǎng)絡(luò)反演濾波[6](network inversion filter，NIF)、網(wǎng)絡(luò)應(yīng)變?yōu)V波[7](network strain filter，NSF)、Kalman- PCA方法[8]、協(xié)方差描述分析[9](covariance descriptor analysis，CDA)、高斯小波變換[10]、Bayesian模型[11]等。此外，在地質(zhì)構(gòu)造運(yùn)動(dòng)活躍地區(qū)，探測(cè)GPS觀測(cè)序列中的瞬態(tài)形變信號(hào)也一直是研究的熱點(diǎn)問(wèn)題[12- 14]。

相對(duì)強(qiáng)弱指數(shù)(relative strength index，RSI)是1978年首次由Wilder提出來(lái)的，最初廣泛用于期貨交易中[15]。RSI對(duì)序列中異常信息特別敏感，因此適用于對(duì)事件的檢測(cè)。本文提出利用RSI對(duì)GPS測(cè)站序列中瞬態(tài)形變進(jìn)行探測(cè)的方法，并利用該方法對(duì)卡斯卡迪亞古陸選取的10個(gè)GPS測(cè)站坐標(biāo)時(shí)間序列進(jìn)行分析，探測(cè)發(fā)生瞬態(tài)形變測(cè)站的時(shí)間尺度、滑移大?。蛔詈罄锰綔y(cè)到的瞬態(tài)形變進(jìn)行建模，減小殘差序列的方差，完善坐標(biāo)時(shí)間序列。

1 RSI基本原理

RSI方法在股票市場(chǎng)中被證明是一種很有價(jià)值的技術(shù)指標(biāo)，其計(jì)算公式為[15]

(1)

(2)

式中，Uj和Dj分別為在指定時(shí)間段n內(nèi)的漲幅與跌幅。采用簡(jiǎn)易移動(dòng)平均法，得到漲幅與跌幅的比率RS，通過(guò)式(2)得到該時(shí)間的RSI，其范圍為0～100。通過(guò)RSI值可以判斷各個(gè)時(shí)間段內(nèi)的整體趨勢(shì)是上升還是下降。如股市RSI分析中普遍選取RSI值上下限的規(guī)則為：當(dāng)n=14時(shí)，范圍為(30,70)；當(dāng)n=20時(shí)，范圍為(40,60)。

然而，GPS數(shù)據(jù)不同于股市數(shù)據(jù)，GPS時(shí)間序列包含時(shí)空噪聲，對(duì)RSI分析產(chǎn)生了不確定性，因此在分析前需對(duì)原始序列進(jìn)行濾波，去除部分共模誤差和高頻噪聲；且RSI值選取范圍也變大，為(20,80)。對(duì)于濾波后的時(shí)間序列，采取n=21 d的簡(jiǎn)易滑動(dòng)平均法計(jì)算各個(gè)時(shí)間段的RSI，并定義RSI值大于50的所有RSI的平均值為上限閾值，小于50的所有RSI的平均值為下限閾值。在計(jì)算出所有測(cè)站各個(gè)方向的RSI值和閾值后，計(jì)算各測(cè)站RSI超出閾值數(shù)，剔除在最大閾值與最小閾值之間的RSI值。本文選取連續(xù)數(shù)超過(guò)7 d的RSI值，將其作為可能發(fā)生瞬態(tài)形變的時(shí)間段，因?yàn)榇蟛糠致剖录r(shí)間跨度是超過(guò)7 d的。

2 基于RSI的瞬態(tài)形變探測(cè)與GPS坐標(biāo)時(shí)間序列改正

2.1 RSI瞬態(tài)形變探測(cè)實(shí)例

數(shù)據(jù)來(lái)源于美國(guó)SOPAC(Scripps orbit and permanent array center)數(shù)據(jù)處理中心,使用GAMIT軟件解算并與QOCA(quasi- observation combination analysis)軟件包結(jié)合處理,選取ITRF2008框架下的60個(gè)GPS測(cè)站。測(cè)站GPS時(shí)間序列經(jīng)過(guò)PCA濾波去除了部分高階誤差，減小了不確定性；通過(guò)最小二乘取出時(shí)間序列的線性項(xiàng)、季節(jié)信號(hào)、同震和震后位移，去除時(shí)間序列中的線性項(xiàng)，并采用三倍中誤差法來(lái)探測(cè)和剔除粗差；對(duì)剔除粗差后的時(shí)間序列進(jìn)行插值補(bǔ)齊，對(duì)于小于3 d的數(shù)據(jù)缺失采用三次樣條插值，對(duì)于缺失較長(zhǎng)的數(shù)據(jù)，選擇直接跳過(guò)該段缺失的數(shù)據(jù)，從下一段連續(xù)的時(shí)間序列繼續(xù)進(jìn)行RSI分析。

本文選取10個(gè)發(fā)生瞬態(tài)形變的測(cè)站進(jìn)行分析，分別為albh、bamf、coup、nano、p397、p403、p405、p425、sc02、sc03測(cè)站。采用RSI方法對(duì)E方向時(shí)間序列的瞬態(tài)形變進(jìn)行研究分析，測(cè)站時(shí)間序列如圖1所示。通過(guò)目測(cè)可以看出，瞬態(tài)形變是往W方向偏移的，且10個(gè)測(cè)站的時(shí)間序列各自在不同的時(shí)段發(fā)生了瞬態(tài)形變，在某個(gè)固定時(shí)段，部分測(cè)站同時(shí)發(fā)生瞬態(tài)形變。

圖2是E方向10個(gè)測(cè)站的RSI值，從圖中可以看出，10個(gè)測(cè)站的RSI值整體比較穩(wěn)定，保持在50左右，說(shuō)明測(cè)站序列整體是比較平穩(wěn)的。然而，每個(gè)測(cè)站部分時(shí)間段的RSI值并不平穩(wěn)，如圖中橢圓標(biāo)記處，這些不平穩(wěn)的時(shí)間段可能就是發(fā)生瞬態(tài)形變的位置。通過(guò)計(jì)算，北方向(N)、西方向(E)和垂直方向(U)3個(gè)方向RSI平均閾值分別為52.8/47.1、53.4/46.4和52.9/47.1，由于E方向偏移最大，RSI的平均閾值也最大，因此利用E方向時(shí)間序列進(jìn)行瞬態(tài)形變探測(cè)也更具代表性。由于E方向各測(cè)站的瞬態(tài)信號(hào)整體往W方向(整體趨勢(shì)下降且RSI值都小于50)偏移，且大部分是同一時(shí)段發(fā)生，而對(duì)于探測(cè)出的極小部分往E方向(整體趨勢(shì)上升且RSI值都大于50)偏移的瞬態(tài)信號(hào)，并沒有作詳細(xì)分析。因此，本文僅對(duì)E方向的測(cè)站序列往W方向發(fā)生的瞬態(tài)形變進(jìn)行探測(cè)研究。

為了獲得發(fā)生瞬態(tài)形變的開始時(shí)間、結(jié)束時(shí)間、持續(xù)時(shí)間和位移，對(duì)圖1中的每個(gè)時(shí)間序列進(jìn)行RSI定量分析，采用簡(jiǎn)單移動(dòng)平均法，計(jì)算連續(xù)21 d的RSI，篩選出至少連續(xù)7 d RSI值小于下限閾值(46.4)的時(shí)間段，認(rèn)為該時(shí)間段內(nèi)發(fā)生了瞬態(tài)形變。圖3為10個(gè)測(cè)站RSI超出閾值數(shù)的統(tǒng)計(jì)結(jié)果，圖中RSI超出值越大，說(shuō)明在此時(shí)間段內(nèi)發(fā)生瞬態(tài)形變的偏移量越大，RSI越密集，說(shuō)明發(fā)生瞬態(tài)形變的時(shí)間段更長(zhǎng)或在此時(shí)間段內(nèi)發(fā)生瞬態(tài)形變的測(cè)站越多。此外，從圖中還可以看出，瞬態(tài)形變的出現(xiàn)具有一定的規(guī)律性，呈似年周期變化。

對(duì)比圖1與圖3，在圖1畫線處的瞬態(tài)形變基本都可以探測(cè)出，表1是對(duì)這10個(gè)測(cè)站發(fā)生瞬態(tài)形變的記錄的統(tǒng)計(jì)分析。從表中也可以看出，瞬態(tài)形變存在12～18個(gè)月的周期性，持續(xù)天數(shù)長(zhǎng)短不同，最大偏移量不等(最大偏移量是該時(shí)段內(nèi)最大最小位移值之差)，同時(shí)段發(fā)生瞬態(tài)形變的測(cè)站數(shù)不同。同時(shí)段發(fā)生瞬態(tài)形變的測(cè)站數(shù)不同，分別為：事件1，albh、bamf和sc02；事件2，albh、bamf、nano和sc02；事件3，albh、bamf、coup、sc02和sc03；事件4，albh、bamf、coup、sc02和sc03；事件5，albh、bamf、coup、p403、sc02和sc03；事件6，albh、bamf、coup、p403、sc02和sc03；事件7，p397和p425；事件8，albh、bamf、coup、p403、p405、p425、sc02和sc03；事件9，albh、coup、p397、p403、p405、p425、sc02和sc03；事件10，albh、bamf、nano、p403、sc02和sc03；事件11，albh、bamf、nano、p403、sc02和sc03；事件12，p397、p405、p425、sc02和sc03。albh、bamf、sc02和sc03瞬態(tài)形變的周期性最明顯，且sc03整體的瞬態(tài)偏移更大，說(shuō)明sc03站受慢滑移事件影響是比較大的。

2.2 時(shí)間序列中瞬態(tài)形變改正

為了利用探測(cè)得到瞬態(tài)形變信息改善GPS坐標(biāo)時(shí)間序列，使用RSI探測(cè)到的瞬態(tài)形變信息對(duì)每個(gè)時(shí)間序列利用偏移改正和不用偏移改正分別進(jìn)行建模，參考模型如下[16]

y(ti)=a+bti+csin(2πti)+dcos(2πti)+

(3)

式中，a為初始位置；b為速率；c、d、e和f分別為年、半年周期項(xiàng)系數(shù)；g為偏移系數(shù)；v為誤差；t為時(shí)間；H為階梯(heaviside step)函數(shù)，用來(lái)描述瞬態(tài)形變的發(fā)生。上文通過(guò)RSI探測(cè)到瞬態(tài)形變發(fā)生的時(shí)間及形變量，結(jié)合探測(cè)信息進(jìn)行建模。圖4為利用RSI方法探測(cè)的瞬態(tài)形變建模后的殘差時(shí)間序列，可以看出改正后的時(shí)間序列一定程度上減小了階躍的影響。使用方差減少率來(lái)評(píng)定時(shí)間序列的改進(jìn)，公式如下

(4)

表1 試驗(yàn)中10個(gè)GPS測(cè)站的慢滑移事件記錄

表2 各測(cè)站E、N、U方向的殘差時(shí)間序列偏移改正后的方差減少比例(%)

從以上分析可以得出，RSI方法是一種有效的探測(cè)方法，它通過(guò)對(duì)一個(gè)區(qū)域各個(gè)測(cè)站序列進(jìn)行RSI分析，討論可能發(fā)生形變的時(shí)刻，并通過(guò)探測(cè)到的形變時(shí)刻，使用模型進(jìn)行改進(jìn)，去除其中的瞬態(tài)形變信號(hào)，從而為序列的速率評(píng)估和噪聲估計(jì)提供便利。

3 結(jié) 語(yǔ)

瞬態(tài)形變?cè)诘厍蛭锢碇杏兄匾饔?，若其沒有被正確探測(cè)與修正，將影響GPS坐標(biāo)時(shí)序的噪聲特性及其測(cè)站的速率估計(jì)，從而會(huì)直接影響GPS坐標(biāo)序列的應(yīng)用水平。本文將RSI方法應(yīng)用到GPS時(shí)間序列地殼瞬態(tài)形變信號(hào)探測(cè)中，實(shí)例應(yīng)用表明，RSI是一種有效的探測(cè)手段，且能給出瞬態(tài)形變發(fā)生的時(shí)間起點(diǎn)、持續(xù)時(shí)間及形變大小等詳細(xì)信息。利用探測(cè)到的瞬態(tài)形變信息對(duì)時(shí)間序列進(jìn)行建模改正，結(jié)果表明，使用偏移改正后的殘差時(shí)間序列方差有明顯減小，有利于進(jìn)一步準(zhǔn)確分析GPS坐標(biāo)時(shí)間序列噪聲特性和進(jìn)行測(cè)站坐標(biāo)速率估計(jì)。