馬 露， 張 帆

(中國煤炭地質(zhì)總局航測(cè)遙感局，西安 710199)

地面沉降對(duì)人類經(jīng)濟(jì)和社會(huì)活動(dòng)造成了嚴(yán)重影響，不僅制約著人類生存區(qū)域的生產(chǎn)經(jīng)營活動(dòng)，也是造成安全事故的重大隱患[1]。合成孔徑雷達(dá)干涉測(cè)量技術(shù)(Interferometric Synthetic Aperture Radar，InSAR)能夠有效地監(jiān)測(cè)地面沉降，與傳統(tǒng)的精密水準(zhǔn)儀監(jiān)測(cè)和GPS監(jiān)測(cè)相結(jié)合能夠極大地提高地面沉降監(jiān)測(cè)的效率，在生產(chǎn)生活中被廣泛應(yīng)用[2-3]。Gabriel等人提出的D-InSAR技術(shù)能夠獲取較短時(shí)間內(nèi)的變形情況及精度較高的形變值[4]，但該項(xiàng)技術(shù)在長時(shí)間連續(xù)監(jiān)測(cè)方面存在不足。為了彌補(bǔ)D-InSAR技術(shù)的不足，F(xiàn)erretti提出了PS-InSAR技術(shù)[5]。PS-InSAR技術(shù)能夠開展長時(shí)間序列的地面沉降監(jiān)測(cè)，在地勢(shì)平坦的城區(qū)取得了較好的應(yīng)用效果，但是在地表起伏較大的自然地貌區(qū)監(jiān)測(cè)結(jié)果很不理想。而Berardino和Lanari等人提出的SBAS-InSAR技術(shù)主要采用短時(shí)空基線的SAR數(shù)據(jù)集形成干涉像對(duì)[6-7]，充分利用像對(duì)間的相干信息，降低相位噪聲對(duì)地形相位的影響，有效獲取地表形變的時(shí)間序列圖以及地表緩慢形變的長時(shí)間變化規(guī)律[8]。SBAS-InSAR技術(shù)實(shí)現(xiàn)了利用有限數(shù)量影像獲取地表毫米級(jí)的沉降量，數(shù)據(jù)處理的運(yùn)算量也隨空間分辨率的降低而減少[9]。

奎達(dá)市是巴基斯坦俾路支省省會(huì)城市，奎達(dá)市存在較嚴(yán)重的地面沉降，2006—2009年平均沉降速率為100mm/a[10]?？_(dá)市以東地區(qū)為山地區(qū)，地形復(fù)雜，構(gòu)造發(fā)育，地質(zhì)災(zāi)害多發(fā)，本次選擇奎達(dá)市及以東地區(qū)采用SBAS-InSAR技術(shù)進(jìn)行地面沉降監(jiān)測(cè)，旨在初步掌握奎達(dá)市及以東地區(qū)沉降區(qū)的位置、最大沉降量和沉降趨勢(shì)，結(jié)合研究區(qū)自然地理概況、地質(zhì)背景等分析沉降發(fā)生的主要原因。

1 短基線集(SBAS-InSAR)技術(shù)原理

SBAS-InSAR技術(shù)廣泛應(yīng)用于長時(shí)間序列地表形變監(jiān)測(cè)領(lǐng)域[11]，基本思路是通過設(shè)定合適的時(shí)間和空間基線閾值得到若干個(gè)小基線集合[12]，應(yīng)用最小二乘法獲得各集合的地表形變序列，采用奇異值分解法(SVD)進(jìn)行求解[13]，該方法不僅能有效避免空間失相干，還能夠提高時(shí)間采樣率[14-18]，主要技術(shù)流程如圖1所示。

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

2 研究區(qū)概況及數(shù)據(jù)處理

2.1 研究區(qū)概況

奎達(dá)市城區(qū)位于研究區(qū)西部，屬?zèng)_積平原區(qū)，地勢(shì)較平坦，海拔在1 600～1 750m，整體呈東南高，西北低。地處俾路支高原，屬亞熱帶沙漠氣候，氣候干旱，耕地靠坎兒井灌溉。該區(qū)靠近阿富汗邊境，是巴基斯坦和阿富汗的貿(mào)易中心，人口密度大。

奎達(dá)市以東地區(qū)為褶皺侵蝕中低山區(qū)，海拔在670～3 300m，山勢(shì)較緩，山脈走向西北—東南、東北—西南以及近南北向均有，以西北—東南為主，山體以直坡為主，局部存在凸坡，坡角多在8°～49°，山間谷地較發(fā)育，多沿河流呈樹枝狀分布(圖2)。

圖2 研究區(qū)位置Figure 2 Location of study area

2.2 數(shù)據(jù)選取

選取覆蓋研究區(qū)的28景Sentinel-1A降軌影像數(shù)據(jù)，數(shù)據(jù)類型為IW_SLC即干涉寬幅模式(TOPS Mode)的斜距單視復(fù)數(shù)數(shù)據(jù)。影像數(shù)據(jù)的成像時(shí)間在2018年1月至2019年12月，空間分辨率為5m×20m，C波段，波長5.6cm，極化方式為VV極化(表1)。此外還用到日本JAXA公司提供的ALOS數(shù)字表面模型“ALOS World 3D-30 m”以及各影像對(duì)應(yīng)的POD精密定軌星歷數(shù)據(jù)(POD Precise Orbit Ephemerides)。

表1 Sentinel-1A數(shù)據(jù)時(shí)相

2.3 數(shù)據(jù)處理流程

利用ENVI5.3軟件中的SARScape模塊，通過設(shè)置時(shí)間和空間基線的閾值，控制生成干涉像對(duì)的數(shù)量[19]。本次研究設(shè)定監(jiān)測(cè)時(shí)間基線閾值為120天，空間基線閾值為最大臨界基線的15%。通過干涉組合共生成115對(duì)干涉像對(duì)(圖3)。

圖3 時(shí)空基線連接圖Figure 3 Space-time baseline connection diagram

通過對(duì)SLC影像配準(zhǔn)、干涉圖生成和去平、自適應(yīng)濾波、相干圖生成和相位解纏獲取一系列解纏圖。解纏相干系數(shù)閾值設(shè)置為0.15，解纏方法采用Minimum Cost Flow，濾波方法為Goldstein。選擇干涉效果較好的1景干涉圖，根據(jù)相干圖選取相干性高的GCP點(diǎn)，為下一步軌道精煉和重去平做準(zhǔn)備。為了降低GCP點(diǎn)的誤差，在谷歌地球上選擇后向散射性高、相干性好的人工建筑，如道路、房屋等，共選取近50個(gè)控制點(diǎn)。經(jīng)過兩次SBAS反演，精確估計(jì)并去除地形殘余相位和大氣效應(yīng)相位，結(jié)合研究區(qū)DEM數(shù)據(jù)進(jìn)行地理編碼后獲取各時(shí)期累積形變圖和平均形變速率圖。

本次監(jiān)測(cè)通過采用SBAS-InSAR技術(shù)對(duì)Sentinel-1數(shù)據(jù)進(jìn)行處理，利用ArcGIS軟件作為輔助，得到了2018年1月至2019年12月兩年內(nèi)的年平均沉降速率圖和累積形變量圖(圖4、圖5)。

圖4 研究區(qū)年平均形變速率Figure 4 Annual average deformation rate in the study area

圖5 研究區(qū)累積形變量Figure 5 Cumulative deformation in the study area

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

在奎達(dá)市城區(qū)存在非常明顯的沉降，尤以城市中心地區(qū)沉降最為集中，從2018年1月16日至2019年12月30日，城區(qū)以外區(qū)域沉降速率大多為3 mm/a ±，平均沉降速率為1.6 mm/a，城市中心沉降速率最大可達(dá)197 mm/a，表明研究區(qū)地面沉降分布很不均衡，空間差異較大。

在奎達(dá)市城區(qū)有A、B、C3個(gè)區(qū)域存在沉降漏斗且沉降量較大，其中A區(qū)域面積最大，是奎達(dá)市主城區(qū)所在地，B、C區(qū)域范圍較小，地表大多為農(nóng)田。重點(diǎn)對(duì)這3個(gè)區(qū)域進(jìn)行分析，從這3個(gè)區(qū)域中分別選擇能夠代表本區(qū)域形變特征的、形變量較大的點(diǎn)進(jìn)行時(shí)間序列分析，A區(qū)域選擇了4個(gè)點(diǎn)，B、C區(qū)域各選擇1個(gè)點(diǎn)。由各形變點(diǎn)的時(shí)間序列分析可以看出：沉降區(qū)A有4個(gè)明顯的沉降中心，A1和A4沉降量>315 mm，A2和A3沉降量>350 mm；沉降區(qū)B最大沉降量約為290mm；沉降區(qū)C的沉降中心沉降量>210 mm(圖6至圖8)。

圖6 A區(qū)域形變時(shí)間序列Figure 6 Deformation time series of A area

圖7 B區(qū)域形變時(shí)間序列Figure 7 Deformation time series of B area

圖8 C區(qū)域形變時(shí)間序列Figure 8 Deformation time series of C area

為了進(jìn)一步研究奎達(dá)市及以東區(qū)域地面沉降成因，詳細(xì)查閱了奎達(dá)市相關(guān)資料，發(fā)現(xiàn)在沉降速率較大的區(qū)域存在明顯的人類活動(dòng)和城市化進(jìn)程。奎達(dá)市人口眾多，氣候干旱，周圍無水源地，只能通過抽取地下水來滿足生活用水需求，耕地大多也是依靠坎兒井來灌溉。因此，可以推測(cè)奎達(dá)市城區(qū)地面沉降是由于城市化的快速發(fā)展，過量開采地下水導(dǎo)致的。

巴基斯坦奎達(dá)市在2006—2009年間，通過GPS測(cè)量沉降速率為100mm/a[10]，而本次監(jiān)測(cè)結(jié)果顯示奎達(dá)市的沉降速率已達(dá)到近200mm/a。

研究區(qū)礦產(chǎn)資源豐富，根據(jù)巴基斯坦1∶100萬地質(zhì)礦產(chǎn)圖及相關(guān)資料，研究區(qū)主要礦產(chǎn)資源有煤礦、鉻鐵礦和銅礦等?？_(dá)市區(qū)以東的山地區(qū)，沉降多分布于研究區(qū)的中部和東南部區(qū)域，累計(jì)形變量大多為60～150mm。通過與高分辨率遙感影像對(duì)比發(fā)現(xiàn)，沉降區(qū)域大多與礦產(chǎn)資源開發(fā)密切相關(guān)，在遙感影像上可見多處采礦集中區(qū)，開采強(qiáng)度較大，個(gè)體企業(yè)的不規(guī)范開采及落后的采礦技術(shù)更加劇了地面塌陷的發(fā)生(圖9)。

圖9 山地區(qū)地面沉降與采礦活動(dòng)的關(guān)系Figure 9 Relationship between land subsidence and mining activities in the mountain area

4 結(jié)語

本文采用短基線SBAS-InSAR技術(shù)對(duì)覆蓋巴基斯坦奎達(dá)市及以東區(qū)域的28景Sentinel-1A降軌數(shù)據(jù)進(jìn)行干涉處理，獲取了研究區(qū)2018年1月至2019年12月地面形變速率與累積沉降量。研究結(jié)果表明：

1)奎達(dá)市城區(qū)存在明顯的地面沉降，分布有3處沉降漏斗，在監(jiān)測(cè)時(shí)段內(nèi)最大沉降量210～390mm，地下水的過量開采是導(dǎo)致地面沉降的主要原因。

2)奎達(dá)市以東的山地區(qū)地面沉降主要發(fā)生在研究區(qū)的東部和中南部，在監(jiān)測(cè)時(shí)段內(nèi)沉降量60～150mm，與煤礦、鉻鐵礦等礦產(chǎn)資源開采密切相關(guān)。

3)通過對(duì)特征點(diǎn)進(jìn)行時(shí)間序列分析，奎達(dá)市地面沉降在監(jiān)測(cè)時(shí)間段內(nèi)基本呈勻速下沉，僅在雨水相對(duì)充足的季節(jié)略有減緩的趨勢(shì)。

4)后續(xù)研究中可收集更長時(shí)間段影像數(shù)據(jù)，增加時(shí)間序列長度，或同時(shí)采用升軌和降軌數(shù)據(jù)進(jìn)行監(jiān)測(cè)，提高監(jiān)測(cè)成果的可靠性。

綜上，奎達(dá)市地面沉降影響范圍廣、沉降速率大，發(fā)展十分迅速，應(yīng)當(dāng)引起相關(guān)部門的高度重視，加強(qiáng)城區(qū)地面沉降監(jiān)測(cè)和防治工作、規(guī)范礦產(chǎn)資源開發(fā)秩序、引進(jìn)先進(jìn)的采礦技術(shù)和設(shè)備。否則，不僅阻礙社會(huì)經(jīng)濟(jì)的可持續(xù)發(fā)展，還會(huì)對(duì)當(dāng)?shù)鼐用竦纳?cái)產(chǎn)安全造成威脅。