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

(1.山東省國土測繪院，山東 濟南 250013; 2.山東省水利科學(xué)研究院，山東 濟南 250014; 3.山東華峰地理信息科技有限公司，山東 濟南 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種成像模式，最大幅寬達400 km，滿足大范圍的數(shù)據(jù)處理。同時歐空局承諾為用戶免費提供哨兵家族數(shù)據(jù)以供科學(xué)研究，大大促進了時序雷達干涉測量(SAR interferometry,InSAR)技術(shù)的發(fā)展和在大區(qū)域的應(yīng)用[1]。

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

黃河三角洲區(qū)域為黃河入?？?，地勢低平，地質(zhì)成因晚，擁有油田、天然氣等豐富礦產(chǎn)資源，由于地質(zhì)原因和人類生產(chǎn)活動，存在地表沉降現(xiàn)象[3]。文獻[4]利用SBAS技術(shù)，采用39景ERS1/2數(shù)據(jù)對面積為5400 km2的東營地區(qū)1992年至2000年間的地表情況進行了監(jiān)測，監(jiān)測出勝利油田中心開采區(qū)最大沉降速率為-33.2 mm/a。文獻[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ū)域進行沉降監(jiān)測研究，并進行了精度評價和成因分析，驗證了利用Sentinel-1數(shù)據(jù)的SABS方法在大區(qū)域地表沉降監(jiān)測的適用性，總結(jié)了黃三角區(qū)域在兩年半時間內(nèi)的地表沉降時空變化的主要特點。

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

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

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

1.2 SBAS技術(shù)介紹

SBAS方法通過增加干涉數(shù)據(jù)的采樣，獲得較高精度的沉降信息和高程[7]，即在一定范圍內(nèi)設(shè)置空間基線和時間基線，將SAR影像組合成若干個小基線集，分別選擇一個主影像進行配準(zhǔn)，利用差分干涉測量的方法獲得每個集合中的差分干涉圖[9]，根據(jù)相關(guān)系數(shù)提取相干性較高的點，應(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個成像時刻的雷達圖像上未知沉降相位點組成的矩陣；δφ為M幅差分干涉圖中觀測相位值矩陣；由于SAR影像由多個小基線集組成，M大于N，系數(shù)矩陣A秩虧。對矩陣A進行奇異值分解(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)

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

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

2.1 區(qū)域概況

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

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

本文研究使用Sentinel-1A IW模式升軌影像進行監(jiān)測，對相鄰2景影像進行拼接，共采用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ù)。

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

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

圖2 SAR影像組合連接

圖3 SRTM DEM

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

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

3 可靠性驗證與分析

3.1 利用CORS站點的精度驗證

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

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

圖5 17組CORS站位置

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

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

3.2 形變特征分析

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

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

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

4 結(jié) 語

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