基于四叉树滤波的InSAR矿区形变监测与分析

杨小鸽，王 华，柏 俊，吴希文，李廷俊，柏玉超，皮廷亮

（1.广东工业大学 土木与交通工程学院，广东 广州 510006；2.广州海川信息科技有限公司，广东 广州 510300；3.云南华联锌铟股份有限公司，云南 文山 663700）

矿区开采导致的地表形变易造成矿区塌陷，破坏矿区建筑设施，影响工作人员的生命安全，因而矿区形变监测是矿区开采的一项重要工作[1-3]。相比传统矿区形变监测手段，InSAR不仅能满足毫米级的精度需求，同时凭借其大尺度和高分辨率的特点，InSAR还能同时获得多个矿区米级空间分辨率的形变图，对于矿区形变监测具有明显的优势[4-6]。时间序列InSAR能有效地监测矿区形变过程，但矿区周围区域植被茂密，容易导致时变去相干，形变结果中会存在大量的离散噪声点。对于露天开挖的矿区，地表属性变化较快，矿区内长期稳定点很少，难以获得长期稳定的永久散射体和分布式散射体。为了获得高分辨率的InSAR形变时间序列，需要对原始短时间基线干涉图的解缠结果进行滤波，从而提高InSAR形变时间序列结果的质量。由于能够在降低数据量的同时保留形变场的主要特征，四叉树采样方法被广泛地用在InSAR形变建模中[7-8]。在四叉树采样思想的基础上，Peng等[9]提出了一种四叉树滤波方法，能够有效消除噪声点，并保留形变特征，在2016年新西兰凯库拉地震的偏移量分析中得到了较好的应用[9-10]。

本文研究的云南某矿区是一个露天开采矿，属亚热带季风气候，年降水量大，矿区周围植被较茂密。矿区范围内有采矿场、选矿车间、排石场3个主要区域,目前主要采用陡帮作业采剥工艺。该矿区形变范围小、形变梯度大，与地震形变的震中区具有相似的特性，而对分辨率的要求比地震更高，在高分辨率下进行相位解缠往往会带来较多的噪声点。图1为20200620-20200702干涉对的解缠相位图，其中黑线表示3个典型剖面的位置(见图4)，黑框表示地面观测点的分布范围(见图5)。从图1可以看出该矿区干涉图受时变去相干影响而含有较多噪声。为了提高InSAR形变结果的质量，本文对该矿区2019年10月至2021年3月的InSAR解缠结果进行四叉树滤波，并通过时间序列分析改正轨道误差、大气延迟误差，最终获得研究区域干净的时间序列形变结果。

图1 20200620-20200702干涉对的解缠相位图Fig.1 Unwrapped phase of the 20200620-20200702 interferogram

图2 滤波效果图和剔除的像素点(a-b分别为粗差阈值为2倍中误差，最小划分窗口为全局64×64的滤波后解缠相位图和剔除的像素点；c-d为粗差阈值为2倍中误差，最小窗口为全局8×8的滤波效果图对应图；e-f为粗差阈值为2倍中误差，近场8×8和远场64×64的滤波效果图对应图；g-h为最小窗口为近场8×8和远场64×64，粗差阈值为1.5倍中误差的滤波效果对应图；i-j为最小窗口为近场8×8和远场64×64，粗差阈值为3倍中误差的滤波效果对应图)Fig.2 Performance of quadtree filtering.(a and b are the filtered unwrapped phase and the removed pixels, respectively, with the gross error threshold of 2 times of the standard deviation and the minimum division window of 64×64; c-d are the corresponding results with the gross error threshold of 2 times of the standard deviation and the window of 8×8; e-f are the corresponding results with the gross error threshold of 2 times of the standard deviation, but with different minimum window size in the near (8×8) and the far field (64×64) ; g-h and i-j are the corresponding results with the gross error threshold of 1.5 and 3 times of the standard deviation, respectively, and the minimum window size as used in e-f.)

图3 矿区时间序列形变结果(远离卫星为正)Fig.3 Time series deformation results of mining area (Positive indicates away from the satellite)

图4 矿区剖面形变时空演化图(远离卫星为正)Fig.4 Spatio-temporalevolution map of deformation along profile(Positive indicates away from the satellite)

图5 矿区形变结果验证图(a)测量机器人监测点位图；(b)InSAR与32个测量机器人点的时间序列形变结果对比分析图，颜色表示差值的大小；(c-e)为P1、P2、P3三个点的InSAR与测量机器人时间序列形变结果对比图Fig.5 Validation between InSAR and total stations (a) Location of 32 total stations in black dots; (b) Scatter plot of the LOS displacement time series from InSAR and total stations, the color-coded circles represent LOS displacement difference of these two means; (c-e) Time series of three points P1, P2 and P3 from InSAR and total stations

1 四叉树滤波

1.1 四叉树滤波原理

四叉树滤波能剔除影像中离散的噪声点。该滤波方法首先将影像进行递归分级，每次将区域分成4个大小相等的子区域，计算子区域内像素值的中误差或减去多项式模型后的中误差，如果标准差超过某一给定的阈值，则继续往下将各子区域进行分割，直到各子区域的中误差满足阈值要求，或者该子区域已划分至预设的最小划分窗口。在四叉树采样中，当各子区域划分达到上述标准后则输出该窗口的平均值(或中值)，而四叉树滤波则是按照一定原则剔除该窗口残差超过一定阈值的像素，并将剩余的所有像素值输出[10]。通常，最小划分窗口和粗差阈值对四叉树滤波的影响较为显著。下文将比较不同的最小划分窗口和粗差阈值对滤波效果的影响。

1.2 最小划分窗口

通常情况下，较大的最小划分窗口有利于剔除离散的噪声。考虑到矿区形变范围小且形变量大，如果对整个解缠图采用同样的最小划分窗口，则会导致矿区形变中心的有效观测数据损失。采用粗差阈值为2倍中误差，全局最小划分窗口为64×64，对图1进行四叉树滤波的效果如图2(a-b)所示，其中图2(a)为滤波后的结果，图2(b)为滤除的点。可见在近场由于矿区形变梯度大，而最小划分窗口过大，在滤波后的解缠图中，矿区中心只留下小部分像素点(见图2(a))，矿区形变中心的大量有效形变点被滤除(见图2(b))。

当最小划分窗口较小时，各子区域划分得比较小，各子区域内的形变梯度相对较小，子区域内像素值的中误差会比大窗口时小，导致近场滤除的点比较少，但此时远场的噪声也会随着被划分至不同子区域，导致各子区域的中误差相对较小，滤波后仍然会存在大量的噪声。同样采用粗差阈值为2倍中误差，图2(c-d)显示了全局最小划分窗口缩小为8×8时对图1进行滤波的效果。由图2(d)可见滤除的点很少，近场的形变点基本保留了，但由图2(c)可见滤波后在远场仍有许多噪声。

由上述比较可知，合理调节最小划分窗口是获得理想滤波效果的一个关键步骤。理论上近场的形变梯度比较大，应采用较小的滤波窗口，而远场形变梯度小，则应采用较大的滤波窗口。将矿区中心2 km范围内的区域作为近场形变中心，将2 km以外区域视为远场，同样采用粗差阈值为2倍标准差，将近场的最小划分窗口设置为8×8，而远场则设置为64×64，滤波效果如图2(e-f)所示。通过对比滤波后的相位解缠图(见图2(e))与原始的相位解缠图(见图1) ，可以看出对近场与远场设置不同的最小划分窗口既能有效剔除噪声点，又不损失近场的形变数据。

1.3 粗差阈值

上节对最小窗口的分析中均采用2倍中误差为粗差阈值，但在实际数据处理中，粗差阈值需根据影像噪声的整体离散程度进行设置，阈值过高会导致影像粗差点滤除不干净，而阈值过低则会将大量非粗差点滤除。将最小划分窗口均设置为远场64×64和近场8×8，图2(e-f)、图2(g-h)、图2(i-j)分别表示粗差阈值为2倍、1.5倍与及3倍中误差的滤波效果图。对比发现，采用1.5倍中误差作为粗差阈值虽然能将几乎所有的噪声点剔除，但是远场与近场的一些非噪声点也被滤除(见图2(g-h))；采用3倍中误差作为粗差阈值虽然能保留形变中心的像素，但是滤波后远近场均存在较多的残余噪声点(见图2(i-j)) ；而采用2倍中误差(见图2(e-f)) 则能在基本保留形变中心有效像素的基础上，远场的滤波效果与采用1.5倍中误差(见图2(g-h)) 阈值滤波效果差异不大。

2 InSAR时间序列解算

2.1 基线选择

本文选取45景2019年10月至2021年3月的Sentinel-1A影像数据。矿区为一露天矿区，矿区中心范围几乎无植被，且在选择时间范围内空间基线均小于150 m，因此在组成干涉对时优先考虑时间基线。由于矿区中心的形变速率快、范围小，且Sentinel-1A为C波段的数据，波长比较短，如果时间基线过长，则会导致干涉图条纹混叠，无法解缠，因此本文的所有时间基线都选取为12 d，即为相邻的图像组成干涉对，一共组成44个干涉对。

2.2 DInSAR数据处理

首先使用ISCE软件[11]进行DInSAR数据处理，干涉图采用Goldstein滤波方法进行滤波[12]，根据形变中心的干涉图条纹密集情况以及远场条纹的清晰情况将Goldstein滤波的强度在0.3～0.8之间进行调整，地形相位去除以及解缠结果的地理编码都采用1角秒的美国SRTM DEM[13]。为使最终方位向与距离向地面分辨率一致，使用方位向与距离向为1:4对数据进行多视处理。干涉图的解缠方法采用枝切法[14]，并通过闭合环自动检测解缠误差[15]。对由于去相干导致的孤岛采用人工桥接相干性良好的区域进行解缠[15]，对形变中心条纹密集区域进行人工检查解缠结果的正确性。

2.3 InSAR时间序列计算

在解缠完成后对解缠结果采用四叉树滤波方法进行滤波，其中近场和远场分别采用2×2和8×8的最小划分窗口，粗差阈值为2倍中误差。接着，将滤波后的解缠图进行地理编码，采用 π −RATE软件计算矿区的时间序列结果[15]。时间序列计算中的轨道误差采用网络法的双线性拟合模型改正。与地形相关的InSAR大气延迟误差改正与轨道误差改正类似，采用线性拟合模型拟合干涉相位与地形的关系，将拟合后的大气延迟误差从解缠后相位中减掉。由于矿区范围比较小，在远场对与地形不相关的大气延迟误差采用空间域低通滤波，获得远场大气延迟部分低通分量，然后将其进行空间插值获得矿区近场的大气延迟误差，然后将插值得到的与地形不相关的大气延迟误差减掉。在时域上，各个干涉图的相干性不一致，且矿区形变的速度快，干涉相位正负交替，在时域上无规律，因此本文对InSAR时间序列不进行时域滤波。

3 InSAR结果与分析

3.1 矿区形变时间序列结果

本文采用四叉树滤波方法对相位解缠结果进行滤波后，采用 π −RATE软件处理得到矿区最终的时间序列形变结果，如图3所示。由图可见从2019年10月至2021年3月，矿区地表形变剧烈，尤其是2020年6月至2020年9月时间范围内更为显著。矿区形变在时域上整体呈波浪状，时沉降时抬升，无明显规律，总体与矿区的开挖与修复工作进度相关。以2019年11月5日为例，从图3可清晰看到，在空间上，矿区形变显著区域大致可分为3个区域，矿区中心区域既有抬升也有沉降。

3.2 矿区形变时空演化分析

矿区在研究时段内发生了剧烈的地表形变，整个矿区间隔12 d的最大沉降增量为67.3 mm，发生在2020年8月31日，最大抬升增量为79.4 mm，发生在2019年11月5日。

本文选取3个形变中心区域作剖面线以进行时域分析，剖面线位置如图1所示，沿剖面线每30 m取一个点。图4为3个形变中心A、B、C区域的剖面时空演化结果，其中上图的折线表示每个时段该区域沿剖面所有点的平均形变量，下图各点的颜色表示不同时段沿剖面所有点的形变值。由这3个剖面图可见，3个形变中心区域的形变增量在时域上呈波浪状，时沉降时抬升，A区域和C区域在时域上均发生了严重的形变，形变起伏波动大于B区域。A区域的形变范围最大，长约2 km，宽约1 km。对比3个剖面图的平均沉降量折线图，在大多数情况下B区域与C区域的形变趋势相同，而与A区域的趋势相反。这是由于在矿区交替出现挖方和填方的施工，导致矿区抬升与沉降交替出现。

A剖面的最大抬升量为52.7 mm，最大平均抬升量为35.0 mm,发生在2019年11月5日，最大沉降量为36.1 mm，最大平均沉降量为27.0 mm，发生在2020年12月29日；B剖面的最大抬升量为20.3 mm，最大平均抬升量为10.2 mm，发生在2020年10月18日，最大沉降量为19.6 mm，最大平均沉降量为13.2 mm，发生在2020年8月31日；C剖面的最大抬升量为36.5 mm，最大平均抬升量为15.2 mm，发生在2020年10月18日,最大沉降量为53.3 mm，最大平均沉降量为31.7 mm，发生在2020年8月31日。

3.3 矿区形变结果验证

该矿区开采区北部布设了一些测量机器人，范围如图1的黑框所示。本文将32个点位41个时段测量机器人的测量结果转换到InSAR视线(line-of-sight，LOS)方向，点位如图5(a)中黑点所示。图5(b)中彩色的圆圈代表这32个点所有时段的结果与InSAR时间序列形变结果的差值。图5(c)～(e)为选取其中3个点P1、P2、P3(点位见图5(a))的测量机器人与InSAR形变结果的对比，两种手段的形变结果趋势一致，大部分时段差异均小于20 mm。对所有点所有时段，两种测量手段形变差值的平均值为0 mm，标准差为10.9 mm，这与测量机器人所能达到的测量精度一致，验证了本文InSAR时间序列形变结果的可靠性。

4 结论

本文将四叉树滤波应用于云南某矿区，并计算了该矿区的地表形变时间序列。研究结果表明，当形变中心的形变梯度较大时，近场与远场采用不同最小划分窗口的四叉树滤波，既能保留近场的形变点，又能去除大量的噪声，从而获得干净的形变结果。通过对该矿区形变时间序列的时空特征分析显示，该矿区在2019年10月至2021年3月地表形变剧烈，形变在时域上总体呈波浪状，时沉降时抬升，在空间上，矿区形变显著区域大致可分为3个区域，并且A区域和B、C区域形变趋势相反，这是由于在矿区开采区与排土区同时存在挖方和填方，导致矿区抬升与沉降交替出现。整个矿区相隔12 d的最大沉降增量达到67.3 mm，最大抬升增量达到79.4 mm。通过对比测量机器人的形变监测结果，验证了本文InSAR结果的可靠性。鉴于InSAR的技术优势和该矿区还有很多区域没有布设地面监测点，建议今后进一步采用InSAR技术加强对矿区的监测，确保矿区安全施工。