黃龍霄，張旭晴，趙強，高明久，程微，安繼魁，吳迪

1.吉林大學 地球探測科學與技術學院,長春 130026；2.遼寧工程勘察設計院，遼寧 錦州 121000

0 引言

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

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

1 SBAS--InSAR技術原理

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

圖1 技術流程Fig.1 Technical process

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

(1)

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

(2)

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

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

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

(3)

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

δφ=Aφ

(4)

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

(5)

當所有干涉對屬于不同基線集時，矩陣A的秩等于N-L+1，出現(xiàn)秩虧，可采用奇異值分解(SVD)的方法計算矩陣A的廣義逆矩陣A+，進而求解方程組(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)中對相位的求解轉化為對平均相位變化速率的求解問題，則待求參數(shù)向量為：

(9)

則有：

δφ=Bv

(10)

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

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

2.1 研究區(qū)概括

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

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

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

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

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

2.3 處理流程

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

圖3 時空連接圖Fig.3 Space-time connection map

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

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

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

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

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

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

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

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

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

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

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

圖10 義縣地質災害詳細調查圖Fig.10 Detailed survey map of geological hazards in Yixian

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

表2 礦產資源匯總

3 討論

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

4 結論

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

(2)義縣的沉降區(qū)與礦產資源的分布有著密切的關系。在礦石資源開采中，土體開挖卸載時開挖面土體向開采面內移動、地下結構整體下沉均可造成地表沉降。

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