基于SBAS-InSAR技術(shù)和Sentinel-1A的城市地表沉降監(jiān)測(cè)

王 磊，連增增,劉 杰

（1.河南省有色金屬地質(zhì)礦產(chǎn)局第一地質(zhì)大隊(duì)，河南 鄭州 450000；2.河南理工大學(xué)測(cè)繪與國(guó)土信息工程學(xué)院，河南 焦作 454000）

地面沉降是由多種因素引起的地面標(biāo)高緩慢降低的地質(zhì)現(xiàn)象，是一種主要由人類(lèi)活動(dòng)引起的環(huán)境地質(zhì)問(wèn)題。隨著城市化進(jìn)程的加快，地面沉降這一地質(zhì)問(wèn)題日益突出，城市基礎(chǔ)設(shè)施遭到破壞，城市規(guī)劃布局、土地有效利用及地下空間開(kāi)發(fā)受到嚴(yán)重影響[1-2]。傳統(tǒng)的地面沉降監(jiān)測(cè)手段有GPS測(cè)量、水準(zhǔn)測(cè)量等，雖然精度高，但具有觀測(cè)周期長(zhǎng)、人力成本高及空間分辨率低等缺點(diǎn)[3]。相對(duì)于傳統(tǒng)大地測(cè)量手段，InSAR技術(shù)具有全天候、全天時(shí)、成本低、覆蓋范圍大、空間分辨率高等優(yōu)勢(shì)[4-5]。D-InSAR技術(shù)是基于InSAR干涉相位獲取地表形變信息的一種新方法，可以探測(cè)亞厘米級(jí)形變，但D-InSAR技術(shù)易受時(shí)空失相關(guān)、大氣延遲等因素影響[6]，進(jìn)而導(dǎo)致形變結(jié)果精度的降低，甚至導(dǎo)致錯(cuò)誤的形變提取結(jié)果，較難完成高精度、長(zhǎng)時(shí)間跨度的地表形變監(jiān)測(cè)[7]。

永久散射體合成孔徑雷達(dá)干涉測(cè)量（persistent scatterer interferometric synthetic apertureradar，PS-In?SAR）技術(shù)[8]能夠通過(guò)對(duì)地面高相干目標(biāo)點(diǎn)進(jìn)行時(shí)序差分干涉相位分析，并剔除低頻大氣的影響，較好地克服失相干與大氣延遲影響問(wèn)題[9]，但要得到比較可靠的監(jiān)測(cè)結(jié)果需要較多的SAR影像數(shù)據(jù)，一般不少于25景[10-11]，一定程度上限制了PS技術(shù)的發(fā)展。小基線集合成孔徑雷達(dá)干涉測(cè)量（small baseline subset inter?ferometric synthetic aperture radar，SBAS-InSAR）技術(shù)[12]將SAR影像數(shù)據(jù)組成若干個(gè)子集，采用最小二乘法求解子集的形變時(shí)間序列，同時(shí)利用奇異值分解法（singularvalue decomposition，SVD）將多個(gè)小基線集聯(lián)合求解，獲取整個(gè)時(shí)間跨度的形變序列[13]，相對(duì)于PS-InSAR技術(shù)，SBAS-InSAR需要的SAR影像數(shù)目較少且獲取非線性形變信息的能力較強(qiáng)[10,14]?，F(xiàn)有成果[15-19]已利用SBAS-InSAR技術(shù)對(duì)城市地表沉降進(jìn)行了監(jiān)測(cè)，均取得了理想成效。但是，由于Sentinel-1A影像數(shù)據(jù)的新穎性，基于SBAS-InSAR技術(shù)和Sen?tinel-1A影像的城市地表沉降監(jiān)測(cè)方面的研究尚顯不足。

本文在對(duì)SBAS-InSAR技術(shù)進(jìn)行分析的基礎(chǔ)上，以位于順義-朝陽(yáng)-通州沉降帶的三河市燕郊鎮(zhèn)為研究區(qū)，對(duì)覆蓋該研究區(qū)2017—2018年的20景Senti?nel-1A數(shù)據(jù)進(jìn)行處理，獲取該區(qū)域地表年平均形變速率及時(shí)序形變信息。同時(shí)，本文還監(jiān)測(cè)了PS-InSAR技術(shù)下的研究區(qū)地表沉降信息，并進(jìn)一步對(duì)比分析了2種監(jiān)測(cè)技術(shù)之間的差異性。最后，本文探討分析了引起研究區(qū)地表沉降的主要成因。

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

設(shè)有按時(shí)間序列t0，…，ti，…，tN獲取到N+1幅單視復(fù)數(shù)SAR影像，將它們以任意影像為主影像進(jìn)行配準(zhǔn)；設(shè)定垂直基線閾值，將垂直基線小于該閾值的SAR影像歸為一組，設(shè)共分L組，對(duì)分組后的影像做差分干涉處理，可得到M幅差分干涉圖，假設(shè)N為奇數(shù)，則差分干涉圖的個(gè)數(shù)M可以表示為：

以t0為初始時(shí)刻，任意時(shí)刻ti(i=1，…，N)相對(duì)于t0時(shí)刻的差分相位φ(ti)為未知參數(shù)，差分干涉相位δφ(tk)(k=1，…，M)為觀測(cè)量。對(duì)于第k(k=1，…，i，…，M)幅差分干涉圖的任意像元(x，r)有

式中，λ為雷達(dá)波長(zhǎng)；d(tA，x，r)和d(tB，x，r)分別為tA和tB時(shí)刻像元(x，r)相對(duì)于初始時(shí)刻t0的雷達(dá)LOS像地表形變，即有d(t0，x，r)=0。為獲取研究區(qū)的時(shí)間形變序列，需要精確預(yù)估出地形相位誤差分量、大氣延遲相位誤差分量以及噪聲相位分量，并將這3個(gè)分量從干涉相位δφk(x，r)中去除。

根據(jù)式（2），在去除各項(xiàng)誤差分量后，M幅差分干涉圖可以得到M個(gè)方程，用矩陣表示為：

式中，A為M×N矩陣，由-1、0、1組成；φΤ=[φ(t1)，…，φ(tN)]為每一景SAR影像中高相干點(diǎn)對(duì)應(yīng)的相位值所組成的向量；δφΤ=[δφ1，…，δφM]為各差分干涉圖對(duì)應(yīng)的解纏相位值所組成的向量。

為求解研究區(qū)域各高相干點(diǎn)的形變速率，可用兩景SAR影像間的平均相位速率代替相位值，則式（3）變?yōu)椋?/p>

式中，B為M×N的系數(shù)矩陣；vT可以表示為：

當(dāng)系數(shù)矩陣B為滿秩（即M≥N）時(shí)，可用最小二乘法則求解出形變速率；當(dāng)M＜N時(shí)，矩陣B出現(xiàn)秩虧，可利用奇異值分解法（SVD）獲取研究區(qū)域的形變速率[16]。獲取形變速率后，依據(jù)研究區(qū)SAR影像時(shí)間區(qū)間求得相應(yīng)時(shí)間段內(nèi)的形變量[15]。

2 研究區(qū)概況與數(shù)據(jù)來(lái)源

2.1 研究區(qū)概況

燕郊鎮(zhèn)位于河北省三河市，如圖1所示，西臨北京市通州區(qū)，東臨天津市，中心地理坐標(biāo)為39°56'N、116°48'E，覆蓋面積約324 km2。區(qū)內(nèi)地勢(shì)平坦，位于潮白河畔，東南部為緩向渤海傾斜的平原。氣候?yàn)榕瘻貛О霛駶?rùn)半干旱季風(fēng)氣候，夏季高溫多雨，冬季寒冷干燥，春、秋短促?，F(xiàn)有資料顯示，研究區(qū)屬于京津冀地區(qū)沉降較為嚴(yán)重的區(qū)域，且與北京五環(huán)外的順義-朝陽(yáng)-通州沉降帶相連接。目前，由于研究區(qū)人口密度大，導(dǎo)致本研究區(qū)相對(duì)于其他區(qū)域的沉降較為顯著[20]。

圖1 研究區(qū)概況圖

2.2 數(shù)據(jù)來(lái)源

本文選取時(shí)間跨度為2017-03-21~2018-10-06的20景覆蓋研究區(qū)域的Sentinel-1A1級(jí)影像數(shù)據(jù)作為研究數(shù)據(jù)，采用波長(zhǎng)為5.6 cm的C波段，軌道方向?yàn)樯?，方位向與距離向分辨率分別為13.985 m與9.318 m，影像中心入射角約為43.8°，極化方式為VV極化。采用美國(guó)宇航局提供的30 m分辨率SRTM1 DEM數(shù)據(jù)去除地形相位影響，同時(shí)利用歐洲空間局發(fā)布的POD精密軌道數(shù)據(jù)去除軌道誤差。

3 數(shù)據(jù)處理與結(jié)果分析

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

本次研究采用SAR-scape系統(tǒng)進(jìn)行數(shù)據(jù)處理，具體數(shù)據(jù)處理過(guò)程如下：

1）數(shù)據(jù)預(yù)處理：將Sentinel-1A原始影像裁剪至本次研究區(qū)域大小，如圖1所示，對(duì)裁剪后的影像進(jìn)行5:1（方位向×距離向）多視處理。

2）干涉處理：為避免空間失相干，設(shè)定時(shí)間基線閾值為180 d、空間基線閾值為150 m進(jìn)行差分干涉處理，由系統(tǒng)自動(dòng)選出日期為2017-11-16獲取的影像為公共主影像，這樣可以保證有足夠多的像對(duì)，在對(duì)像對(duì)編輯時(shí)可以將相干性低的像對(duì)進(jìn)行去除。將其余所有影像進(jìn)行配準(zhǔn)并重采樣至公共主影像，并選用3D解纏，共生成92對(duì)小基線差分干涉像對(duì)，像對(duì)組合方式如圖2所示，圖2中黃顏色點(diǎn)代表公共主影像，線段長(zhǎng)度代表干涉基線長(zhǎng)度，統(tǒng)計(jì)組合后干涉相對(duì)最大時(shí)間基線為180 d，最大空間基線為110 m。

圖2 時(shí)空基線組合

3）軌道精煉及重去平：對(duì)恒定相位進(jìn)行估算去除，在此步驟進(jìn)行之前需要對(duì)上一步生成的干涉結(jié)果進(jìn)行質(zhì)量檢查，若在干涉圖、解纏結(jié)果圖以及相干性圖中出現(xiàn)質(zhì)量較差（失相干現(xiàn)象嚴(yán)重）的結(jié)果，則考慮將該像對(duì)做刪除處理。利用GCP控制點(diǎn)對(duì)所有像對(duì)進(jìn)行重去平，本次研究共選取了30個(gè)GCP點(diǎn)。

4）SBAS第一次反演：設(shè)置相干系數(shù)閾值為0.35，估算形變速率和殘余地形，同時(shí)進(jìn)行二次解纏對(duì)干涉圖進(jìn)行優(yōu)化。

5）SBAS第二次反演：分別設(shè)置大氣高通濾波、大氣低通濾波2個(gè)參數(shù)進(jìn)行大氣濾波，估算并去除大氣相位，計(jì)算時(shí)間序列位移信息。

6）地理編碼：依據(jù)振幅與相位的穩(wěn)定性篩選研究區(qū)內(nèi)的高相干點(diǎn)，通過(guò)奇異值分解（SVD）法求解高相干目標(biāo)點(diǎn)的沉降速率，將第二次反演結(jié)果從SAR坐標(biāo)系轉(zhuǎn)換到地理坐標(biāo)系，最后得到研究區(qū)的LOS向形變速率及時(shí)序形變信息。

3.2 結(jié)果分析

1）研究區(qū)形變速率。由SBAS-InSAR技術(shù)獲取到研究區(qū)2017-03-21~2018-10-06的地表年平均沉降速率信息，如圖3所示，研究區(qū)下沉速率較大區(qū)域主要集中在西部與中部（圖中所示多邊形A區(qū)和B區(qū)）。整個(gè)研究區(qū)的年平均沉降速率范圍達(dá)到-62~26 mm/a，不均勻沉降較為明顯，平均沉降速率超過(guò)20 mm/a，同時(shí)存在的沉降漏斗平均沉降速率超過(guò)36 mm/a。

圖3 SBAS技術(shù)研究區(qū)LOS向平均沉降速率

2）研究區(qū)時(shí)序形變。以2017-03-21為起始時(shí)間（參考形變量為零）獲得研究區(qū)的時(shí)序形變信息，如圖4所示，至2017-09-17研究區(qū)中部（圖中深藍(lán)色區(qū)域）已經(jīng)發(fā)生了較大沉降，最大沉降達(dá)到34 mm；至2018-03-28研究區(qū)西部發(fā)生較大沉降，最大沉降達(dá)到45 mm，研究區(qū)沉降范圍基本穩(wěn)定；至2018-10-06，研究區(qū)中部與西部最大沉降均達(dá)到74 mm。

圖4 SBAS技術(shù)研究區(qū)LOS向時(shí)序沉降變化

為驗(yàn)證SBAS-InSAR監(jiān)測(cè)結(jié)果的準(zhǔn)確性，同時(shí)采用PS-InSAR技術(shù)進(jìn)行地表形變監(jiān)測(cè)，監(jiān)測(cè)結(jié)果如圖5所示。由圖5可知，沉降同樣主要發(fā)生在研究區(qū)的中部與西部（圖中所示多邊形A區(qū)和B區(qū)），最大沉降速率達(dá)到24 mm/a。

圖5 PS技術(shù)研究區(qū)LOS向平均沉降速率

為分析2種技術(shù)在地表形變監(jiān)測(cè)結(jié)果的差異性，在研究區(qū)的集中沉降區(qū)域提取同名特征點(diǎn)，其中A5特征點(diǎn)位于橋面，其余所有特征點(diǎn)均位于建筑物表面，點(diǎn)位分布如圖6所示。通過(guò)放大圖6可以發(fā)現(xiàn)，SBAS技術(shù)在監(jiān)測(cè)地表形變細(xì)節(jié)上更具優(yōu)勢(shì)，在A2、A3、B1點(diǎn)所在的特征點(diǎn)區(qū)域均形成明顯的下沉盆地，且下沉范圍連續(xù)集中；而在PS監(jiān)測(cè)結(jié)果中，A2、A3、B1點(diǎn)所處區(qū)域沒(méi)有監(jiān)測(cè)到下沉盆地形態(tài)。此外，整體而言，相較于SBAS技術(shù)的監(jiān)測(cè)結(jié)果，通過(guò)PS技術(shù)所監(jiān)測(cè)到的高相干點(diǎn)出現(xiàn)了大面積缺失。

圖6 同名特征點(diǎn)分布圖

為定量分析SBAS與PS的監(jiān)測(cè)結(jié)果，分別從研究區(qū)的A區(qū)和B區(qū)提取特征點(diǎn)的時(shí)序形變信息，提取結(jié)果如圖7所示。從圖7可以看出，SBAS獲取的沉降值整體大于PS獲取的沉降值。盡管兩組曲線的走向是一致的，但是利用SBAS技術(shù)監(jiān)測(cè)的特征點(diǎn)時(shí)序曲線較為平滑，而PS技術(shù)監(jiān)測(cè)的特征點(diǎn)時(shí)序曲線出現(xiàn)較大幅度波動(dòng)，且在A3、A6點(diǎn)無(wú)法獲取有效監(jiān)測(cè)值。

圖7 SBAS與PS技術(shù)結(jié)果特征點(diǎn)時(shí)序累計(jì)沉降量

通過(guò)計(jì)算各特征點(diǎn)至2018-10-06的累計(jì)沉降值可知（表1），位于橋面的A5特征點(diǎn)累計(jì)沉降值為18 mm，兩者沉降值差為零；A2與B1特征點(diǎn)沉降值差較大，原因是這2個(gè)特征點(diǎn)在SBAS結(jié)果中均出現(xiàn)在下沉較為集中區(qū)域，而在PS結(jié)果中此區(qū)域無(wú)高相干點(diǎn)被探測(cè)，且其周邊也出現(xiàn)大面積缺失值；對(duì)其余各特征點(diǎn)沉降值差求取均方根誤差為5.6 mm，說(shuō)明SBAS-InSAR方法在監(jiān)測(cè)地表形變是可靠的。

表1 各同名特征點(diǎn)最大沉降值差/mm

通過(guò)研究結(jié)果的交叉對(duì)比可知，SBAS技術(shù)和PS技術(shù)獲取的沉降結(jié)果在監(jiān)測(cè)細(xì)節(jié)上存在輕微差異，但是總體來(lái)講，2種方法獲取的形變區(qū)域具有較高的一致性，這也證明了該研究區(qū)域在2017年3月至2018年10月期間存在地面沉降問(wèn)題，且平均沉降速率均達(dá)到20 mm/a。導(dǎo)致2種方法監(jiān)測(cè)下沉信息有所差異的原因是由于SBAS是非線性沉降模型，PS是線性沉降模型，兩者所描述的沉降過(guò)程不完全一樣，存在局部一致或局部不一致的現(xiàn)象。本次研究使用20景Senti?nel-1A數(shù)據(jù)，數(shù)據(jù)量相對(duì)較少，在同等數(shù)據(jù)量的情況下，SBAS可以形成高密度的干涉像對(duì)組合，而PS方法過(guò)于依賴于單一公共主影像的質(zhì)量，且在識(shí)別高相干點(diǎn)時(shí)對(duì)影像數(shù)據(jù)量有一定要求，因此在數(shù)據(jù)量有限時(shí)，SBAS-InSAR方法監(jiān)測(cè)地表形變可以獲得更為可靠的結(jié)果[21-22]。

3.3 形變成因分析

從SBAS和PS技術(shù)提取的結(jié)果可以看出，沉降區(qū)域位于研究區(qū)中部與西部，造成沉降的原因有：①該區(qū)域主要分布于潮白河、溫榆河和泃河流域的沖積、洪積扇平原上，多條河流橫跨研究區(qū)，屬于沖積扇的中下部，土質(zhì)主要由細(xì)沙黏土構(gòu)成[22]，同時(shí)地下水過(guò)度開(kāi)采會(huì)導(dǎo)致地面沉降；②近幾年北京地區(qū)的城市化向東擴(kuò)張依舊呈現(xiàn)逐年增加趨勢(shì)，本研究區(qū)與北京市沉降區(qū)逐漸連成一片[20]；隨著一些重工業(yè)的東移，地表構(gòu)筑物的建設(shè)以及地下工程建設(shè)打破了原有地下巖土體應(yīng)力平衡，且這些建設(shè)工程在空間分布不均勻，造成的區(qū)域性集中沉降分布不均勻，形變量級(jí)存在明顯差異，造成該區(qū)域出現(xiàn)較為顯著的地表沉降現(xiàn)象，形成明顯的下沉盆地。

4 結(jié)語(yǔ)

本文利用SBAS-InSAR技術(shù)對(duì)20景Sentinel-1A數(shù)據(jù)進(jìn)行處理，得到了研究區(qū)2017~2018年的沉降速率及時(shí)序形變信息，發(fā)現(xiàn)下沉速率較大區(qū)域主要集中在研究區(qū)中部與西部，平均沉降速率超過(guò)20 mm/a，部分集中沉降區(qū)域平均沉降速率超過(guò)36 mm/a，累計(jì)沉降達(dá)到74 mm。同時(shí)利用PS-InSAR做結(jié)果交叉驗(yàn)證，發(fā)現(xiàn)SBAS方法和PS方法獲取的形變區(qū)域具有較高的一致性，但SBAS方法在監(jiān)測(cè)地表細(xì)節(jié)更占優(yōu)勢(shì)，在數(shù)據(jù)量有限時(shí)，SBAS-InSAR監(jiān)測(cè)地表形變可以獲得更為可靠的結(jié)果。研究區(qū)沉降主要受到其地質(zhì)構(gòu)造、地下水開(kāi)采及北京地區(qū)的城市化向東擴(kuò)張有關(guān)。