陳志達(dá)，龐校光

（1.紹興市上虞區(qū)自然資源監(jiān)測(cè)中心，浙江 紹興 312399；2.中國科學(xué)院精密測(cè)量科學(xué)與技術(shù)創(chuàng)新研究院 大地測(cè)量與地球動(dòng)力學(xué)國家重點(diǎn)實(shí)驗(yàn)室，湖北 武漢 430077；3.中國科學(xué)院大學(xué)，北京 100049）

InSAR技術(shù)是20世紀(jì)后期發(fā)展起來的衛(wèi)星遙感新方法，可以獲取高精度地形和地表形變信息[1-2]。與常規(guī)地面水準(zhǔn)測(cè)量、全球定位系統(tǒng)（GPS）、光學(xué)遙感等測(cè)繪手段相比，InSAR技術(shù)進(jìn)行地面形變監(jiān)測(cè)時(shí)優(yōu)勢(shì)明顯[3-5]，且已在地表形變[6-8]、礦區(qū)[9-10]、滑坡[11]及地震形變[12]等研究與監(jiān)測(cè)中得到了廣泛的應(yīng)用。

礦區(qū)開采引起的地面沉降是非常復(fù)雜的過程，其危害十分嚴(yán)重，研究礦區(qū)塌陷具有重要的意義[13-14]。2011年，Jiang[15]等利用時(shí)序InSAR對(duì)礦區(qū)煤火引起的地表形變進(jìn)行了監(jiān)測(cè)，結(jié)果顯示時(shí)序InSAR技術(shù)對(duì)礦區(qū)煤火引起的地表形變具有可靠的應(yīng)用潛力。2014年，Dong[16]等利用差分干涉測(cè)量短基線集時(shí)序分析技術(shù)（SBAS-InSAR）對(duì)淮南地下煤礦開采引起的地面沉降進(jìn)行時(shí)空分析，證明了L波段InSAR可用于監(jiān)測(cè)礦區(qū)地下開采活動(dòng)引起的地面沉降。同年，Zhao[17]等通過InSAR技術(shù)對(duì)大同不同采礦沉降階段及時(shí)間演變進(jìn)行分析，發(fā)現(xiàn)1992—2011年間約80%煤礦發(fā)生形變。2015年，LI[3]等對(duì)徐州礦區(qū)進(jìn)行實(shí)驗(yàn)，再次驗(yàn)證了In-SAR技術(shù)用于礦區(qū)形變監(jiān)測(cè)結(jié)果的可靠性。2019年，朱建軍[18]利用InSAR技術(shù)求取了煤礦開采沉陷的三維形變場(chǎng)，并進(jìn)行了實(shí)驗(yàn)?zāi)M。

山西省大同市渾源縣煤炭資源豐富，從20世紀(jì)70年代以來，煤礦無序開采導(dǎo)致周邊地區(qū)大范圍地面沉降，嚴(yán)重破壞了周圍的生態(tài)環(huán)境，直接威脅著當(dāng)?shù)鼐用竦纳柏?cái)產(chǎn)安全。為此，本文利用覆蓋大同市渾源縣46景SAR影像，采用時(shí)序InSAR技術(shù)進(jìn)行分析處理，獲取了2007—2013年間渾源縣平均年形變速率結(jié)果和時(shí)序結(jié)果，并根據(jù)實(shí)地調(diào)查資料分析了地面沉降時(shí)空變化特征及其與采空區(qū)分布的關(guān)系。

1 研究區(qū)概況及實(shí)驗(yàn)數(shù)據(jù)

1.1 研究區(qū)概況

渾源縣（39.67°～39.83°N、113.72°～113.92°E）地處黃土高原的邊緣地帶，位于山西省東北部大同盆地東南邊緣，東臨廣靈縣，南鄰繁峙縣，西毗應(yīng)縣，北接云州區(qū)。渾源縣地處桑干河支流渾河中上游，全縣南北長(zhǎng)57 km，東西寬約52 km，總面積1 968 km2。渾源縣地形復(fù)雜，東、南、北三面環(huán)山，境內(nèi)起伏較大，海拔最高的臥羊場(chǎng)山主峰海拔2 333 m，最低的渾河西部河床海拔1 026 m，相對(duì)高差1 307 m。境內(nèi)山脈較多，海拔均在2 000 m以上，山勢(shì)陡峭，溝壑縱橫[19]。

渾源縣礦產(chǎn)資源豐富，儲(chǔ)量大，包含煤炭、花崗巖、鐵礦等，已探明的煤礦儲(chǔ)量約17億t。自上世紀(jì)70年代以來，粗放式煤礦開采破壞了當(dāng)?shù)氐纳鷳B(tài)環(huán)境，威脅著人民的生命財(cái)產(chǎn)安全[17,20]。本文選取整個(gè)渾源縣為研究區(qū)域，圖1為渾源縣概況及數(shù)據(jù)覆蓋情況。

圖1 研究區(qū)域和數(shù)據(jù)覆蓋情況

1.2 實(shí)驗(yàn)數(shù)據(jù)

ALOS-1于2006-01-24發(fā)射升空，是日本新一代的對(duì)地觀測(cè)衛(wèi)星。其攜帶的L波段SAR適用于地表破壞嚴(yán)重的煤礦開采區(qū)域，可有效提取渾源縣礦區(qū)地表形變范圍和形變量信息。RADARSAT-2是一顆搭載C波段傳感器的高分辨率商用雷達(dá)衛(wèi)星，于2007-12-14發(fā)射升空，是RADARSAT-1的后續(xù)衛(wèi)星，在新的圖像獲取能力及性能方面，有長(zhǎng)足進(jìn)展。該衛(wèi)星具有最高1 m高分辨率成像能力，多種極化方式。

本文 獲取2007年6月—2011年1月36景ALOS PALSAR影 像 和2012年1月—2013年4月10景RADARSAT-2影像數(shù)據(jù)，以獲取渾源縣長(zhǎng)期地表高程變化速率及時(shí)間序列。ALOS/PALSAR數(shù)據(jù)分辨率為10 m，入射角為34.3°。RADARSAT-2數(shù)據(jù)分辨率為30 m，入射角為21.25°。DEM數(shù)據(jù)采用NASA提供的SRTM，其空間分辨率為90 m，用以去除干涉對(duì)中的地形相位。此外，還搜集5大煤礦區(qū)域、99座小煤礦坐標(biāo)以及煤礦的地質(zhì)報(bào)告等煤礦開采資料，用于InSAR監(jiān)測(cè)結(jié)果的驗(yàn)證與分析。通過獲取渾源縣地面形變范圍和形變量信息，結(jié)合煤礦開采存檔、野外實(shí)地調(diào)查等輔助數(shù)據(jù)，進(jìn)行相互驗(yàn)證。

2 SBAS-InSAR數(shù)據(jù)處理

2.1 SBAS-InSAR基本原理

SBAS-InSAR方法原理是將D-InSAR形變結(jié)果作為觀測(cè)值，通過SVD方法將多個(gè)小基線集聯(lián)合進(jìn)行最小二乘解求解，獲得的形變時(shí)間序列可有效反演出地表形變范圍和形變量信息。對(duì)某一區(qū)域不同時(shí)間獲取N+1幅SAR影像，選擇主影像后將其他影像配準(zhǔn)至主影像，為避免時(shí)空失相干造成觀測(cè)誤差及影響監(jiān)測(cè)結(jié)果可靠性，設(shè)置時(shí)空基線閾值以保證干涉對(duì)的有效性，后利用外部DEM去除參考橢球相位與地形相位，對(duì)形變相位進(jìn)行解纏及地理編碼得到地理坐標(biāo)系下研究區(qū)域形變結(jié)果[1-2,21]。兩景SAR影像干涉相位可用φint表示為：

式中，φflat為參考橢球相位；φtopo為地形相位；φatm為大氣延遲相位；φorbit為軌道誤差；φnoise為噪聲相位；φdef為地表形變相位。對(duì)去除各項(xiàng)誤差后的差分干涉對(duì)地表形變相位進(jìn)行解纏后，可獲取視線向的位移量研究區(qū)域各像元形變量。

2.2 數(shù)據(jù)處理

本文對(duì)ALOS/PALSAR與RADARSAT-2影像數(shù)據(jù)分別處理，干涉對(duì)主輔SAR數(shù)據(jù)精配準(zhǔn)精度均控制在0.1個(gè)像元內(nèi)，空間基線均為最大基線的45%，最大基線分別為1 984.26 m和-271.49 m，時(shí)間基線分別為365 d和180 d，最終分別生成70和24個(gè)干涉對(duì)。對(duì)各干涉對(duì)分別進(jìn)行精密軌道數(shù)據(jù)優(yōu)化空間基線，消除軌道誤差，去除地平效應(yīng)，利用外部SRTM-DEM去除地形相位，干涉圖濾波，相位解纏，消除大氣相位延遲等處理得到地表形變相位，對(duì)獲取的形變干涉圖采用多基線方法獲取時(shí)間序列[22-23]，流程圖如圖2所示。

圖2 SBAS-InSAR方法流程

3 結(jié)果與討論

3.1 地表形變結(jié)果

利用覆蓋渾源縣的ALOS/PALSAR和RADARSAT-2衛(wèi)星雷達(dá)遙感數(shù)據(jù)，采用SBAS-InSAR技術(shù)獲取了渾源縣地表形變結(jié)果，包括渾源縣年均形變速率與重點(diǎn)區(qū)域的累積形變量。根據(jù)InSAR地表形變結(jié)果，結(jié)合99座小煤礦范圍、5大煤礦位置和野外地面調(diào)查資料，圈定了49個(gè)采空區(qū)分布范圍，如圖3所示。針對(duì)49個(gè)采空區(qū)獲取了采空區(qū)歷年累計(jì)區(qū)域平均形變量統(tǒng)計(jì)結(jié)果，其時(shí)序形變過程如圖4所示。

圖4 渾源縣49個(gè)形變區(qū)域中心2007-2013年間時(shí)序結(jié)果

從時(shí)序圖上可知各采空區(qū)7 a間地表高程持續(xù)下降，最大累計(jì)形變達(dá)-211.21 cm。其中15、17、18、19、20等區(qū)域平均形變量在2010年前后發(fā)生了劇烈變化，從圖3可知，5個(gè)區(qū)域位于黃花灘鄉(xiāng)-大磁窯鎮(zhèn)-青瓷窯鄉(xiāng)交界處，屬渾源縣東部沉降群。上述5個(gè)區(qū)域自2010年后煤礦開采規(guī)模與速度均有所加劇，其中20號(hào)區(qū)域平均形變量為當(dāng)年渾源縣最為嚴(yán)重區(qū)域。2012年地面沉降進(jìn)一步加劇，區(qū)域平均沉降量超過-50 cm的有3個(gè)區(qū)域，分別為20號(hào)的-62.01 cm、36號(hào)的-60.34 cm、41號(hào)的-56.32 cm。2013年平均沉降量超過50 cm的區(qū)域?yàn)?3號(hào)和20號(hào)，其值分別為-53.33 cm、-52.45 cm。

圖3 渾源縣年均沉陷速率及沉陷區(qū)分布范圍

為驗(yàn)證InSAR監(jiān)測(cè)結(jié)果的可靠性，結(jié)合實(shí)地調(diào)查結(jié)果與煤礦開采記錄資料等輔助數(shù)據(jù)對(duì)結(jié)果進(jìn)行驗(yàn)證。實(shí)地調(diào)查了沙圪坨鎮(zhèn)西坪村、大仁莊鄉(xiāng)吳家圪坨村、蔡村鎮(zhèn)白道村、青磁窯鄉(xiāng)青磁窯村，如圖5所示，結(jié)合形變時(shí)間序列結(jié)果可知7 a間累計(jì)形變量分別達(dá)-90.0 cm、-141.22 cm、-125.51 cm、-75.77 cm。受煤礦開采影響，煤礦區(qū)域周圍居民遭遇不同程度的地面沉降危害，村莊房屋建筑出現(xiàn)不同程度受損情況。

圖5 部分村莊房屋損毀情況

3.2 采空區(qū)形變分布特征

從圖3可知，渾源縣7 a間在東部、北部、西北部分別形成了三大沉降漏斗群，其中東部范圍最廣，包括大磁窯鎮(zhèn)-黃花灘鄉(xiāng)-青瓷窯鄉(xiāng)北部-千佛嶺鄉(xiāng)北部-大仁莊鄉(xiāng)中部等廣闊區(qū)域，該區(qū)域分布著果子園、古瓷窯、白莊等眾多煤礦。對(duì)各區(qū)域時(shí)序形變結(jié)果中歷年不同形變量比例進(jìn)行統(tǒng)計(jì)得到圖6，從圖中可知，隨著時(shí)間變化，渾源縣49個(gè)區(qū)域中歷年形變量大的比重不斷增加。2007—2013年間，沉降量＞-10 cm的區(qū)域從75.5%下降至2%，年沉降量超過-50 cm區(qū)域從無到有，其中2012年共3個(gè)，占6.1%；2013年2個(gè)，占4.1%。但從2007年看形變量較2008年大，分析其原因可能是當(dāng)年數(shù)據(jù)從2007年6月起，僅半年數(shù)據(jù)，獲取的形變量受季節(jié)性變化影響較大，難以準(zhǔn)確反映出全年的形變規(guī)律；2013年最大沉降值較2012年更小，可能有兩個(gè)原因，一是2013年采用的數(shù)據(jù)不足半年，由于數(shù)據(jù)量不足及季節(jié)性影響可能無法完整反映出渾源縣2013年形變結(jié)果；另一個(gè)原因是2012年政府對(duì)煤礦開采加以控制，地面形變情況在2013年得到改善，針對(duì)原因不確定性，后續(xù)采用完整數(shù)據(jù)驗(yàn)證具體影響因素。

圖6 渾源縣歷年不同形變量所占比例

根據(jù)時(shí)序InSAR結(jié)果可知渾源縣地表形變?cè)诳臻g上呈散點(diǎn)式漏斗形狀，面積小而多，均呈持續(xù)下沉狀態(tài)。結(jié)合時(shí)序InSAR結(jié)果獲取了渾源縣采空區(qū)中心隨時(shí)間變化的空間演變結(jié)果（圖7）。從圖上可知渾源縣北部與西部在2007—2013年間中心形變發(fā)生明顯形變，其中部分采空區(qū)中心形變具有明顯階段性特征，以北部41號(hào)采空區(qū)為例，該區(qū)形變以2011年為分界線；2007—2011年間呈穩(wěn)步線性趨勢(shì)變化，發(fā)生突然加速。渾源縣東部在2007—2010年間采空區(qū)中心形變較為平穩(wěn)，無明顯突變，但2011—2013年該區(qū)域地表形變異常，通過與時(shí)序結(jié)果對(duì)比發(fā)現(xiàn)在2010年后地表發(fā)生明顯變化的15、17、18、19、20號(hào)采空區(qū)均位于東部漏斗群，7 a間最大累計(jì)形變超過-2.3 m。

圖7 渾源縣采空區(qū)中心歷年累計(jì)變化量

2007—2013年間，渾源縣三大沉降漏斗群均出現(xiàn)不同程度形變。在時(shí)間緯度，49個(gè)采空區(qū)均呈持續(xù)的、大量級(jí)的形變，最大累計(jì)沉降超過-2.3 m；在空間緯度，采空區(qū)分布于11個(gè)鄉(xiāng)鎮(zhèn)，其中以東部的黃花灘鄉(xiāng)、大磁窯鎮(zhèn)、青瓷窯鄉(xiāng)、千佛嶺鄉(xiāng)、大仁莊鄉(xiāng)最為嚴(yán)重。7 a間渾源縣地面沉降無論從沉降量級(jí)還是沉降范圍均發(fā)生了巨大的變化。

3.3 討論

為分析時(shí)序InSAR監(jiān)測(cè)結(jié)果可靠性，本文搜集了渾源縣采煤沉陷區(qū)村莊治理、大同市原煤產(chǎn)量與渾源縣煤炭銷量[24]等資料。首先對(duì)涉及采煤沉陷區(qū)治理的村莊分析，發(fā)現(xiàn)村莊主要分布在大仁莊鄉(xiāng)、大磁窯鎮(zhèn)、東坊城鄉(xiāng)、黃花灘鄉(xiāng)、青瓷窯鄉(xiāng)、千佛嶺鄉(xiāng)、蔡村鎮(zhèn)、沙圪坨鎮(zhèn)等鄉(xiāng)鎮(zhèn)，分布于渾源縣東部、北部、西北部，其空間分布上與根據(jù)時(shí)序InSAR結(jié)果圈定的49個(gè)采空區(qū)空間位置分布相一致。說明煤礦開采是渾源縣出現(xiàn)形變漏斗的主要原因。

此外，通過對(duì)大同市原煤產(chǎn)量與渾源縣煤炭銷量進(jìn)行分析發(fā)現(xiàn)，自2008年起，當(dāng)?shù)孛旱V開采呈現(xiàn)爆炸式增長(zhǎng)。大同市原煤銷量從2007年的7 394萬t增長(zhǎng)至2010年的9 403萬t，2015年達(dá)122 59萬t；渾源縣煤炭銷量從2008年的15.46萬t，增長(zhǎng)至2011年1 060萬t，年均增長(zhǎng)幅度達(dá)309%。當(dāng)?shù)孛旱V的產(chǎn)量與銷量變化情況與部分采空區(qū)形變量在2010年前后發(fā)生突變趨勢(shì)相吻合（圖8）。經(jīng)調(diào)查發(fā)現(xiàn)，自2008年起，煤炭?jī)r(jià)格飛速上漲，煤礦開采日益嚴(yán)重，尤其開采難度較低的露天煤礦已初具規(guī)模，使得多數(shù)區(qū)域年形變量開始呈指數(shù)加速增長(zhǎng)狀態(tài)。

圖8 大同市原煤產(chǎn)量與渾源縣煤炭銷量及沉降中心年均形變時(shí)序圖

4 結(jié)論

本文分析了2007年6月—2011年1月36景ALOS PALSAR影像和2012年1月—2013年4月10景RADARSAT-2影像數(shù)據(jù)，獲取了渾源地區(qū)地表形變速率結(jié)果和時(shí)序形變結(jié)果。通過對(duì)結(jié)果分析可知：

1）根據(jù)時(shí)序InSAR地表形變結(jié)果，結(jié)合99座小煤礦范圍、5大煤礦位置和野外地面調(diào)查資料，圈定49個(gè)采空區(qū)范圍及各采空區(qū)歷年累計(jì)區(qū)域平均形變量統(tǒng)計(jì)結(jié)果。7 a間渾源縣地表呈持續(xù)虧損狀態(tài)，累計(jì)形變量不斷加劇，且在部分區(qū)域2010年后形變呈突然加劇狀態(tài)。

2）2007—2013年間，渾源縣出現(xiàn)嚴(yán)重的地表形變，在東部、北部、西北部分別形成了三大沉降漏斗群，其中范圍最廣的東部區(qū)域分布著果子園、古瓷窯、白莊等眾多煤礦。隨著時(shí)間推移，渾源縣49個(gè)區(qū)域中歷年形變量大的比重不斷增加，沉降量級(jí)與沉降范圍均呈惡化趨勢(shì)。

3）根據(jù)時(shí)序InSAR結(jié)果圈定的形變區(qū)空間分布與采煤沉陷區(qū)開采量一致性較高，表明煤礦開采是造成地表沉陷加劇的主要原因。