基于RSI的GPS坐標時間序列瞬態(tài)地殼形變信號探測

成國輝，羅 勇，劉 斌，匡翠林

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

瞬態(tài)形變是由地殼構(gòu)造運動(如震后地殼松弛形變、板塊邊界走滑斷層、消減帶慢地震等)導(dǎo)致的位移變化[1]。GPS觀測精度突破毫米級使得對坐標序列中的瞬態(tài)形變信息進行探測成為可能。GPS坐標時間序列中瞬態(tài)形變探測在早期主要是基于物理模型來反演坐標序列中的瞬態(tài)變形信號[2- 4]，該類方法實現(xiàn)起來稍復(fù)雜，且適用于監(jiān)測網(wǎng)較小的情況，在斷層幾何模型或格林函數(shù)相關(guān)的先驗信息不準確的情況下，給形變信息的提取帶來較大的偏差。南加州地震中心(southern California earthquake center,SCEC)于2010年提出基于SCIGN(sdouthern California integrated GPS network)監(jiān)測網(wǎng)開展瞬態(tài)形變信息探測的研究[5]，各學(xué)者陸續(xù)提出并發(fā)展了一些探測瞬態(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)造運動活躍地區(qū)，探測GPS觀測序列中的瞬態(tài)形變信號也一直是研究的熱點問題[12- 14]。

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

1 RSI基本原理

RSI方法在股票市場中被證明是一種很有價值的技術(shù)指標，其計算公式為[15]

(1)

(2)

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

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

2 基于RSI的瞬態(tài)形變探測與GPS坐標時間序列改正

2.1 RSI瞬態(tài)形變探測實例

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

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

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

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

對比圖1與圖3，在圖1畫線處的瞬態(tài)形變基本都可以探測出，表1是對這10個測站發(fā)生瞬態(tài)形變的記錄的統(tǒng)計分析。從表中也可以看出，瞬態(tài)形變存在12～18個月的周期性，持續(xù)天數(shù)長短不同，最大偏移量不等(最大偏移量是該時段內(nèi)最大最小位移值之差)，同時段發(fā)生瞬態(tài)形變的測站數(shù)不同。同時段發(fā)生瞬態(tài)形變的測站數(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)偏移更大，說明sc03站受慢滑移事件影響是比較大的。

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

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

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

(3)

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

(4)

表1 試驗中10個GPS測站的慢滑移事件記錄

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

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

3 結(jié) 語

瞬態(tài)形變在地球物理中有著重要作用，若其沒有被正確探測與修正，將影響GPS坐標時序的噪聲特性及其測站的速率估計，從而會直接影響GPS坐標序列的應(yīng)用水平。本文將RSI方法應(yīng)用到GPS時間序列地殼瞬態(tài)形變信號探測中，實例應(yīng)用表明，RSI是一種有效的探測手段，且能給出瞬態(tài)形變發(fā)生的時間起點、持續(xù)時間及形變大小等詳細信息。利用探測到的瞬態(tài)形變信息對時間序列進行建模改正，結(jié)果表明，使用偏移改正后的殘差時間序列方差有明顯減小，有利于進一步準確分析GPS坐標時間序列噪聲特性和進行測站坐標速率估計。