適用于露天礦時序形變監測的優化DS-InSAR技術

李鳴庚 張書畢 高延東 李世金 張艷鎖 賈義琨 孔令辰

(1.中國礦業大學環境與測繪學院,江蘇 徐州 221116;2.自然資源部國土環境與災害監測重點實驗室,江蘇 徐州 221116)

露天開采是現階段我國資源獲取的主要手段之一[1],然而隨著露天礦開采深度的增加,工作面滑坡、排土場邊坡失穩、礦區地表沉陷等重大災害對露天礦安全生產造成的威脅也隨之增加[2-3]。遼寧省鞍山市是我國東北老工業基地重要的鐵礦資源城市,其周邊的西鞍山、東鞍山、大孤山、齊大山、眼前山、鞍千礦等十幾個露天礦區不可避免地面臨由于長期大規模露天開采引發的地表形變問題[4-5]。因此,對鞍山市周邊露天礦區進行地表形變監測對于城市健康發展具有重要意義[6]。

傳統地表形變監測方法如水準測量、全站儀、GPS (Global Position System)測量等,雖然精度較高,但也存在耗時長、成本高、且難以實現長時序、大范圍的地表形變監測等不足[7]。時序InSAR技術是在合成孔徑雷達干涉測量(Interferometric Synthetic Aperture Radar,InSAR)技術的基礎上發展而來,如小基線集(Small base lines subset InSAR,SBAS-InSAR)技術與永久散射體 (Permanent Scatters InSAR,PS-InSAR)技術,已在大范圍形變監測中得到了廣泛應用,具有高效率、低成本、長時序、大范圍等優勢[8]。時序In-SAR技術雖然在城市地區具有較高的監測精度與點位密度,但在裸地、植被覆蓋區與非城區則存在因相干性較低引起的監測點密度不足的問題[9]。在此基礎上,學者們提出了分布式目標InSAR (Distributed Scatterers InSAR,DS-InSAR)技術,可有效提高低相干區域的監測點密度[10],在長時序、大范圍的地表形變監測方面具有明顯優勢。

近年來,不少學者將DS-InSAR技術應用于大范圍地表形變監測中,取得了較為理想的結果。王波等[11]利用DS-InSAR技術獲取了西部礦區時序形變結果,并將DS-InSAR監測結果與水準數據進行對比,相關性系數為0.97,證明該技術具有較高的可靠性。李柱等[12]通過DS-InSAR技術對烏達煤田火區進行了長時序地表形變監測分析,通過與TCP-InSAR監測數據的對比分析,反映出DS-InSAR技術可靠性較好。吳昊等[13]通過DS-InSAR技術對巴西Brumadinho礦區尾礦庫進行了形變監測,結果表明:DSInSAR獲取的尾礦庫監測點密度為PS-InSAR的2.2倍,兩者獲取的形變速率具有較強的相關性。王新玲等[14]通過DS-InSAR技術對鄆城礦區進行了地表形變監測,結果表明:DS-InSAR技術可顯著提升裸地和稀疏低矮植被區的監測點數量,獲取的形變場更為完整。趙立峰等[15]針對常規InSAR技術在礦區受失相干影響嚴重、監測點數量低的問題,通過DS-InSAR技術提取了張雙樓煤礦長時序地表形變結果,并與傳統SBAS-InSAR技術進行對比分析,驗證了DS-InSAR監測精度。楊嘉威等[16]基于DS-InSAR技術對礦區鐵路線進行了形變監測,通過相位優化顯著提升了高相干像元數量與解算質量。

上述成果分析表明:雖然DS-InSAR技術在大范圍地表形變監測方面已有較多應用,但研究區以煤礦或露天礦尾礦庫、排土場居多,對于該技術在露天礦形變監測方面的應用研究相對較少。為此,本研究提出了一種基于FaSHPS同質像元提取與復相干矩陣特征值分解相位優化的DS-InSAR時序地表形變監測技術,獲取了2020年1月—2021年8月遼寧省鞍山地區齊大山礦、大孤山礦、鞍千礦3個礦區的時序地表形變特征,結合常規時序InSAR技術監測結果對比驗證了DS-InSAR技術的可靠性,豐富了DSInSAR技術在露天礦形變監測方面的研究,提高了DS-InSAR技術露天礦形變監測的精度,對城市周邊大范圍露天礦形變分析具有一定的參考意義。

1 研究區概況與數據來源

1.1 研究區概況

鞍山市位于遼寧省中部,是我國重要的鋼鐵生產基地,其周邊主要礦區分布如圖1所示。由圖1可知:東鞍山、齊大山、大孤山、眼前山等大型露天礦采場從東北至東南向南環繞鞍山市城區,最遠距市中心12 km,最近距市中心區僅2 km,形成一個寬約6 km、長為40 km的礦區生產活動帶[17]。鞍千礦業鐵礦位于胡家廟村和肖家堡村之間,據統計原有19個小型鐵礦采場和2家小型硅石礦采場,占地面積41.18 km2[18]。東鞍山、大孤山、齊大山各礦區占地面積均在1 000 km2以上,且距離市區較近,東鞍山鐵礦距離市中心僅7 km,距離最近的建成居民區僅2 km。由此可見,露天礦的安全監測對于鞍山市健康發展具有極其重要的意義。

1.2 數據源

為對鞍山市周邊礦區進行大范圍、長時序的地表形變監測,本研究選取歐空局2014年發射的Sentinel-1A衛星、單視復數(Single Look Complex,SLC)、升軌SAR影像,時間跨度為2020年1月9日—2021年8月31日,共50景,覆蓋范圍為41.013°～41.203°N,123.000°～123.224°E,如圖1所示。數字高程模型(Digital Elevation Model,DEM)采用SRTM (Shuttle Radar Topography Mission) DEM,分辨率為30 m。為削弱大氣誤差,本研究引入了合成孔徑雷達通用性大氣改正在線服務(Generic Atmospheric Correction Online Service for InSAR,GACOS)的天頂對流層延遲(Zenith Troposphere Delay,ZTD )產品[18]。另外,試驗還使用了Sentinel-1A POD精密定軌星歷(Precise Orbit Ephemerides,POD)數據,定位精度優于5 cm。

圖1 研究區地理位置Fig.1 Location of the study area

2 試驗方法

為實現對鞍山地區齊大山礦、大孤山礦、鞍千礦3個礦區大范圍、高精度的時序地表形變監測,本研究采用了一種基于FaSHPS同質像元提取與復相干矩陣特征值分解相位優化的DS-InSAR技術。該技術與傳統時序InSAR技術的區別在于增加了同質像元提取與相位優化2個步驟,能夠有效解決傳統時序InSAR技術在低相干區域監測點數量不足且空間分布不均勻的問題,在露天礦區長時序、大范圍地表形變監測方面具有明顯優勢。

2.1 FaSHPS同質像元提取

本研究采用快速同質像元識別算法FaSHPS(Fast

Statistically

homogeneouspixels

selection,FaSHPS)進行同質像元提取。該算法基于SAR影像任一像元的振幅均值服從高斯分布的假設,通過采用置信區間估計實現快速高效的同質像元識別,有效解決傳統的基于假設檢驗的同質像元識別算法效率不高的問題,且該算法提取的同質像元間異質性較低,在露天礦等高相干點分布不均且數量不足的區域具有明顯優勢[19]。設有N景已配準的SAR影像,對于影像上任一像元I,其時間維度上的振幅A(I)為

式中,振幅A(I)的均值(I)服從高斯分布,可由式(2)獲取(I)的區間估計:

式中,P{·}為概率;E(·)為期望;Var(·)為方差;Z1-α/2表示標準正態分布中置信度為α時所對應的分位點。

單視復數SAR影像上的任一像元在時間維度上服從瑞利分布,則變異系數Cv為一常數,計算公式為

假設在均質地區像元I后向散射穩定,則式(2)可轉化為

通過確定置信區間α,若待估計像元I的振幅平均值(I)在目標像元置信區間內,則為同質像元。

2.2 復相干矩陣分解法相位優化

為提高獲取的同質像元密度及可靠性,本研究采用復相干矩陣分解法進行相位優化。由于分布式目標在SAR影像上具有相似的后向散射性,因此復相干矩陣可通過特征值分解實現相位優化。該方法在運算過程中無需進行迭代,相位優化效率較高,可獲取高質量的優化相位[20]。分布式目標的相干性矩陣可表示為

式中,向量y=[y1,y2,…,yN]為分布式目標的同質像元在N景影像上的復數觀測量經過時間維振幅歸一化后的復數觀測值;NSHPs為分布式目標同質點數量;Ω為分布式目標同質點集;yH表示向量y的Hermitian共軛轉置。

由于分布式目標的相干性矩陣為半正定Hermitian矩陣,因此可根據特征值分解原理將式(5)轉化為

式中,λi為特征值;μi為與λi對應的特征向量;μ1為優化后的相位;為μ1的Hermitian共軛轉置;signal為主散射體信號;noise為噪聲信號。

完成同質像元提取與相位優化后,可利用干涉相位在經過相位優化前后的擬合度進行相位優化質量評價。通過設置相干性閾值篩選分布式目標中具有高相干性與高信噪比的部分,將篩選出的高相干點融合PS點進行時序InSAR處理,實現對研究區高精度形變監測。

本研究方法的數據處理流程如圖2所示。為保證運算效率,本研究對SLC采用2×1多視處理。

圖2 DS-InSAR數據處理流程Fig.2 Flow of DS-InSAR data processing

3 監測結果與分析

3.1 監測結果分析

利用DS-InSAR技術,基于2020年1月—2021年8月共50景Sentinel-1A影像,通過如圖3所示的時空基線組合,獲取鞍山市周邊礦區年形變速率變化特征。在缺少水準數據的情況下,采用交叉驗證[21]方法將DS-InSAR技術與常用的StaMPS-InSAR技術獲取的年沉降速率值進行對比分析,以驗證DS-In-SAR技術的可靠性。本研究以StaMPS監測點為基準,通過經緯度搜索最鄰近DS點為同名點,共獲取了16 890對同名點,占StaMPS總監測點數量的73.53%。通過對同名點的形變速率值進行相關性分析,獲取二者形變速率相關性與形變速率差值直方圖。

圖3 時空基線分布Fig.3 Layout of spatial-temporal baseline

鞍山市周邊礦區形變速率分布如圖4所示。由圖4可知:DS-InSAR與StaMPS監測結果具有較好的空間一致性。在研究區內, StaMPS技術共選出22 971個監測點;DS-InSAR技術共選出28 804個監測點,為StaMPS技術的1.253倍。DS-InSAR技術在尾礦庫、丘陵等低相干區域的點位數量明顯提升,且DS點在具有較高點位密度的同時,在空間分布上也較為均勻,在一定程度上提高了形變區的監測精度,提供的形變信息更加準確、完整。

圖4 鞍山市周邊礦區衛星視線向年平均沉降速率Fig.4 Annual average subsidence rates of LOS direction in Anshan City

StaMPS-InSAR與DS-InSAR監測結果交叉驗證獲取的形變速率相關圖及差值分布如圖5所示。

由圖5(a)可知:StaMPS技術與DS-InSAR技術獲取的同名點地表形變速率具有較高的相關性,通過Pearson函數計算出的兩者相關系數R為0.855 6。由圖5(b)可知:2種技術獲取的同名點形變速率差值較小,標準差為13.91 mm/a。在缺乏實測水準數據的情況下,通過DS-InSAR技術與StaMPS技術監測結果的交叉驗證,可有效反映出DS-InSAR技術的可靠性。

3.2 鞍山市周邊礦區大范圍地表形變分析

為系統研究鞍山市周邊礦區地表形變特征,本研究基于DS-InSAR技術,對鞍山市周邊礦區進行了長時序的地表形變監測,獲取了各時段相對于2020年1月9日的時序累積形變量,如圖6所示。從時間序列上分析,鞍山市周邊礦區自2020年1月9日起,出現明顯的、連續的地表形變現象。2020年1月9日—2020年9月17日,各形變區的形變范圍呈緩慢擴大趨勢,形成了A～E5個主要形變區。2020年11月4日—2021年8月31日,各形變區的擴大趨勢逐漸減緩,形變區內的累積形變量呈增加趨勢。從空間分布上分析,在監測時段內,鞍山市周邊礦區地表形變現象較為嚴重,出現了A～E5個主要形變區。其中,形變區A、B位于齊大山礦范圍內,形變區C、D主要位于鞍千礦范圍內,形變區E則位于大孤山尾礦庫以西、露天礦坑以東的區域。齊大山礦區形變范圍與形變量均較高,形變區A與形變區B呈自西北到東南的連續分布趨勢,與齊大山礦床的NW—SE走向一致。鞍千礦范圍內存在多個形變中心,主要形變集中于形變區C,形變區D存在多個形變中心,形變量均較小。形變區E位于大孤山范圍內,形變現象較為明顯,且大孤山露天礦坑西側也存在少量的形變現象。

圖6 鞍山市周邊礦區時序形變累積結果Fig.6 Time series deformation accumulation results of surrounding mining areas in Anshan City

為更詳細地分析鞍山市周邊礦區地表形變,結合圖4的DS-InSAR監測結果對各形變區進行了定量分析。綜合圖4、圖6可知:A區最高形變速率為-260.43 mm/a,B區為-258.49 mm/a;A區的形變中心存在因形變梯度較大導致失相干嚴重而出現DS點數量不足的現象;C區最高沉降速率為-229.05 mm/a,D區為-175.92 mm/a;E區最高沉降速率為-312.18 mm/a,高于上述各區。

3.3 鞍山市周邊礦區時序地表形變分析

為分析鞍山市周邊礦區時序形變特征,在如圖6所示的A～E5個主要形變區內,根據DS-InSAR監測結果選取了Q1～Q8共8個監測點,位置分布如圖7所示。提取上述8個監測點對應的時序形變結果,繪制時序形變曲線。同時選取Q4、Q6、Q7、Q84個DS監測點對應的同名StaMPS監測點Q4′、Q6′、Q7′、Q8′繪制時序形變曲線進行對比分析,以驗證DS-InSAR技術的可靠性。最終獲取的時序形變曲線如圖8所示。由圖8(c)、圖8(d)可知:DS監測點與對應的同名StaMPS監測點在時序形變趨勢上具有較高的一致性,但整體上形變量存在部分差異?？紤]到DS監測點與StaMPS監測點間存在一定距離與誤差的影響,認為形變量的差異處于合理范圍內,可以使用選取的DS監測點進行時序形變分析。

圖7 監測點位置Fig.7 Location of the monitoring points

圖8 各監測點的時序形變曲線及StaMPS對比結果Fig.8 Time series deformation curves and stamps comparison results of each monitoring point

圖8(a)為監測點Q1～Q3時序形變結果。監測點Q1形變速率最高,為-260.43 mm/a,累積沉降量為397.48 mm,且Q1點周圍出現因相干性較低而引起的DS點缺失現象,說明A區存在更為嚴重的沉降趨勢。Q2點沉降速率為-254.66 mm/a,累積沉降量為-407.35 mm。Q3點沉降速率為-258.49 mm/a,累積沉降量為-426.99 mm/a。3個監測點均位于齊大山礦采區內,自2020年1月起呈連續沉降趨勢。Q1點自2020年1月開始呈現較高的沉降速率,在2020年6月—2020年8月沉降速率出現減緩,2020年8月后沉降速率再次提升至較高水平。Q2、Q3點位于齊大山礦東南部,地處同一采區,整體形變趨勢上具有一定的相似性。Q2點的形變速率與形變量均低于Q3點,且Q3點在2020年全年均保持較高的沉降速率,進入2021年則明顯放緩。推測這種沉降速率的變化可能與監測時段內齊大山鐵礦的擴幫生產[22]規劃有關。

圖8(b)中,Q4點位于C區,為鞍千礦主要形變區。Q5、Q6位于D區。Q4點形變速率為-229.53 mm/a,累積形變量為-315.86 mm。Q5、Q6點的形變速率分別為-175.92 mm/a、-166.76 mm/a。累積形變量分別為-240.12 mm、-256.03 mm。分析圖8(b)可知:2020年1—8月,監測點Q4～Q6的形變速率均保持較低水平,在2020年9月后形變速率出現明顯提升,且呈連續的沉降趨勢。Q4點的沉降速率與沉降量均高于Q5與Q6點。Q5、Q6點在形變速率的變化趨勢與累積形變量方面具有一定的一致性。鞍千礦區內小型采區較多,因此沉降中心數量相對較多、空間分布上較為分散。C區推測可能為監測時段內鞍千礦的主要采區,沉降面積、沉降速率與沉降量均高于D區。

圖8(c)為監測點Q7、Q8時序形變結果。Q7、Q8點主要處于大孤山礦區內,Q7點位于大孤山露天礦坑以東、尾礦庫以西,形變速率為-312.18 mm/a,累積沉降量為-493.43 mm。Q8點位于大孤山露天礦坑西北側,形變速率為-110.61 mm/a,累積形變量為-169.88 mm。Q7點自2020年1月即呈連續沉降趨勢。除了在2020年10月、2020年12月、2021年5月等少數時段沉降速率有所降低外,其余時段內均保持較高的沉降速率。Q8點沉降速率與沉降量均較低,推測是由自然原因引發的地表沉降。2020年,大孤山鐵礦正處于露天開采轉為地下開采的過渡期[23],原本作為主要開采區的大孤山露天礦坑內部沉降速率相對較小,而礦坑西部、尾礦庫東部區域發生明顯的地表沉降,這種變化可能是由于大孤山礦由露天開采轉為地下開采引起的。

3.4 典型露天礦區形變特征分析

本研究選取A～E5個形變區中監測點位密度較高、分布較為均勻的C區作為典型露天礦區進行形變特征分析。為定量分析C區的形變特征,如圖9所示繪制了兩條剖面線x、y,并按照圖示方向在兩條剖面線上各選取了13個監測點,編號分別為X1,X2,…,X13,Y1,Y2,…,Y13繪制了形變區剖面時序下沉曲線,如圖10所示。

圖9 剖面線與監測點分布示意Fig.9 Schematic of profile lines and monitoring points

圖10 形變區剖面線時序下沉曲線Fig.10 Time series subsidence curves of profile lines in deformation area

由圖10(a)可知:2020年1月9日—2020年7月7日,C區的主要沉降集中在位于形變區西側的X3點附近,其余監測點整體呈下沉趨勢,但未出現明顯的沉降中心。2020年7月7日后,X3發生持續沉降的同時,位于形變區中部的監測點X6、X7點開始出現明顯的下沉趨勢,2020年12月10日前,X6點沉降量略高于X7點。2021年1月15日后,沉降中心開始向沉降區中部轉移,X6、X7點附近地區的沉降量較X3點更為明顯,取代X3點成為新的沉降中心。

由圖10(b)可知:2020年1月9日—2020年7月7日,位于C區南部的監測點Y10沉降現象明顯。8月12日后,位于沉降區中部的Y6、Y7點開始發生明顯的沉降現象,在2020年9月29日后沉降量超過Y10點。Y10點在2020年12月10日后形變速率降低,Y6、Y7點則持續保持較高的沉降速率,成為新的沉降中心。綜合剖面線x、y的沉降中心變化情況,可以推測2020年1月9日—7月7日,C區的主要開采區可能位于露天礦西南側。7月7日后,露天礦中心地帶的開采強度逐漸增大,而西南側的開采強度開始減小。進入2021年后,礦區西南側的開采強度明顯減小,主要開采區轉至礦區中部,成為新的沉降中心。

4 結 論

為實現對露天礦區的高精度、大范圍、長時序的地表形變監測,以鞍山市周邊礦區為例,研究了一種適用于露天礦時序形變監測的優化DS-InSAR技術。通過該技術對2020年1月9日—2021年8月31日共50景Sentinel-1A SAR影像進行處理,獲取了研究區的時序地表形變特征,并結合監測結果與相關資料進行了地表形變分析。得出如下結論:

(1)與傳統的StaMPS-InSAR技術相比,本研究采用的優化DS-InSAR技術在通過FaSHPs算法提取的研究區內部分布式目標的基礎上,采用復相干矩陣特征值分解法實現相位優化,再通過時序InSAR處理獲取露天礦區時序形變特征。該方法有效提高了獲取的同質像元密度及可靠性,在一定程度上解決了傳統時序InSAR技術在露天礦區因失相干嚴重而產生的監測點數量不足、空間分布不均勻的問題,實現了對露天礦區大范圍、長時序、高精度的地表形變監測。

(2)研究時段內,通過DS-InSAR技術監測到大孤山礦、齊大山礦與鞍千礦共存在5處沉降量級與形變范圍不同的形變區,最大形變速率為-312.18 mm/a。結合相關文獻資料分析表明,本研究采用的DS-InSAR技術監測到的形變區位置與范圍較為準確,獲取的時序地表形變監測結果具有較高的精度和可靠性。該方法可為露天礦區地表形變監測及其形變機理研究提供參考。

(3)本研究僅通過DS-InSAR技術對鞍山市周邊礦區進行地表形變監測,結合傳統的StaMPS-InSAR監測結果交叉驗證其可靠性,并通過相關資料推測各礦區開采活動,缺少具體實測數據資料與礦區開采資料進行輔助驗證。后續研究中可結合各露天礦具體開采量、開采時段、開采計劃及實測水準數據等資料與DS-InSAR時序監測結果進行對比,進一步分析DS-InSAR技術在露天礦地表形變監測中的適用性與可靠性。