基于合成孔径雷达差分干涉测量技术的矿区地面形变分析①

张文展

(新乡航空工业(集团)有限公司,河南 新乡453049)

0 引 言

20世纪60年代末,合成孔径雷达干涉测量技术(InSAR)作为合成孔径雷达(SAR)技术的一个新的应用领域发展起来,尢其是利用合成孔径雷达的相位信息提取地表的三维信息和高程变化信息的一项技术。由其拓展的合成孔径雷达差分干涉测量技术(D-InSAR)具有全天候、全天时、覆盖范围广等特点,被广泛应用于火山、地震形变、滑坡等监测[1-4],但是,将该技术应用于煤矿开采地面形变监测还处于初步研究阶段[5],对其进行研究具有重要的理论和实际意义。利用双轨法D-InSAR处理ENVISAT卫星获取4景覆盖济宁某矿区的先进的合成孔径雷达(ASAR)数据,分别获取了济宁某矿区的地面形变相位图、差分干涉图、沉降分布以及沉降面积等,并对该矿区地面形变进行了研究分析。

1 D-InSAR测量基本原理

图1是差分干涉测量的成像几何示意图[6]。A1和A2是卫星两次同一地区成像的位(即天线的位置),R1和R2表示目标点P到天线A1和A2距离,其中,R2=R1+ΔR.天线A1和A2对目标点P的测量相位差为[7]

式中,λ为雷达信号的波长。

由图1中的几何关系模型和余弦定理可得:

图1 差分干涉测量成像几何示意图

式中:B为天线A1和 A2之间的距离,称为基线距;θ为雷达入射角;α为水平线与基线的夹角。通常 ΔR≪R1,则式(2)可化简为

式中,B∥是基线距沿雷达视线方向的分量。把式(3)代入式(1)可得:

在地表发生形变后,进行观测获取的相位差不但与地形有关,而且与沿雷达视线方向的形变量ΔRd有关。此时对目标点P的测量相位差为

式中,B′∥表示地形未发生形变时,获取的干涉图的基线距在雷达视线方向上的分量。

由式(4)和式(5),可得雷达视线方向变量Δ Rd所引起的相位为

式中左边各量可由干涉纹图的相位和轨道参数计算得到,进而可确定图像上每点的视线向形变量ΔRd.

2 数据源以及数据处理的流程

2.1 数据源

收集了4景覆盖济宁某矿区的ENVISAT ASAR数据,这4景ASAR数据组合的干涉像对的具体参数如表1所示。本文所研究煤矿,矿井范围南北长4～4.5 km,东西宽4～6 km。矿区内地形平坦,地势东北略高,西南稍低,地面标高为+37.04～+41.28 m,平均为+38 m,自然地形坡度为万分之七。气候温和,四季分明,气候宜人;交通方便。ENVISAT卫星[8]是由欧洲空间局在2002年3月1日通过阿里亚纳5号火箭发射的一颗先进的太阳同步轨道地球环境监测卫星,该卫星上配有的ASAR,是目前先进的合成孔径雷达,具有多极化、多入射角、大宽度等特性。

表1 ASAR像对的基线列表

主图像选取考虑到时间基线和空间基线对图像相干性的影响以及成像期间植被的影响以及形成结果的可靠性等,利用D-InSAR技术中的双轨差分干涉法对表1中的影像数据对进行处理。采用双轨差分的优点是外部数字高程模型(DEM)和满足干涉的影像对容易获取,缺点是如果DEM精度不高,会对数据处理结果造成严重影响。本文选用分辨率为90 m的雷达地形测绘(SRTM)[9-10]数据作为外部DEM.

2.2 数据处理流程

采用双轨差分干涉测量法处理数据,该方法利用已知DEM反演的干涉相位进行二次差分处理,即从干涉相位中减去地形相位,得到只表征地表形变信息的相位,从而可计算地表形变量。其具体数据处理流程如图2所示[11]。

图2 双轨法差分干涉测量处理流程

3 D-InSAR结果及其矿区地面形变分析

3.1 D-InSAR处理结果

1)差分干涉图

干涉图是像对的对应点复共轭相乘得到的,其相位等于原始两SAR图像的相位差,反映了不同观测位置和天线到对应的分辨单元中心的路程差的函数。在干涉图中,2π的相位变化(图中表现为一个红蓝绿颜色周期)称为一个干涉条纹,由于复数对其相位的周期性,干涉得出的不是直接两次测量的相位差原值,而是其被周期折叠后的主值,即真实相位 2π的模。研究数据通过双轨法 D-In-SAR处理后获得增强差分干涉图和相干图如图3所示。

从图3中可以看出,2004年 12月26日和2005年1月30日两幅影像之间的相干性很差,只有极小的区域相干性可以,在差分干涉图中,除去地形信息,也难以明确的看出干涉条纹。但2005年1月30日和2005年3月6日,两幅影像之间的相干性较好,在差分干涉图中可以可出明显的沉降变化区域。

2)地面沉降图

对差分干涉图进行相位解缠,并进行相高转换后,即可得到沿雷达视线方向的形变图,即为沉降图。

在图4的沉降图可以看到所研究区域的某些区域形成了明显的“漏斗”沉降;某些区域发生了一些微小的沉降可能是由于开采地下水或数据处理中不可避免的误差造成;有些区域发生了上升现象,这可能与开采矿物有关,也可能与数据处理中有关参数设置不当有关。

3.2 矿区地面形变分析

根据矿区范围,从整体上分析该研究区域不同时间段的沉降量、沉降分布及其沉降面积统计。

1)2004年12月26日-2005年1月30日期间沉降分布及其沉降面积估算

2)2005年1月30日-2005年3月6日期间沉降分布及其沉降面积估算

从图5和图7中,可以看出整个矿区有不同沉降值的地面面积沉降,在这些地面沉降中,一些沉降程度不太明显,这些沉降可能是由于地下水开采过度所引起的,也可能是数据处理中的误差所造成的。图5和图7的左上方,都存在着由于煤矿开采而造成的明显沉降,分布区域也大体相同。从图6和图8中,可以看出地面形变分布以及沉降程度有所不同,不同沉降值所对应的沉降面积有所不同。

图4 双轨法D-InSAR沉降图

图8 2005年1月30日-2005年3月6日沉降面积曲线图

3.3 精度分析

为了验证D-InSAR监测矿区地面沉降的精度,收集了水准监测手段获取2004.11.20-2005.12.31时间段内所研究矿区地面沉降监测值,与DInSAR监测结果相比较。表2和图9给出了DInSAR监测结果与水准结果的比较情况。

图9 水准结果与D-InSAR监测结果比较

从表2和图9中,可以看出水准监测的沉降量结果和D-InSAR监测的沉降结果存在着一定差异,两者之间的相符程度不太高。这可能是由于造成这种现象的原因有:①成像分辨率和波长的制约;②空间上不一致;③时间上不一致;④ASAR数据的质量和数量:ASAR成像波段容易受到植被影响,所研究矿区植被分布比较多;缺少2004.11时间段的ASAR数据,造成了时间“漏洞”;⑤双轨法D-InSAR数据处理中,引入外部DEM精度不高。

4 结 论

采用双轨法D-InSAR处理覆盖济北某矿区的ENVISATASAR数据,提取了增强差分干涉图、相干图、沉降图,并从矿区地面沉降分析和精度分析等方面对D-InSAR监测矿区地面沉降结果进行了分析讨论,可以得出:

1)D-InSAR技术可以从整体上得到矿区沉降的整体分布、不同时间段的沉降量以及估算不同沉降值的沉降面积。

2)由于受到ASAR分辨率和成像波长等因素制约,造成D-InSAR技术监测矿区地面沉降结果与水准监测沉降结果之间存在一定差异。

[1] 路 旭.用InSA R作地面沉降监测的实验研究[J].大地测量与地球动力学,2002,22(4):66-70.

[2] Chris R.Tectonics of the Italian earthquake[EB/OL].[2009-04-06].http://scienceblogs.com./highlyalloch tho nous/2009/04/tectonics of the italicmearth.php.

[3] 龙四春,李 陶.D-InSAR中参考DEM误差与轨道误差对相位贡献的灵敏度研究[J].遥感信息,2009(1):3-6.

[4] Didier M.The displacement field of the Landers earthquake mapped by radar interferometry[J].Nature 1993(364):138-142.

[5] Howard.On the derivation of coseismic displacement fields usingdifferential radarinterferometry:the Lander earthquake[J].Journal of Geophysical Research Solid Earth,1994,99(B10):617-634.

[6] 邹积亭,刘远明.D-InSAR在地面形变监测中的应用研究[J].北京测绘,2009(2):19-20.

[7] 王 超,张 红,刘 智.星载合成孔径雷达干涉测量[M].北京:科学出版社,2002:59-60.

[8] Ramon F H.Radar interferometry data interpretation and error analysis[M].England:Kluwer Academic Publishers,2002:315-340.

[9] 郭华东,王长林.全天候全天时三维航天遥感技术-介绍航天飞机雷达地形测图计划[J].遥感信息,2000(1):47-48.

[10] 邵裕森,戴先中.过程控制工程[M].北京:机械工业出版社,2006:139-175.

[11] 周金国.InSAR数据处理方法与应用研究[D].重庆:重庆大学,2008:47-50.