高 建 陳炳乾 馮遵德

（1.江蘇師范大學(xué)地理測繪與城鄉(xiāng)規(guī)劃學(xué)院，江蘇徐州221116；2.煤炭資源清潔利用與礦山環(huán)境保護湖南省重點實驗室，湖南湘潭411201）

合成孔徑雷達干涉測量（Interferometric Synthetic Aperture Radar，InSAR）是近年來發(fā)展起來的一種新型大地測量手段，該技術(shù)能夠進行大范圍、高精度、低成本和長周期地表沉降監(jiān)測［1］，很大程度上彌補了傳統(tǒng)監(jiān)測方法的不足，為礦區(qū)地表變形監(jiān)測提供了新方法。尤其是高級InSAR 技術(shù)（如小基線集技術(shù)、永久散射體技術(shù)等）的出現(xiàn)，進一步克服了傳統(tǒng)In-SAR技術(shù)時空基線的限制，并在一定程度上消除了大氣影響，提高了InSAR技術(shù)監(jiān)測礦區(qū)形變的精度和可靠性，使InSAR技術(shù)成為礦區(qū)變形監(jiān)測的重要手段之一［2-7］。然而，當(dāng)?shù)V區(qū)地表發(fā)生大梯度形變時（如地表塌陷），其形變梯度超出InSAR 最大可檢測形變梯度（MDDG）時，對應(yīng)的干涉條紋會發(fā)生混疊［8-9］，此時若采用常規(guī)解纏算法則無法恢復(fù)礦區(qū)的真實形變相位。為有效判斷InSAR 方法是否適用于地表變形區(qū)域探測，減少不必要的人力、物力與財力浪費，一些學(xué)者相繼提出了InSAR 形變梯度檢測模型。其中，Massonnet 等［10］于1998 年首次根據(jù)干涉測量的基本條件提出了MDDG 的概念并給出了InSAR 的MDDG理論公式，即為雷達波長的與像素大小的比值［10］。然而在實際InSAR 測量中，由于受時空去相干、大氣、噪聲等因素的影響，真實InSAR能夠檢測的MDDG 值遠小于理論值。因此，Baran 等［11］于2005 年將相干性引入到MDDG 函數(shù)模型中，提出了基于In-SAR 相干性的MDDG 函數(shù)模型。在此基礎(chǔ)上，蔣彌等［12］對Baran 模型進行改進，引入多視視數(shù)作為自變量，基于C 波段的SAR 影像，提出了融合InSAR 相干性和視數(shù)的MDDG函數(shù)模型。之后，王琪潔等［13］基于L波段的PALSAR數(shù)據(jù)建立了形式類似的函數(shù)模型。

以上函數(shù)模型為合理選擇變形監(jiān)測對象提供了有益參考，然而現(xiàn)有模型在實際應(yīng)用中需要提前獲取SAR影像計算相干性或相干性+視數(shù)來判斷InSAR技術(shù)應(yīng)用的適宜性，并且在應(yīng)用于礦區(qū)開采沉陷時未能考慮礦區(qū)開采的實際情況和不同的地質(zhì)采礦條件對開采沉陷的影響。本研究從煤礦開采引起的地表形變出發(fā)，結(jié)合開采沉陷的地質(zhì)采礦影響因素，提出一種適用于礦區(qū)開采沉陷的InSAR 形變梯度函數(shù)模型。該模型無需事先獲取SAR 影像計算相干性，只需利用地質(zhì)采礦條件參數(shù)即可判斷不同波長、不同分辨率和不同入射角的InSAR 數(shù)據(jù)的最大可檢測值，其可靠性通過實際礦區(qū)的SAR 影像試驗進行了驗證。

1 影響開采沉陷的地質(zhì)采礦因素

大量工程實踐表明，礦區(qū)地表沉陷是自然地質(zhì)與采礦因素綜合作用的結(jié)果。自然地質(zhì)因素包括煤層埋藏的幾何條件（如煤層厚度、傾角和埋深等），構(gòu)造因素（如褶皺程度和斷層密度等）和巖石物理、力學(xué)和化學(xué)性質(zhì)（如巖石硬度系數(shù)和水文地質(zhì)性能等）［14］。以上因素屬于不可控因素，因此本研究在后續(xù)建模過程中不予考慮。

影響地表沉陷的另一類主要因素是采礦技術(shù)因素，主要包括煤層開采深度、開采厚度、頂板管理方法和采空區(qū)尺寸等。現(xiàn)有研究表明：開采厚度越大，垮落帶、導(dǎo)水?dāng)嗔褞Ц叨仍酱?，地表移動和變形值也越大；同時隨著開采深度的增加，地表各種移動變形值減小，地表移動范圍擴大，移動盆地更平緩［15］。地表移動和變形既與開采厚度成正比，又與開采深度成反比，因此開采深厚比可作為衡量開采條件對地表沉陷影響的估計指標(biāo)［16］。

頂板管理方法決定了覆巖及地表的破壞程度、移動量大小。頂板管理方法有全部垮落法、充填法和煤柱支撐法等［17］。目前，全部垮落法是國內(nèi)礦區(qū)普遍采用的頂板管理方法，因此假定頂板管理均采用全部垮落法，采空區(qū)尺寸影響地表的充分采動程度。當(dāng)采動程度較小時，隨著工作面尺寸的增加，地表移動盆地范圍內(nèi)各種移動變形值均增大，當(dāng)達到充分采動后，工作面尺寸大小對地表移動和變形值大小無影響［18］。工作面尺寸大小對地表移動和變形值的影響規(guī)律可以用下沉系數(shù)來表征［14］。隨著工作面尺寸的增加，采動程度逐漸變大，同時下沉系數(shù)逐漸增加，當(dāng)工作面達到一定尺寸時，地表達到充分采動后，下沉系數(shù)趨于定值。因此，本研究將下沉系數(shù)作為后續(xù)建模的另一參數(shù)。

2 研究方法

本研究方法分為數(shù)據(jù)預(yù)處理過程和建模過程。數(shù)據(jù)預(yù)處理的主要目的是確定MDDG 與下沉系數(shù)之間的關(guān)系，生成初始模型，并探討SAR 影像參數(shù)對兩者關(guān)系的影響；建模過程是基于預(yù)處理中得到的結(jié)論來分析初始模型與開采深厚比的關(guān)系，從而將開采深厚比引入函數(shù)模型，最后選取合適的擬合函數(shù)得到最終的函數(shù)模型。由于數(shù)據(jù)預(yù)處理過程和后續(xù)建模過程的數(shù)據(jù)處理步驟一致，因此本節(jié)將著重闡述數(shù)據(jù)預(yù)處理步驟，如圖1所示。

為確定下沉系數(shù)與MDDG 的關(guān)系，并研究SAR影像參數(shù)對兩者關(guān)系的影響，在開采深厚比不變（所選開采深厚比為150）的前提下，選取不同的SAR 影像參數(shù)（波長、分辨率和入射角）和下沉系數(shù)進行分析。具體步驟如下：

（1）設(shè)置不同的下沉系數(shù)，利用概率積分法［19］模擬生成不同量級的形變場。

（2）基于步驟（1）中不同量級的形變場，設(shè)置SAR影像的波長、分辨率和入射角，并結(jié)合雷達干涉測量基本原理模擬得到對應(yīng)的差分干涉相位。

（3）采用最小費用流方法對步驟（2）得到的模擬差分干涉相位進行相位解纏。以概率積分法得到的形變相位作為參考，對解纏后的相位進行準(zhǔn)確性驗證，并提取解纏正確的相鄰兩像元以計算形變梯度。形變梯度計算公式為

式中，Gradient 為O、P 兩點間的形變梯度；ΔO和ΔP分別為兩相鄰像元O 點和P 點的下沉值；D 為像元間距離。

（4）統(tǒng)計分析步驟（3）獲取的可檢測形變梯度，得到不同SAR 影像參數(shù)（不同波長、不同分辨率）下下沉系數(shù)對應(yīng)的MDDG 值。通過線性擬合分析下沉系數(shù)與MDDG值的關(guān)系，確定臨界下沉系數(shù)。

（5）將步驟（4）得到的MDDG 值與臨界下沉系數(shù)進行第一次回歸分析，建立初始函數(shù)模型。

（6）預(yù)處理過程確定了臨界下沉系數(shù)與MDDG的關(guān)系，排除了SAR 影像參數(shù)對初始模型的影響，即模型無需引入SAR參數(shù)變量。因此，建模過程中的模擬數(shù)據(jù)只需以開采深厚比為變量，數(shù)據(jù)模擬與預(yù)處理過程一致。以步驟（5）得到的初始函數(shù)模型為基礎(chǔ)，探尋開采深厚比與初始模型的關(guān)系，并采用經(jīng)驗統(tǒng)計和二次回歸分析方法，確定最終的MDDG函數(shù)模型。

3 試 驗

3.1 數(shù)據(jù)預(yù)處理

本研究所有的模擬數(shù)據(jù)均在長度400 m、寬度200 m 的模擬工作面上獲得，并且為保證礦區(qū)變形影響范圍能夠被完整地模擬，設(shè)置SAR 影像區(qū)域為1 000 m×1 000 m。采用概率積分法進行形變場模擬，選取的參數(shù)如表1所示。

?

3.1.1 數(shù)據(jù)模擬與分析

現(xiàn)有研究表明［17］：當(dāng)開采深厚比小于30時，地表極易產(chǎn)生地裂縫、塌陷坑破壞，該地質(zhì)采礦條件下的形變梯度較大，地表移動和變形在空間和時間上呈現(xiàn)不連續(xù)性，且地表變形不具有明顯的規(guī)律性；當(dāng)開深厚比大于200 時，地表變形也無明顯規(guī)律。因此，本研究在保證表1中參數(shù)不變的前提下，固定開采深厚比為150，下沉系數(shù)在[ ]0.01,0.9 區(qū)間內(nèi)變化，同時下沉系數(shù)以0.01為增量模擬生成90個形變場。

為確保所建立的函數(shù)模型適用于不同波長、不同分辨率和不同入射角的SAR 影像，模擬試驗在將形變值轉(zhuǎn)化為差分干涉相位時，分別選用了常用的C、L 和X 波段對應(yīng)的不同分辨率、不同入射角等參數(shù)，具體參數(shù)取值如表2所示。

?

由表2 可知：3 種不同波長、不同分辨率和不同入射角條件下能夠?qū)?yīng)生成18 種不同SAR 影像參數(shù)下的形變相位，因此預(yù)處理階段模擬的90 個形變場可對應(yīng)生成1 620 幅差分干涉相位圖。將差分干涉圖經(jīng)最小費用流解纏得到解纏相位圖，以模擬變形轉(zhuǎn)化的形變相位為基準(zhǔn)值，規(guī)定基準(zhǔn)值的為限差，篩選解纏正確的像元并求取像元間的MDDG值。利用經(jīng)驗統(tǒng)計方法分析MDDG 值和下沉系數(shù)的數(shù)據(jù)發(fā)現(xiàn)：表2 中每一組參數(shù)對應(yīng)的90 幅差分干涉圖（下沉系數(shù)變化區(qū)間為[ 0.0 1,0.9]），其求得的像元間MDDG 值與下沉系數(shù)之間的規(guī)律具有相似性，具 體如圖2所示。

由圖2可知：各個子圖雖然選取的SAR影像參數(shù)不同，但其MDDG 值與下沉系數(shù)的關(guān)系在整體趨勢上保持一致，即前期MDDG 值隨著下沉系數(shù)呈線性增長，當(dāng)MDDG 增大到某一值后，MDDG 將基本保持不變。這是因為：在前期下沉系數(shù)較小，地表最大形變值沒有超過當(dāng)前影像參數(shù)下的MDDG 值，InSAR 技術(shù)能夠準(zhǔn)確地監(jiān)測礦區(qū)形變；隨著下沉系數(shù)增大到某一臨界值，地表最大形變梯度將等于該影像參數(shù)下的MDDG 值，這個臨界值可稱之為臨界下沉系數(shù)，即圖中的轉(zhuǎn)點處；當(dāng)下沉系數(shù)超過臨界下沉系數(shù)后，解纏方法只能恢復(fù)MDDG 范圍以內(nèi)的像元，此時MDDG值將不隨下沉系數(shù)發(fā)生變化。

通過以上分析發(fā)現(xiàn)，InSAR 技術(shù)MDDG 值的求取關(guān)鍵在于臨界下沉系數(shù)的確定。從模擬數(shù)據(jù)中提取了與波段、分辨率和入射角相對應(yīng)的臨界下沉系數(shù)及其最大可檢測形變梯度，如表3所示。

由表3 可知：不同分辨率和波段條件下，臨界下沉系數(shù)和MDDG 之間具有明顯的正相關(guān)關(guān)系。因此，后續(xù)回歸分析只需提取出臨界下沉系數(shù)及其模擬的差分干涉圖，即可構(gòu)建函數(shù)模型。

3.1.2 初始模型的生成

為研究MDDG 和臨界下沉系數(shù)之間的定量關(guān)系，并確定波長、分辨率和入射角對兩者關(guān)系的影響，對表3 中的6 組數(shù)據(jù)進行了擬合分析，結(jié)果如圖3所示。

?

由圖3 分析可知：MDDG 與臨界下沉系數(shù)之間存在明顯的線性關(guān)系，波長、分辨率和入射角不影響MDDG與臨界下沉系數(shù)之間的關(guān)系。

因此，利用表3 中的數(shù)據(jù)，對MDDG 與臨界下沉系數(shù)進行了線性回歸分析，得到的初始函數(shù)模型為

式中，Gmax為最大可檢測的形變梯度MDDG 值；qc為臨界下沉系數(shù)。

式（2）函數(shù)模型的均方根誤差為9.877×10-5，說明擬合程度極為顯著。

3.2 模型構(gòu)建

為建立開采深厚比和臨界下沉系數(shù)與MDDG 之間的函數(shù)模型，仍以表1 中參數(shù)不變?yōu)榍疤?，下沉系?shù)以0.01 為增量在[ ]0.01,0.9 區(qū)間內(nèi)變化，開采深厚比以5 為增量在[ 30 ,60 ]區(qū)間內(nèi)變化，共模擬生成了630個不同量級的形變場。

考慮到回歸分析的精度問題和數(shù)據(jù)量冗余造成運算效率不高的問題，在該階段的形變相位轉(zhuǎn)化過程中去掉入射角38.7°的影像參數(shù)，其他參數(shù)仍與表2中一致。最終，形變場轉(zhuǎn)化為形變相位的過程相比預(yù)處理過程減少了即一個形變場對應(yīng)9 幅差分干涉相位圖），則630個形變場共生成5 670幅差分干涉圖。圖4展示了不同下沉系數(shù)、采深采厚比對應(yīng)的部分差分干涉圖。其中，圖4（a）、圖4（d）、圖4（g）影像參數(shù)均為C 波段，分辨率為5 m×5 m；圖4（b）、圖4（e）、圖4（h）的影像參數(shù)均為X 波段，分辨率為1 m×1 m；圖4（c）、圖4（f）、圖4（i）的影像參數(shù)均為L 波段，分辨率為5 m×5 m。所有影像的入射角θ均為23°。

本階段對于干涉圖處理、模型參數(shù)提取和擬合均與預(yù)處理過程一致，不再贅述。假設(shè)模型為

式中，a為模型系數(shù)；b為模型常數(shù)。

按開采深厚比的變化，一共擬合得到了7個初始函數(shù)模型，模型參數(shù)統(tǒng)計如表4 所示。分析表4可知：當(dāng)開采深厚比在[ 30 ,60 ]區(qū)間內(nèi)變化時，常數(shù)b較小，基本在0.000 2左右變化；結(jié)合式（3）可知，開采深厚比為150時，常數(shù)b量級縮小到原來的1/10。因此，開采深厚比對MDDG 的影響主要集中在系數(shù)a，常數(shù)b基本可以忽略。故而可將系數(shù)a視為有關(guān)開采深厚比的函數(shù)，即為。為分析系數(shù)a與開采深厚比的具體關(guān)系，根據(jù)表4 數(shù)據(jù)繪制了如圖5 所示的曲線圖。

?

由圖5 可知：隨著開采深厚比的增大，系數(shù)a 逐漸減小，系數(shù)a 與開采深厚比成反比。由于參數(shù)b 與之間并無明顯的線性函數(shù)關(guān)系，且其數(shù)值在[ 0.0 00 2,0.000 4 ]區(qū)間內(nèi)小范圍變化，因此可進一步分析其變化情況。綜上分析，為探究MDDG 與臨界下沉系數(shù)和開采深厚比之間的關(guān)系，其模型可假定為

利用上述模型對模擬數(shù)據(jù)中提取的臨界下沉系數(shù)和相應(yīng)的開采深厚比進行了回歸分析，結(jié)果如圖6所示。由此可得到最終的函數(shù)模型為

式中，模型對應(yīng)的均方根誤差為0.000 079 38，說明該模型的準(zhǔn)確性顯著，符合要求。

4 實例驗證

試驗區(qū)域所在礦區(qū)的煤層平均厚度為6.94 m，開采深度為235 m，煤層傾角1°～3°，屬于近水平煤層，工作面開采時間為2012 年10 月—2013 年3 月，研究所用的SAR 數(shù)據(jù)為TerraSAR-X 的High Resolution SpotLight 模式，方位向分辨率為0.85 m，距離向分辨率為0.91 m，時間跨度為2012 年11 月21 日—2013 年4 月2 日，共13 景數(shù)據(jù)。為獲取地表充分采動后的最大下沉值，同時減少時間失相關(guān)對SAR 影像的影響，采用文獻［20］方法，按時間先后順序依次對13 景數(shù)據(jù)進行了干涉組合，生成了12組干涉對，之后對所有組合的干涉對進行疊加獲取其總形變量，最終得到的垂直位移如圖7 所示。此外，選用地表觀測站的GPS 數(shù)據(jù)作為參考形變值檢驗所建立模型預(yù)計結(jié)果的準(zhǔn)確性。

由圖7 可知：礦區(qū)邊緣位置條紋清晰，但在工作面中心由于變形較大發(fā)生干涉條紋混疊，條紋模糊無法分辨，該區(qū)域無法反映礦區(qū)真實形變信息。根據(jù)臨界下沉系數(shù)定義并結(jié)合研究區(qū)工作面參數(shù)（開采深度235 m，煤層厚度6.94 m），通過模擬獲得該工作面的下臨界沉系數(shù)為0.14。利用式（5）求得在該時間段內(nèi)InSAR 的MDDG 值為0.008 56。由于觀測線的控制點以20 m 為間隔進行布設(shè)，因此根據(jù)模型分析結(jié)果可知，InSAR 技術(shù)在相鄰控制點間能夠檢測的形變量不超過0.171 2 m。

為驗證上述分析結(jié)論，提取的工作面地表控制點的GPS 點間監(jiān)測差值如圖8（a）所示。由圖8（a）可知：在模型預(yù)計結(jié)果內(nèi)點間形變梯度有5個部分。第1 部 分 包 括 TQ1～TQ2、TQ2～TQ3、TQ3～TQ4、TQ4～TQ5；第2 部分包括TQ6～TQ7、TQ7～TQ8；第3 部分在礦區(qū)中心區(qū)域，包括TQ12～TQ13、TQ13～14；第4部 分 包 括TQ19～TQ20、TQ20～TQ21；第5 部 分 包 括TQ22～TQ25。根據(jù)模型定義，當(dāng)控制點兩側(cè)形變梯度小于MDDG時，該位置的變形能被InSAR技術(shù)有效監(jiān) 測。故 分 別 選 取TQ1～TQ3、TQ7、TQ13、TQ20、T23～T25作為模型驗證對象。為此，分別提取的地表控制點的InSAR和GPS監(jiān)測值如圖8（b）所示。

圖8（b）顯示：在控制點TQ1～TQ3 部分，InSAR 與GPS 監(jiān)測結(jié)果相近，最大誤差不超過0.083 m；控制點TQ7 處監(jiān)測失真，誤差達到0.272 m；控制點TQ13 處監(jiān)測失真，誤差超過3 m；控制點TQ20 處監(jiān)測有效，誤差為0.021 m；在控制點TQ23～TQ25 部分，GPS 和InSAR監(jiān)測值基本吻合，相差最大不超過0.019 m。

由此可以看出，模型對連續(xù)點的檢測能力較強，如：TQ1～TQ3這3個相鄰控制點的點間形變梯度都在檢測區(qū)域內(nèi)，表明檢測結(jié)果可靠。而當(dāng)孤點處于檢測區(qū)域內(nèi)時，檢測精度受到影響，如TQ7、TQ13 點都是這種情況。但有一個特例TQ20 點，其檢測精度達到0.021 m。為分析其原因，將TQ7、TQ13 點與TQ20點的相鄰點進行對比，發(fā)現(xiàn)TQ21～TQ22 點間的形變差值接近0.171 2 m，誤差為9.8 mm，可稱為近檢測點；對于TQ7和TQ13兩點，其相鄰點的點間形變差值與0.171 2 m 的誤差均達到分米級以上，說明這兩點均為遠檢測點。因此可以判斷：若孤點的相鄰控制點為近檢測點時，模型檢測精度也會提高。綜合上述分析可知，本研究模型對連續(xù)點或相鄰控制點為近檢測點的孤點有較高的檢測能力。

此外，進一步利用蔣彌等［12］提出的形變梯度模型與本研究模型進行對比試驗。采用文獻［12］中的計算方法得到的礦區(qū)形變梯度為0.014 29，然后將SAR 影像多視比例設(shè)置為1∶1，利用生成的相干圖計算出礦區(qū)相干均值為0.385。根據(jù)蔣彌模型［12］，該礦區(qū)在相干系數(shù)—形變梯度坐標(biāo)系中的位置即為圖9中五角星符號所在區(qū)域。圖9 中陰影區(qū)域即為模型所求的可檢測區(qū)，顯然礦區(qū)的相干性與形變梯度無法采用InSAR技術(shù)進行有效檢測。

通過與蔣彌模型［12］的對比試驗可以發(fā)現(xiàn)，本研究模型克服了SAR 影像源和相干性的限制，能夠在影像獲取前提前計算出InSAR技術(shù)的礦區(qū)可檢測值。同時，該模型精度優(yōu)于蔣彌模型［12］，能夠細化到判斷任一控制點的可檢測性。

5 結(jié) 論

（1）利用數(shù)值模擬與回歸分析方法，分析了地質(zhì)采礦因素與InSAR 技術(shù)可檢測的最大形變梯度之間的關(guān)系，建立了適用于開采沉陷的InSAR最大可檢測形變梯度模型。利用真實礦區(qū)的SAR 影像數(shù)據(jù)，將所提出的模型與已有模型進行了對比試驗。結(jié)果顯示：該模型無需提前獲取SAR 影像，也無需計算SAR影像的相干性，只需利用開采深厚比和下沉系數(shù)兩個參數(shù)即可預(yù)計不同開采深厚比、不同下沉系數(shù)條件下的不同波長、不同分辨率和不同入射角的InSAR數(shù)據(jù)的最大可檢測值，此外該模型具有較高的預(yù)計精度，較強的細化水平，對于促進InSAR 技術(shù)在礦區(qū)開采沉陷監(jiān)測中的應(yīng)用有一定的參考價值。

（2）所建立的函數(shù)模型是在保證礦區(qū)煤層傾角為0的條件下獲得的，且該模型僅考慮開采穩(wěn)定后的變形情況，沒有考慮動態(tài)開采過程中InSAR技術(shù)的可適用性，這將對該模型的應(yīng)用造成一定的限制。因此，不同傾角煤層下的InSAR可檢測函數(shù)模型的建立將是下一步的研究方向。