王 磊，連增增,劉 杰

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2.1 研究區(qū)概況

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

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

2.2 數(shù)據(jù)來源

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

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

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

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

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

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

圖2 時空基線組合

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

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

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

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

3.2 結(jié)果分析

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

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

2）研究區(qū)時序形變。以2017-03-21為起始時間（參考形變量為零）獲得研究區(qū)的時序形變信息，如圖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向時序沉降變化

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

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

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

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

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

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

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

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

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

3.3 形變成因分析

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

4 結(jié)語

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