基于SBAS時序分析方法的斷裂帶形變監(jiān)測

羅 想，王 萍，張艷梅，張 蕎，陳 慧

（1.國家測繪地理信息局第三航測遙感院，四川 成都 610100）

海原活動斷裂帶位于青藏高原東北邊緣，東起寧夏哨口西至甘肅興泉堡，全長約為240 km，是中國大陸西部重要的活動地塊邊界構(gòu)造帶和強震活動帶。1920年的海原8.5級大地震導(dǎo)致該斷裂帶10.5 m的左旋走滑位移和幾乎整個斷裂帶的貫通性破裂。近一個世紀(jì)過去了，斷裂帶現(xiàn)今的活動性如何，地震危險性多大一直是人們關(guān)注的熱點[1-2]。

合成孔徑雷達干涉測量（InSAR）技術(shù)作為一種空對地大范圍觀測的有效手段，利用SAR影像中攜帶的相位信息獲取地表高程信息及其變化，有效克服了精密水準(zhǔn)測量、靜態(tài)或動態(tài)GPS測量等常規(guī)測量方法工作量大、成本高、效率低等問題[2]。短基線集干涉測量（SBAS-InSAR）技術(shù)是2002年由Berardino等針對傳統(tǒng)D-InSAR中時間、空間失相關(guān)和大氣效應(yīng)等問題而提出的形變時序分析方法，極大地提高了InSAR監(jiān)測精度，可應(yīng)用于地殼長時間緩慢形變監(jiān)測[3]。

本文基于20景Sentinel-1數(shù)據(jù)， 利用SBAS-InSAR技術(shù)，以祁連山海原斷裂帶為研究區(qū)，分析了海原斷裂帶近年來的地殼形變規(guī)律和特點。

1 數(shù)據(jù)源與技術(shù)原理

1.1 Sentinel-1數(shù)據(jù)

Sentinel-1是一個全天時、全天候、極地軌道雷達成像系統(tǒng)，由歐洲航天局“哥白尼全球監(jiān)測計劃”中的兩顆對地觀測衛(wèi)星組成，主要用于城市規(guī)劃與環(huán)境制圖、地表形變監(jiān)測、海洋環(huán)境監(jiān)控與海冰監(jiān)測、土壤保護與水體管理、森林火災(zāi)監(jiān)測等。Sentinel-1 A、B兩星分別于2014年4月3日和2016年4月26日發(fā)射升空，單顆衛(wèi)星的重訪周期為12 d，雙星系統(tǒng)重訪周期為6 d[4]。

Sentinel-1參數(shù)設(shè)計延續(xù)了ERS1/2和ENVISAT的基本特征，搭載有C波段合成孔徑雷達傳感器，工作頻率為5.4GHz；包括條帶成像（SM）、干涉寬幅（IW）、超寬幅（EW）和波浪（WV）4種工作模式[4-5]，不同成像模式參數(shù)見表1。

表1 Sentinel-1主要參數(shù)指標(biāo)

歐空局對外免費開放SM、IW和EW三個模式的Level-0、Level-1級別產(chǎn)品，以及各模式下的Level-2 級別產(chǎn)品。Level-0級產(chǎn)品為壓縮、未聚焦的原始數(shù)據(jù)，包括回波數(shù)據(jù)、校準(zhǔn)參數(shù)、衛(wèi)星姿態(tài)和軌道信息等；Level-1級產(chǎn)品分為單視復(fù)數(shù)影像（SLC）和地距影像（GRD）兩類，SLC包括聚焦的SAR數(shù)據(jù)、基于衛(wèi)星軌道和姿態(tài)的斜距模式地理參考信息，GRD包括多視處理、采用WGS84橢球投影至地距的聚焦數(shù)據(jù)；Level-2 級產(chǎn)品為海洋產(chǎn)品，主要包括海洋風(fēng)場、海洋膨脹光譜、地表徑向速度等信息[5]。

1.2 SBAS-InSAR技術(shù)原理

SBAS-InSAR技術(shù)原理是將SAR影像分成若干個較短時空基線子集，再利用像對差分干涉和奇異值分解求得整個時間序列的形變[3]。

假設(shè)有時間序列t0,…,tN的N+1幅同軌SAR影像，以任意影像為主影像進行配準(zhǔn)，通過設(shè)定時間、空間基線閾值得到一個或若干個基線集合，再分別對每個集合內(nèi)的像對進行差分干涉處理，得到M幅差分干涉圖，且有：

對于時刻ta和tb獲取的兩幅SAR影像生成的第 j (j=1,…,M)幅差干涉圖，去除平地相位和地形相位，若ta時刻早于tb時刻，則任意像元(x,r)對應(yīng)的解纏相位為：

式中，λ為波長；d(ta,x,r)和d(tb,x,r)分別為時刻i和j相對于初始時刻t0的視線方向累積形變量。但式（2）未考慮殘余地形相位、大氣延遲相位、時空失相干等噪聲相位的影響。

因有M幅差分干涉圖，故可得到M個方程組，用矩陣表示該方程組為：

式中，A為M×N的近似關(guān)聯(lián)矩陣，每一行對應(yīng)一個解纏差分干涉相位，每一列對應(yīng)相位未知參數(shù)，當(dāng)M≥N，可用最小二乘法求解：

當(dāng)SAR影像被分為多個子集時，矩陣A秩虧，ATA為奇異矩陣，假設(shè)有L個基線子集，那么A的秩為N-L+1，則方程組有無窮多個解。Berardino等提出奇異值分解（SVD）方法，將相鄰時間獲取的平均相位速率作為未知量求解，通過矩陣A的廣義逆來求得方程最小范數(shù)下的最小二乘解，從而得到時間序列間的形變速率和累計形變量[6-8]。

2 海原斷裂帶形變監(jiān)測

2.1 研究區(qū)與實驗數(shù)據(jù)

研究區(qū)位于寧夏海原縣境內(nèi)的祁連山海原斷裂帶（圖1），區(qū)域內(nèi)主要為黃土質(zhì)山區(qū)，植被稀少，氣候干燥。本文收集的數(shù)據(jù)為2015年3月14日～2016年7 月30日20景Sentinel-1 A IW模式的覆蓋整個研究區(qū)域的SLC數(shù)據(jù)，極化方式為VV，分辨率為5 m×20 m，具體見表2；同時還收集了對應(yīng)的AUX_POEORB精密軌道參數(shù)文件以及NASA獲取的空間分辨率為90 m的SRTM4參考DEM數(shù)據(jù)。

圖1 研究區(qū)范圍和地震構(gòu)造環(huán)境

2.2 數(shù)據(jù)處理

本文以ENVI SARScape為軟件平臺，對研究區(qū)域20景Sentinel-1影像進行SBAS-InSAR處理，主要包括差分干涉像對的生成、相干目標(biāo)選擇、差分干涉圖解纏、軌道精煉及重去平、形變速率估算、地形殘余相位和大氣相位去除等步驟。

表2 實驗SAR影像與參數(shù)

1）干涉連接圖生成。由于垂直基線都遠遠小于臨界基線（約5 440 m），這里不用設(shè)置基線閾值，設(shè)置時間基線為200 d。以20150922這一期為超級主影像，對SAR影像進行組合，形成短基線集，得到117個干涉像對。

2）干涉生成到相位解纏，主要包括干涉圖生成、去平地效應(yīng)、干涉圖濾波、相干系數(shù)計算、相位解纏。本文采用Goldstein濾波方法以及基于Delaunay三角網(wǎng)的最小費用流法解纏，解纏相干系數(shù)閾值設(shè)為0.35，得到的部分結(jié)果如圖2所示。通過相干性和解纏結(jié)果的對比篩選，剔除了23對相干性差、解纏效果不理想的像對。

3）軌道精煉和重去平。選取25個覆蓋所有像對、遠離形變區(qū)域且位于有效解纏結(jié)果上的GCP點，對所有像對進行軌道精煉和相位偏移計算，消除可能的斜坡相位。

4）形變速率和高程系數(shù)估算。第一次估算位移速率和殘余地形，用以對合成的干涉圖進行去平，重新作相位解纏和精煉，生成更優(yōu)化的結(jié)果；基于三次方模型計算得到所有像對的形變（形變、加速度、加速度變化）和高程（校正值和新的DEM）。

圖2 相干系數(shù)圖、去平后干涉圖和相位解纏圖

5）大氣相位和地形殘余相位去除。在第一次得到的形變速率基礎(chǔ)上進行定制的大氣濾波，從而估計和去除大氣相位，得到更加純凈的時間序列上的最終位移結(jié)果。

2.3 結(jié)果分析

20景Sentinel-1A數(shù)據(jù)經(jīng)過SBAS-InSAR技術(shù)處理后，得到研究區(qū)平均地表形變結(jié)果如圖3所示。正、負(fù)值分別表示地表抬升和沉降，監(jiān)測分析表明：

圖3 平均形變速率圖

1）2015年3月～2016年7月，隨時間推移，研究區(qū)地表形變量在不斷積累，但形變量整體較小，1.3 a形變累積量集中在-10～+5 mm，其中最大形變累積量位于海原縣第一中學(xué)附近，達28 mm。進一步分析發(fā)現(xiàn)，沉降漏斗區(qū)域恰好是近年海原縣建設(shè)發(fā)展最為劇烈的地區(qū)之一，大量的人工建筑和人類活動是城市地表下沉的主要原因。

2）研究區(qū)平均形變速率具有良好的趨勢性，呈明顯的跨斷層形變梯度。斷裂帶北盤形變速率整體為正值，呈抬升趨勢，斷裂帶南盤形變速率為負(fù)值，呈沉降趨勢，南北兩盤相對形變速率約為6 mm/a。假設(shè)研究區(qū)目標(biāo)只做水平運動，沒有垂直運動，則根據(jù)海原斷裂帶南北兩盤形變位移差異可判斷海原斷裂帶呈左旋走滑趨勢。這一研究結(jié)果與基于地形地貌特征和測年的地質(zhì)學(xué)方法得到的海原斷裂帶呈3～8 mm/a左旋走滑速率，基于GPS等大地測量方法得到的海原斷裂帶呈3～10 mm/a等研究結(jié)果一致[9-10]。

3 結(jié) 語

本文基于SBAS-InSAR技術(shù)，利用20景Sentinel-1數(shù)據(jù)對海原斷裂帶2015年3月～2016年7月的地表形變進行了監(jiān)測，獲得了時間序列差分干涉圖，反演、解析了形變量與形變速率，并利用已有研究和監(jiān)測資料佐證了研究結(jié)果的可靠性。本文為海原斷裂帶地區(qū)地震預(yù)測提供了可靠的斷裂活動數(shù)據(jù)，同時有助于Sentinel-1數(shù)據(jù)和SBAS-InSAR技術(shù)在活動斷裂帶形變、城市地表沉降以及地理國情監(jiān)測中推廣應(yīng)用。

[1] 崔篤信,胡亞軒,王文萍,等.海原斷裂帶庫侖應(yīng)力積累[J].地球科學(xué)(中國地質(zhì)大學(xué)學(xué)報),2009,34(4):641-650

[2] 楊成生,張勤,趙超英,等.短基線集InSAR技術(shù)用于大同盆地地面沉降、地裂縫及斷裂活動監(jiān)測[J].武漢大學(xué)學(xué)報(信息科學(xué)版),2014,39(8):945-950

[3] 胡樂銀,張景發(fā),商曉青.SBAS-InSAR技術(shù)原理及其在地殼形變監(jiān)測中的應(yīng)用[J].地殼構(gòu)造與地殼應(yīng)力文集,2010(22):82-89

[4] 李永生,馮萬鵬,張景發(fā),等.2014年美國加州納帕Mw6.1地震斷層參數(shù)的Sentinel-1A InSAR反演[J].地球物理學(xué)報,2015,58(7):2 339-2 349

[5] 楊魁,楊建兵,江冰茹.Sentinel-1衛(wèi)星綜述[J].城市勘測,2015(2):24-27

[6] 羅鋮.基于SBAS-InSAR的西安地表沉降監(jiān)測[D].西安:長安大學(xué),2012

[7] 王如意.基于SBAS-InSAR的高分辨率地面沉降監(jiān)測技術(shù)研究[D].北京:中國地質(zhì)大學(xué),2015

[8] 李永生,張景發(fā),李振洪.利用短基線集干涉測量時序分析方法監(jiān)測北京市地面沉降[J].武漢大學(xué)學(xué)報(信息科學(xué)版),2013,38(11):1 374-1 377

[9] 屈春燕,單新建,張國宏,等.時序INSAR斷層活動性觀測研究進展及若干問題探討[J].地震地質(zhì),2014,36(3):731-748

[10] 屈春燕,單新建,宋小剛,等.基于PSInSAR技術(shù)的海原斷裂帶地殼形變初步研究[J].地球物理學(xué)報,2011,54(4):984-993