差分干涉技术在变形观测中的应用
2013-04-23 来源：中国鸣网
引言
地面沉降监测是对地观测的方法之一，经典的方法是采用精密水准测量，观测一、二等水准网，然后计算平差，得到微小地面沉降的变化值。近代GPS技术开始融入到地面沉降监测工作中。然而，受到若干条件的限制，GPS测量的数据密度往往不够，而且目前其空间域的分辨率也难以满足高空间分辨率形变监测的要求。
星载雷达差分干涉测量（D-InSAR）技术是20世纪90年代发展起来的1项对地观测新技术，其是利用雷达图像的相位差信息，精确测量出图像上每个点的三维位置，用于提取地面目标点微小地形变化的技术，精度能达到厘米级甚至毫米级，干涉图质量较好时，可以从干涉图识别出毫米级的形变。
D-InSAR技术的特点是面积大、快速、准确、成本低，其实用性大大优于常规地面沉降监测技术。目前，SAR数据已覆盖了全球大部分区域，在无常规测量数据区域的地震监测中发挥着显著优势。笔者在深入研究D-InSAR技术的基础上，对在地面沉降监测应用中的研究，解决了D-InSAR技术中影像处理面临的一些关键问题，得到最终的地表形变结果。
1、D-InSAR的三轨法
通常，进行地面沉降差分干涉，需要获取被观测地区沉降前后的2幅影像，通过比较得出形变量。
在2幅影像干涉图中，同时包含有相位信息和地形信息。欲得到形变量，需除去其中的地形信息，因此，必须先太清楚地面的DEM值，被称为“二轨法”。
在难以获得地面DEM值时，可采用“三轨法”。三轨法需要用到同一区域的3幅影像，分别是形变前的2幅和形变后的1幅。然后，利用这3幅影像可生成2幅干涉纹图，其中，1幅呈现形变前的地形信息，另1幅呈现形变后的地形信息。同样的，对2幅干涉图进行匹配、重采样、去除平地效应和相位解缠，最后将2幅影像进行差分，便可得到地表沿视线方向的地形变化。
这方法的优点是无需地面信息，如DEM，对于一些无地形数据的变化监测尤为重要;缺点是相位解缠较复杂，其质量将影响到最终的形变结果。
2、D-InSAR技术在巴姆地区的应用
2.1 被观测数据的选取
针对不同的干涉应用需选取合适的SAR干涉影像对。伊朗东南部城市巴姆地区曾于2003年发生了强度为6.6级的大地震，为了利用D-InSAR技术获取巴姆地震引起的地面形变，欧空局（ESA）提供了7幅ENVISAT ASAR数据。其中，4幅为降轨数据, 巴姆古城位于图像中部，Track号为120;3幅为升轨数据, 巴姆古城位于图像的近距端，Track号为385。
进行干涉处理，要求像对必须相干，因此，选取像对时必须考虑2个因素：a）临界基线距的限制。过短的基线将无法监测出形变变化量，过长的基线将导致2幅影像间的相关系数过低而无法形成干涉;b）存在时间间隔的影像。对于地形变监测为目的的干涉测量，其时间基线的选取要尽可能地短，以减小时间的去相关影响。
利用Eoli-sa软件进行数据选取时，根据被观测区域的经纬度和时间，查找合适的干涉数据。其中，Orbit确定轨道，景号（Frame）确定雷达照射区域位置，航迹号（track）确定成像时间。用于干涉测量的SAR数据必须具有相同的 track和 frame。
从中看出，形变前的影像9192和6687以及形变后的影像9693具有相同的Track号。6687的Frame为3020，即中心结点号为3020，该幅影像的结点跨度3011～3029。因为，ENVISAT卫星的轨道被分割成7200个结点，通过与影像中心相近的400个结点来标识SAR的Frame号。因此，2个相邻的Frame用18个结点来标识。同理，9192影像的结点跨度3004～3021，即这2幅影像有共同的结点，即有共同区域，进而得出9693与9192也有共同的区域。
2.2 实验过程
笔者采取三轨法进行差分处理，分别利用9192和6687构成的像对，9192和9693构成的像对，得出地震形变前后的地形，然后，再对这2个像对影像进行差分，获得地面点高程的形变信息。差分的关键步骤包括干涉影像的配准、干涉图的滤波、去除平地效应、相位解缠，最后把相位信息转换到高程。
2.2.1 SAR复数影像配准
在雷达干涉处理中，配准的目的是将2幅影像上的对应点处理成地面上同名地物点。先采用匹配方法找到2幅影像上的同名点，继而确定2点的偏移量，再对其中1幅影像进行重采样，使2幅影像上的同名点一一对应。
干涉雷达复图像的匹配，分粗匹配和精匹配两步进行[2]。确定图像搜索窗内按行列不同的整像元偏移量为粗匹配，计算复图像对之间的匹配质量评价指标,得到精度在像元级的匹配;精匹配通过对复图像作亚像元插值，再计算匹配指标。理论上，匹配度需要达到子像素级（1/10像素）[3]。匹配策略如下：
a）粗匹配采用基于卫星轨道参数的方法。由于卫星飞行的轨道参数是可确定的，且较准确和稳定，因此是最基本的粗匹配方法。其优点是能自动完成偏移量的估算，缺点是匹配精度较低。
式中，Pm（x,y）和Ps（x,y）对应主影像和辅影像的行列值。粗定位就是利用卫星轨道参数计算同名点对的偏移量offset（x,y）;
b）像元级匹配。干涉雷达测量中复图像的像元级配准通常采用基于窗口的配准技术,利用主图像和辅图像在空间域或频率域进行配准。笔者采用的匹配相似性评价参数是较基本相关系数法，它也是其他很多匹配算法的基础，且计算起较简单，不会出现太大误差;
c）亚像元级匹配。亚像元级匹配是在像元级匹配的基础上进行的。影像对（9192和9693）配准的统计分析见表1。
随意选取10个精配准点，检查其偏移量及精匹配的结果。匹配完成后，2幅影像上的同名点仅存在较小差异，用这些点作为控制点对，用二次多项式拟合确定2幅影像中任何1个匹配位置的相对偏移量;
2.2.2 干涉图的生成和处理
计算干涉系数由2幅影像上对应点共轭相乘得到，是作为评价干涉图质量的依据。干涉系数产生的相干图能直观地显示干涉图的质量，因此，相干图也可作为解缠质量的依据。通过处理得到，相干图上相干系数较大的区域是干涉效果较好的区域，干涉条纹较为明显，解缠也相对容易，反之，相干系数小或者干涉条纹不清楚的地方，一般不易解缠成功。
从中看出，震前震后的2幅干涉图中垂直基线都很大，在干涉图相位中的山区部分即右下角，基线去相关现象较严重，致使相干性很差，表现为黑色。
解缠前还要进行滤波，滤波的目的不仅为提高信噪比，且能为解缠计算提供较好的数据，在去除干涉
图噪声的基础上能保持边缘信息，防止产生误差。笔者采用Goldstein滤波算法，去除噪声后条纹清晰，很好地保留了地形信息。图4为主、辅影像干涉图去除平地效应的结果，为进行dstein滤波的结果。
将震前后去除平地效应的地形相减，得到形变的地形。这时，显示的还是周期性的条纹，表示地形的相位还是缠绕的，没有得到真实地形，需要通过解缠来得到地形的真实高程。
2.2.3 形变地形的相位解缠
相位解缠的结果直接影像地面沉降量的精度,是InSAR数据处理过程中的关键。相位解缠的基本思想从数学方面讲，假设给定1个二维相位矩阵Ψi,j，为了获得该矩阵的真实值，需要对每个点（i,j）加上或减去2π的整数倍，得到1个连续的函数Φi,j [6]，即
其中，ki,j为整数，-π≤Ψi,j<π，i ∈[0,M-1]， j ∈[0,N-1]。图6给出的是形变干涉的结果。
进行相位解缠计算时必须兼顾到一致性和精确性。一致性是指解缠后任意2点之间的相位差与2点之间路径无关，精确性则是指解缠后的相位要能忠实地恢复原始相位函数。将解缠相位梯度和缠绕相梯度之间的差异最小化，即最小成本网络流法，见图7。
3、实验结果
从干涉图可以看出,相隔175天和相隔35天的干涉图其相关性依然保持很好,这是因为巴姆地区气候干旱,植被稀少，保持了很好的相干性。
从中可以看到，条纹变化较大的为形变区。根据条纹的分布也可得出形变结果。ASAR雷达采用的是C波段微波成像，波长λ为5.666 cm，可利用式（3）计算出，
每个周期的条纹变化代表了2.8 cm的形变量。从图6中还可以看到，干涉条纹的颜色周期是不同的。根据GMT画出的图8、图9表示形变的三维状况。其中，北边区域为下沉区，其形变的最大值为-17.8 cm，南边区域为上升区，其形变最大值达到+29.9 cm。
4、结语
使用多级配准方法，利用三轨法差分干涉测量技术对ENVISAT ASAR雷达影像进行处理，可以使配准结果达到1/10像素;通过对干涉图进行去除平地效应、滤波等处理，分析获取影像期间地物的变化情况，针对被观测地区的具体条件，选择最优算法，可获取地震。
