一种像素偏移量追踪监测及校正方法

技术领域

本发明属于合成孔径雷达（SAR）像素偏移量追踪（offset-tracking）变形监测领域，涉及一种像素偏移量追踪监测及校正方法。

背景技术

地震、滑坡、采矿、冰川等造成的地表形变往往能达到几米甚至几十米，这种大量级形变会导致基于SAR影像相位获取形变信息的干涉测量InSAR技术的干涉相位完全失相干，最终导致利用InSAR技术获取的形变值比真实形变值严重偏小，甚至无法获取准确的形变信息。在此情况下，相关学者提出使用基于SAR影像强度信息的像素偏移量追踪技术来获取大量级形变。基于SAR影像强度信息的像素偏移量追踪方法来源于图像匹配的归一化互相关算法，该方法利用两幅影像的强度信息，通过归一化的相关性度量公式来计算两幅影像的相似度，并利用如下公式计算两幅影像的互相关系数

其中

偏移量追踪技术是利用配准后的两幅SAR影像的强度信息计算两幅影像的互相关系数矩阵，并找出互相关系数矩阵峰值（通常是互相关系数矩阵的最大值）的位置，互相关系数矩阵的最大值位置与影像中心在方位向和距离向的距离就是两幅影像在方位向和距离向的偏移量，最后将偏移量转换成形变。互相关系数矩阵直接关系到两幅影像的偏移量，而互相关系数矩阵是利用影像的强度信息通过数学计算直接得到，因此该方法的监测精度直接受影像的强度信息和互相关窗口的大小的影响。

作者在研究过程中发现，互相关系数矩阵可能会存在多个峰值，且利用最大的峰值计算得到的偏移量与实际形变差距较大。原因可能是由于季节变化或人类活动等原因导致两幅影像的强度信息发生剧烈变化，造成两幅影像的同名点误匹配，使互相关系数矩阵出现多个峰值，最终导致互相关系数矩阵的最大峰值位置与影像中心的距离偏大，使监测结果发生异常。

发明内容

本发明的目的是针对上述基于SAR影像强度信息的偏移量追踪技术存在的问题，提出一种像素偏移量追踪监测及校正方法，它能够去除因地表强度发生剧烈变化导致偏移量追踪结果异常，并有效解决因互相关系数矩阵存在多个峰值导致监测结果不准确的问题，使偏移量追踪技术能够获取准确的地表形变信息。

一种像素偏移量追踪监测及校正方法，包括如下步骤：

S1 将研究区域同一SAR卫星拍摄的两幅SAR影像进行配准并生成主影像和辅影像，利用经验方法预计该研究区域内最大的距离向偏移量r和方位向偏移量a；

S2 基于互相关算法对主影像和辅影像进行偏移量追踪，并基于距离向偏移量r和方位向偏移量a的约束获取某一像元p(i,j)的互相关系数矩阵C1，计算获得p(i,j)在距离向的偏移量；

S3 重复S2直至完成对主影像和辅影像的偏移量追踪，获取主影像和辅影像的距离向偏移量矩阵RA；

S4基于下式计算像元p(i,j)的距离向偏移量阈值Tr：

如果

获取像元p(i,j)对应的矩阵C1，找出矩阵C1中所有的峰值，保留大于预设上限值的峰值用于计算新的距离向偏移量，如果新的距离向偏移量Tr差值的绝对值仍大于或等于所述预设阈值，则将像元p(i,j)对应的各影像块中像元强度最大值使用其他统计值替代后重新提取矩阵C1，计算像元(i,j)的距离向偏移量；

重复S4直至新的距离向偏移量与距离向偏移量阈值Tr差值的绝对值小于所述预设阈值。

作为一种优选的实施方式，所述S2中，基于距离向偏移量r和方位向偏移量a的约束获取某一像元p(i,j)的互相关系数矩阵，计算获得p(i,j)在距离向的偏移量的方式为：

先基于互相关算法获取某一像元p(i,j)的互相关系数矩阵C，再从C中提取新的矩阵：

其中N为互相关系数矩阵的列数，M为互相关系数矩阵的行数，n为过采样系数；

以矩阵C1峰值所在的位置与矩阵C1中心在距离向的距离作为像元p(i,j)在距离向的偏移量。

作为一种优选的实施方式，所述预设阈值的大小根据偏移量追踪方法的精度确定。

作为一种优选的实施方式，所述预设阈值的大小为偏移量追踪方法的精度的3倍数值。

作为一种优选的实施方式，所述预设上限值为矩阵C1的均值。

作为一种优选的实施方式，利用各影像块中的像元强度均值替代影像块中的像元强度最大值。

作为一种优选的实施方式，所述S4中，如果保留的峰值只有1个，则直接使用该峰值计算新的距离向偏移量，重复S4直至新的距离向偏移量与距离向偏移量阈值Tr差值的绝对值小于所述预设阈值；如果重复次数超过预设重复上限值后，新的距离向偏移量与距离向偏移量阈值Tr差值的绝对值仍大于或等于所述预设阈值，则令该像素点的距离向偏移量为0。

作为一种优选的实施方式，所述S4中，如果保留的峰值超过1个，则分别计算各峰值与C1矩阵中心在距离向的距离offset

偏移量追踪的互相关系数矩阵可能会存在多个峰值，因峰值识别错误会导致监测结果与实际情况相差较大，此外地表强度发生剧烈变化也会导致其监测结果异常。本发明先使用经验知识确定研究区域内的偏移量阈值，通过对互相关系数矩阵进行处理保证获取的偏移量不会出现与实际情况相差较大的情况，再利用峰值函数找出互相关系数矩阵内大于互相关系数均值的峰值及其位置，如果峰值只有一个，则对用于计算该像元偏移量的影像块强度进行优化，利用优化后的影像块再次进行互相关计算取准确的形变信息，如果峰值不止一个，则找出与该处形变阈值最近的峰值及其位置，这种情况下无需再次进行互相关计算，在找出准确的偏移量的情况下能够节约大量的计算时间，大大提高了运算效率，并能够准确识别互相关系数矩阵峰值、去除因强度发生剧烈变化导致偏移量追踪的异常结果，在提高偏移量追踪监测方法的精度的同时提高运算效率，为其更好的应用于大量级形变监测提供技术支撑。

附图说明

图1是覆盖研究区域的两幅SAR影像强度图。

图2是用于计算像元p(620,320)偏移量的主辅影像块强度。

图3中（a）是像元p(620,320)的互相关系数矩阵C；（b）是设定偏移量阈值后的像元p(620,320)的互相关系数矩阵C1。

图4中（a）是使用两幅强度图进行偏移量追踪获取的距离向形变图；（b）是设定偏移量阈值后获取的距离向形变图。

图5中（a）是用于计算像元p(640,530)偏移量的主影像块强度；（b）是用于计算像元p(640,530)偏移量的辅影像块强度；（c）是像元p(640,530)的互相关系数矩阵C；（d）是像元p(640,530)的互相关系数矩阵C1。

图6是使用本发明的方法获取的距离向偏移量。

具体实施方式

下面结合附图对本发明的具体实施方式做出详细的说明。

实施案例选取陕西省榆林市大柳塔矿区52304工作面的竖直向变形监测作为实验对象，52304工作面的煤层开采厚度为6.45m，开采时间为2011年11月1日至2013年3月25日。

1）选取覆盖研究区域的2景TerraSAR-X影像为实验影像，影像距离向像元尺寸为0.91m，影像入射角为42.43°，影像获取时间为2013年1月26日、2013年2月6日，matlabR2020b为本发明影像处理软件；对2景TerraSAR-X影像进行裁剪并配准，配准后的两幅影像强度图如图1所示，影像大小为900*600像元，m=900，n=600。从图1中可以看出，黑色方框内的强度发生了剧烈变化。

2）利用GPS监测得到的研究区域内在2012年11月10至2013年4月10日期间最大的沉降值

式中，

3）将偏移量追踪技术的互相关窗口大小设置为64像元*64像元，对于像元p(i,j)，A中用于计算像元p(i,j)偏移量的影像块A1=A(i-32:i+32,j-32:j+32)，B中用于计算像元p(i,j)偏移量的影像块B1=B(i-32:i+32,j-32:j+32)，当i=620、j=320时，该点在研究区域中的位置如图1中黑色圆点所示，该点处于形变区域边缘，其形变应属于小量级形变，图2为用于计算像元p(620,320)偏移量的影像块强度；将过采样系数设置为8，对A1和B1进行互相关计算获取A1和B1的互相关系数矩阵C（大小为520*520，如图3中（a）所示），找出矩阵C峰值所在的位置（如（a）中黑色方框所示），矩阵C峰值所在的位置与矩阵C中心在距离向和方位向的距离分别为1.6875像元和15.4375像元，利用矩阵C使用原始方法计算得到像元p(620,320)的距离向和方位向偏移量为1.6875像元和15.4375像元，显然与实际情况不符。提取新的互相关系数矩阵C1=C(260-a*8: 260+a*8, 260-r*8: 260+r*8)，如图3中（b）所示，找出矩阵C1峰值所在的位置，矩阵C1峰值所在的位置与矩阵C1中心在距离向和方位向之间的距离为0.0625像元和0.3125像元，即像元p(620,320)在距离向和方位向的偏移量为0.0625像元和0.3125像元，这样计算得到的距离向和方位向偏移量就不会超过研究区域内距离向和方位向偏移量的两个阈值r和a。如果i≤32或i>m-32或j≤32或j>n-32，则直接令像元p(i,j)在距离向和方位向的偏移量均为零。按照上述方法对A和B的每一个像元进行偏移量追踪处理即可获取强度图A和B的距离向偏移量矩阵RA（大小为m*n）。图4中（a）是使用原始偏移量追踪方法获取的距离向形变图，（b）是使用设定偏移量阈值后获取的距离向形变图，从图中可以看出，使用原始偏移量追踪方法获取的距离向形变存在大量的异常值，而使用设定偏移量阈值的偏移量追踪方法能够大大减少距离向形变的异常值。

4）令i=33，j=33，此时i的取值范围为33~m-32的正整数，j的取值范围为33~n-32的正整数，对于像元p(i,j)，按如下公式计算该像元的距离向偏移量阈值Tr：

如果|Tr-RA(i,j)|<0.3，说明像元p(i,j)的距离向偏移量是准确的，不需要校正，则计算下一个像元p(i,j+1)的Tr，如果j>n-32，则计算像元p(i+1,33)的Tr；如果|Tr-RA(i,j)|≥0.3，说明像元p(i,j)的距离向偏移量不准确，需要校正，则执行如下步骤：

①提取A中像元p(i,j)用于偏移量追踪的影像块A1=A(i-32:i+32,j-32:j+32)，提取B中像元p(i,j)用于偏移量追踪的影像块B1=B(i-32:i+32,j-32:j+32)；

②计算A1和B1的互相关系数矩阵C，并提取新的互相关系数矩阵C1=C(260-a*8:260+a*8, 260-r*8: 260+r*8)，计算C1的均值mean_c；

③令temp=0，使用imregionalmax函数找出C1所有的峰值及所有峰值的位置，并去除峰值小于或等于mean_c的峰值，只保留峰值大于mean_c的峰值P及其位置，计算保留下来的峰值的位置与C1矩阵中心的在距离向的距离offset1，如果保留下来的峰值只有一个，则执行步骤④-⑤；如果保留下来的峰值不只一个，则执行步骤⑥；

④如果|Tr-offset1|<0.3，令RA(i,j)=offset1，执行步骤5）；如果|Tr-offset1|>=0.3，执行步骤⑤；

⑤如果temp>5，令RA(i,j)=0，执行步骤5）；如果temp<=5，令temp=temp+1，计算A1和B1的均值mean_A和mean_B，找出A1和B1最大值max_A和max_B，令max_A所在的像元的强度值等于mean_A，令max_B所在的像元的强度值等于mean_B，形成新的A1和B1；计算新的A1和B1的互相关系数矩阵C，并提取新的互相关系数矩阵C1=C(260-a*8: 260+a*8, 260-r*8:260+r*8)，计算出矩阵C1峰值与矩阵C1中心在距离向的距离offset1，执行步骤④；

⑥计算offset1中的各偏移量与Tr的差值Diff，Diff=|Tr-offset1|，找出Diff的最小值Rmin及Rmin对应的偏移量offset2，如果Rmin<0.3，令RA(i,j)=offset2，执行步骤5）；如果Rmin≥0.3，则执行步骤④-⑤；

以像元p(640,530)为例对步骤4）做进一步说明：RA(639,529)=0.0625，RA(639,530)=0.0625，RA(639,531)=0.0625，RA(640,529)=0.0625，RA(640,530)=3.1875，计算像元p(640,530)的距离向偏移量阈值Tr=0.0625，此时|Tr-RA(640, 530)|≥0.3，说明像元p(640,530)的距离向偏移量不准确，需要校正。

执行步骤①，提取A中像元p(640,530)用于偏移量追踪的影像块A1=A(608:672,498:562)，提取B中像元p(640,530)用于偏移量追踪的影像块B1=B(608:672,498:562)，影像块A1、B1的强度如图5中（a）和（b）所示。

执行步骤②，计算A1和B1的互相关系数矩阵C（520*520）（如图5中（c）所示），从图5的（c）中可以看出，矩阵C中最大峰值的位置距离矩阵中心较远，说明使用矩阵C获取的该像元的偏移量严重偏大，不符合实际情况；提取新的互相关系数矩阵C1=C(228: 292, 228:292)（如图5d所示），从图5的（d）中可以看出C1中存在多个峰值，计算C1的均值mean_c=0.0033。

执行步骤③，令temp=0，使用imregionalmax函数找出C1所有的峰值及所有峰值的位置，并去除峰值小于或等于mean_c的峰值，只保留峰值大于mean_c的峰值P及其位置，计算保留下来的9个峰值的位置与C1矩阵中心的在距离向的距离offset1=[3.1875; 0.1875;0.9375; -4.0625; 3.9375; 2.0625; -0.8125; 3.8125; -2.3125]。

执行步骤⑥，计算offset1中的各偏移量与Tr的差值Diff=[3.125; 0.125;0.875; 4.125; 3.875; 2; 0.875; 3.75; 2.375]，找出Diff的最小值Rmin=0.125，Rmin对应的偏移量offset2=0.1875，Rmin<0.3，令RA(640,530)=offset2，然后对下一个像元p(640,531)进行上述处理。

5）接下来对像元p(i,j+1)执行步骤4），直至完成像元p(m-32,n-32)的偏移量异常校正，至此可去除使用偏移量追踪技术获取的强度图A和B的监测结果因互相关系数矩阵存在多个峰值或因地表强度发生剧烈变化导致偏移量追踪结果的异常。校正后的距离向的形变图如图6所示。