为了利用星载合成孔径雷达差分干涉测量技术获取巴姆地震引起的形变场,欧空局(ESA)提供了7幅ENVISAT ASAR数据,7幅ASAR数据均为SLC产品,Image模式,Swath 2,极化方式为VV极化。其中4幅为降轨数据,Track号为120,巴姆古城位于图像中部; 3幅为升轨数据,Track号为385,巴姆古城位于图像的近距端。为了更全面地监测巴姆古城的地震情况,我们选取Track120的四景数据进行差分干涉测量处理(图1)。

图1 Track120四景数据所组成的干涉对

(垂直基线距由Descw软件计算出来,与实际的垂直基线距还存在偏差)

Track120共有震前数据两景,震后数据两景,总共可以组成6个干涉对。为了求得巴姆地震的形变场,根据雷达数据和巴姆地区的情况,差分干涉测量采用三轨法。根据三轨差分对空间基线和时间基线的要求,选取轨道号为9192的数据为主影像,轨道号为10194的数据为从影像,轨道号为6687的数据为辅助影像。

表1 Track120差分干涉基本参数信息

第一步是对9192-6687组成的topo-pair干涉对进行处理,首先需要对两景影像进行精确配准。配准过程首先采用基于荷兰Delft大学对地空间研究所(DEOS)的精密轨道数据的粗略配准,可以得出景中心的初始偏移量为(18,212),再采用基于相关系数的精确配准方法(WANG等,2005),通过最小二乘法进行多项式拟合可以得到几何偏移的拟合参数(表2)。

表2 Topo pair干涉对配准拟合参数

配准完成后就可以生成干涉图,多视处理为2:10(即距离向采样因子为2,方位向为10),生成的干涉图大小为2 583列×2 046行。

再采用基于轨道参数的方法去除平地效应。为了消除相位噪声,还需要对干涉条纹图进行滤波处理,滤波方法采用Goldstein滤波法(Goldstein,1998)。生成的干涉条纹图如图2a所示。同时进行精密基线估计,采用的方法是基于Delft研究所的精密轨道数据的基线估计方法,之后对滤波后的干涉条纹图进行相位解缠,相位解缠方法采用Minimum Cost Flow的方法(Eineder,1998)。

第二步采用上述相同的处理步骤对9192-10194组成的defo-pair干涉对进行处理,精确配准的拟合参数见表3。

表3 Defo-pair干涉对配准拟合参数

生成的干涉条纹图如图2b所示,然后用上面相同的方法进行相位解缠。生成的干涉图大小为 2 583列×2 064行。

第三步就是进行差分干涉处理,由于两个干涉对共用同一幅影像,所以干涉对之间不需要再进行配准,通过对基线进行比例标定,然后从defo-pair的相位中减去topo-pair的地形相位,再用公式(2)进行一定的转换,最终就可以得到雷达视线方向上的形变图(图2c)。

图2 Track120差分干涉结果(a)震前干涉条纹图(9192-6687);(b)同震干涉条纹图(9192-10194);(c)地表形变图

由于地形干涉对(9192-6687)的垂直基线较长(-470～-480 m),在去除了平地效应后,在干涉条纹图中显示出了与当地地形起伏相似的干涉条纹。

而形变干涉对(9192-10194)的垂直基线很小(-2～+2 m),在去除了平地效应后,干涉相位主要是由地形高程和地表形变引起的。但因形变干涉对的垂直基线很小,地形高程引起的相位由于与垂直基线成正比(Hanssen,2001),所以在干涉图中也就变得很小,从而使地表形变引起的相位在干涉图中占主导地位,故在图2b中显示出了与地表形变图相似的干涉条纹。

图2c就是所得到的雷达视线方向上的形变图,也就是我们求得的最终的地震同震形变场。在形变图中,一个颜色周期表示2π的相位变化,代表实际的地表形变为2.8 cm,即一个色周代表2.8 cm的形变量。从图中可以清晰地看到,由于地震的影响,在地表形成了一个蝴蝶形状的地表形变场,其中形变主要发生在东边的南北两个“花瓣”上。由颜色周期变化的不同,我们发现一个“花瓣”是上升,而另外一个“花瓣”却是下沉的。进一步分析发现,南边的“花瓣”是上升的,在雷达视线方向上形变量的最大值为+29.9 cm,北边的“花瓣”是下沉的,在雷达视线方向上的形变量的最大值为-17.8 cm。

为了更详细地分析地表形变的情况,我们取形变量最大处所在的行进行剖面图显示。图3a是地表下沉最大值所在的第897行的剖面图,图3b是地表上升最大值所在的第1 213行的剖面图。其中水平方向表示列像元值,竖直方向表示D-InSAR技术监测到的地表形变量,单位为米。图3c是地表下沉的最大值所在的第1 029列的剖面图,图3d是地表上升最大值所在的第1 113列的剖面图。其中水平方向表示行像元值,竖直方向表示D-InSAR技术监测到的地表形变量,单位为米。通过剖面图,我们可以清楚地看到地表形变的情况,包括形变的最大值、形变发生的距离等。

图3 形变最大处剖面图(a)897行的剖面图;(b)1 213行的剖面图;(c)1 029列的剖面图;(d)1 113列的剖面图

为了更进一步地分析地震对地表形变造成的影响,我们选取了地表形变最大值所在位置的部分区域进行进一步的分析,所选的区域大小为1 000行*1 000列,完全覆盖蝴蝶形状的地表形变场。

为了更直观地显示巴姆古城的地表形变量,我们将D-InSAR技术监测到的雷达视线方向上形变量相同的点相连,生成了等值线图。生成的等值线主要是+1、+5、+10、+15、+20、+25、-2、-5、-10、-15 cm。

图4是根据D-InSAR差分干涉测量得到的地表形变的等值线图。图中等值线上的数字表示地表形变量,正值表示上升,负值表示下降,单位为cm。图中的蓝色圆点为地表形变量的最大值所在的位置。

图4 巴姆地震地表形变的等值线图

表4 地面沉降面积统计

通过图4及表4,我们可以很清楚地看到地震在不同区域造成的地面沉降的大小和范围。

在InSAR干涉测量中,相干图表示的是两幅影像的相关性,也就是指的相关系数,相干图是最常用的相位质量图,也是最直观的质量评价图(王超等,2002)。相干图中相关系数的高低表明图像在不同区域的相干性。相干图在InSAR中起着重要的作用,一方面,它可以作为相位解缠时的重要依据,通过对相干图进行分析,可以确定在哪些区域相位是一致的,在哪些区域相位是不一致的。根据相干图来确定相位解缠的策略,可以指导解缠路径的设置,甚至可以作为解缠过程中的权值。另一方面,相干图的变化也表征了在图像获取期间地物的变化情况,所以相干图也可以用于地物的分类等研究。我们可以根据相干图的变化来识别一些特殊的现象。

下面,我们重点分析巴姆地震的相干图(图5)。

由于在震前相干图和震后相干图中垂直基线都很大,都在400 m以上,在震前、震后两幅相干图中的山区,基线去相关现象特别严重,致使在山区很多地方的相干性很差,表现为黑色(图5a,

图5 相干图的应用(a)震前相干图(9192-6687);(b)震后相干图(9693-10194); (c)同震相干图(9192-10194);(d)同震相干图(局部放大)

b的左下角区域)。而在同震相干图中,其垂直基线很小(-2～+2 m),基线去相关现象并不是很突出,也仅仅只是在某些山区相干性不好,如图5c的左下角部分区域。

由于巴姆地区气候干旱,植被稀少,巴姆古城主要是由土所筑成的砖墙建造而成的,正常情况下,在InSAR影像中会保持很高的相干性。而巴姆地震造成了巴姆古城区的房屋大面积倒塌,我们很想知道巴姆地震造成的破坏情况,然而由于SAR幅度影像中斑点噪声的大量存在,房屋的大面积倒塌无法直接反映在雷达影像上,无法直接通过雷达幅度影像对其进行分析。然而可以考察InSAR相干图,在InSAR相干图上反映的是地物的相干性,由于房屋的大面积倒塌,会使得InSAR的相干性大大降低,甚至在某些地区的相干性会很低。我们通过对震前、震后、同震相干图的对比,就可以检测出地震造成房屋大面积破坏的区域的大致位置和分布情况。

通过对比震前、震后、同震相干图,在去除了基线去相关的影响后,我们可以看到在同震相干图中,巴姆古城附近明显有一块连续区域的相干性很差,如图5c中的黑色矩形框所示,这就是地震造成破坏最严重区域的位置以及大致的范围。图5d是图5c黑色方框内的局部放大图。

通过对图5d中相干系数小于0.5的像元进行统计,在相干图中有24 346个点小于0.5,再根据像元大小,我们可以得到巴姆地震破坏最严重区域的面积大约是39 km²。