利用Sentinel-1 SAR數(shù)據(jù)及SBAS技術(shù)的大區(qū)域地表形變監(jiān)測

張金盈，崔 靚，劉增珉，王新田，林 琳，徐鳳玲

(1.山東省國土測繪院，山東 濟(jì)南 250013; 2.山東省水利科學(xué)研究院，山東 濟(jì)南 250014; 3.山東華峰地理信息科技有限公司，山東 濟(jì)南 250101)

星載SAR系統(tǒng)對全球地表形變監(jiān)測具有重要作用。近幾年廣泛應(yīng)用于城市、礦區(qū)、滑坡監(jiān)測的SAR數(shù)據(jù)有如下衛(wèi)星數(shù)據(jù)：2014年日本發(fā)射的ALOS衛(wèi)星，幅寬最大為70 km；2007年德國發(fā)射的TerraSAR衛(wèi)星，提高了干涉數(shù)據(jù)的相干性，幅寬最大為30 km；2007年加拿大發(fā)射的radarsat衛(wèi)星，提供高分辨率全極化影像，幅寬最大為150 km。這幾種數(shù)據(jù)的共性為覆蓋范圍較小，在大范圍監(jiān)測中面臨數(shù)據(jù)拼接等復(fù)雜問題。2014年歐空局發(fā)射Sentinel-1衛(wèi)星，該衛(wèi)星具有4種成像模式，最大幅寬達(dá)400 km，滿足大范圍的數(shù)據(jù)處理。同時歐空局承諾為用戶免費(fèi)提供哨兵家族數(shù)據(jù)以供科學(xué)研究，大大促進(jìn)了時序雷達(dá)干涉測量(SAR interferometry,InSAR)技術(shù)的發(fā)展和在大區(qū)域的應(yīng)用[1]。

相比永久散射體合成孔徑雷達(dá)干涉測量(PS-InSAR技術(shù))，小基線集(small baseline subsets InSAR，SBAS InSAR)技術(shù)是在限制時空基線閾值的前提下，對獲取的SAR影像集進(jìn)行自由干涉組合，從而得到一系列短基線差分干涉相位圖集。這些差分干涉圖能夠較好地克服時空失相關(guān)的影響[2]，可以顯著增加相干影像對的數(shù)目，一定程度上解決了大范圍時間、空間失相關(guān)和大氣效應(yīng)對InSAR的影響。

黃河三角洲區(qū)域?yàn)辄S河入海口，地勢低平，地質(zhì)成因晚，擁有油田、天然氣等豐富礦產(chǎn)資源，由于地質(zhì)原因和人類生產(chǎn)活動，存在地表沉降現(xiàn)象[3]。文獻(xiàn)[4]利用SBAS技術(shù)，采用39景ERS1/2數(shù)據(jù)對面積為5400 km2的東營地區(qū)1992年至2000年間的地表情況進(jìn)行了監(jiān)測，監(jiān)測出勝利油田中心開采區(qū)最大沉降速率為-33.2 mm/a。文獻(xiàn)[5]利用28景ALOS數(shù)據(jù)，監(jiān)測東營市勝利油田2007—2011年間的地表沉降，分區(qū)域分析了地表沉降及成因。大面積區(qū)域性地面沉降調(diào)查存在數(shù)據(jù)相干性問題、低相干區(qū)沉降信息難提取、大氣波動影響等問題[6]。本文利用Sentinel-1數(shù)據(jù)，采用SBAS方法對覆蓋山東近1/4范圍的6.3×104km2的黃河三角洲區(qū)域進(jìn)行沉降監(jiān)測研究，并進(jìn)行了精度評價(jià)和成因分析，驗(yàn)證了利用Sentinel-1數(shù)據(jù)的SABS方法在大區(qū)域地表沉降監(jiān)測的適用性，總結(jié)了黃三角區(qū)域在兩年半時間內(nèi)的地表沉降時空變化的主要特點(diǎn)。

1 Sentinel-1數(shù)據(jù)及SBAS技術(shù)介紹

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

Sentinel-1衛(wèi)星是歐洲航天局哥白尼計(jì)劃(GMES)中的地球觀測衛(wèi)星，由兩顆衛(wèi)星組成，載有C波段合成孔徑雷達(dá)，可提供連續(xù)圖像(白天、夜晚和各種天氣)。該衛(wèi)星具有超高的輻射分辨率，能有效提升雷達(dá)圖像參數(shù)反演精度，能實(shí)現(xiàn)全球80%的覆蓋，單個衛(wèi)星在中國境內(nèi)12 d的重訪周期且實(shí)現(xiàn)80%覆蓋。兩個衛(wèi)星星座可縮短至6 d。其干涉寬測繪帶模式(interferometric wide swath，IW)帶寬達(dá)250 km，單景影像可覆蓋面積為4.5×104m2，空間分辨率為5 m×20 m，兩景影像即可實(shí)現(xiàn)黃河三角洲區(qū)域全覆蓋。

1.2 SBAS技術(shù)介紹

SBAS方法通過增加干涉數(shù)據(jù)的采樣，獲得較高精度的沉降信息和高程[7]，即在一定范圍內(nèi)設(shè)置空間基線和時間基線，將SAR影像組合成若干個小基線集，分別選擇一個主影像進(jìn)行配準(zhǔn)，利用差分干涉測量的方法獲得每個集合中的差分干涉圖[9]，根據(jù)相關(guān)系數(shù)提取相干性較高的點(diǎn)，應(yīng)用奇異值分解法(singular value decomposition，SVD)將多個小基線集聯(lián)合起來求解，最終獲得時間序列沉降圖及沉降速率[8]。

假設(shè)不考慮由于地形、大氣延遲、噪聲等造成的相位誤差，可建立如下模型

Aφ=δφ

(1)

式中，A為M×N稀疏系數(shù)矩陣，每一行對應(yīng)一幅差分干涉圖，每一列對應(yīng)不同采集日期的同一幅影像，主影像的系數(shù)為1，輔影像的系數(shù)為-1，其他為0；φ為N+1個成像時刻的雷達(dá)圖像上未知沉降相位點(diǎn)組成的矩陣；δφ為M幅差分干涉圖中觀測相位值矩陣；由于SAR影像由多個小基線集組成，M大于N，系數(shù)矩陣A秩虧。對矩陣A進(jìn)行奇異值分解(SVD)

A=USVT

(2)

式中，U為M×N正交矩陣，由AAT的特征向量ui組成；V為N×M的正交矩陣，由ATA的特征向量vi組成；S為一個M×N的對角矩陣，對角線元素是AAT的前R個非零特征值，后面的M-R個特征值為0，A+為A的偽逆矩陣

(3)

經(jīng)過奇異值分解，可以得到沉降相位φ的最小范數(shù)意義上的最小二乘解

(4)

解算出線性沉降速率和高程誤差，去除線性部分，得到殘余相位。對殘余相位在空間上進(jìn)行低通濾波，在時間上進(jìn)行高通濾波，分離出大氣相位，就可以得到非線性沉降。將非線性沉降疊加到線性沉降中，就可以求出總的地表沉降信息[10-11]。

2 試驗(yàn)處理流程及結(jié)果

2.1 區(qū)域概況

研究區(qū)域是黃河攜帶泥沙在入海口淤積形成的沖積平原，地勢低平。近年來，由于人類日益頻繁的經(jīng)濟(jì)活動和黃河帶來的巨量沉積物的自重固結(jié)影響，黃河三角洲區(qū)域出現(xiàn)明顯的地表沉降。根據(jù)已有研究，1992—2000年該地區(qū)最大沉降速率達(dá)-33.2 mm/a[4]；2015年5月—2016年4月期間，該地區(qū)最大年平均沉降速率為224.6 mm/a[12]。

2.2 SBAS技術(shù)處理流程

本文研究使用Sentinel-1A IW模式升軌影像進(jìn)行監(jiān)測，對相鄰2景影像進(jìn)行拼接，共采用32期64景Sentinel-1A影像，時間跨度為2015年7月—2017年8月，極化方式為VV極化，C波段。處理范圍如圖1所示，該測區(qū)覆蓋東營、濰坊、淄博、濱州等地，陸地面積約為6.3×104km2。軌道數(shù)據(jù)使用與Sentinel-1A數(shù)據(jù)間隔21 d分發(fā)的精密軌道數(shù)據(jù)，為了精細(xì)化模擬地形相位，采用30 m分辨率的SRTM DEM數(shù)據(jù)。

驗(yàn)證數(shù)據(jù)采用均勻覆蓋山東全省的SDCORS基準(zhǔn)網(wǎng)和2015年12月—2016年12月逐年觀測的GPS數(shù)據(jù)點(diǎn)。

根據(jù)SBAS網(wǎng)絡(luò)構(gòu)建原則，通過設(shè)定不同的時間閾值及空間閾值，優(yōu)化影像組合方式，以提高干涉圖的相干性。根據(jù)干涉的效果和垂直基線大小，選取基線小于100 m時間基線小于120 d的194對干涉圖，如圖2所示。采用分辨率為30 m的SRTM DEM生成模擬地形相位如圖3所示，然后對干涉圖進(jìn)行差分干涉處理[13-18]。

圖2 SAR影像組合連接

圖3 SRTM DEM

首先提取相干系數(shù)大于0.4的點(diǎn)為永久散射體點(diǎn)，然后利用Delaunay不規(guī)則三角網(wǎng)建立候選點(diǎn)之間的連接關(guān)系，對永久散射體點(diǎn)進(jìn)行解纏和地理編碼，最后利用矩陣奇異值分解的方法估算沉降速率和高程系數(shù)，對結(jié)果進(jìn)行時間域的高通濾波，對每幅干涉圖進(jìn)行空間域的低通濾波，從而獲取大氣延遲相位，將其從原始差分干涉圖中去除，得到去除了大氣影響的沉降結(jié)果。黃三角區(qū)域監(jiān)測時段內(nèi)形變速率如圖4所示。

圖4 SBAS方法監(jiān)測黃河三角洲地區(qū)形變速率

3 可靠性驗(yàn)證與分析

3.1 利用CORS站點(diǎn)的精度驗(yàn)證

為驗(yàn)證該方法反演地表沉降的精度，利用收集到的監(jiān)測區(qū)域2016年3月和2017年3月兩期CORS數(shù)據(jù)點(diǎn)沉降觀測數(shù)據(jù)進(jìn)行驗(yàn)證。

根據(jù)驗(yàn)證點(diǎn)的地理坐標(biāo)，采用最近點(diǎn)法搜索站點(diǎn)5個像元(150 m)內(nèi)的相干目標(biāo)點(diǎn)，共搜索到符合條件的17組CORS驗(yàn)證數(shù)據(jù)如圖5所示。17組驗(yàn)證數(shù)據(jù)見表1，分別對每組的相干點(diǎn)沉降觀測值和驗(yàn)證點(diǎn)沉降測量值進(jìn)行時間序列的曲線擬合，部分驗(yàn)證點(diǎn)時間序列沉降圖如圖6所示。SBAS技術(shù)觀測的相干點(diǎn)變化趨勢與CORS站及GPS點(diǎn)解算的沉降變化趨勢基本一致，解算中誤差為7.8 mm。

圖5 17組CORS站位置

圖6 部分SBAS-InSAR觀測值和CORS站觀測值曲線擬合

表1 17組CORS站驗(yàn)證數(shù)據(jù)組中誤差 mm

3.2 形變特征分析

圖7為研究區(qū)域提取出的沉降漏斗，地面沉降嚴(yán)重區(qū)域主要分布在黃河下游(A、B區(qū)域)，平均沉降速率在50～152 mm/a。A區(qū)域主要呈現(xiàn)6個沉降明顯的漏斗，面積為1500 km2，最大沉降速率達(dá)152 mm/a。實(shí)地調(diào)查發(fā)現(xiàn)該區(qū)域分布石油化工、鋼板廠等多家高耗水企業(yè)。同時通過文獻(xiàn)收集[19]，發(fā)現(xiàn)該沉降區(qū)域與地下水降落漏斗區(qū)域基本重合。

圖7 黃三角地區(qū)平均沉降速率點(diǎn)分布

B區(qū)域分布在黃三角沿海區(qū)域，存在3個典型沉降漏斗，面積為500 km2，最大沉降速率達(dá)111.72 mm/a。通過現(xiàn)場調(diào)查B地區(qū)，主要分布鹽場和溴素等工廠，開采地下鹵水造成一定程度沉降，同時該區(qū)域地層沉積較晚，地層巖性以粉土為主，因此土的先期固結(jié)較低，在自重壓力下長期固結(jié)壓密，在一定程度上也會造成地面沉降。

4 結(jié) 語

本文重點(diǎn)介紹了采用SBAS技術(shù)，利用Sentinel-1A影像獲取6.3×104km2黃河三角洲區(qū)域地表沉降的方法。與CORS站點(diǎn)測量結(jié)果相比，中誤差為6.9 mm，表明利用Sentinel-1數(shù)據(jù)結(jié)合SBAS技術(shù)可以進(jìn)行大范圍地表沉降監(jiān)測。研究結(jié)果發(fā)現(xiàn)，黃河入海口的平原區(qū)域存在多個降落漏斗，面積約2000 km2。最大沉降達(dá)152 mm/a。根據(jù)文獻(xiàn)調(diào)查及地質(zhì)資料分析，分析地表下沉主要是地下水及沿海鹵水長期超量開采且得不到補(bǔ)償導(dǎo)致。