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

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

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

瞬態(tài)形變是由地殼構(gòu)造運(yùn)動(如震后地殼松弛形變、板塊邊界走滑斷層、消減帶慢地震等)導(dǎo)致的位移變化[1]。GPS觀測精度突破毫米級使得對坐標(biāo)序列中的瞬態(tài)形變信息進(jìn)行探測成為可能。GPS坐標(biāo)時(shí)間序列中瞬態(tài)形變探測在早期主要是基于物理模型來反演坐標(biāo)序列中的瞬態(tài)變形信號[2- 4]，該類方法實(shí)現(xiàn)起來稍復(fù)雜，且適用于監(jiān)測網(wǎng)較小的情況，在斷層幾何模型或格林函數(shù)相關(guān)的先驗(yàn)信息不準(zhǔ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)造運(yùn)動活躍地區(qū)，探測GPS觀測序列中的瞬態(tài)形變信號也一直是研究的熱點(diǎn)問題[12- 14]。

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

1 RSI基本原理

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

(1)

(2)

式中，Uj和Dj分別為在指定時(shí)間段n內(nèi)的漲幅與跌幅。采用簡易移動平均法，得到漲幅與跌幅的比率RS，通過式(2)得到該時(shí)間的RSI，其范圍為0～100。通過RSI值可以判斷各個時(shí)間段內(nèi)的整體趨勢是上升還是下降。如股市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í)空噪聲，對RSI分析產(chǎn)生了不確定性，因此在分析前需對原始序列進(jìn)行濾波，去除部分共模誤差和高頻噪聲；且RSI值選取范圍也變大，為(20,80)。對于濾波后的時(shí)間序列，采取n=21 d的簡易滑動平均法計(jì)算各個時(shí)間段的RSI，并定義RSI值大于50的所有RSI的平均值為上限閾值，小于50的所有RSI的平均值為下限閾值。在計(jì)算出所有測站各個方向的RSI值和閾值后，計(jì)算各測站RSI超出閾值數(shù)，剔除在最大閾值與最小閾值之間的RSI值。本文選取連續(xù)數(shù)超過7 d的RSI值，將其作為可能發(fā)生瞬態(tài)形變的時(shí)間段，因?yàn)榇蟛糠致剖录r(shí)間跨度是超過7 d的。

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

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

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

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

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

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

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

為了利用探測得到瞬態(tài)形變信息改善GPS坐標(biāo)時(shí)間序列，使用RSI探測到的瞬態(tài)形變信息對每個時(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ù)，用來描述瞬態(tài)形變的發(fā)生。上文通過RSI探測到瞬態(tài)形變發(fā)生的時(shí)間及形變量，結(jié)合探測信息進(jìn)行建模。圖4為利用RSI方法探測的瞬態(tài)形變建模后的殘差時(shí)間序列，可以看出改正后的時(shí)間序列一定程度上減小了階躍的影響。使用方差減少率來評定時(shí)間序列的改進(jìn)，公式如下

(4)

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

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

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

3 結(jié) 語

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