基于D-InSAR 技術(shù)的煤礦區(qū)地面沉降監(jiān)測(cè)研究*

宋繼德，邸志眾

(1．陜西省地質(zhì)礦產(chǎn)勘查開發(fā)局 第二綜合物探大隊(duì)，陜西 西安 710016;2．北京中勘邁普科技有限公司，北京 100086)

0 引言

InSAR(SAR Interferometry 或Interferometric SAR)，即合成孔徑雷達(dá)干涉測(cè)量，是指利用SAR 數(shù)據(jù)中的相位信息進(jìn)行干涉測(cè)量處理，結(jié)合雷達(dá)參數(shù)和衛(wèi)星位置信息反演地表三維及其變化信息的遙感手段［1，2］。D-InSAR 即差分合成孔徑雷達(dá)干涉測(cè)量，是指對(duì)干涉相位進(jìn)行差分處理提取變化信息的干涉測(cè)量手段，是InSAR 技術(shù)的延伸。InSAR 技術(shù)已被廣泛用于大區(qū)域地形測(cè)圖和地表信息提取。常見的InSAR 方式有順軌(Along track)，交軌(Across track)和重復(fù)軌(Repeat pass)。由于地表形變監(jiān)測(cè)是測(cè)量一定時(shí)間間隔內(nèi)的地表形變量，因而重復(fù)軌是目前應(yīng)用的主要方式，具體是指利用同一地區(qū)不同時(shí)刻獲取的SAR 數(shù)據(jù)進(jìn)行差分干涉處理以提取地表形變信息。由于InSAR具有全天候、全天時(shí)、覆蓋面廣和高精度獲取地表信息的能力，因而在地球科學(xué)的諸多領(lǐng)域如地震、火山、礦區(qū)塌陷、地面沉降、滑坡、冰川等災(zāi)害中逐漸得到了廣泛的應(yīng)用［1，5］。

本文以山東省重點(diǎn)地區(qū)礦山開采的遙感動(dòng)態(tài)監(jiān)測(cè)為研究?jī)?nèi)容，展開濟(jì)寧、兗州地區(qū)地面沉降InSAR 監(jiān)測(cè)試驗(yàn)研究，圍繞濟(jì)寧和兗州地區(qū)煤礦開采、地下水開采等多種因素作用下地面沉降(陷)的分布與發(fā)展?fàn)顩r，介紹以D-InSAR 技術(shù)為主要手段的地面沉降調(diào)查與監(jiān)測(cè)工作的研究情況與成果。

1 差分合成孔徑雷達(dá)干涉測(cè)量原理與數(shù)據(jù)處理方法

1．1 基本原理

重復(fù)軌干涉測(cè)量的幾何原理［2］，如圖1 所示。圖中O1、O2分別為兩次獲取SAR 信號(hào)的天線位置，B 為兩天線間空間距離，稱為空間基線。

圖1 重復(fù)軌D-InSAR 幾何示意圖Fig．1 Geometric schematic diagram of repeat pass D-InSAR

設(shè)基線與水平方向的夾角為α，H 為平臺(tái)高度，地面點(diǎn)P到兩天線間的距離分別為R1和R2，θ 為第一副天線的側(cè)視角，P 點(diǎn)高程為h。將基線沿雷達(dá)視線向分解，垂直于視線向分量為B⊥，稱為垂直基線。雷達(dá)對(duì)地物觀測(cè)時(shí)，接收到地物的后向散射回波信號(hào)可以表示為:

兩次觀測(cè)所獲取的SAR 影像分別稱為主影像(Master)和輔影像(Slave)，接收到P 點(diǎn)的雷達(dá)信號(hào)分別為:

雷達(dá)后向散射回波信號(hào)中的相位信息包含兩部分信息，一是與往返路徑相關(guān)的相位，另一部分是目標(biāo)后向散射相位，可以表示為:

兩次獲取的相位分別為:

式中:ΔR 為路程差，與雷達(dá)與目標(biāo)的距離相關(guān)。若兩次雷達(dá)獲取的相位中的散射相位保持穩(wěn)定，即:

則經(jīng)過共軛相乘處理后得到的干涉相位為:

干涉相位φint的貢獻(xiàn)值(Contributions)由以下5 個(gè)因素構(gòu)成:平地相位δφflat;地形相位δφtopo;兩次觀測(cè)期間觀測(cè)目標(biāo)沿雷達(dá)視線向(LOS)移動(dòng)引起的相位變化δφmov;電磁波傳輸過程中大氣波動(dòng)引起的相位延遲δφatm;其他噪聲相位δφnoise。則干涉相位可以表示為:

1．2 數(shù)據(jù)處理方法

根據(jù)式(8)，要得到反映地形信息的干涉相位δφtopo，應(yīng)當(dāng)首先從干涉相位和中消除平地相位的影響，其次從中減去其他各項(xiàng)得到地形相位，再進(jìn)行相位解纏和相位到高程的轉(zhuǎn)換，將生成的高程信息轉(zhuǎn)換到特定的坐標(biāo)系統(tǒng)下，生成DEM 數(shù)據(jù)［1，6］。由于地形形變相位、熱噪聲相位、大氣波動(dòng)相位以及軌道偏差引起的相位波動(dòng)可以通過合理的處理方法來消除或者降低。對(duì)于ERS-1/2，其特殊的“串行模式”(Tandem Mode)SAR 數(shù)據(jù)是生成DEM 的最佳選擇［6］。因而，利用InSAR 獲取DEM 數(shù)據(jù)的處理流程，見圖2。

圖2 InSAR 數(shù)據(jù)處理流程Fig．2 Data processing flow of InSAR

1．3 D-InSAR 地表形變監(jiān)測(cè)技術(shù)及應(yīng)用

隨著D-InSAR 技術(shù)研究的不斷深入，利用其進(jìn)行地表形變信息獲取的數(shù)據(jù)處理方法也在不斷進(jìn)步。針對(duì)不同形變體的變形特征，產(chǎn)生了不同的監(jiān)測(cè)方法。按照監(jiān)測(cè)目標(biāo)空間分布的不同和數(shù)據(jù)處理方式的差別，可將InSAR 信息提取方式分為兩類:基于點(diǎn)目標(biāo)的D-InSAR 方法和基于連續(xù)目標(biāo)的D-InSAR 監(jiān)測(cè)方法。

基于單個(gè)點(diǎn)目標(biāo)的方法以角反射器方法為代表，主要針對(duì)的是形變區(qū)域無相干目標(biāo)或存在極少數(shù)相干目標(biāo)的情況，監(jiān)測(cè)區(qū)域較小，適合于局部變形監(jiān)測(cè)，如滑坡體、危巖活動(dòng)監(jiān)測(cè)等［1，7，8］。

連續(xù)目標(biāo)方法指的是在空間大面積范圍內(nèi)，長(zhǎng)時(shí)間條件下依然存在較多相干目標(biāo)的情況［1］。常規(guī)差分干涉測(cè)量方法難以解決多時(shí)相動(dòng)態(tài)監(jiān)測(cè)問題。緩慢連續(xù)變形體的移動(dòng)監(jiān)測(cè)依賴于長(zhǎng)時(shí)間連續(xù)監(jiān)測(cè)，對(duì)于重復(fù)軌監(jiān)測(cè)，長(zhǎng)時(shí)間周期會(huì)引起失相干，這種情況下，對(duì)于大區(qū)域連續(xù)形變體的移動(dòng)監(jiān)測(cè)則比較困難。由于監(jiān)測(cè)對(duì)象多是與人們生產(chǎn)和生活緊密相關(guān)的城市或城鎮(zhèn)，這些區(qū)域在長(zhǎng)時(shí)間間隔下依然保持了很好的相干性，因而，可以通過多時(shí)相分析，獲取這些區(qū)域的地表形變信息。由此，則產(chǎn)生了累積紋圖法(Stacking Interferograms)方法和永久性散射體(Permanent Scatterers)技術(shù)。這兩種方法通過長(zhǎng)時(shí)間序列分析，彌補(bǔ)了常規(guī)重復(fù)軌干涉測(cè)量的不足，提高了時(shí)間采樣率，能更準(zhǔn)確的獲取形變隨時(shí)間變化而產(chǎn)生的形變?？梢詫⑿巫兲卣鞣譃榫€性和非線性分別進(jìn)行處理［3，4，8］。但由于永久性散射體技術(shù)對(duì)SAR 數(shù)據(jù)數(shù)目要求較高，一般情況下難以完成大面積形變。

2 研究區(qū)域地表形變InSAR 監(jiān)測(cè)的技術(shù)方法

山東省兗州地區(qū)為典型的煤礦開采區(qū)，開采沉陷誘發(fā)的礦區(qū)環(huán)境破壞不僅出現(xiàn)在井田開采上覆地區(qū)，同時(shí)也影響到了城市和城鎮(zhèn)的地面穩(wěn)定性，因而礦山沉陷InSAR 監(jiān)測(cè)的首要任務(wù)是明確2004 ～2010 年開采沉陷的變化特征，進(jìn)而選擇合理的數(shù)據(jù)處理方法和相關(guān)技術(shù)參數(shù)。

針對(duì)緩慢變形和快速沉降兩種不同的形變類型，采用常規(guī)D-InSAR 技術(shù)與永久散射體干涉測(cè)量相結(jié)合的工作方法、利用永久性散射體干涉測(cè)量(PSI)技術(shù)相較其他技術(shù)，在剔除地形誤差和大氣波動(dòng)對(duì)相位信號(hào)的影響方面的優(yōu)勢(shì)，且目前可獲得的SAR 數(shù)據(jù)在精度和分辨率方面也足以滿足PSI 對(duì)數(shù)據(jù)的要求。

2．1 永久散射體干涉測(cè)量(PSInSAR)技術(shù)

永久散射體干涉測(cè)量(Permanent Scatterers Interferometry SAR，PSInSAR)的核心思想是對(duì)永久散射體干涉相位進(jìn)行時(shí)間序列分析，根據(jù)各相位分量的時(shí)空特征，估算大氣波動(dòng)，DEM 誤差以及噪聲等，將其從差分干涉相位中逐個(gè)分離，最終獲取每個(gè)PS 的線性和非線形形變速率、大氣延遲量(Atmosphere Phase Screen)以及DEM 誤差。經(jīng)PS 方法處理，獲取的年度形變率的精度可以達(dá)到毫米級(jí)。其基本思想是基于大量的SAR 數(shù)據(jù)(20-30 景甚至更多)，從中篩選出具有穩(wěn)定散射特性的相干點(diǎn)目標(biāo)，構(gòu)成離散點(diǎn)觀測(cè)網(wǎng)絡(luò)(較之常規(guī)的變形監(jiān)測(cè)網(wǎng)密度更高)，通過分析PS 點(diǎn)目標(biāo)相位變化獲取地表形變狀況。由于將永久散射體作為觀測(cè)對(duì)象，降低了空間基線對(duì)相干性的影響，即使在臨界基線的條件下，仍然可以通過分析PS 差分干涉相位的變化反演形變信息、提取地表形變速率和累積形變量［8］(見圖3)。

PS 技術(shù)一般采用線性形變模型提取點(diǎn)目標(biāo)對(duì)應(yīng)的形變量，非線性模擬的處理過程復(fù)雜，而且非常耗時(shí)，限制了其用于大面積的形變測(cè)量，但隨著處理技術(shù)的進(jìn)步，處理時(shí)間將逐步縮小，處理的范圍也可再進(jìn)一步擴(kuò)大。

PSInSAR 對(duì)數(shù)據(jù)量要求較高，只有SAR 影像個(gè)數(shù)達(dá)到一定的程度才能篩選出在整個(gè)時(shí)間跨度內(nèi)具有穩(wěn)定信號(hào)的PS 點(diǎn)，對(duì)于相干目標(biāo)較多的地區(qū)，如城區(qū)等，由于地物在長(zhǎng)時(shí)間間隔下保持了較高的相干性，數(shù)據(jù)量較少的情況下也可進(jìn)行處理。一般情況下，至少必須滿足5 個(gè)/km2才能完成相位解纏。

圖3 PSInSAR 數(shù)據(jù)處理的基本方法Fig．3 The basic method of data processing of PSInSAR

PS 方法在處理過程中通過分析每個(gè)點(diǎn)目標(biāo)對(duì)應(yīng)的相位，將其中包含的因地表形變，高程以及大氣波動(dòng)引起的相位分量加以區(qū)分，最終得到各個(gè)點(diǎn)對(duì)應(yīng)的形變相位。通過分析PS 數(shù)據(jù)集，估算每個(gè)干涉紋圖中由大氣波動(dòng)引起的相位分量，此外可以確定出每一點(diǎn)精確的高程信息。完成上述過程后，在干涉處理時(shí)將二者對(duì)相位的貢獻(xiàn)量加以去除。經(jīng)PS 方法處理，獲取的年度形變率的精度可以達(dá)到±0．1 mm/a。

2．2 可獲取的主要雷達(dá)數(shù)據(jù)及參數(shù)分析

當(dāng)前用于InSAR 技術(shù)研究和應(yīng)用較多的是ERS -1/2、ENVISAT、RADARSAT -1/2 以及ALOS PALSAR 獲取的SAR 數(shù)據(jù)，其空間分辨多為20 ～30 m，可歸為中等分辨率雷達(dá)數(shù)據(jù)。與之相比較而言，近幾年發(fā)射的雷達(dá)衛(wèi)星以高分辨率為主，包括Cosmo- Skymed 及TerraSAR - X，其最高分辨率可以達(dá)到1 ～3 m。表1 所列為當(dāng)前幾種雷達(dá)傳感器所提供的SAR 數(shù)據(jù)。

表1 SAR 數(shù)據(jù)及相關(guān)參數(shù)Tab．1 SAR data and relevant parameters

從雷達(dá)波段而言，目前主要的中等分辨率數(shù)據(jù)為C 波段數(shù)據(jù)(ENVISAT)。由于這些衛(wèi)星具有較高的定位精度(可以達(dá)到5 cm)，空間基線解算結(jié)果可靠，數(shù)據(jù)相干性較好，因此，在多數(shù)地表形變監(jiān)測(cè)中均可達(dá)到較好的監(jiān)測(cè)結(jié)果。高分辨率雷達(dá)衛(wèi)星TerraSAR-X 重復(fù)周期為11 d，采用X 波段，對(duì)地表形變等變化信息更為敏感，可以應(yīng)用于形變速率極快的監(jiān)測(cè)，礦山開采沉陷等非常適用。Cosmo-Skymed 與TerraSAR -X 相似，也屬于高分辨率雷達(dá)數(shù)據(jù)，最高分辨率達(dá)1 m，掃描模式獲取數(shù)據(jù)為3 m地面分辨率。在空間基線不大的條件下，L 波段的SAR 數(shù)據(jù)能夠在很長(zhǎng)的時(shí)間內(nèi)保持較高的相干性，能滿足變化幅度較大的形變監(jiān)測(cè)。特別對(duì)于快速的地表移動(dòng)監(jiān)測(cè)，應(yīng)用L 波段數(shù)據(jù)也較為合適。同時(shí)，由于這些數(shù)據(jù)分辨率的提高，空間采樣單元的個(gè)數(shù)將會(huì)增大，能更為準(zhǔn)確的描述地表形狀。因此，在很大程度上滿足了從短時(shí)間，變化速率高以及長(zhǎng)時(shí)間，變化幅度大的形變監(jiān)測(cè)。

2．3 工作方法與技術(shù)

針對(duì)兗州地區(qū)地表形變調(diào)查和監(jiān)測(cè)的工作重點(diǎn)，主要目標(biāo)為查明工作區(qū)地面沉降(陷)的歷史狀況和現(xiàn)狀，調(diào)查兗州、濟(jì)寧和鄒城下屬的多個(gè)礦區(qū)開采沉陷的基本特征和分布范圍，開展兗州市地面沉降InSAR 監(jiān)測(cè)。因而，采用常規(guī)D-InSAR 技術(shù)與永久散射體干涉測(cè)量相結(jié)合的工作方法進(jìn)行。

2．3．1 工作區(qū)SAR 數(shù)據(jù)獲取

采用的雷達(dá)數(shù)據(jù)為歐空局自2004 年11 月至2010 年10 月間接收的ENVISAT SAR 數(shù)據(jù)，其有效工作波段是C 波段，數(shù)據(jù)獲取周期為35 d，空間分辨率為20 m，是目前常用且應(yīng)用效果顯著的SAR 數(shù)據(jù)。獲取ASAR 數(shù)據(jù)共計(jì)24 景，如表2 所示。覆蓋范圍，見圖4。

表2 ENVISAT 衛(wèi)星數(shù)據(jù)列表Tab．2 List of ENVISAT satellite data

圖4 ENVISAT 衛(wèi)星數(shù)據(jù)覆蓋范圍Fig．4 Coverage of ENVISAT satellite data

區(qū)域性地面沉降監(jiān)測(cè)實(shí)際使用的數(shù)據(jù)從2004 年12 月到2010 年10 月計(jì)20 景，根據(jù)區(qū)域地面沉降長(zhǎng)時(shí)序監(jiān)測(cè)需求形成干涉像對(duì)。

2．3．2 數(shù)據(jù)處理和分析

數(shù)據(jù)處理和分析系統(tǒng)以專業(yè)InSAR 數(shù)據(jù)處理軟件，如GAMMA、Doris 和SARSCAPE 為主，選擇利用ArcGIS 平臺(tái)為數(shù)據(jù)后處理系統(tǒng)，針對(duì)工作區(qū)環(huán)境特點(diǎn)及礦區(qū)開采沉陷關(guān)鍵In-SAR 處理分析技術(shù)，分別利用常規(guī)D-InSAR 技術(shù)和永久散射體InSAR 技術(shù)分別針對(duì)快速沉陷和緩慢沉降開展監(jiān)測(cè)。InSAR 數(shù)據(jù)處理過程主要包括最優(yōu)干涉數(shù)據(jù)選擇、干涉像對(duì)的基線估算、平地效應(yīng)去除、相干性分析、大氣效應(yīng)分析、數(shù)據(jù)配準(zhǔn)、離散點(diǎn)相位解纏，時(shí)間序列分析方法等。

2．3．3 地表形變信息提取

采用永久散射體InSAR 時(shí)間序列分析技術(shù)提取地表形變信息，并結(jié)合若干實(shí)地觀測(cè)結(jié)果獲取不同時(shí)間間隔的地面沉降速率圖。針對(duì)礦山開采沉陷的發(fā)生特征，利用常規(guī)D-InSAR 技術(shù)開展數(shù)據(jù)處理，城市地表緩慢變化和活動(dòng)斷裂移動(dòng)監(jiān)測(cè)利用PSInSAR 技術(shù)開展數(shù)據(jù)處理與分析。

2．3．4 監(jiān)測(cè)結(jié)果后處理

在多期地面沉降信息分析與綜合的基礎(chǔ)上編制工作區(qū)地面沉降(陷)速率圖和累積沉降量圖。

3 區(qū)域地面沉降與開采沉陷InSAR 監(jiān)測(cè)與分析

3．1 煤礦區(qū)開采沉陷歷史狀況InSAR 動(dòng)態(tài)監(jiān)測(cè)

考慮到礦山塌陷的大變形過程對(duì)相干性的影響，以及季節(jié)性植被覆蓋變化，選擇冬季干涉像對(duì)，且間隔時(shí)間較短的數(shù)據(jù)進(jìn)行D-InSAR 數(shù)據(jù)處理，共獲取8 個(gè)能清晰反應(yīng)礦山沉陷的干涉紋圖，用于反映工作區(qū)礦區(qū)開采沉陷2011 年之前歷史狀況的動(dòng)態(tài)過程。

全區(qū)總的開采沉陷區(qū)主要分布于兗州、曲阜和鄒城沿市泗河沿岸以及濟(jì)寧市區(qū)周邊。這些地區(qū)不同時(shí)段內(nèi)差分干涉條紋密集，表明了各個(gè)煤礦在監(jiān)測(cè)時(shí)段內(nèi)開采沉陷的發(fā)生位置和范圍。在2004 年11 月到2005 年2 月3 個(gè)月的時(shí)間內(nèi)，共有4 個(gè)差分干涉紋圖(圖5a、b、c、d)顯示了該時(shí)間段的塌陷變化過程，各個(gè)采區(qū)的干涉條紋密度隨著時(shí)間間隔的增加而整體上增大。由041121 -050306 差分干涉紋圖(圖5b)可知兗州市城區(qū)東部與曲阜交界(序號(hào)2)的開采沉陷漏斗干涉條紋達(dá)5 個(gè)色周，沉陷量達(dá)到15 cm，表明在105 d內(nèi)由于煤礦開采引起了最大15 cm的沉降，沉降呈現(xiàn)出盆地形狀的分布特征，與傳統(tǒng)開采沉陷的沉陷盆地分布相同。位于濟(jì)寧南張鎮(zhèn)西北(序號(hào)1)的沉陷漏斗在041121 -050130(圖5a)的干涉紋圖超過2 個(gè)色周，沉陷量超過6 cm。050103 和050306 兩景SAR 數(shù)據(jù)相干性好，差分干涉紋圖(圖5d)各沉陷中心清晰，可以明顯看到該沉陷漏斗中心向北偏移。

2005 年12 月到2006 年2 月有1 個(gè)相干性好的干涉紋圖(圖5e)，該干涉紋圖表現(xiàn)了該時(shí)間段內(nèi)各沉陷中心的沉陷狀況，曲阜與兗州在泗河交界處的沉陷漏斗干涉條紋有5 -6 個(gè)色周，中心沉陷量達(dá)15 ～18 cm;濟(jì)寧南張鎮(zhèn)-劉堤頭間(序號(hào)3)的沉陷漏斗條紋達(dá)4 個(gè)色周，中心沉陷量超過12 cm。2009 年12 月到2010 年10 月相干性良好的3 個(gè)干涉紋圖(圖5f、g、h)清晰地表現(xiàn)了該時(shí)間段內(nèi)的開采沉陷狀況。由091220 -100124(圖5f)和100124 -100228(圖5g)干涉紋圖可知，位于曲阜與兗州在泗河交界處張家村煤礦一帶(序號(hào)4)的沉陷范圍在35 d的時(shí)間間隔內(nèi)達(dá)到3 個(gè)色周的沉降變化，表明在091220 -100228共70 d的時(shí)間內(nèi)沉陷量達(dá)到近18 cm，屬于快速沉降階段。濟(jì)寧主城區(qū)西北前王村南沉陷中心(序號(hào)5)在091220 -100124 和100124 -100228 干涉紋圖中分別達(dá)到了6 個(gè)和3 個(gè)以上色周，共計(jì)約10 個(gè)條紋，表明最大沉陷量超過30 cm。圖5h 為直接利用091220 和100228 兩景數(shù)據(jù)干涉后的條紋圖，反映了105 d內(nèi)整個(gè)沉降區(qū)的變化狀況，是圖5f 和圖5g 反映的沉降量的總和。

總體上，干涉圖的空間位置和分布形態(tài)一致，但局部干涉圖的條紋清晰程度降低。這種結(jié)果產(chǎn)生的原因是由于時(shí)間間隔進(jìn)一步增大，相干性下降引起的。同時(shí)由于沉陷量增大，局部出現(xiàn)非相干的移動(dòng)變化，這種快速變化限制了利用InSAR 干涉圖之間確定沉降量的精準(zhǔn)程度。同時(shí)，位于樓莊和梁家營(yíng)一帶出現(xiàn)了局部回彈，干涉條紋與其它沉降引起的形態(tài)相反，變化量為1個(gè)色周，影響范圍約2 ～3 km。

3．2 全區(qū)2004 ～2010 年區(qū)域性地面沉降InSAR 連續(xù)監(jiān)測(cè)

利用長(zhǎng)時(shí)間序列進(jìn)行PSInSAR 處理，提取相干目標(biāo)的沉降速率。研究區(qū)內(nèi)的濟(jì)寧、兗州等城市，地面沉降主要分為緩慢沉降以及煤礦開采沉陷。前者為緩慢微小變化，主要發(fā)生于濟(jì)寧市城區(qū)以及城區(qū)周邊，而后者則多發(fā)生于煤礦區(qū)，與采礦條件、煤層厚度與賦存條件，開采方式等因素相關(guān)。因而，在煤礦區(qū)開采沉陷動(dòng)態(tài)監(jiān)測(cè)的基礎(chǔ)上，實(shí)現(xiàn)全區(qū)地面沉降的連續(xù)監(jiān)測(cè)。

圖5 工作區(qū)煤礦開采沉陷InSAR 動(dòng)態(tài)監(jiān)測(cè)結(jié)果圖Fig．5 Dynamic monitoring results for applying InSAR in subsidence area caused by coal mining

利用2004 年12 月～2006 年02 月和2008 年9 月～2010 年10 月兩個(gè)時(shí)間段的ENVISAT ASAR 數(shù)據(jù)進(jìn)行PSInSAR 時(shí)序分析，得到兩個(gè)時(shí)間段的平均沉降速率，如圖6、圖7 所示。監(jiān)測(cè)結(jié)果表明:工作區(qū)總的沉降范圍與礦區(qū)分布一致，分別位于濟(jì)寧市任城區(qū)北部、兗州市東南部、曲阜市西部以及鄒城市西北部;工作區(qū)內(nèi)有兩大明顯的沉降區(qū)域，兗州市東南、曲阜市西部以及鄒城市西北部沉降連接成片，沉降速率大于20 mm/a的區(qū)域超過180 km2，另一個(gè)明顯的沉降區(qū)域?yàn)闈?jì)寧市任城區(qū)北部。

2008 年09 月～2010 年10 月InSAR 監(jiān)測(cè)地面沉降結(jié)果(見圖7)表明兩大沉降區(qū)沉降范圍有所擴(kuò)大，兗州—曲阜—鄒城沉降區(qū)內(nèi)沉降中心沉降速率有所減緩，但西部整體沉降加劇，東部有所減緩;濟(jì)寧任城區(qū)沉降區(qū)與汶上縣、嘉祥縣交界處沉降明顯減緩，主城區(qū)西北郊南張鎮(zhèn)一帶有所加劇，北部岱莊—栗鄉(xiāng)一帶有所減緩，但任城區(qū)東部和南部出現(xiàn)明顯沉降中心。

3．3 全區(qū)累積地面沉降量分析

通過2004 年12 月至2010 年10 月共24 景ASAR 數(shù)據(jù)的PSInSAR 時(shí)序分析處理，得到6 年間工作區(qū)累積地面沉降量，如圖8 所示。由InSAR 監(jiān)測(cè)累積地面沉降量所示可知，在近6 年的時(shí)間內(nèi)，全區(qū)最大累積沉降量超過250 mm，兩個(gè)沉降區(qū)分別為濟(jì)寧市北部煤礦區(qū)和曲阜、兗州交界沿泗河一線。

圖6 2004．12 ～2006．02 InSAR 監(jiān)測(cè)地面沉降速率圖Fig．6 Land subsidence rate by applying InSAR from December 2004 to February 2006

濟(jì)寧、兗州轄區(qū)約4 060 km2的面積，去除大面積水域共約3 922 km2，通過對(duì)累積地面沉降量空間插值分析統(tǒng)計(jì)可得各沉降范圍影響面積，如表3 所示。

圖7 2008．09 ～2010．10 InSAR 監(jiān)測(cè)地面沉降速率圖Fig．7 Land subsidence rate by applying InSAR from September 2008 to October 2010

圖8 2004．12 ～2010．10 InSAR 監(jiān)測(cè)累積地面沉降量Fig．8 The cumulative land subsidence amount by applying In-SAR from December 2004 to October 2010

表3 累積地面沉降影響范圍統(tǒng)計(jì)表Tab．3 The statistical table of influence range about cumulative land subsidence

對(duì)濟(jì)寧、兗州各沉降范圍影響面積進(jìn)行單獨(dú)統(tǒng)計(jì)，各地區(qū)主要沉降主要分布，如表4、表5 所示。

表4 濟(jì)寧市累積地面沉降影響范圍統(tǒng)計(jì)Tab．4 The statistical table of influence range about cumulative land subsidence in Jining city

表5 兗州市累積地面沉降影響范圍統(tǒng)計(jì)Tab．5 The statistical table of influence range about cumulative land subsidence in Yanzhou city

3．4 礦區(qū)范圍與沉降范圍比較

通過對(duì)累積地面沉降量空間插值分析得到該區(qū)域累積沉降等值線圖，將濟(jì)寧和兗州地區(qū)煤礦區(qū)礦權(quán)范圍數(shù)據(jù)與實(shí)際開采沉陷范圍進(jìn)行疊加，得到如圖9 所示的累積沉降范圍與開采范圍。圖中左下角深灰色范圍為主要煤礦區(qū)礦權(quán)劃定范圍，比較發(fā)現(xiàn)，兩個(gè)主要沉降區(qū)分別對(duì)應(yīng)著兗州和濟(jì)寧地區(qū)兩個(gè)主要礦區(qū)范圍，其它無礦權(quán)范圍地區(qū)均為較小沉降區(qū)。對(duì)比沉降量的大小可知，大于60 mm沉降量的地區(qū)為煤礦開采區(qū)、整個(gè)沉降量為60 mm等值線所包圍的區(qū)域均為煤礦礦權(quán)覆蓋范圍;東山王樓煤礦沉降突出，中心區(qū)大于80 mm，影響范圍較廣。

圖9 工作區(qū)累積地面沉降范圍與礦區(qū)開采范圍Fig．9 The mining range and the cumulative land subsidence range in study area

4 結(jié)論

本次研究工作針對(duì)兗州和鄒城地區(qū)區(qū)域性地面沉降和煤礦開采區(qū)快速沉陷兩種類型地表形變，利用2004 ～2010 年多期ENVISAT ASAR 雷達(dá)數(shù)據(jù)，應(yīng)用PSInSAR 處理技術(shù)與常規(guī)D-In-SAR 差分干涉處理方法，獲取了全區(qū)2 個(gè)不同時(shí)間段的區(qū)域性地面沉降速率和6 年時(shí)間內(nèi)全區(qū)累積沉降量結(jié)果，并得到了多個(gè)時(shí)段礦區(qū)動(dòng)態(tài)沉陷監(jiān)測(cè)成果。監(jiān)測(cè)結(jié)果查明了濟(jì)寧、兗州地區(qū)地面沉降分布狀況，統(tǒng)計(jì)得到了全區(qū)累積沉降量和影響地區(qū)范圍，全面掌握了該區(qū)域目前地面沉降的總體分布態(tài)勢(shì)和局部沉降特征。區(qū)域性地面沉降結(jié)果表明，濟(jì)寧、兗州的地面沉降主要與礦山開采相關(guān)，城市活動(dòng)等其他因素帶來的地面沉降甚微;工作區(qū)內(nèi)地面沉降顯著區(qū)有兩個(gè)，一個(gè)為兗州市東南、曲阜市西部以及鄒城市西北部沉降區(qū)，另一個(gè)明顯的沉降區(qū)域?yàn)闈?jì)寧市任城區(qū)北部。比較兩期監(jiān)測(cè)成果，雖然最大沉降速率有所減小，但沉降嚴(yán)重區(qū)域(年沉降速率＞20 mm/a)有擴(kuò)大的趨勢(shì)。

煤礦開采沉陷是上述區(qū)域沉降的主要因素，由于煤礦開采和開采過程的輸排水等，改變了工作區(qū)局部地質(zhì)環(huán)境，出現(xiàn)地面破壞。這種破壞與大規(guī)模開采地下水引發(fā)的區(qū)域性大范圍沉降漏斗不同，從空間分布上沉降中心比較集中，沉降速率比較大，局部甚至出現(xiàn)地裂縫和塌陷。因而，從監(jiān)測(cè)手段上InSAR 技術(shù)實(shí)現(xiàn)面上控制和城市全覆蓋，對(duì)于單個(gè)點(diǎn)上的連續(xù)快速沉降，還需進(jìn)一步利用其它地面觀測(cè)的方式開展工作。

為了保證本研究工作的連續(xù)性和發(fā)展性，下一步工作中擬采用RADARSAT-2 衛(wèi)星獲取的精細(xì)模式SAR 數(shù)據(jù)開展相應(yīng)的監(jiān)測(cè)工作。該數(shù)據(jù)分辨率為8 m，優(yōu)于ENVISAT 數(shù)據(jù)，覆蓋范圍為50 ×50 km2，獲取周期為24 d，應(yīng)用PSInSAR時(shí)序分析方法，連續(xù)獲取數(shù)據(jù)既可滿足礦區(qū)動(dòng)態(tài)監(jiān)測(cè)，也可滿足項(xiàng)目全區(qū)緩慢沉降監(jiān)測(cè)的需要。

［1］ 劉國(guó)祥．利用雷達(dá)干涉技術(shù)監(jiān)測(cè)區(qū)域地表形變［M］．北京:測(cè)繪出版社，2006:148．

［2］ 舒寧．雷達(dá)影像干涉測(cè)量原理［M］．武漢:武漢大學(xué)出版社，2003:99．

［3］ 程滔，單新建． CR、PS 干涉測(cè)量聯(lián)合解算算法研究［J］． 地震，2007，27(2):64 -71．

［4］ 范景輝，李梅，郭小方，等．基于PSInSAR 方法和ASAR 數(shù)據(jù)監(jiān)測(cè)天津地面沉降的試驗(yàn)研究［J］．國(guó)土資源遙感，2007(4):23 -27．

［5］ 張潔．干涉雷達(dá)數(shù)據(jù)處理技術(shù)與應(yīng)用［D］． 武漢:中國(guó)地質(zhì)大學(xué)(武漢)，2005．

［6］ 韓春林．SAR 成像技術(shù)和SAR 圖像處理算法研究［D］．成都:電子科技大學(xué)，2004．

［7］ 楊成生． 基于D-InSAR 技術(shù)的煤礦沉陷監(jiān)測(cè)［D］． 西安:長(zhǎng)安大學(xué)，2008．

［8］ 董玉森．干涉雷達(dá)中的永久散射體和交叉干涉技術(shù)的研究與應(yīng)用［D］．武漢:中國(guó)地質(zhì)大學(xué)(武漢)，2007．