甘肅金川礦區地表三維形變InSAR監測

楊旺， 何毅， 張立峰, 王文輝， 陳有東， 陳毅

(1.蘭州交通大學測繪與地理信息學院,蘭州 730070； 2.地理國情監測技術應用國家地方聯合工程研究中心,蘭州 730070； 3.甘肅省地理國情監測工程實驗室,蘭州 730070)

0 引言

礦產資源在經濟發展中起到了極大的推動作用，其中，礦業占了國內生產總值的5.5%，礦產資源的開發利用已成為助推中國經濟發展的重要動力[1]。但是，礦產資源的大規模開發不可避免地導致地表形變、塌陷和其他潛在的災害隱患。中國地質環境監測院調查結果顯示，2015年末，中國礦山因地下開采形成的采空區體積達到12.8億m3，引發了礦區坍塌、塌陷等多種地質災害，造成大量的財產損失和人員傷亡[2]。因此，有必要監測礦區地表形變，并進行多維分析。

目前，已有眾多學者利用GPS和水準測量等方式對礦區的地表形變進行了監測[3-5]，但這些方法存在監測范圍小、效率低、成本高和監測點密度低的缺點[6]。合成孔徑雷達干涉測量(interferometric synthetic aperture Radar,InSAR)是一種利用相位差來獲取地表形變信息的高精度測量技術[7]，該方法監測點的密度高、范圍大且效率高[8]。諸多學者的研究表明，InSAR技術是監測礦區地表形變的有效手段[9-12]，例如Li等[13]使用TS-InSAR(time series InSAR)技術結合建筑物破壞水平經驗模型對淮南礦區的居民住宅進行了穩定性評價； Ma等[14]利用SBAS-InSAR(small baseline subset InSAR)技術對神東煤田的開采形變狀況進行了調查和分析。總結上述研究并結合文獻綜述獲知，雷達數據的合理選擇對礦區地表形變的監測至關重要。目前主流的雷達衛星有C波段的Sentinel-1和GF-3，L波段的ALOS，X波段的TerraSAR-X，COSMO-SkyMed和TanDEM-X等。由于Sentinel-1A數據開放，且具備影像幅寬大、重返周期短、成像質量高等優勢，在礦區地表形變監測中得到了廣泛應用[15-16]。然而，常規InSAR技術只能獲取視線向(line-of-sight,LOS)的形變，聯合多個軌道模式的影像不僅可以求解地表的三維形變，而且可以緩解由雷達側視引起的透視收縮、疊掩和陰影[17]。

金川礦區是我國最大的鎳礦開采基地之一，其中龍首礦作為金川鎳礦的主力礦山，由于連續開采，2016年出現閉環形裂縫，就此停止采礦。2018年首次采用無底柱分段崩落法繼續開采，該方法具有機械化程度高、開采強度大以及采礦成本相對較低等特點，但引起的地表形變量級較大[18-19]。地表形變是龍首礦嚴重的次生災害，對礦山安全生產活動造成了巨大威脅。因此，有必要探究復采以來龍首礦區地表形變的多維時空特征。基于此，丁寧等[20]利用SBAS-InSAR技術對龍首礦露天礦坑2009—2011年和2015—2016年的地表形變進行了監測，分析了三維形變特征，但該研究存在監測時段不連續以及最新研究資料欠缺的疏漏，并缺乏對金川礦區整體形變的監測分析； 李強[18]聯合SBAS-InSAR、無人機傾斜攝影測量和水準測量對龍首礦工業化試驗區的地表形變特征進行了研究。仍存在監測范圍小、數據一致性差、形變機理揭示不夠精細的缺陷[20]。因此，本文基于升降軌模式下3個軌道(升軌128、降軌33和135)的Sentinel-1A數據，利用SBAS-InSAR技術并結合先驗條件的最小二乘迭代法，獲取了金川礦區管轄范圍內所有生產區2018—2020年間的地表三維形變時空特征，并就其形變機理展開了全面剖析。

1 研究區概況及數據源

1.1 研究區概況

金川礦區位于龍首山斷裂帶，北接阿拉善臺塊，南鄰北祁連山縫合帶[20](圖1)。該礦是金昌市最早發現并開采利用的老礦山[21]，擁有我國最大的鎳儲蓄量，也是我國鉑金屬的主要產地。該礦體呈NW-SE走向分布，長約6 000 m，寬約300 m，中間部分厚度約1 000 m[22]。礦區被一系列SW-NE走向的沖斷層切割，將侵入礦體分成若干段，甘肅省地質調查院將其劃分為Ⅰ，Ⅱ，Ⅲ，Ⅳ礦區[23]，其中Ⅰ和Ⅲ礦區由金川公司龍首礦開采和管理，Ⅱ礦區由西二礦開采和管理，Ⅳ礦區由東三礦開采和管理[19](圖1(b))。金川礦區土地覆蓋以沙化土地為主，包含砂、黏土和礫石等，植被覆蓋率較低，后向散射強度較高，開采區地形起伏程度相對較小，適用 InSAR進行地表形變監測。

(a) 金昌市地形特征及影像覆蓋范圍(b) 金川礦區采礦區分布

1.2 數據源

Sentinel-1A衛星是歐洲航天局和歐盟委員會推出的哥白尼計劃中的首個空間對地觀測衛星，搭載C波段雷達傳感器[24]，重訪周期為12 d，包含4種工作模式： 條帶模式(strip map,SM)，干涉寬幅模式(interferometric wide,IW)，超寬幅模式(extra-wide swath,EW)以及波束采集模式(wave mode acquisitions,WM)[25]。其中IW模式采用TOPSAR(Terrain observation with progressive scans SAR)技術成像，幅寬為240 km,空間分辨率為5 m×20 m，影像質量較高，噪聲更易處理。此外，相同極化(VV)比交叉極化對地物的介電常數和表面粗糙度更敏感[26]，因此，本研究選擇IW模式、VV極化的雷達影像。

本研究選用3個軌道(升軌128，降軌33和135)的Sentinel-1A數據反演金川礦區時序三維形變，下載獲取了覆蓋研究區128軌道(2018年1月2日—2020年10月6日)81景、33軌道(2018年 1月7日—2020年10月23日)81景和135軌道(2018年3月27日—2020年10月11日)71景影像，影像參數如表1所示[19]。此外，本文利用數字高程模型(digital elevation model，DEM)去除地形相位并對形變結果進行地理編碼，DEM選用美國國家航空航天局、美國國家測繪局和德國、意大利航天局聯合繪制的SRTM(shuttle Radar topography mission)地形數據，空間分辨率為90 m[26]。為了驗證本文形變監測結果的有效性和準確性，選取部分水準測量數據[18]與研究區的時序InSAR結果做了檢驗。

2 研究方法

2.1 SBAS-InSAR技術

SBAS-InSAR是由Berardino等[25]和Lanari等[27]提出的時間序列InSAR分析技術，該技術通過將獲取的數據進行分組配對，得到一系列短的空間基線差分干涉圖集，利用奇異值分解法(singular value decomposition, SVD)將獨立的數據集進行連接并生成相干目標的平均形變速率和時間序列累計沉降值[28]。SBAS-InSAR技術的優勢是充分利用了多時相的SAR數據，提高了形變測量的時間采樣率和空間覆蓋率[28]。本文數據處理中，SBAS-InSAR流程如圖2所示。具體流程如下:

1)干涉對組合和連接圖生成。將獲取的數據按時間和空間閾值進行分組配對，綜合考慮時相數、影像質量和計算機的處理能力設置閾值，本研究中3個軌道時間基線閾值均設為240 d，空間基線閾值設為45%。超級主影像由軟件自動選取，原則為時間基線和最小，影像日期分別為： 2018年7月13日(升軌128)、2018年12月28(降軌33)、2019年3月22日(降軌135)，生成時空基線圖如圖3所示。

圖2 SBAS-InSAR處理流程

(a) 128軌道時間基線 (b) 128軌道空間基線

(c) 33軌道時間基線 (d) 33軌道空間基線

(e) 135軌道時間基線 (f) 135軌道空間基線

2)干涉工作流處理，包括干涉圖生成、去平、濾波和相位解纏。首先，利用復數像對共軛相乘得到干涉圖，為了提高原始影像的信噪比和相干性，多視比設置為4∶1； 其次，利用參考DEM去除地形相位； 接著選用Goldstein自適應濾波算法進行干涉相位濾波； 最后，利用最小費用流解纏方法獲得差分干涉序列的解纏相位圖，解纏相干性閾值設置為0.2。

3)軌道精煉和重去平。各軌道數據分別選取30個地面控制點(ground control points，GCP)，且滿足GCP相干性良好、遠離形變區，用殘余相位優化方法來估算和去除殘余的固定相位和解纏后還存在的相位坡道。

4)SBAS-InSAR反演。首先，通過選擇和細化GCP來估算相位坡道[29]； 其次，利用線性模型估計位移速率和殘余地形； 最后，采用SVD方法獲取時間序列地表形變。在估計平均形變速率后，通過利用時間域高通濾波和空間域低通濾波去除大氣相位延遲。

5)反演結果進行地理編碼。DEM作為坐標參考，生成坐標系(GCS-WGS-84)下的LOS形變和垂直向形變。

2.2 三維形變構建方法

2.2.1 InSAR觀測方程

SAR衛星多采用近極地軌道，雷達波束照射到同一地面區域時飛行狀態有2種: 一種是由南向北，稱為升軌; 一種是由北向南，稱為降軌(圖4 (a))。由于升降軌使用同一衛星平臺，成像時各項物理參數幾乎未發生變化，所以影像質量一致性較好。升降軌成像最大的區別就是雷達波束入射角和衛星飛行的方位角不同，利用這一特征，可以構建形變監測結果的三維分解模型[30]。InSAR觀測結果是將真實形變投影到LOS方向上，LOS向的形變又可以分解到垂直、東西、南北向上，其幾何關系如圖4(b)表示，其關系式為：

DLOS=DUcosθ-DNsinθcos(α-3π/2)-DEsinθsin(α-3π/2)，

(1)

式中：DLOS為真實形變投影到LOS向形變值；DU，DN和DE分別表示垂直向、南北向和東西向形變(約定垂直向上、北向、東向為正)；θ為SAR傳感器入射角;α為衛星飛行方位角。

(a) 衛星飛行軌道示意圖 (b) 地表形變三維幾何關系

2.2.2 結合先驗條件的最小二乘迭代法

結合先驗條件的最小二乘迭代法[30]是利用已知條件或者假設可信已知條件作為多余觀測，迭代使用最小二乘法求解矛盾方程獲得變量近似解的間接求解方法。Sentinel-1A升降軌道均為近極地軌道，三維形變分解中南北向對形變信息極不敏感[31]，則可以將南北向位移值向極小值逼近作為初始先驗條件，結合式(1)進行最小二乘解算，依據LOS形變中三維分量的貢獻程度和理論敏感度的高低，依次作為先驗約束條件迭代使用最小二乘求解獲得三維分解值，具體步驟如下：

1)獲得多個軌道時序InSAR結果之后，根據式(1)，聯立三維求解方程組，下標表示3個不同軌道，矩陣形式表示為：

(2)

首先設DN為0，增加多余觀測，此時未知數為DE和DU，而3個軌道可構成3個觀測方程，利用最小二乘法可以求得DE和DU，由于DU分解敏感度較高，認為DU數值是可靠的。此步驟將DN設為0，因為南北向分解敏感度較低，設為0可作為條件約束，但不表示南北向位移為0，在前兩步約束求解之后可依據DE和DU求解DN的近似值。

2) 獲得DU之后，將其代入式(2)構成多余觀測，2個未知數為DN和DE，再次運用最小二乘法求解DN和DE，為了提高DU求解精度，只將DE作為可靠值代入下一步。

3) 設式(1)—(2)中求得的DU和DE視為真值，此時未知數只有一個即DN。3個條帶的方程均可求出DN，此時構成2個多余觀測，再次運用最小二乘求解得到DN。

3 金川礦區地表三維形變特征

3.1 地表三維形變速率特征

利用SBAS-InSAR和三維形變解算后，得到三維形變速率(圖5)。圖5(a)為垂直向形變速率，圖中龍首礦、西二礦、東三礦開采區出現了3個明顯的形變區A，B和C，從西北至東南方向依次分布，其中A區形成一個長軸約為2.6 km，短軸約為0.9 km，面積約為1.1 km2的橢圓狀形變區，垂直向最大形變速率為-132.4 mm/a； B區形成一個形變面積約為1.8 km2的圓形形變區，垂直向最大形變速率為-162.8 mm/a； C區形成一個面積約為0.71 km2的半圓形形變區，垂直向最大形變速率為-52 mm/a。圖5(b)為東西向形變速率，圖中3個形變區的西部均向東發生少量形變，最大形變速率為20.5 mm/a； 東部向西收縮，最大形變速率約為-23.8 mm/a。圖5(c)為南北向形變速率，3個沉降區西部均向南發生形變，最大形變速率約為-32.4 mm/a； 東部均向北發生形變，最大形變速率約為59.2 mm/a。綜上，B形變區整體形變速率最大，A區次之，C區由于東三礦建成投產時間較短，開采力度較小，整體形變速率最小。

(c) 南北向形變速率(d) 三維分解幾何條件

3.2 地表三維形變時序累計形變量特征

為分析礦區形變過程，解算2018年1月—2020年10月6個時段的三維累積形變值(圖6)，便可知曉該過程。從圖6中可以看出，3個形變區的形變范圍逐漸擴大，相對于2018年1月，垂直向累積形變值在逐漸增大，研究時段內的最大累積形變值為-408.9 m。此外，東西向和南北向的累積形變值也在逐漸增大，其中，東西向形變范圍為-85.4～152.8 mm，南北向形變范圍為-60.6～54.0 mm。金川礦區垂直向的形變趨勢主要表現為地面沉降，西二礦開采區的累積形變值最大，龍首礦次之，東三礦開采區的累積形變值最小。水平位移趨勢均表現為形變區西北部和東南部向中部收縮。

(a) 2018.1—2018.9(垂直向)(b) 2018.1—2019.3(垂直向)(c) 2018.1—2019.8(垂直向)

(d) 2018.1—2020.1(垂直向)(e) 2018.1—2020.6(垂直向)(f) 2018.1—2020.10(垂直向)

(g) 2018.1—2018.9(東西向)(h) 2018.1—2019.3(東西向)(i) 2018.1—2019.8(東西向)

(j) 2018.1—2020.1(東西向)(k) 2018.1—2020.6(東西向)(l) 2018.1—2020.10(東西向)

(m) 2018.1—2018.9(南北向)(n) 2018.1—2019.3(南北向)(o) 2018.1—2019.8(南北向)

(p) 2018.1—2020.1(南北向)(q) 2018.1—2020.6(南北向)(r) 2018.1—2020.10(南北向)

為進一步分析形變區的時序特征，求取圖5(a)中剖面線a，b，c和d在2018年1月—2020年10月期間6個時段的剖面圖(圖7)。

(a) 剖面線a (b) 剖面線b

(c) 剖面線c (d) 剖面線d

剖面線a的坐標方位角為124°，全長8 233 m，呈NW-SE走向同時穿過A，B和C形變較為顯著的區域，在1 796 m，4 321 m和6 282 m處出現3個沉降漏斗； 剖面線b，c和d坐標方位角均為214°，全長分別為2 703 m，3 920 m以及3 274 m，在b，c和d上1 451 m，1 952 m和1 528 m處分別出現沉降漏斗。由剖面線趨勢可以發現，3個形變區中沿NW-SE向的剖面形變曲線基本符合Peck公式[32]，NE-SW向的剖面圖顯示： 形變區A，B和C的東北邊界出現驟降趨勢，可能原因是受到自西向南開采活動的影響且受到斷層的控制。圖7中剖面線橫坐標取特定值時，各時相累積形變曲線的間距近似相等，可推斷3個形變區在近勻速沉降。為驗證該論斷，在形變區均勻抽樣選取測試點，位置分布如圖5(a)所示。繪制測試點累積形變曲線見圖8，結果顯示，50個測試點的形變曲線表現為均近直線形變趨勢。因此，可認為金川礦區的3個形變區，整體在勻速沉降，形變速率伴隨空間位置的不同而發生變化。

圖8 測試點形變曲線

4 討論

4.1 多軌道相關性分析

圖9給出了3個軌道垂直向形變速率的相關性密度散點圖。從圖9中可以看出，升軌128與降軌33和135軌道散點圖的相關系數分別為0.68和0.81，由于衛星飛行的方向和視角不同，雖然導致了散點圖左右兩側的偏離度增加，但卻形成了三維形變分解的幾何條件。由于衛星飛行方向視角接近，降軌33和135軌道的密度散點圖表現出較高的相關性，相關系數為0.96，相關直線的斜率接近0.5。綜上，3個不同軌道既構成了立體觀測條件，其監測結果又體現出較好的相關性，說明本文SBAS-InSAR處理結果可信度較好，三維形變分解方法較為合理。為了進一步證明時序InSAR的監測性能，選取實測數據對本文的監測結果進行驗證。

(a) 33軌道與128軌道(b) 135軌道與128軌道(c) 135軌道與33軌道

4.2 三維分解結果對比驗證

根據2019年6月和7月的水準實測數據[18]，選擇點位清晰、完整性較好的水準點1和2位置如圖5(a)所示，用于驗證本文三維形變的結果。由于水準實測數據記錄了所在點位的當月垂直向累積形變值，因此僅對垂直向分解結果進行驗證。提取水準點位置20 m半徑圓形緩沖區內同時期的SBAS-InSAR時序垂向形變矢量監測點，求其平均，將結果與水準監測值進行比較。從圖10可以看出，2個水準點的分解結果與水準測量結果的趨勢基本相同，最大誤差的絕對值分別為1.3 mm和0.7 mm。說明，SBAS-InSAR的監測結果與水準測量結果具有較高的一致性，證明本文實驗參數的設置較為合理，表明本文中的三維解算結果總體上是可靠的。

(a) 水準點1對比驗證(b) 水準點2對比驗證

圖10 垂向分解結果與水準數據的驗證

上述通過水準數據驗證了基于先驗條件的最小二乘迭代法進行地表形變三維分解的可靠性，為凸顯本文方法的優勢，將本文方法與常用三維分解算法進行對比分析。

聯合GPS數據和InSAR監測數據的三維形變求解方法可以高效高精度地求解地表三維形變，但是這種方法對數據有較高的要求，主要為： 時間域統一、空間域統一、數據融合權重的統一[33]。本文方法采用基于同源數據多視角觀測的形變三維分解策略，擺脫了對數據統一性的要求，并且在精度控制上不斷進步[34-36]。直接求解法和基于先驗條件的最小二乘迭代法均為基于同源數據多視角觀測的形變三維分解策略，直接求解法是對式(1)直接聯立方程組進行求解。將利用直接求解法計算獲得的三維形變速率與本文方法計算的結果進行對比分析，圖11為直接求解法和基于先驗條件的最小二乘法求解的2019年三維形變速率。從圖11中可以看出，直接求解法解算的三維形變速率東西向和南北向均沒本文方法效果好，尤其南北向的三維形變速率精度較差，這可能是直接求解法中由于InSAR觀測對南北向形變極不敏感，解算結果必然存在誤差放大現象，誤差最大可達數十厘米[37]，獲取的形變并不能代表真實的運動情況，因此，本文方法嘗試在有效降低形變解算誤差的同時，重建了更為真實可靠的三維形變場。此外，本文將圖10中水準數據形變量換算至形變速率，并與直接求解法和本文方法進行垂直向對比驗證(表2)，結果顯示： 直接求解法解算結果整體偏小，最大誤差達到33.8 mm/a； 本文方法中最大誤差為6.2 mm/a,結果較為理想，證明了本文方法的可靠性和可行性。依據精準的三維形變監測結果，可以預防并防治相關災害的發生。然而，形變機理的科學分析，是減緩形變演化行之有效的方法，鑒于此，就金川礦區2018年初—2020年末形變發生的機理進行討論。

(a) 本文方法(垂直向)(b) 本文方法(東西向)(c) 本文方法(南北向)

(d) 直接求解法(垂直向)(e) 直接求解法(東西向)(f) 直接求解法(南北向)

表2 垂直向形變速率對比驗證

4.3 金川礦區地表形變機理

在采礦區發生地表形變的關鍵因素一般涉及到機械化開采導致的地表應力變化、斷層的空間分布和巖性的結構強度[19]，因此，將從地表形變與人機采動、研究區斷層分布和巖性構造3方面分析礦區地表形變機理。

1)礦區地表形變與人機采動的關系。2018年復采后，采用無底柱分段崩落法采礦[19]，礦物從地下被開采出以后，礦房內部和周圍的應力狀態失穩，為了重新達到應力平衡，在采礦區上方覆蓋的巖層以及地表地物重力的作用下，巖體的變形和破壞主要表現為坍塌、斷裂、褶皺等形式，繼而引發地表形變[37]，形成以礦脈為中心的沉降漏斗。隨著開采工作的有序進行，金川礦區的形變趨于穩定。因此，2018—2020年金川礦區地表形變的原因之一可能是人機采動引起。

2)礦區地表形變與斷層的關系。如圖12所示，金川礦區北部被F1斷裂帶分割，下盤巖系主要包括侏羅紀—新近紀陸相碎屑沉積巖，上盤主要包括元古代變質巖和古生代沉積巖[23]，巖層曾被多代超基性巖脈侵入，形成了龍首山礦體[38-39]。隨著采礦進行，受F1斷層控制，上盤未發生明顯形變； 與礦體相間分布的F16，F17，F23與F1斷層相交，傾角變化范圍為60°～80°，將巖體分割為4個礦區，金川礦區形變區域也相應地表現出間斷性差異，龍首礦區北鄰F8斷層和F6斷層，形變區在此斷層邊界上出現間斷，斷層北部形變量整體減小，在剖面線b中1 000 m附近形變量驟增處即為F6斷層切割； 西南部受到F3斷層控制出現近似弧形形變邊界，整個開采區在斷層的作用下形變僅局限在近似橢圓范圍內。西二礦開采區除北部分布有F16斷層外周圍再無主要斷層存在，故開采區以礦脈為中心不斷向周圍擴張形成近似圓形形變區，但北部仍然受到斷層限制擴張范圍有限，且預測在開采工作不斷進行過程中除北部受到斷層控制外，形變區還會向外繼續延伸。東三礦開采區受到F17和F23斷層的控制形變區范圍基本局限在2個斷層中間形成呈半圓形分布。通過以上分析可推斷，金川礦區斷層的分布控制該區域地表形變。

圖12 金川礦區地質構造與累積形變疊加簡圖[23]

3)礦區地表形變與巖性的關系。金川礦區巖性復雜，區域差異性較大，龍首礦和西二礦開采區東北部主要為大理巖巖組、片麻巖巖組，主要由白云石和方解石組成，含有少量滑石、石英、方美石和云母等組成，顆粒均勻細小，質地均勻，硬度、結構強度較小，在應力和侵蝕作用下易于發生形變[28]； 西南部主要為片麻巖和復片麻巖巖組，主要由長石、石英和云母等組成，此外還包括角閃石、輝石等，顆粒大小不均，質地不均，結構強度整體較強。除巖層結構外，受大陸性溫帶干旱氣候和地表巖性的影響，金川礦區地表形成較厚的一層砂、礫石、黏土和黃土混合松散堆積層，該層整體密度較低，質松軟、結構強度較差，堆積層長期受到雨水侵蝕和滲透，降低了黃土、黏土、砂的黏聚強度和基質吸力，破壞了堆積層的穩定狀態，砂土覆蓋層軟化，濕陷性增加[29]，致使礦區人工建筑基礎發生坍塌，進而導致覆蓋層表面發生形變。金川礦區地表形變會隨著人工開采繼續進行，且西部開采區的巖體穩定性大于東部開采區[39]，隨著鎂鐵超鎂鐵侵入巖和超基性巖的多期侵入，節理裂隙發育破碎，東部開采區的沉降也將會更加明顯。因此，推斷巖性特征也是金川礦區地表形變的控制性因素。

5 結論

本文采用SBAS-InSAR技術分別利用3個軌道(升軌128、降軌33和135)的Sentinel-1A數據獲取金川礦區地表形變速率和形變量，并利用結合先驗條件的最小二乘迭代法解算礦區三維地表形變，對3個形變區空間變化規律和地表形變機理進行了分析。主要結論如下：

1)通過3個軌道垂直向累計形變值相關性分析實現解算結果的內部互檢，并使用水準實測數據對三維形變分解后垂直向累計形變值進行驗證，保證了地面形變監測的可靠性。

2)金川礦區出現3個形變區，其中B區形變范圍、整體形變速率和形變量均為最大； A區次之； C區最小。形變方式以沉降為主，最大年平均垂直向形變速率達到-162.8 mm/a，最大垂直向累計形變量達到-408.9 mm。此外，形變區也同時發生了水平面上的位移，主要表現為3個形變區西北部和東南部向形變區中部收縮。

3)地質斷層和巖性構造是礦區形變的控制性因素，F1斷層限制了礦區地表形變向東北向延伸，F16，F17和F23斷層的存在使得開采區形變出現間斷性差異。人機采動是礦區地表形變的主要因素，無底柱分段崩落法開采造成了開采區在30個月內近似均勻形變。