周子琪 周世健 歐陽雙艷
(1. 東華理工大學 測繪工程學院, 江西 南昌 330013; 2. 南昌航空大學 環(huán)境科學與工程學院, 江西 南昌 330063)
改革開放至今,我國國民經濟高速發(fā)展的同時,對礦產資源的需求量也逐漸增大。各類礦產資源的大力開發(fā)對工業(yè)發(fā)展和社會經濟有著良好的促進作用[1],但開采過程可能會對礦區(qū)的地質結構帶來較大影響,甚至產生一系列地質災害問題,最易產生塌陷、地表沉降等現象,更為嚴重會導致地下水資源污染、巖體變形誘發(fā)崩塌、滑坡以及泥石流等地質災害的發(fā)生,嚴重影響當地地質結構和周圍居民的人身財產安全[2]。因此,為了保障礦區(qū)的可持續(xù)開采以及周邊居民的正常生活,除了選擇合理的開采方式以及減少礦區(qū)的開采量外,對礦區(qū)進行沉降監(jiān)測分析也顯得尤為重要。如今,通過技術手段獲取地面沉降信息及其分布特征可以及時有效采取對應措施,在一定程度上可以降低發(fā)生地質災害的風險并開展防治工作[3-5]。
在近幾十年大面積的沉降監(jiān)測研究當中,沉降監(jiān)測的方法也發(fā)生了變化,用于地表觀測的方式、手段也趨于多樣化,以往傳統(tǒng)觀測的方法包含:精密水準測量、三角高程測量、靜態(tài)全球導航衛(wèi)星系統(tǒng)(Global Navigation Satellite System,GNSS)測量和動態(tài)GNSS監(jiān)測等[6],這類傳統(tǒng)的監(jiān)測方式雖有較高的測量精度和可靠性,但對于大面積、需要長期動態(tài)監(jiān)測的區(qū)域而言,這些散點式的監(jiān)測方式不能很好地滿足實際的工程應用需要。為了滿足大面積的沉降監(jiān)測需要,近年來,合成孔徑雷達干涉測量(Interferometric Synthetic Rader,InSAR)的快速發(fā)展彌補了傳統(tǒng)變形監(jiān)測對于大面積監(jiān)測的不足。InSAR作為新興遙感技術,其中合成孔徑雷達差分干涉測量(Diffe-rential Interferometric Synthetic Aperture Rader,D-InSAR)技術具有連續(xù)覆蓋地表、高精度、自動化程度高的優(yōu)點[7-8],已成為獲取表面變形信息的有效手段,其精度可以達到厘米甚至毫米級別[9-11]。1989年,Gabrie等首次應用D-InSAR技術獲取了加州ImperialValley灌溉區(qū)的影像數據,首次說明D-InSAR技術可以探測厘米級的形變信息[12]。但由于D-InSAR技術易受時空失相干性以及大氣延遲效應的影響,應用的領域受限[13-15]。為了對D-InSAR技術進行優(yōu)化,2002年BERARDINO等提出了小基線集技術(SBAS-InSAR),該技術能夠針對時間跨度長、規(guī)模大的區(qū)域進行分析[16]。2011年江利明獲取了烏達的34景影像,其研究結果顯示,利用SBAS-InSAR技術處理出來的結果與GPS(Global Positioning System)結果一致[17],證明廣域差分增強系統(tǒng)(Satellite-Based Augmentation System,SBAS)-InSAR技術在地面沉降監(jiān)測應用的可靠性以及與其他技術監(jiān)測結果的一致性[17-18]。本文利用SBAS-InSAR技術,對德興礦區(qū) 2018年12月28日 至 2019年 12月24日 期間的30景SAR影像進行處理,得到該區(qū)域的沉降速率和沉降分布,對沉降嚴重的區(qū)域進行分析,此研究結果可以為該區(qū)域地質災害的防治問題提供一定幫助。
德興礦區(qū)位于江西省上饒市德興境內,距德興市22 km,其礦區(qū)位置位于117°43′40″E,29°01′26″N處。此次實驗研究區(qū)的地勢呈東南部偏高、西北部偏低,德興礦區(qū)地理位置屬于江南丘陵地區(qū),其整體地貌起伏較大,海拔高度約為65~500 m,礦區(qū)內地面切割較為強烈,山坡較陡,坡度在30°左右。德興礦區(qū)內礦產資源種類豐富,存礦量較大,其中大型和特大型銅金礦高度集中,產業(yè)特色明顯,該礦區(qū)是我國的礦業(yè)經濟中最為發(fā)達與發(fā)展?jié)摿ψ顬樾酆竦牡貐^(qū)之一,同時也是國家級重要有色貴金屬能源基地之一。德興礦區(qū)周圍分布有祝家社區(qū)、富家塢社區(qū)、詹家塢社區(qū)以及梨園社區(qū),礦區(qū)的開采在一定程度上間接影響了當地居民的生活環(huán)境。同時,德興礦區(qū)在長期的挖掘開采過程中,積累了一定程度的礦區(qū)地質問題,這些問題會給礦區(qū)從業(yè)者和周圍村民帶來潛在的安全隱患。
本文實驗數據選取覆蓋德興市礦區(qū)的30景升軌Sentinel-1A數據,Sentinel-1A數據,其時間跨度為2018年12月29日至2019年12月24日,影像時間間隔為12 d。本實驗使用的Sentinel-1A數據均為C波段,波長為5.66 cm,入射角大約為39.3°。采用由美國航空航天局(National Aeronautics and Space Administration, NASA)提供的30 m分辨率航天飛機雷達地形測繪使命(Shuttle Radar Topography Mission,SRTM)1數據作為外部參考數字高程模型(Digital Elevation Model,DEM)以消除地形誤差對實驗的影響。Sentinel-1A數據參數見表1。
表1 Sentinel-1A數據參數
SBAS-InSAR是以傳統(tǒng)D-InSAR技術為基礎再進行研究,用來獲取工作區(qū)地表形變的時間序列圖。該技術選取N+1幅SAR影像,通過設置合適的時間基線和空間基線生成N個小基線集差分干涉對,去除地形相位后生成差分干涉圖并且進行相位解纏,將常規(guī)D-InSAR監(jiān)測的觀測結果用作單個觀測值,采用奇異值分解(Singular Value Decomposition,SVD)的方法將每一組小基線數據集連接起來,解決時間上采樣過于稀疏的問題,又結合穩(wěn)定散射體的干涉相位信息,得到更高的空間分辨率,然后根據最小二乘法求出高精度的變形時間序列。
本文實驗假定在感興趣區(qū)域的有序時間(t1,t2,…tN)內,采集研究區(qū)內同軌道上的一組N+1景SAR圖像,根據干涉條件,將輔影像配準到主影像上,在配準完成后根據連接圖生成的干涉圖進行復相位運算。通過設置合理的時間和空間基線生成M個差分干涉圖。其中M需要滿足以下條件[19]
(1)
在干涉步驟中,在tA和tB時刻(tA>tB)獲取的兩幅SAR影像作為主輔影像生成的第j幅差分干涉圖,第j幅差分干涉圖中任何像元的干涉相位可寫成
(2)
式中,λ為雷達波長;d(tA,x,r)和d(tB,x,r)是在tA和tB時刻相對于參考時間t0的雷達視線方向的累積形變;Δφjtop是由所參考的DEM不精確而引起的地形相位誤差;Δφjatm為大氣延遲相位誤差;Δφnoise是噪聲影響的相位誤差。
為去除地形相位誤差、大氣延遲誤差以及噪聲相位誤差,通過簡化式(2),干涉相位可以表示為
(3)
為了獲得含有物理意義的沉降序列,將式(3)中相位與時間之比來表示獲取影像時間中的平均相位速度
(4)
其中,第j幅干涉圖的相位值可以表示為
(5)
用SBAS-InSAR算法計算時間序列的變形可表示為
AV=δφ
(6)
其中,AV為M×N矩陣;δφ是代表干涉相位值的向量。最后使用SVD方法用于獲得最終形變速率。
SBAS-InSAR技術處理流程圖如圖1所示。
圖1 SBAS-InSAR技術處理流程圖
實驗使用GAMMA軟件進行SBAS-InSAR處理。GAMMA軟件可以完整地處理雷達信號,其中的功能模塊包括:從合成孔徑雷達(Synthetic Aperture Radar,SAR)原始信號處理到單視復數(Single Look Complex,SLC)成像、單視以及多視處理、基于雷達信號濾波、正射糾正和配準、DEM提取(干涉)、形變分析(差分干涉、點目標干涉分析)、土地利用等,同時也可以處理各種星載、機載及地基雷達的觀測數據。
3.2.1影像裁剪
影像覆蓋范圍廣,數據量大,為了提高數據的處理效率和監(jiān)測精度,根據研究區(qū)域的經緯度,裁剪出研究區(qū)域的影像數據。
3.2.2差分干涉對組合
實驗選取2019年6月15日的SAR影像作為公共主影像,將輔影像配準到主影像上,使配準精度達到0.001像元。將多視系數設置為4∶1,對干涉對進行多視處理。本實驗的時間基線和空間基線閾值分別為60 d和200 m,最終構成112對干涉對用于時序分析。圖2為干涉對時空基線分布圖。
圖2 干涉對時空基線分布
3.2.3干涉圖的生成及相位解纏
根據干涉對時空基線分布圖對影像進行干涉處理,利用外部DEM數據去除地形相位的影響,生成差分干涉圖,對差分干涉圖進行濾波增強處理,采用最小費用流的方法進行相位解纏,得到研究區(qū)域相位的真實值。通過查看每組生成的相干系數圖、濾波后的干涉圖和解纏圖,移除干涉質量較差的干涉對。
3.2.4研究區(qū)平均沉降速率的獲取
在具有高相干性的像素點上建立模型,通過 SVD 法反演估算形變速率,并去除殘余地形,計算相位殘差,進行殘差分離,最終得到研究區(qū)域內的地表沉降速率以及地表沉降分布情況。
采用SBAS-InSAR技術對30景升軌Sentinel-1A數據進行處理,其結果表明,在該研究區(qū)域內有三處明顯沉降,分布在研究區(qū)西南部、東南部以及北部。
圖3為研究區(qū)域圖,如圖所示,將這三個沉降區(qū)域分別命名為A、B、C。
圖3 研究區(qū)域圖
圖4為A區(qū)域沉降速率圖,從實驗結果可知:研究區(qū)域西南面的A區(qū)域沉降最為嚴重,最大年平均沉降速率為-490 mm/a。經核查該區(qū)域采礦活動頻繁,地表有多處塌陷且出現積水。為了進一步研究監(jiān)測期間礦區(qū)的累積沉降量,在沉降嚴重的A區(qū)域選取三個特征點進行時序分析,分別為點A#1、A#2和A#3,對特征點進行時間序列分析得到時間序列形變圖,通過結果可知,在2018年12月29日~2019年12月24日內,隨著時間的推移,該區(qū)域形變呈現連續(xù)下降的趨勢。從圖5可知,位于該區(qū)域沉降中心的采樣點A#3,該點的最大累積沉降量高達440 mm,該點存在連續(xù)下沉的趨勢,且沉降速率較快。A#1和A#2兩點選取在沉降區(qū)邊緣,A#1和A#2兩個特征點的沉降量分別為-184 mm、-131 mm。沉降區(qū)中心特征點的累積沉降量明顯高于沉降區(qū)邊緣,該區(qū)域沉降形成明顯的漏斗狀。沉降區(qū)域主要集中在祝家社區(qū)一帶,該區(qū)域的沉降會對村民的生活造成影響,相關部門應引起高度重視。
圖4 A區(qū)域沉降速率圖
圖5 A區(qū)域特征點時間序列形變圖
沉降區(qū)域B位于研究區(qū)東南處,從圖6可知其最大年平均沉降速率可達-390 mm/a,在B區(qū)域選取三個特征點B#1、B#2和B#3,對其進行時間序列分析得到特征點時間序列形變圖,從圖7可以看到沉降中心特征點B#3的最大累積沉降量達到358 mm,在整個研究的時間跨度之內一直處于持續(xù)下沉的狀態(tài),并且未來還有進一步下沉的趨勢。B#1和B#2的累積沉降量分別為-120和-118 mm。沉降中心的沉降速率向沉降邊緣逐漸降低,特征點的累積形變量從沉降中心往沉降邊緣逐漸減少。
圖6 B區(qū)域沉降速率圖
圖7 B區(qū)域特征點時間序列形變圖
沉降區(qū)域C位于研究區(qū)北處,從圖8可知,其最大年平均沉降速率可達233 mm/a,在C區(qū)域選取兩個特征點C#1、C#2和C#3,對其進行時間序列分析得到時間序列形變圖(圖9),從圖9可以看出,特征點C#1沉降量最大,其累積沉降
圖8 C區(qū)域沉降速率圖
圖9 C區(qū)域特征點時間序列形變圖
量達到-148 mm。由于該區(qū)域靠近水體,該地的沉降不僅因為頻繁的采礦活動,也和水體周圍的地質有著密不可分的關系。
本文采用2018—2019年間30景C波段Sentinel-1A升軌雷達數據為實驗樣本進行監(jiān)測,結合SBAS-InSAR技術對江西德興礦區(qū)進行地表沉降監(jiān)測分析,反演出研究區(qū)域真實地表形變,進一步分析該區(qū)域沉降速率及嚴重區(qū)域特征點的累積形變量,得出以下結論:
(1)在2018年12月29日~ 2019年 12月24日研究期間內,研究區(qū)域有多處沉降,分布在研究區(qū)西南部、東南部、北部。西北部未見明顯沉降。
(2)研究區(qū)中最為嚴重的沉降區(qū)域集中在礦區(qū)西南部,年平均沉降速率最大達到490 mm/a,在沉降最嚴重的區(qū)域中心選取特征點進行分析,該點的累積沉降量最大可達440 mm,該特征點未來還有進一步下沉的趨勢。
(3)祝家社區(qū)位于研究區(qū)德興銅礦內,礦區(qū)的過度開采嚴重影響了當地村民的生活,導致該區(qū)域地面塌陷、農田減少、地下水資源污染,嚴重影響當地居民的生命財產安全。今后需把該區(qū)域作為重點監(jiān)測區(qū)域并對其進行持續(xù)監(jiān)測,其應有效合理制定開采計劃,防止礦區(qū)沉降問題失控、周圍居民生活得不到保障等情況出現。
綜上所述,通過SBAS-InSAR技術可以對礦區(qū)開采過程中產生的地表沉降進行有效監(jiān)測,可以快速、精準地捕捉興趣區(qū)內沉降區(qū)域的具體位置和大致范圍,結果可知其沉降范圍與礦區(qū)的開采過程大致相同。對沉降結果分析有助于防治此類災害的發(fā)生,若持續(xù)開采可能會造成各類地質災害問題。在礦山持續(xù)開采的過程中,針對地表沉降嚴重而易引發(fā)地質災害等相關問題,對沉降數據進行準確預測就顯得尤為重要,如何對沉降的時間序列數據進行有效精準預測則是下一步研究的重點。