基于SAR偏移量追蹤技術的露天礦大梯度形變監測

孫 軍，張 錦

(太原理工大學 礦業工程學院，太原 030024)

差分雷達干涉技術(differential interferometric synthetic aperture radar，DInSAR)和多時相InSAR技術(multi-temporal InSAR，MT-InSAR)在監測煤礦開采區大梯度沉陷時受到一定限制，表現為干涉結果中存在空白區域即失相干區域，SAR圖像偏移量追蹤(offsets tracking，OT)方法可以探測到大梯度形變，這種方法不需要相位解纏，因此受研究區域相位失相干影響較小。偏移追蹤技術最初應用于冰川移動監測，冰川移動具有移動速度快且不相干等特點。與差分干涉技術相比，強度相關偏移解算精度較低，但是可以估計出冰川速度矢量的大小和方向[1]。礦區形變監測中可以首先利用單視復數(single look complex，SLC)影像對煤礦周邊區域進行干涉分析，而失相干區域主要位于開采區，使用偏移追蹤技術測得其沉陷量，結果明顯優于差分干涉，偏移量追蹤方法可以彌補DInSAR技術監測大梯度形變信息的不足，兩種方法結合可獲取采區的大量級形變[2]。

偏移量追蹤法可分為兩類，強度追蹤法和相干性追蹤法。強度追蹤法利用SAR(synthetic aperture radar)影像的強度信息即振幅，計算影像塊間的歸一化互相關值(normalized cross correlation，NCC)，互相關值最大的位置即為最佳匹配位置，最后估計出斜距向和方位向位移[3]。相干性追蹤法也稱條紋可見算法或相干性優化過程，匹配窗口內的影像塊共軛相乘形成干涉條紋，隨后搜索相關性峰值所在位置并估計偏移量[4]。ERTEN et al[5]在強度追蹤法中引入振幅影像的斑點結構信息，提出最大似然追蹤法，與傳統方法計算結果對比分析表明最大似然法受匹配窗口尺寸影響較小，而傳統方法在大梯度形變區域易產生誤匹配結果。WANG et al[6]歸納了偏移計算中匹配窗口均勻分布的缺點：高相干區域會導致偏移值的高估；強散射點則會造成鄰近匹配結果出現斑塊狀現象。并提出通過散射性強的點狀目標計算偏移，計算結果表明基于強散射點計算偏移提高了匹配像素對的相干性和結果的可靠性。自適應互相關(adaptive cross correlation，ACC)方法中將影像匹配窗口劃分為較小窗口計算互相關值，這種方法可以提高偏移估計精度[7]，但較小的匹配窗口易受噪聲影響且無法探測出大量級形變。CAI et al[8]分析滑坡及其周邊區域的特征，從規則匹配窗口中提取出與移動區域邊界相近的自適應窗口，只有自適應窗口中的像素參與互相關計算，這種方法可以有效降低未發生移動或者不相關區域的灰度值對相關值計算結果的影響，與通過規則窗口計算相關性的方法對比表明：自適應匹配窗口下位移估計結果能較好地保持穩定。上述偏移量追蹤方法提高匹配準確性的方式可以歸納為三種：1) 通過SAR影像斑點結構信息，根據最大似然原則計算匹配影像的概率，從而獲得最佳匹配位置；2) 與永久散射體相似，通過提取高相干目標點可以估算出偏移；3) 調整窗口尺寸，選擇規則或不規則窗口，提高匹配影像間的互相關性。

本文以露天煤礦為研究區域，首先分析了哨兵應用平臺(sentinel application platform，SNAP)偏移追蹤工具中不同參數對偏移結果的影響，其次計算了研究區域的強度互相關的分布，針對互相關值較低的區域，通過在偏移追蹤工具源代碼中增加基于紋理特征的偏移追蹤方法，將兩種方法相結合進一步分析了露天礦的大梯度形變。

1 研究區概況和數據源

1.1 研究區概況

平朔礦區地處山西省西北部朔州市平魯區，地理坐標范圍為東經112°10′-113°30′，北緯39°23′-39°37′，地貌類型為黃土低山丘陵，海拔高度1.3 km～1.4 km.礦區東西寬18 km，南北長21 km，總面積380 km2，主要分布有三座露天煤礦和四座井工礦，其中安家嶺礦(AJL)、安太堡礦(ATB)和東露天礦(DLT)為露天煤礦，如圖1所示。

圖1 研究區域Fig.1 Overview of the study area

本文選取安太堡露天煤礦作為研究區域，安太堡露天礦東西長5.7 km，南北長3.1 km，該露天礦的開采方向為自西向東，分布有排土場、開采區和剝離區。目前安太堡露天礦主要剝離位置在礦坑東北部，排土作業區位于礦坑的西北側，各作業區之間分布有運輸道路。

1.2 數據源

Sentinel-1A IW SLC影像用于提取偏移量。SLC是1級單視復數影像標準產品,空間分辨率為2.7 m×22 m至3.5 m×22 m，像素大小為2.3 m×14.1 m.標準SNAP8.0軟件在POT計算中可使用地距檢測(ground range detected，GRD)影像，GRD是聚焦的SAR數據經過去除相位信息、多視處理并依據WGS84橢球模型投影到地距向的1級產品，投影之后影像像素近似方形，斑點噪聲減弱同時也使空間分辨率降低，像素大小約為10 m，距離向與方位向分辨率為20 m×22 m[9].利用SLC數據可得到視線方向(line-of-sight，LOS)向位移，與GRD數據計算的地距向位移相比可較好地分析礦區地面沉降。

2 改進的偏移量追蹤方法

SNAP8.0中偏移追蹤工具主要處理流程包括：初始化參數、創建主影像GCP、確定副影像GCP、計算偏移量、計算偏移速率、偏移值平滑和補洞處理。偏移追蹤算法主要是基于SAR影像的強度信息計算匹配窗口內的NCC值，由于不同區域后向散射性差異明顯，未發生地形變化或僅存在少量形變的區域的影像塊間相關性高；地形發生變化區域的影像塊相關性降低。另一方面，強度影像中存在明顯的斑點結構，斑點結構信息反映了地面不同目標的散射特性和地形起伏。改進的偏移追蹤方法在匹配過程中引入紋理特征。紋理特征能夠較好地反映地形信息，對SNAP8.0中偏移追蹤工具的修改包括以下3部分：

1) 確定互相關值較低的區域。首先設置一定尺寸的匹配窗口，根據窗口大小在主、副影像中均勻分布影像控制點。在匹配窗口移動前，利用窗口中主、副影像的強度值計算互相關值。遍歷所有影像控制點得到整幅影像的相關性值。

2) 計算匹配窗口內的紋理特征值并構建紋理特征向量。根據研究區域設置互相關閾值，SAR影像強度互相關值較高的區域可認為不存在大量級沉降或地形變化；互相關值較小的區域可認為存在地形變化，通過匹配窗口逐像素移動計算窗口內的紋理特征值。

利用灰度共生矩陣(gray level co-occurance matrix，GLCM)統計像素灰度對組合頻率。由于影像強度值分布于零至幾百之間，因此首先選擇累積頻率2%和98%的灰度值作為原圖像灰度閾值，將圖像灰度等級壓縮為8級，在匹配窗口內分別計算相鄰像素點之間4種空間關系對應的GLCM.根據HARALICK et al[10]提出的14種紋理選擇其中對比度、同質度、角二階矩、平均值、標準差、熵共6種紋理特征描述子，計算每個GLCM對應的紋理特征值。不同紋理特征值分別組成特征向量：

(1)

(2)

(3)

式中：α(fi)表示在遍歷過程中同一紋理特征的標準差。

通過遍歷所有位置并記錄每個位置計算的Df，選擇Df最小值所在位置作為最佳匹配位置并估計偏移量和偏移速率。

3 結果與分析

3.1 偏移量追蹤參數分析

Sentinel-1A SLC影像預處理步驟包括：子帶分離、軌道數據更新、熱噪聲去除、圖像校準、Deburst、DEM輔助配準和多視處理。偏移計算前需要選擇合適的參數值。文中選擇四個參數進行分析，包括匹配窗口、過采樣因子、互相關閾值和速率閾值，通過控制單一變量對各參數計算結果進行統計分析，GCP格網間距、偏移值平滑和孔洞填充半徑均采用默認值，文中設置的4個參數值如表1所示。

表1 參數分析中設置的參數值Table 1 Parameter values set in parameter analysis

經過偏移追蹤工具處理后共輸出6 878個像素點，分別對設置的4種參數計算結果的最大值(max)、最小值(min)、平均值(mean)、標準差(std)和變異系數(Coefficient of Variation，CV)進行統計分析，其中變異系數為標準差與平均值之比，其倒數可表示信噪比。

過采樣因子影響可探測的最小位移。由表2可知過采樣因子增大時，會影響局部區域的計算結果，表現為最大值和最小值隨著采樣因子的增大而增大，但是平均值和標準差變化小于0.10 m/d，變異系數的變化小于0.10，表明過采樣因子對總體偏移量計算影響較小，可根據精度要求和計算效率選擇合適的過采樣值。

表2 不同過采樣因子計算結果的統計值Table 2 Statistics of calculation results for different oversampling factors

匹配窗口用于在搜索區域內按順序移動并計算互相關值。較小的匹配窗口對局部噪聲更加敏感，因而計算結果易受強散射體影響；過大的匹配窗口可探測到大尺度位移，同時局部的形變現象也被掩蓋。由于不同匹配窗口計算的結果數值差異較大，于是使用ln(x)函數對統計結果進行縮放得到了圖2.由圖2可知,隨著匹配窗口的增大，偏移速率最大值迅速增大，從而導致整體計算結果的平均值、標準差增大，變異系數逐漸減小并趨于平穩，表明增大匹配窗口可以抑制局部噪聲，提高信噪比；最小值變化較小也反映出偏移追蹤結果因不同目標的后向散射強度的互相關性而異。窗口大小為128像素、256像素和512像素時對應的標準差分別為3.87 m/d、設置速率閾值可以降低誤匹配結果，大于速率閾值的像素點將會被剔除，并通過空間插值的方法重新估計其偏移速率。與匹配窗口統計結果類似，由圖3可知隨著速率閾值的增大，最大值、平均值和標準差都呈現出不同的增長，圖3中各曲線的增長趨勢可以分為三部分，即速率閾值在小于1.5 m/d時，除了最小值外，各項統計值均呈現快速增長趨勢，過小的閾值也會抑制偏移結果；速率閾值在1.5 m/d～5 m/d時，各項統計值變化相對平穩，起伏不大；當速率閾值大于5 m/d時，起伏明顯，且設置的閾值大于9 m/d時，各項統計值波動較大。5 m/d作為速率閾值可以減小對偏移結果的過度限制，也能剔除誤匹配點，因此本文將其作為速率閾值。

圖2 不同的匹配窗口計算結果統計值Fig.2 Statistics of calculation results for different registration windows

12.07 m/d、37.92 m/d，對于本研究區域這三種窗口不建議用于偏移量計算。考慮到多視處理后SLC影像像素大小為13.98 m，大小為64像素的匹配窗口內像素空間相關性較小，因此研究中選擇的匹配窗口大小為32像素。

圖3 不同的速率閾值計算結果統計值Fig.3 Statistics of calculation results for different velocity thresholds

互相關閾值用于匹配窗口移動至某一位置時進行判斷，從表3可以看出互相關閾值對偏移結果影響有限，各項統計值變化不大，參考文獻[12]和[13]中相關性閾值設置為0.2和0.25，本文中計算了每個GCP點所在位置的互相關值，整個露天礦區的互相關平均值為0.63，礦坑內部的互相關值范圍為-0.01～0.2，為保證結果的可靠性，后續計算中設置互相關閾值為0.3.

表3 不同的互相關閾值計算結果統計值Table 3 Statistics of calculation results for different cross-correlation thresholds

3.2 露天礦地面沉降分析

通過上述對偏移追蹤參數的分析，選擇過采樣因子為16，匹配窗口大小為32像素，速率閾值為5.00 m/d，互相關閾值為0.3，將兩種偏移追蹤方法相結合，對覆蓋ATB露天礦的連續四景升軌的SLC影像進行處理和分析。圖4中(a)、(b)、(c)是經過計算得到的偏移速率；(d)、(e)、(f)為距離向偏移量，(g)、(h)、(i)為方位向的偏移量。偏移速率最大值為4.04 m/d，主要位于ATB礦開采區，該區域距離向和方位向偏移也達到最大，分別為-40.28 m～-18.87m、-32.28 m～16.55 m；同時剝離區存在明顯的沉降區域，距離向和方位向偏移量分別為-14.98 m～-1.18 m和-18.87 m～-4.42 m；排土場偏移量較小，距離向和方位向偏移量為-8.90 m～9.20 m、-0.20 m～11.4 m.

圖4 ATB礦偏移量和偏移速率Fig.4 ATB mine offset and velocity

進一步利用偏移追蹤方法對同一時期的降軌數據進行了處理，并用于解算ATB礦的二維形變量。由于升軌和降軌數據采集時間并不完全一致，文中采用了SAMSONOV et al[14]提出的補償因子的方法對偏移計算結果進行了補償，同時偏移計算的點位也不一致，文中通過反距離加權法將降軌數據計算的偏移量進行了重采樣，保證偏移計算結果中像素坐標一致。通過升軌和降軌偏移追蹤結果解算了ATB露天礦的垂直向和東西向的形變量。圖5為垂直向的形變量，A～F是在剝離、開采和排土區識別的大梯度形變的區域。

圖5 ATB礦垂直向形變量Fig.5 Vertical deformation of ATB mine

圖6為A～F區域的垂直向形變量直方圖，各形變區域的形變量集中在-30.00 m～20.00 m之間；A、D、E、F區域同時存在沉降與上升，B、C區域以上升為主；形變量最大值出現在F、B區域，F區域位于剝離區，其沉降量最大值為-42.41 m，B區域位于排土場東南側，形變量為31.87 m.

圖6 各形變區域的形變量 Fig.6 Deformation of each deformation area

由升、降軌解算的東西向的形變量如圖7所示，東西向形變較為明顯的區域與垂直向分布基本一致，A、D、E區域的形變量范圍為-29.40 m～39.52 m，B、C、F區域形變量為負值，表現為向西移動；E、F區域東西向的形變量達到最大，分別為42.80 m，-50.41 m.

圖7 ATB礦東西向形變量Fig.7 East-West deformation of ATB mine

進一步與遙感影像相結合分析了像素點的偏移方向和偏移速率，圖8中箭頭的方向表示偏移方向，箭頭的尺寸表示偏移速率的大小。由圖8可知，偏移速率較大的區域位于開挖或排土作業較為頻繁的區域，如A～E區域均存在排土和剝離作業；像素點的偏移方向具有局部特性，如形變區F像素點向西北方向偏移，該區域主要存在邊坡，主要由邊坡的沉降引起。由于高強度的采動導致地形發生變化，使得地面目標的后向散射強度也發生改變，因此偏移追蹤結果存在一定誤差。

圖8 ATB礦像素偏移速率和偏移方向Fig.8 ATB Offset tracking velocity and direction

4 結論

以露天煤礦為研究區域，運用SAR偏移量追蹤技術計算了研究區域的大梯度形變，主要結論有：

1) 匹配窗口是影響偏移追蹤計算結果的主要因素，較大的匹配窗口可抑制噪聲，反之則易受噪聲影響；本文計算結果也發現偏移速率閾值會影響被剔除點的區域和插值結果， 文中通過對比不同參數的計算結果選擇合適的速率閾值，但是較可靠的速率閾值仍需通過實地調查得到。

2) 本文偏移追蹤計算的露天煤礦形變量達到分米級以上，盡管利用DEM輔助配準法可以降低軌道誤差、地形對偏移結果的貢獻，但露天煤礦局部區域的偏移量仍然達到幾十米，其原因是該區域存在采掘、排土作業，同時哨兵數據經過多視處理像素大小為13.9 m，因此偏移追蹤結果也與影像分辨率有關。

3) 與井工開采引起的地面形變相比，露天煤礦形變局部特性明顯且較為復雜，需要結合形變模型進一步研究；露天煤礦區的邊端幫、運輸線路、邊坡等散射性較強且為線狀目標，基于規則窗口的偏移追蹤法并不完全適用于不規則形變區域的監測。