馬道鳴，魯?shù)顒偅煳幕?，張雪麗，蘇永江

(河南省航空物探遙感中心，河南 鄭州 453000)

0 前 言

礦區(qū)地表形變監(jiān)測(cè)一直是礦山安全生產(chǎn)的基礎(chǔ)保障，也是嚴(yán)重制約礦山升級(jí)發(fā)展的重要因素，是在人類采礦活動(dòng)和自然地質(zhì)條件變化的雙重作用下的巖石應(yīng)力變化，進(jìn)而使得礦區(qū)深部采空區(qū)上部區(qū)域及鄰近區(qū)域的巖礦石發(fā)生不均一的形變，甚至塌陷等，最終威脅著資源的安全開(kāi)發(fā)活動(dòng)，同時(shí)造成大量的礦石損耗等[1]。因此，提高礦區(qū)地表形變監(jiān)測(cè)效果不僅是有效防治礦山災(zāi)害事故的主要途徑，也是提高資源利用效率的基礎(chǔ)。傳統(tǒng)的GPS、全站儀等監(jiān)測(cè)僅可獲得監(jiān)測(cè)點(diǎn)在某個(gè)時(shí)間點(diǎn)的沉降量，而無(wú)法獲得監(jiān)測(cè)區(qū)域整體的形變特征，這就限制了對(duì)監(jiān)測(cè)區(qū)域整體形變規(guī)律的研究，進(jìn)而造成形變特征分析出現(xiàn)偏差的局面，所以選擇監(jiān)測(cè)方法是提高監(jiān)測(cè)質(zhì)量的前提[2]。合成孔徑雷達(dá)差分干涉測(cè)量技術(shù)(D-InSAR)是在合成孔徑雷達(dá)干涉測(cè)量技術(shù)(InSAR)上發(fā)展起來(lái)的新型監(jiān)測(cè)技術(shù)，該監(jiān)測(cè)方法充分利用了高精度DEM成果及干涉圖等基礎(chǔ)土建，能夠更精確地反映礦區(qū)地表微小的形變量，故常用于較精細(xì)及緩慢形變監(jiān)測(cè)領(lǐng)域中，如滑坡形變監(jiān)測(cè)[3-4]、凍土區(qū)公路地表形變監(jiān)測(cè)[5-6]、城市地表形變監(jiān)測(cè)[7-8]、鐵路地表形變監(jiān)測(cè)[9]、火山附近地表形變監(jiān)測(cè)[10]、地震監(jiān)測(cè)領(lǐng)域[11-13]以及礦區(qū)地表形變監(jiān)測(cè)等領(lǐng)域，在上述形變監(jiān)測(cè)中取得了良好的應(yīng)用效果。鑒于此，本文以金銅礦山的多幅ENVISAT ASAR數(shù)據(jù)為基礎(chǔ)，試圖采用D-InSAR技術(shù)分析該礦區(qū)地表形變規(guī)律，綜合研究礦區(qū)地表形變發(fā)展規(guī)律，為礦山資源的安全開(kāi)發(fā)利用提供參考依據(jù)。

1 D-InSAR的技術(shù)流程

礦區(qū)地表形變是人類采礦活動(dòng)及地質(zhì)條件變化的雙重作用下誘發(fā)引起的地殼淺表層高程變化以及位移的一種表現(xiàn)形式，D-InSAR技術(shù)就是利用了不同時(shí)段高分辨率影像數(shù)據(jù)中的微小形變特征而獲取地表形變信息的。因此，D-InSAR技術(shù)在礦區(qū)地表形變監(jiān)測(cè)中的應(yīng)用就是利用不同時(shí)段獲得的ENVISAT ASAR影像數(shù)據(jù)，通過(guò)影像數(shù)據(jù)的進(jìn)一步解纏處理等，就可獲得該時(shí)段內(nèi)礦區(qū)地表的微小形變特征。影像數(shù)據(jù)的處理方法有二軌法、三軌法和四軌法等幾種類型，其中后2種數(shù)據(jù)處理方法中引入了大氣誤差，容易造成監(jiān)測(cè)精度降低。二軌法是最為常見(jiàn)的和常用的數(shù)據(jù)處理方法，本文也采用二軌法進(jìn)行處理跨越形變期不同時(shí)段的ENVISAT ASAR影像數(shù)據(jù)，進(jìn)而借助高精度DEM模型進(jìn)行反演，將地形相位等剔除，就可獲得監(jiān)測(cè)區(qū)域內(nèi)任何一點(diǎn)的形變量數(shù)據(jù)。綜上所述，本文采取的二軌法在數(shù)據(jù)處理方面的技術(shù)流程見(jiàn)圖1。

圖1 D-InSAR技術(shù)的二軌法數(shù)據(jù)處理流程

2 礦區(qū)基本概況

研究區(qū)為一金銅多金屬礦床，礦區(qū)面積為2.5 km2，礦區(qū)邊界為一不規(guī)則的多邊形。礦區(qū)地形高差變化較大，屬于高山峽谷地貌與丘陵山區(qū)地貌的結(jié)合部位，導(dǎo)致地形地貌較為復(fù)雜。礦床成因類型為矽卡巖型金銅多金屬礦床，礦山已有的開(kāi)采方式為井下硐采，且局域北西區(qū)域礦體埋深較淺，南東區(qū)域礦體埋深較深。該礦山為一開(kāi)采礦山，現(xiàn)狀已有6處采空區(qū)，采空區(qū)上方區(qū)域的地表已出現(xiàn)不同程度的地表沉降問(wèn)題，甚至在某些采空區(qū)上方已出現(xiàn)細(xì)微的地裂縫。因此，為了了解該礦區(qū)地表形變規(guī)律以及指導(dǎo)資源開(kāi)發(fā)，本文以ENVISAT ASAR影像數(shù)據(jù)為基礎(chǔ)，通過(guò)D-InSAR技術(shù)研究礦區(qū)地表形變發(fā)展規(guī)律，為礦山安全生產(chǎn)提供參考依據(jù)。

3 D-InSAR技術(shù)在礦區(qū)地表形變監(jiān)測(cè)中的應(yīng)用

3.1 配準(zhǔn)影像數(shù)據(jù)及估算基線

前人的大量研究資料[3-9]顯示，二軌法的影像數(shù)據(jù)配準(zhǔn)方法主要包括2種類型：一是以精密星歷數(shù)據(jù)以及相干系數(shù)配準(zhǔn)影像數(shù)據(jù)；二是以地面控制點(diǎn)為基準(zhǔn)點(diǎn)進(jìn)行影像數(shù)據(jù)的配準(zhǔn)。本次所使用的影像數(shù)據(jù)來(lái)源于COSMO-SkyMed高分辨率雷達(dá)衛(wèi)星，相對(duì)而言該衛(wèi)星的影像數(shù)據(jù)精度較低，故本次影像數(shù)據(jù)的配準(zhǔn)方式選擇以精密星歷數(shù)據(jù)以及相干系數(shù)配準(zhǔn)，可以通過(guò)該種配準(zhǔn)方式有效地提高監(jiān)測(cè)精度。在使用第一種配準(zhǔn)方式過(guò)程中一般采用DORIS軌道數(shù)據(jù)對(duì)相應(yīng)的ENVISAT ASAR影像數(shù)據(jù)進(jìn)行基線的估算處理，在完成上述操作的基礎(chǔ)上對(duì)礦區(qū)內(nèi)地形地貌相對(duì)較為平緩的區(qū)域以條紋率進(jìn)行精估計(jì)[2]，就可獲得精度較高的配準(zhǔn)后的影像數(shù)據(jù)。影像數(shù)據(jù)精配準(zhǔn)后結(jié)合礦區(qū)高精度的DEM數(shù)據(jù)對(duì)影像數(shù)據(jù)進(jìn)行相位模型處理，將影像數(shù)據(jù)中包含的地形信息、形變信息等轉(zhuǎn)換為相位信息，為后期進(jìn)一步展開(kāi)影像數(shù)據(jù)的相位差分干涉處理以及去地形相位解纏等處理做準(zhǔn)備。

3.2 差分處理及相位濾波

在完成ENVISAT ASAR影像數(shù)據(jù)配準(zhǔn)和基線估算后形成了包含地形信息和形變信息的相位信息，即相位干涉圖，此時(shí)的相位干涉圖中包含了地形信息和形變信息。因此，需要將相位干涉圖中的地形信息剔除，就可獲得去地形干涉相位，該過(guò)程就是干涉影像數(shù)據(jù)的差分處理[13]。去地形干涉相位中不僅包含了大量的微小形變信息，而且包含了其他的噪聲信息，如噪聲相位、非線性相位等，此時(shí)需要進(jìn)行噪聲相位和非線性相位的剔除處理，即差分干涉相位的濾波處理。

3.3 相位解纏及地理編碼

完成差分干涉相位濾波處理后可進(jìn)行相位的解纏過(guò)程，該過(guò)程也是獲取礦區(qū)地表形變信息的基礎(chǔ)。本文相位的解纏采用最小費(fèi)用流法，該方法能夠降低影像數(shù)據(jù)中微小形變模糊的問(wèn)題[5]，進(jìn)而獲得與線性變化地形相互一一對(duì)應(yīng)的相位數(shù)據(jù)，在GAMMA軟件中數(shù)據(jù)配準(zhǔn)換算關(guān)系等解算出不同解纏點(diǎn)的形變量信息[12]。采集形變信息的數(shù)據(jù)并編輯地理編碼等信息，并將其統(tǒng)一至同一坐標(biāo)系統(tǒng)下，進(jìn)而通過(guò)圖件整飾處理等就可獲得礦區(qū)地表形變圖等成果性圖件，能夠?qū)⒌V區(qū)地表形變直觀地表達(dá)出來(lái)(見(jiàn)圖2)。

圖2 研究區(qū)地表沉降形變分布

4 礦區(qū)地表沉降形變分析

4.1 沉降形變監(jiān)測(cè)結(jié)果分析

在完成相位解纏后就可在MAPGIS軟件等中采集并編輯形變信息以及地理編碼等，再通過(guò)圖件整飾就可生成礦區(qū)地表形變分布圖(見(jiàn)圖2)，圖件中包含了礦山中的形變區(qū)域、分布位置、形變范圍、沉降量等信息。由圖2可知：研究區(qū)地表沉降范圍共識(shí)別出6個(gè)區(qū)域，地表形變范圍與深部采空區(qū)具有良好的空間吻合關(guān)系；同時(shí)，礦區(qū)地表形變面積以及沉降量大小與深部資源的開(kāi)發(fā)利用關(guān)系密切，其中采空區(qū)形成時(shí)間長(zhǎng)且面積較大，地表形變面積較大且沉降量較大，如圖2南東側(cè)和北西側(cè)的采空區(qū)所對(duì)應(yīng)的地表形變區(qū)，即在垂向上地表形變分布區(qū)域深部采空區(qū)具有良好的對(duì)應(yīng)關(guān)系；地表形變分布區(qū)的空間展布形態(tài)與深部采空區(qū)的展布方向和形態(tài)基本一致，形變分布區(qū)的長(zhǎng)軸方向與采空區(qū)的長(zhǎng)軸方向一致。由礦山沉降量形變圖可知：礦區(qū)內(nèi)地表形變可劃分為3個(gè)等級(jí)(-25～0 mm、-50～-25 mm和-100～-50 mm)，其中礦區(qū)地表最大沉降量位于礦區(qū)南東側(cè)的采空區(qū)上部，即礦區(qū)主礦體采空區(qū)上部，最大沉降量可達(dá)92.5 mm；礦區(qū)內(nèi)最小沉降形變區(qū)分布在礦區(qū)南西側(cè)，其沉降量一般小于25 mm；礦區(qū)內(nèi)總形變區(qū)域可達(dá)0.78 km2，其中形變量大于-50 mm的面積約為0.07 km2，形變量介于-50～-25 mm的面積為0.31 km2，形變量介于-25～0 mm的面積為0.40 km2。

為了解礦區(qū)地表形變?cè)诖瓜蜃呦虻淖兓?guī)律，在形變最發(fā)育的范圍內(nèi)對(duì)垂直形變主體走向進(jìn)行了剖面研究(見(jiàn)圖3)。由圖3可知：研究區(qū)礦區(qū)地表形變?cè)谄拭嫔暇哂写笾聻椤癡”形變化的規(guī)律，其最大沉降量位于采空區(qū)中心部位靠近南東一側(cè)，與采空區(qū)的分布范圍基本吻合；同時(shí)，沉降量從采空區(qū)中心部位向兩側(cè)具有大致對(duì)稱分布的特征，且形變量向遠(yuǎn)離采空區(qū)范圍逐漸降低。

圖3 研究區(qū)A-B剖面沉降特征

4.2 精度對(duì)比分析

為了對(duì)比分析D-InSAR技術(shù)在礦區(qū)地表形變監(jiān)測(cè)中的精度問(wèn)題，本文利用礦區(qū)形變區(qū)內(nèi)已有的控制點(diǎn)坐標(biāo)信息等，使用傳統(tǒng)的全站儀測(cè)量方法對(duì)D-InSAR技術(shù)的精度進(jìn)行驗(yàn)證。本文累計(jì)收集了10個(gè)控制點(diǎn)的坐標(biāo)信息，通過(guò)開(kāi)采前后的三維坐標(biāo)對(duì)比，獲得了控制點(diǎn)的形變量數(shù)據(jù)，進(jìn)而與D-InSAR形成的沉降形變圖中的沉降量對(duì)比分析。本次為了獲得更加準(zhǔn)確的D-InSAR沉降形變數(shù)據(jù)，將控制點(diǎn)周邊的約30個(gè)高相干點(diǎn)數(shù)據(jù)進(jìn)行加權(quán)平均，進(jìn)而作為該控制點(diǎn)的D-InSAR方法的數(shù)據(jù)，最終的數(shù)據(jù)統(tǒng)計(jì)見(jiàn)表1。由表1可知：本文所獲取的D-InSAR加權(quán)平均形變量數(shù)據(jù)與全站儀測(cè)量所獲形變量的變化趨勢(shì)是一致的，二者的誤差百分比均小于5%，說(shuō)明本文所獲礦區(qū)地表形變量是可靠的，其精度能夠滿足礦山形變監(jiān)測(cè)應(yīng)用。因此，可以將合成孔徑雷達(dá)差分干涉測(cè)量技術(shù)(D-InSAR)應(yīng)用至礦區(qū)地表形變監(jiān)測(cè)應(yīng)用中，其監(jiān)測(cè)效果良好。

表1 D-InSAR技術(shù)監(jiān)測(cè)結(jié)果與全站儀監(jiān)測(cè)結(jié)果對(duì)比

5 結(jié) 論

①D-InSAR技術(shù)在礦區(qū)地表形變監(jiān)測(cè)中具有良好的應(yīng)用優(yōu)勢(shì)，該技術(shù)方法能夠獲取礦區(qū)微小的形變信息，對(duì)全面了解礦區(qū)地表形變特征以及綜合研究形變規(guī)律，進(jìn)而預(yù)測(cè)形變發(fā)展趨勢(shì)等意義重大。

②研究區(qū)地表形變范圍與礦區(qū)深部采空區(qū)具有良好的空間吻合關(guān)系，形變長(zhǎng)軸與礦區(qū)采空區(qū)的長(zhǎng)軸方向基本一致，沉降量從采空區(qū)中心向兩側(cè)呈大致對(duì)稱狀分布，且向遠(yuǎn)離采空區(qū)中心向兩側(cè)形變量逐漸減小，總體上具有“V”形的變化特征。

③礦區(qū)地表最大沉降量可達(dá)92.5 mm，最小沉降形變區(qū)分布在礦區(qū)南西側(cè)，其沉降量一般小于25 mm；礦區(qū)內(nèi)總形變區(qū)域可達(dá)0.78 km2，其中形變量大于-50 mm的面積約為0.07 km2，形變量介于-50～-25 mm的面積為0.31 km2，形變量介于-25～0 mm的面積為0.40 km2。