鄆城礦區(qū)地表形變時空演變DS-InSAR監(jiān)測

王新玲，杜玉玲，蔣金雄，趙 峰，閆世勇

（1.青島濱海學院 機電工程學院，山東 青島266555；2.中國礦業(yè)大學 自然資源部國土環(huán)境與災(zāi)害監(jiān)測重點實驗室，江蘇 徐州221116）

我國是世界上煤炭資源最豐富的國家之一，隨著煤炭資源的持續(xù)開采，礦區(qū)地面沉降日益嚴重，導(dǎo)致沉降區(qū)內(nèi)地面損毀、莊稼及道路被淹、橋體開裂等危害發(fā)生，給當?shù)厝嗣竦纳a(chǎn)、生活帶來了極大的安全隱患[1-2]。因此，及時掌握礦區(qū)地表形變及其動態(tài)變化尤為重要。

傳統(tǒng)的監(jiān)測技術(shù)如精密水準測量和GNSS等方法耗時耗力、空間采樣低，難以滿足礦區(qū)大范圍長時序地表沉降監(jiān)測及其穩(wěn)定性評估的需求[3-4]。隨著新型空間對地觀測技術(shù)的不斷發(fā)展建設(shè)，合成孔徑差分干涉測量（D-InSAR）技術(shù)因其測量精度高、觀測范圍大、不受天氣限制等技術(shù)優(yōu)勢，在火山[5]、地震[6]、冰川[7]、滑坡[8]和礦區(qū)沉降[9]等領(lǐng)域得到了廣泛的應(yīng)用，但D-InSAR技術(shù)易受各類失相干及大氣相位延遲等因素影響，不適用于形變速度緩慢、量級較小的地表形變監(jiān)測研究。時序InSAR分析技術(shù)[10-11]則在一定程度上克服了上述不足，為廣域小量級形變的監(jiān)測提供了可能的方法。近年來，為進一步提高時序形變估計質(zhì)量，融合DS和PS目標的時序InSAR技術(shù)受到廣泛的關(guān)注[12-14]，圍繞同質(zhì)像元的識別[12,15]和相位優(yōu)化[14]等重要方面開展了相關(guān)研究，并逐漸應(yīng)用于多個領(lǐng)域的地表形變監(jiān)測[11-19]，但目前針對大采深煤礦引起地表形變的應(yīng)用較少，針對性研究有所不足。為此在已有研究的基礎(chǔ)上，采用基于相干矩陣特征值分解的DS-InSAR技術(shù)，開展鄆城礦區(qū)地面沉降及其時空演變的高精度監(jiān)測分析，實現(xiàn)對礦區(qū)地面采空塌陷區(qū)域的穩(wěn)定性評估，為塌陷區(qū)域的恢復(fù)治理提供必要的技術(shù)與數(shù)據(jù)支撐。

1 研究區(qū)概況與數(shù)據(jù)源

1.1 礦區(qū)概況

鄆城礦區(qū)位于山東省西南部、菏澤市東北部，面積約624.6 km2，屬于溫帶半濕潤季風區(qū)海洋-大陸性氣候，四季分明，春旱多風，夏熱多雨。該地區(qū)大地構(gòu)造屬中朝陸臺山東臺背斜的西南邊緣凹陷帶，處于“魯西斷塊”中的“魯西南塊陷”上，地表全為新生界覆蓋，沒有基巖出露。鄆城地區(qū)煤炭礦產(chǎn)資源豐富，主要有李樓煤礦、彭莊煤礦、郭屯煤礦、趙樓煤礦等4個生產(chǎn)礦井，地表覆蓋以耕地、城鎮(zhèn)工礦和水域水利設(shè)施為主，礦區(qū)開采沉陷對耕地及交通、水利等基礎(chǔ)設(shè)施造成了一定的破壞[20]，嚴重威脅了當?shù)厝嗣竦纳敭a(chǎn)安全，因此對該區(qū)域地表形變進行快速高效監(jiān)測顯得至關(guān)重要。鄆城礦區(qū)位置示意圖如圖1。

1.2 數(shù)據(jù)源

選取了2018年1月3日至2019年11月30日期間Sentinel-1衛(wèi)星獲取的55景TOPS模式SAR影像，并以2019年3月11日的SAR數(shù)據(jù)為主影像進行數(shù)據(jù)的準確配準。歐洲航空局提供的精密軌道數(shù)據(jù)用于干涉處理中的軌道精估計和平地效應(yīng)去除。此外，采用30 m分辨率的SRTM DEM數(shù)據(jù)進行地理編碼和地形相位補償。

圖1 鄆城礦區(qū)位置示意圖Fig.1 Location distribution map of Yuncheng Mining Area

2 研究方法

不同于傳統(tǒng)基于高相干點目標的時序InSAR技術(shù)，融合分布式目標的DS-InSAR技術(shù)能夠提取到DS點和PS點目標進行形變解算，從而提升時序InSAR測量中相干目標的空間分布密度，其中同質(zhì)像元（SHPs）識別[21]和相位優(yōu)化估計[22]是DS-InSAR的2個核心內(nèi)容。融合分布目標的時序InSAR技術(shù)主要流程如圖2。

圖2 DS-InSAR數(shù)據(jù)處理流程圖Fig.2 DS-InSAR data processing flowchart

2.1 同質(zhì)像元識別

同質(zhì)像元的概念由Ferretti首次提出[23]，同質(zhì)像元的識別是在假設(shè)檢驗的基礎(chǔ)上，識別與樣本數(shù)據(jù)具有相同/相似散射特性的分布式目標像元作為同質(zhì)像元的集合。目前基于后向散射信息的參數(shù)假設(shè)檢驗是同質(zhì)像元的識別常用方法，其在已知樣本分布函數(shù)的前提下，判斷樣本數(shù)據(jù)的總體分布統(tǒng)計量是否存在顯著性差異來進行同質(zhì)像元的識別，選擇雙總體T檢驗進行同質(zhì)像元的識別。

2.2 相位優(yōu)化估計

針對分布式目標多視干涉相位，需要采取一定的相位優(yōu)化算法在滿足相位一致性的條件下構(gòu)建出1組單一主影像優(yōu)化相位值，獲取經(jīng)自適應(yīng)多視降噪處理后的1組最佳擬合相位，以盡可能地減弱分布式目標失相干現(xiàn)象的影響。所選用的特征值分解方法是通過從相干矩陣中有效地分離、估計出多元散射機制并得到所對應(yīng)的相位分量，實現(xiàn)對單一主影像干涉相位的優(yōu)化估計。為評估DS優(yōu)化相位的質(zhì)量，通過擬合比較原相干矩陣中自適應(yīng)多視干涉相位和優(yōu)化相位得到的干涉相位，即用時間相干性衡量γDS分布式目標受時間去相干的影響程度，同時也可用于估計相位優(yōu)化的質(zhì)量。

式中：γDS為分布式目標的時間相干性；φ^om、φ^on為分布式目標在單一主影像下的優(yōu)化相位；φmn為自適應(yīng)多視干涉相位；m、n、o為影像序列號。

3 結(jié)果分析

3.1 DS-InSAR監(jiān)測結(jié)果

采用融合分布式目標的DS-InSAR時序分析方法能夠在研究區(qū)內(nèi)裸地和農(nóng)田區(qū)域選取到更多的觀測點，從而更好地達到表征礦區(qū)形變空間分布特征的目的，鄆城礦區(qū)平均沉降速率圖如圖3。地面沉降速率等值線圖如圖4。

圖3 鄆城礦區(qū)平均沉降速率圖Fig.3 Average subsidence rate chart of Yuncheng Mining Area

2018—2019年鄆城礦區(qū)內(nèi)4個礦井均有沉降發(fā)生，其中李樓煤礦和郭屯煤礦沉降較為嚴重，影響范圍最廣，產(chǎn)生了1個較大的沉降漏斗，形變中心沉降速率達到了100 mm/a；彭莊煤礦和趙樓煤礦沉降量相對較小，最大沉降速率在70 mm/a。年形變速率對比可知，變形均集中在李樓煤礦和趙屯煤礦，相較于2018年1個大的沉降漏斗，2019年的形變在南北方向上表現(xiàn)出3個形變漏斗，形變中心速率分別可達80、80、70 mm/a。整體監(jiān)測結(jié)果表明，鄆城礦區(qū)煤礦開采導(dǎo)致地表發(fā)生了不同程度的沉降且分布不均。

3.2 重點礦區(qū)形變分析

為進一步分析煤炭開采引起的地面沉降，結(jié)合礦區(qū)信息，繪制了地面沉降等值線圖，并針對研究區(qū)域內(nèi)郭屯煤礦和彭莊煤礦2個典型礦井，進行了時間和空間維度的精細化解譯分析。

圖4 地面沉降速率等值線圖Fig.4 Contour map of land subsidence rate

郭屯煤礦北臨李樓煤礦、東臨彭莊煤礦、南臨趙樓煤礦，研究時段內(nèi)該煤礦的開采工作面為一采區(qū)1325工作面和四采區(qū)的4303和4302工作面，根據(jù)2018年和2019年的沉降速率等值線圖可知，2018年礦區(qū)中部和北部連接形成1個較大范圍的沉降區(qū)域，其最大沉降量約為100 mm/a，而在2019年，之前的大區(qū)域沉降表現(xiàn)為3塊小的沉降漏斗，其中礦區(qū)北部與李樓煤礦形成1個完整的沉降漏斗，且沉降中心發(fā)生了一定的偏移，最大形變量為80 mm/a；中部地區(qū)的形變主要集中在一采區(qū)和四采區(qū)工作面附近，與實際開采情況相一致；還有一處沉降集中在工業(yè)廣場附近，最大形變量從2018年的90 mm/a變?yōu)?019年的70 mm/a?？梢?，郭屯煤礦的地表形變主要與煤礦開采相關(guān)。

彭莊煤礦位于鄆城礦區(qū)的東側(cè)，在2018-2019年間，該煤礦對東二采區(qū)的3303和3304工作面進行開采，期間形成了4個小的沉降漏斗。由平均沉降速率等值線圖可知，2018年該煤礦的沉降主要發(fā)生在礦區(qū)西部，東部邊緣僅有小范圍沉降，而2019年東西區(qū)域內(nèi)的沉降區(qū)域趨于融合形成1個大的沉降漏斗。監(jiān)測結(jié)果表明，彭莊煤礦左側(cè)3處沉降漏斗分布與先前研究發(fā)現(xiàn)的采空塌陷區(qū)分布基本一致[24]，可見老采空區(qū)依舊有持續(xù)沉降現(xiàn)象的發(fā)生但相對有減弱趨勢，而新的采掘活動也逐步形成了新的沉降漏斗。

3.3 累計沉降量分析

為進一步說明各礦區(qū)的形變趨勢，在礦區(qū)顯著沉降區(qū)內(nèi)部選取特征點，（圖4中點1～點6），繪制的礦區(qū)典型沉降中心時序累積沉降圖如圖5。

圖5 礦區(qū)典型沉降中心時序累積沉降圖Fig.5 Time series cumulative subscidence diagram of typical subsidence center in mining area

由圖5可知，李樓煤礦（點1）、彭莊煤礦（點2）和郭屯煤礦（點3）的形變中心均表現(xiàn)為持續(xù)下沉特征，其中，點3和點4處的形變量級和形變趨勢趨于一致，且形變量均達到了130 mm；點5處的形變相較于點3和點4較為緩和，而礦區(qū)東側(cè)（梁寶寺煤礦，點6）的形變中心在2018年雖然有一定波動但總體保持穩(wěn)定，在2019年則呈現(xiàn)出下沉趨勢。資料顯示李樓煤礦和郭屯煤礦的生產(chǎn)能力均為240萬t/a，彭莊煤礦的生產(chǎn)能力為110萬t/a，可見，點3、點4與點5處的不同形變速率可能與開采量有關(guān)。而點6處沉降中心位于梁寶寺煤礦，開采塌陷未沉穩(wěn)[25]，且根據(jù)該地的夏季炎熱多雨氣候特點，沉陷區(qū)形變監(jiān)測結(jié)果會在夏季產(chǎn)生較大的形變波動。

4結(jié)語

1）基于DS-InSAR技術(shù)和Sentinel-1雷達影像獲取了鄆城礦區(qū)2018—2019年的地表沉降分布，并對典型區(qū)域的地表形變及其時空演變規(guī)律進行了重點分析，結(jié)果表明鄆城礦區(qū)存在多個明顯的沉降中心，且最大形變可達145 mm，鄆城礦區(qū)形變主要體現(xiàn)在李樓煤礦和郭屯煤礦，彭莊煤礦和趙樓煤礦形變量級相對較小。

2）與傳統(tǒng)的時序InSAR方法相比，DS-InSAR技術(shù)通過同質(zhì)像元識別、相位優(yōu)化等步驟能夠在裸地和稀疏低矮植被等區(qū)域有效識別出分布式目標，提高了測量點數(shù)量和空間分布密度，從而準確、全面地反映了礦區(qū)開采沉降的影響范圍和形變區(qū)域的空間分布特征，為礦區(qū)穩(wěn)定性分析和土地合理利用等提供了良好的技術(shù)和數(shù)據(jù)支撐，有助于預(yù)防和減少礦山建設(shè)及生產(chǎn)活動造成的礦山地質(zhì)環(huán)境及地質(zhì)災(zāi)害。