孫 軍，張 錦

(太原理工大學(xué) 礦業(yè)工程學(xué)院，太原 030024)

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

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

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

1 研究區(qū)概況和數(shù)據(jù)源

1.1 研究區(qū)概況

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

圖1 研究區(qū)域Fig.1 Overview of the study area

本文選取安太堡露天煤礦作為研究區(qū)域，安太堡露天礦東西長5.7 km，南北長3.1 km，該露天礦的開采方向為自西向東，分布有排土場、開采區(qū)和剝離區(qū)。目前安太堡露天礦主要剝離位置在礦坑?xùn)|北部，排土作業(yè)區(qū)位于礦坑的西北側(cè)，各作業(yè)區(qū)之間分布有運(yùn)輸?shù)缆贰?/p>

1.2 數(shù)據(jù)源

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

2 改進(jìn)的偏移量追蹤方法

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

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

2) 計算匹配窗口內(nèi)的紋理特征值并構(gòu)建紋理特征向量。根據(jù)研究區(qū)域設(shè)置互相關(guān)閾值，SAR影像強(qiáng)度互相關(guān)值較高的區(qū)域可認(rèn)為不存在大量級沉降或地形變化；互相關(guān)值較小的區(qū)域可認(rèn)為存在地形變化，通過匹配窗口逐像素移動計算窗口內(nèi)的紋理特征值。

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

(1)

(2)

(3)

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

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

3 結(jié)果與分析

3.1 偏移量追蹤參數(shù)分析

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

表1 參數(shù)分析中設(shè)置的參數(shù)值Table 1 Parameter values set in parameter analysis

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

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

表2 不同過采樣因子計算結(jié)果的統(tǒng)計值Table 2 Statistics of calculation results for different oversampling factors

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

圖2 不同的匹配窗口計算結(jié)果統(tǒng)計值Fig.2 Statistics of calculation results for different registration windows

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

圖3 不同的速率閾值計算結(jié)果統(tǒng)計值Fig.3 Statistics of calculation results for different velocity thresholds

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

表3 不同的互相關(guān)閾值計算結(jié)果統(tǒng)計值Table 3 Statistics of calculation results for different cross-correlation thresholds

3.2 露天礦地面沉降分析

通過上述對偏移追蹤參數(shù)的分析，選擇過采樣因子為16，匹配窗口大小為32像素，速率閾值為5.00 m/d，互相關(guān)閾值為0.3，將兩種偏移追蹤方法相結(jié)合，對覆蓋ATB露天礦的連續(xù)四景升軌的SLC影像進(jìn)行處理和分析。圖4中(a)、(b)、(c)是經(jīng)過計算得到的偏移速率；(d)、(e)、(f)為距離向偏移量，(g)、(h)、(i)為方位向的偏移量。偏移速率最大值為4.04 m/d，主要位于ATB礦開采區(qū)，該區(qū)域距離向和方位向偏移也達(dá)到最大，分別為-40.28 m～-18.87m、-32.28 m～16.55 m；同時剝離區(qū)存在明顯的沉降區(qū)域，距離向和方位向偏移量分別為-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

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

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

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

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

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

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

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

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

4 結(jié)論

以露天煤礦為研究區(qū)域，運(yùn)用SAR偏移量追蹤技術(shù)計算了研究區(qū)域的大梯度形變，主要結(jié)論有：

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

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

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