基于SBAS
--InSAR的遼寧義縣中東部區(qū)域地表沉降監(jiān)測(cè)

黃龍霄，張旭晴，趙強(qiáng)，高明久，程微，安繼魁，吳迪

1.吉林大學(xué) 地球探測(cè)科學(xué)與技術(shù)學(xué)院,長(zhǎng)春 130026；2.遼寧工程勘察設(shè)計(jì)院，遼寧 錦州 121000

0 引言

地面沉降又稱地面下陷，是由于地下松散地層固結(jié)壓縮，導(dǎo)致地殼表面標(biāo)高降低的一種局部的工程地質(zhì)現(xiàn)象[1]。合成孔徑雷達(dá)干涉測(cè)量技術(shù)(Interferometric Synthetic Aperture Radar，InSAR)是近20年發(fā)展起來的一種先進(jìn)的遙感技術(shù)，補(bǔ)充了已有的地面沉降監(jiān)測(cè)方法:精密水準(zhǔn)測(cè)量和GPS監(jiān)測(cè)[2]。Gabriel等人將D--InSAR技術(shù)應(yīng)用在地面沉降監(jiān)測(cè)領(lǐng)域[3]，D--InSAR技術(shù)雖然可以獲取密度較大的形變值[4]，但它著重研究短時(shí)間間隔內(nèi)的單一形變情況。為了彌補(bǔ)D--InSAR技術(shù)的不足，學(xué)者Ferretti提出PS--InSAR技術(shù)[5]。PS--InSAR 技術(shù)只關(guān)注穩(wěn)定性較高點(diǎn)，可以進(jìn)行長(zhǎng)時(shí)間序列的地面沉降監(jiān)測(cè)，但在地表情況變化較大的一些非城市區(qū)域使用PS--InSAR技術(shù)時(shí)很難獲取高密度的沉降監(jiān)測(cè)結(jié)果。因此，Berardino P et al.提出 SBAS--InSAR 技術(shù)[6]，該技術(shù)利用時(shí)空基線較短的SAR數(shù)據(jù)集形成干涉對(duì)[7]，充分利用短時(shí)空基線影像的相干信息，有效地抑制相位噪聲對(duì)地形相位的影響，從而獲取長(zhǎng)時(shí)間緩慢地表形變的演變規(guī)律[8]。SBAS技術(shù)不僅可以減輕時(shí)空失相關(guān)、大氣延遲的影響，還可以應(yīng)用有限數(shù)量的影像得到毫米級(jí)的時(shí)序沉降量[9]，空間分辨率的降低也相應(yīng)的減少了數(shù)據(jù)處理的運(yùn)算量。

由于遼寧義縣中東部地區(qū)地形較復(fù)雜，為了得到好的監(jiān)測(cè)結(jié)果，綜上考慮，采用SBAS--InSAR技術(shù)對(duì)該區(qū)域進(jìn)行地面沉降監(jiān)測(cè)研究。

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

SBAS技術(shù)的主要思想為通過設(shè)置時(shí)間和空間基線閾值得到L個(gè)小基線集合，每個(gè)集合之內(nèi)，干涉對(duì)基線較小，而集合之間，干涉對(duì)基數(shù)較大，共包括M幅差分干涉圖，篩選出基線較短的干涉對(duì)組成不同的集合。集合與集合之間的基線較長(zhǎng)，可通過使用最小二乘法得到每個(gè)集合的形變序列，再通過奇異值分解法將每個(gè)集合聯(lián)合求解，這種方法可以在有限的數(shù)據(jù)基礎(chǔ)上得到高度精度反演，同時(shí)不要求主影像必須相同[10--13](圖1)。

圖1 技術(shù)流程Fig.1 Technical process

選取同區(qū)域的N+1幅SAR影像，獲取時(shí)間依次為t0,t1,t2,…,tn,根據(jù)設(shè)置的干涉條件可以組合得到M幅干涉圖，其中M滿足下式：

(1)

現(xiàn)選取tA和tB兩時(shí)刻獲得的SAR影像，然后進(jìn)行差分干涉處理，形成第i個(gè)差分干涉對(duì)，則差分干涉相位δφi可由下式表示：

(2)

式中：d(tA)和d(tB)表示tA和tB時(shí)刻相對(duì)于參考時(shí)刻t0(d(t0)≡0)的累積形變信息，φ(tA)和φ(tB)表示相應(yīng)相位值。

現(xiàn)用向量φT=[φ(t1),…,φ(tN)]表示N個(gè)相位圖，向量δφT=[δφ1,…,δφM]x表示M個(gè)干涉相位圖，如果主圖像先于輔圖像的獲得時(shí)間，則有：

δφj=φtSj-φtMmj,j=1,…,M

(3)

式中：Mm、S分別表示干涉像對(duì)中主、輔圖像。所有參與生成干涉圖的SAR影像可以在不同的短基線集中。現(xiàn)定義一個(gè)含有M個(gè)方程，N個(gè)未知參數(shù)的方程組，其表達(dá)式為：

δφ=Aφ

(4)

式中：矩陣A[M×N]中每一行對(duì)應(yīng)一個(gè)干涉對(duì)，每一列對(duì)應(yīng)一副SAR影像，A[j,Sj]=1和A[j,Mmj]=-1，其余元素為0。矩陣A是一個(gè)由干涉圖組合方式?jīng)Q定的近似關(guān)聯(lián)矩陣。當(dāng)所有干涉對(duì)屬于同一個(gè)基線集時(shí)，A是一個(gè)列滿秩矩陣。當(dāng)M=N時(shí)，方程組(8)有唯一解；當(dāng)M>N時(shí)，方程組是超定方程，此時(shí)可使用最小二乘法求此方程組的唯一解：

(5)

當(dāng)所有干涉對(duì)屬于不同基線集時(shí)，矩陣A的秩等于N-L+1，出現(xiàn)秩虧，可采用奇異值分解(SVD)的方法計(jì)算矩陣A的廣義逆矩陣A+，進(jìn)而求解方程組(4)的最小范數(shù)解，其過程如下：

A+=US+VT

(6)

式中：正交矩陣U[M×N]，其前N列是ATA特征向量,稱為矩陣A的左奇異矩陣，正交矩陣V[M×N]，其所有列ATA稱為矩陣A的右奇異矩陣。S為M×N矩陣。則：

A+VS+UT

(7)

(8)

式中：ui和vi分別表示正交矩陣U和V的列向量。

將方程組(4)中對(duì)相位的求解轉(zhuǎn)化為對(duì)平均相位變化速率的求解問題，則待求參數(shù)向量為：

(9)

則有：

δφ=Bv

(10)

式中：B[M×N]，其中,B[i,j]=tj+1-tj,(Sj+1≤j≤Mmi,i=1,…,M)其他元素值為0。對(duì)B進(jìn)行奇異值分解，即可解出各時(shí)間段平均速度v。

2 地面沉降監(jiān)測(cè)

2.1 研究區(qū)概括

遼寧義縣東部為剝蝕構(gòu)造中低山區(qū)，由變質(zhì)巖和花崗巖侵入體組成，山勢(shì)較陡峭，海拔在400～800 m之間，山脈走向由東北向西南延展，山體多為直坡，局部為凸坡，坡度角為30°～50°，山間谷地，多呈樹枝狀分布。義縣的中部為沖積平原區(qū)，分布于大凌河、細(xì)河?xùn)|岸；上部為亞砂土覆蓋；下部為砂礫石層；底部混土，地面較平，微傾向河床。中部區(qū)域部分巖層性質(zhì)為白堊系：灰白色砂礫巖夾砂巖和頁(yè)巖及上部夾薄層煤，中部賦存煤層(圖2)。

圖2 研究區(qū)域Fig.2 Research area

2.2 監(jiān)測(cè)數(shù)據(jù)

Sentinel--1數(shù)據(jù)具有較好的空間基線控制，并且重返周期較短，適合地表形變監(jiān)測(cè)[13]。由于遼寧地處中國(guó)東北部，冬季多雪，而C波段影像受降雪影響會(huì)出現(xiàn)嚴(yán)重的失相干現(xiàn)象，導(dǎo)致成果獲取的相干點(diǎn)數(shù)量降低，誤差增大，所以在本次監(jiān)測(cè)中未采用冬季影像。最終選取由歐空局提供的IW寬條干涉的降軌SAR影像，C波段，分辨率為5 m×20 m，極化方式為VV，2017年—2018年2月到11月的18景影像數(shù)據(jù)，具體時(shí)間見表1。同時(shí)使用Sentinel--1衛(wèi)星精密軌道數(shù)據(jù)(POD精密定軌星歷數(shù)據(jù))和由地理空間數(shù)據(jù)云將ASTER GDEM(V1)數(shù)據(jù)進(jìn)行加工得來的GDEMDEM 30M 分辨率數(shù)字高程數(shù)據(jù)作為輔助數(shù)據(jù)。

表1 Sentiel--1數(shù)據(jù)

2.3 處理流程

利用SAR Scape 軟件，通過設(shè)置時(shí)間基線、空間基閾值控制生成干涉對(duì)的數(shù)量。本次監(jiān)測(cè)設(shè)置270 d的時(shí)間基線和最大45%的臨界空間基線閾值，生成88對(duì)像對(duì)，生成的時(shí)空連接如圖3所示。

圖3 時(shí)空連接圖Fig.3 Space-time connection map

然后基于干涉對(duì)進(jìn)行SLC影像配準(zhǔn)，干涉圖的生成、去平，經(jīng)過自適應(yīng)濾波、相干性生成及相位解纏獲取一系列相位圖。解纏相干系數(shù)閾值為0.2，解纏方法為Delaunay MCF，濾波方法為Goldstein。然后再選擇其中1景干涉效果中等的干涉圖，根據(jù)相干性圖選取相干性高、相位好的GCP點(diǎn)，為軌道精煉和重去平做準(zhǔn)備。選取近100個(gè)控制點(diǎn)，通過多次篩選將各個(gè)點(diǎn)的誤差降到了1.5以下。然后經(jīng)過兩次SBAS反演，精確估計(jì)且去除地形殘余相位、大氣效應(yīng)相位，最后結(jié)合DEM數(shù)據(jù)進(jìn)行地理編碼后獲取WGS--84坐標(biāo)系下的各期累積形變圖和平均形變速率圖[1,10]。

2.4 監(jiān)測(cè)結(jié)果與分析

本次監(jiān)測(cè)使用SBAS--InSAR技術(shù)通過對(duì)Sentinel--1數(shù)據(jù)的處理，結(jié)合ArcGIS軟件，得到2017—2018年成果沉降速率渲染圖(圖4)、沉降速率統(tǒng)計(jì)圖(圖5)和2017—2018年成果沉降速率等值線圖(圖6)。

圖4 2017年—2018年成果沉降速率渲染圖Fig.4 Rending map of results settlement rate from 2017 to 2018

圖5 沉降速率統(tǒng)計(jì)圖Fig.5 Statistical map of settlement rate

由圖4可知，在研究區(qū)的左下角部分即義縣縣城周邊存在較明顯的沉降，義縣中心地區(qū)沉降最為集中，由圖5可知，2017 年 2 月 21日—2018年 10月26日期間，研究區(qū)范圍內(nèi)大部分PS點(diǎn)的沉降速率為5 mm/a±，平均沉降速率為3.5 mm/a，最大沉降速率為 98 mm/a，說明研究區(qū)地面沉降空間分布有較大差異，地面沉降不均衡。

使用克里金法進(jìn)行插值，繪制出該區(qū)域的形變量等值線圖(圖6)，疊加到谷歌地圖上，便于沉降區(qū)定位；發(fā)現(xiàn)在A、B、C 3個(gè)小區(qū)域存在沉降漏斗且沉降量較大，對(duì)這3個(gè)小區(qū)域進(jìn)行重點(diǎn)分析，從這3個(gè)區(qū)域中分別選擇一個(gè)能夠代表本區(qū)域形變特征的、形變量較大的 PS 點(diǎn)進(jìn)行時(shí)間序列分析。 時(shí)間序列分析是處理變形觀測(cè)數(shù)據(jù)的一種有效方法，通過時(shí)間序列分析可以發(fā)現(xiàn)各形變點(diǎn)的時(shí)間變化規(guī)律，可以對(duì)沉降區(qū)域進(jìn)行短期的預(yù)測(cè)與分析[14]。由圖7～9各曲線可得，沉降區(qū)A的沉降中心沉降量>-180 mm；沉降區(qū)B最大沉降量大約為259 mm；沉降區(qū)C有3個(gè)較明顯的沉降中心，沉降量>-275 mm。

圖6 2017年—2018年成果沉降速率等值線圖Fig.6 Contour map of results settlement rate from 2017 to 2018

圖7 A區(qū)PS點(diǎn)形變時(shí)間序列Fig.7 Deformation time series of PS point in region A

圖8 B區(qū)PS點(diǎn)形變時(shí)間序列Fig.8 Deformation time series of PS point in region B

圖9 C區(qū)PS點(diǎn)形變時(shí)間序列Fig.9 Deformation time series of PS point in region C

圖10 義縣地質(zhì)災(zāi)害詳細(xì)調(diào)查圖Fig.10 Detailed survey map of geological hazards in Yixian

為了檢驗(yàn)結(jié)果的可靠性，將圖10與監(jiān)測(cè)結(jié)果對(duì)比，發(fā)現(xiàn)義縣地面塌陷情況與所得結(jié)果基本吻合。地面沉降主要是由地殼構(gòu)造活動(dòng)、礦產(chǎn)資源開采、地下水過量采集與城市建筑荷載增加等因素造成的。為分析義縣地面沉降成因，通過查閱相關(guān)資料發(fā)現(xiàn)：研究區(qū)內(nèi)礦產(chǎn)資源較豐富，有煤礦21個(gè)、黏土礦10個(gè)、采石廠11個(gè)、鐵礦2個(gè)和硅石礦1個(gè)(表2)。義縣礦業(yè)開采強(qiáng)烈，尤其是對(duì)地質(zhì)環(huán)境影響較大且易發(fā)生災(zāi)害的煤礦、花崗巖礦的開采，其中對(duì)煤礦、黏土礦的開采會(huì)造成地面塌陷、地裂縫災(zāi)害。據(jù)統(tǒng)計(jì)，截至2004年義縣中部地區(qū)的聚糧屯煤礦和九道嶺煤礦采空區(qū)，使2 km2多耕地成為沉陷區(qū)。一些個(gè)體企業(yè)無規(guī)則開采對(duì)山體破壞嚴(yán)重，極容易發(fā)生崩塌、滑坡等地質(zhì)災(zāi)害。沉降區(qū)的分布與義縣的礦區(qū)分布較為一致，且處于采空塌陷易發(fā)區(qū)，表明研究區(qū)地面沉降主要是由人類工程經(jīng)濟(jì)活動(dòng)導(dǎo)致的，其中礦產(chǎn)資源開采是最大的影響因素。

表2 礦產(chǎn)資源匯總

3 討論

本次監(jiān)測(cè)改進(jìn)之處：①數(shù)據(jù)處理過程中為了能夠選出形變量更小、穩(wěn)定性更好、誤差更小的 GCP 點(diǎn)，通過控制每個(gè)GCP點(diǎn)的誤差值，使結(jié)果更加的精確，減弱了人為因素對(duì)軌道精煉誤差的影響。②監(jiān)測(cè)結(jié)果分析采用了等值線分析法、時(shí)間序列分析法和統(tǒng)計(jì)學(xué)法，從不同的角度研究了義縣地面沉降情況，成功地獲取研究區(qū)的平均沉降速率和形變序列并結(jié)合野外實(shí)地驗(yàn)證，確定了主要沉降區(qū)的位置、最大沉降量以及沉降趨勢(shì)。

4 結(jié)論

(1)監(jiān)測(cè)發(fā)現(xiàn)義縣地表沉降分布不均勻，在義縣中部地區(qū)存在多處較明顯的沉降區(qū)，如聚糧屯鄉(xiāng)和前楊鄉(xiāng)附近。在監(jiān)測(cè)時(shí)段內(nèi)研究區(qū)沉降量在-275～5 mm范圍內(nèi)。

(2)義縣的沉降區(qū)與礦產(chǎn)資源的分布有著密切的關(guān)系。在礦石資源開采中，土體開挖卸載時(shí)開挖面土體向開采面內(nèi)移動(dòng)、地下結(jié)構(gòu)整體下沉均可造成地表沉降。

(3)義縣地表沉降趨勢(shì)顯示下降，時(shí)間序列與季節(jié)變化存在一定的相關(guān)性，雨水充沛期和冬季降雪期，沉降速率都會(huì)出現(xiàn)減緩或與總趨勢(shì)相反的小幅度變化。