差分雷达干涉技术(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)
式中：fi表示特征向量；i=1,2,3,…,6表示不同的纹理特征描述子；
表示同一描述子但不同空间关系对应的特征值。
3) 计算纹理特征向量间的距离。由匹配窗口内的主、副影像计算的纹理特征向量
和
根据SMITS et al[11]使用的方法计算纹理特征向量之间的距离Df：
(2)
式中：di表示纹理特征向量
和
之间的距离，计算方式如下式所示：
(3)
式中：α(fi)表示在遍历过程中同一纹理特征的标准差。
通过遍历所有位置并记录每个位置计算的Df，选择Df最小值所在位置作为最佳匹配位置并估计偏移量和偏移速率。
3 结果与分析
3.1 偏移量追踪参数分析
Sentinel-1A SLC影像预处理步骤包括：子带分离、轨道数据更新、热噪声去除、图像校准、Deburst、DEM辅助配准和多视处理。偏移计算前需要选择合适的参数值。文中选择四个参数进行分析，包括匹配窗口、过采样因子、互相关阈值和速率阈值，通过控制单一变量对各参数计算结果进行统计分析，GCP格网间距、偏移值平滑和孔洞填充半径均采用默认值，文中设置的4个参数值如表1所示。
表1 参数分析中设置的参数值
Table 1 Parameter values set in parameter analysis
参数参数值过采样因子2、4、8、16、32、64匹配窗口尺寸/pixels32、64、128、256、512速率阈值/(m·d-1)0.05、0.1、0.5、1.5、3、5、7、9互相关阈值0.1、0.2、0.3、0.4、0.5、0.6
经过偏移追踪工具处理后共输出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
过采样值最大值/(m·d-1)最小值/(m·d-1)平均值/(m·d-1)标准差/(m·d-1)变异系数221.130.003.653.911.07421.440.003.623.901.08821.410.103.653.901.071620.850.013.623.861.073220.860.023.613.871.076423.130.023.603.851.06
匹配窗口用于在搜索区域内按顺序移动并计算互相关值。较小的匹配窗口对局部噪声更加敏感，因而计算结果易受强散射体影响；过大的匹配窗口可探测到大尺度位移，同时局部的形变现象也被掩盖。由于不同匹配窗口计算的结果数值差异较大，于是使用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
互相关阈值最大值/(m·d-1)最小值/(m·d-1)平均值/(m·d-1)标准差/(m·d-1)变异系数0.10.230.010.050.030.530.20.230.000.050.030.530.30.240.000.050.030.540.40.220.010.040.020.530.50.260.010.050.040.52
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) 与井工开采引起的地面形变相比，露天煤矿形变局部特性明显且较为复杂，需要结合形变模型进一步研究；露天煤矿区的边端帮、运输线路、边坡等散射性较强且为线状目标，基于规则窗口的偏移追踪法并不完全适用于不规则形变区域的监测。
[1] STROZZI T,WEGMULLER U,WERNER C,et al.Potential of a C-band SAR mission with 12-day repeat cycle to derive Ice surface velocity with interferometry and offset tracking[C]∥IEEE International Geoscience and Remote Sensing Symposium,2007:4229-4232.
[2] 徐小波,马超,单新建,等.联合DInSAR与Offset-tracking技术监测高强度采区开采沉陷的方法[J].地球信息科学学报,2020,22(12):2425-2435.
XU X B,MA C,SHAN X J,et al.Monitoring ground subsidence in high-intensity mining area by integrating DInSAR and offset-tracking technology[J].Journal of Geo-information Science,2020,22(12):2425-2435.
[3] STROZZI T,LUCKMAN A,MURRAY T,et al.Glacier motion estimation using SAR offset-tracking procedures[J].IEEE Transactions on Geoscience and Remote Sensing,2002,40(11):2384-2391.
[4] PATTN F.DInSAR and coherence tracking applied to glaciology:the example of shirase glacier[M].John Wiley & Sons,Ltd,1999.
[5] ERTEN E,REIGBER A,HELLWICH O,et al.Glacier velocity monitoring by maximum likelihood texture tracking[J].IEEE Transactions on Geoscience and Remote Sensing,2009,47(2):394-405.
[6] WANG T,JONSSON S.Improved SAR amplitude image offset measurements for deriving three-dimensional coseismic displacements[J].IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing,2015,8(7):3271-3278.
[7] LIU J,ISKANFDE M.Adaptive cross correlation for imaging displacements in soil[J].Journal of Computing in Civil Engineering,2004,18(1):46-57.
[8] CAI J,WANG C,MAO X,et al.An adaptive offset tracking method with SAR images for landslide displacement monitoring[J].Remote Sensing,2017,9(8):830.
[9] European Space Agency.Sentinel-1 product types and processing levels[EB/OL].https://sentinel.esa.int/web/sentinel/user-guides/sentinel-1-sar/product-types-processing-levels/level-1.
[10] HANMUGAM K,DINSTEIN I.Textural features for image classification[J].Studies in Media and Communication,1973,SMC-3(6):610-621.
[11] SMITS P C,ANNONI A.Toward specification-driven change detection[J].IEEE Transactions on Geoscience and Remote Sensing,2000,38(3):1484-1488.
[12] POLCARI M,TOLOMEI C,BIGNAMI C,et al.SAR and optical data comparison for detecting co-seismic slip and induced phenomena during the 2018 Mw 7.5 Sulawesi Earthquake[J].Sensors,2019,19(18)：3976.
[13] SUN L,JANPETER M.Evaluation of the use of sub-pixel offset tracking techniques to monitor landslides in densely vegetated steeply sloped areas[J].Remote Sensing,2016,8(8):659.
[14] SAMSONOV S V,D’OREYE N.Multidimensional small baseline subset (MSBAS) for two-dimensional deformation analysis:case study mexico city[J].Canadian Journal of Remote Sensing,2017,43(4):318-329.
