陶秋香,刘国林,刘伟科
山东科技大学 测绘科学与工程学院,青岛 266510
矿区地面沉降严重危害着地面建筑设施和自然环境,影响着人类的生存环境、生命财产安全和当地经济发展.对矿区地面沉降进行有效监测,深入研究地面沉降的形成机理与变化规律,对合理开采地下矿产资源,控制矿区地面沉降相当重要[1-3].传统大地水准测量、静态GPS测量或动态GPS监测矿区地面沉降时得到的都是一个点的沉降信息,随着测绘手段的发展和地面沉降监测应用范围的进一步扩展,其缺点也越来越突出.星载合成孔径雷达干涉测量 (Interferometric Synthetic Aperture Radar,InSAR)技术是近年来快速发展起来的空间对地观测新技术,其差分干涉测量模式(Differential InSAR,D-InSAR)在矿区地面沉降监测的应用国内外均取得若干成功案例,监测精度已达到毫米级.如澳大利亚新南威尔士大学Ge等在2001年就已经将GPS和InSAR技术结合起来进行矿区地面沉降监测,精确获取了研究区煤矿开采造成的地面沉降信息[4-8].国内一些学校和科研院所如中国矿业大学(北京)和山东科技大学也相继提出应用DInSAR技术进行煤矿区地表演变规律及开采沉陷监测实验研究,并开展了理论探索、数据处理等研究[9-10].
国内外现有D-InSAR监测矿区地面沉降的研究主要集中在利用单一SAR数据,对其模型和算法进行改进,提高矿区地面沉降监测精度.而对所用SAR数据因成像参数的不同必然会引发不同的监测结果和精度的研究却很少[11-14].但在实际应用中,SAR成像的波长、入射角、地面分辨率等参数严重影响着SAR差分干涉测量监测地面沉降的能力和精度,合适SAR数据源的选取是D-InSAR处理能否成功的最为关键的因素之一.本文以目前矿区地面沉降监测中最常用的ALOS PALSAR和ENVISAT ASAR数据为例,通过理论推导和矿区实际沉降差分相位模拟,对比分析了L和C波段雷达干涉数据监测到的最大沉降梯度、最大沉降量、保相能力、对微小沉降的敏感程度、监测的沉降范围大小等;以济宁某矿区为实验区,精化双轨差分干涉的数据处理流程和方法,分别选取成像时间最相近的ALOS PALSAR影像对和ENVISAT ASAR影像对进行对比实验,研究L和C波段雷达干涉数据在矿区地面沉降监测中的实际应用情况,得出一些有益的结论.研究结果对如何为指定矿区进行地面沉降监测选取合适的SAR数据、提高监测精度具有重要的参考价值和实用意义.
以ALOS PALSAR和ENVISAT ASAR影像数据为例,对L和C波段的雷达干涉数据进行简单介绍.ALOS(Advanced Land Observation Satellite,ALOS)是由日本于2006年1月发射的一颗先进的太阳同步极轨地球环境监测卫星.ALOS上搭载的L波段合成孔径雷达PALSAR(Phased Array type L-band Synthetic Aperture Radar,PALSAR)具有高分辨率、扫描式合成孔径雷达和3种极化观测模式,能获取比普通SAR更宽的地面幅宽,主要应用于自然灾害监测、资源环境调查、雷达遥感教学与科研等领域.ENVISAT卫星是由欧空局于2002年3月发射的一颗太阳同步极地轨道卫星,可以提供关于大气、海洋、陆地和冰的测量信息,支持地球科学研究,可以对环境、气候变化进行监测.ENVISAT上搭载的ASAR工作于C波段,具有7种入射角度和5种工作模式.其先进性表现在雷达数据获取的覆盖范围、入射角范围、极化方式和操作模式的灵活性等方面[15-16].表1列出了本文所用L波段ALOS PALSAR和C波段ENVISAT ASAR数据的部分重要参数.
表1 ALOS PALSAR和ENVISAT ASAR数据部分重要参数Table 1 Parameters of ALOS PALSAR and ENVISAT ASAR images
1998年 Massonnet和Feigl[17-18]就提出了差分InSAR监测到的最大形变梯度理论,即
其中,dx为InSAR可监测的最大形变梯度,λ为雷达入射波波长,ps为像元大小(SAR传感器的地面分辨率).若差分干涉图相邻像元的形变梯度小于dx,则可认为其形变是连续形变,反之,则认为其形变是一种空间不连续形变.
干涉图中两个相邻像元之间的相位差如果小于1/4个波长,则认为其相位是连续的[19-20].假设研究区地面沉降相位是连续的,则InSAR监测到的整个下沉盆地沿视线方向的最大形变量可以表示为[20]:
其中,ΔRmaxLOS表示InSAR监测的整个下沉盆地中视线方向(Line Of Sight,LOS)上的最大形变量,r为下沉盆地的主要影响半径.从而可以得出InSAR监测的垂直向上的最大沉降量为[20]:
其中,ΔRmax表示垂直方向上的最大沉降量,θ为雷达入射角.
利用公式(1)可以计算出L波段的 ALOS PALSAR和C波段ENVISAT ASAR数据监测地面沉降的最大沉降梯度分别为7.9×10-3和9.3×10-4.这样计算出来的最大沉降梯度是一个无量纲的值.而当假定下沉盆地的主要影响半径为150m时,利用公式(2)和(3)可以计算出二者能监测到的垂直向上的最大沉降量分别为756mm和76mm.由此可以看出,由于C波段数据波长较短、地面分辨率较低,能够监测到的最大沉降值和最大沉降梯度都远远小于L波段雷达干涉数据.因此,在沉降量和沉降梯度较大的区域(如下沉盆地中心),由于波长和地面分辨率的原因,C波段雷达干涉数据就会出现相位不连续或混叠的现象,难以探测出准确的沉降值.而L波段SAR数据由于波长较长,地面分辨率高,则能够较好地避免这种现象的出现,监测到较大的沉降梯度和沉降量.
公式(2)和(3)是在假定研究区的沉降连续的情况下得出的,若沉降不连续,就无法利用这两个公式计算最大沉降量了.利用公式(1)—(3)计算雷达数据监测到的最大沉降梯度和沉降量均未考虑各种去相关因素如时空去相关、热噪声去相关、大气延迟等的影响.因此,在实际应用中,L波段的 ALOS PALSAR和C波段ENVISAT ASAR数据可监测的沉降梯度和沉降量要小于其理论值[21-22].
不同类型雷达干涉数据对地面沉降监测的能力因其星载SAR系统所用的波长、地面分辨率、入射角以及重复周期等参数而不同.为进一步形象说明L和C波段雷达干涉数据矿区地面沉降的监测能力,本文利用矿区实际沉降,考虑成像的波长、入射角和地面分辨率等,模拟了矿区地面沉降的差分干涉相位.图1为矿区实际的地面沉降,其垂直向上最大沉降值为69cm,图2和图3分别是由图1所示的矿区地面沉降模拟出的L和C波段差分干涉图的绝对相位和缠绕相位.
模拟相位均未考虑去相干因素、大气延迟以及轨道误差对形变信息的影响.从图1—3可以看出:
(1)比较图2a和图3a,由于L和C波段SAR干涉数据的波长、入射角和地面分辨率不同,对于同一矿区地面沉降而言,其在差分干涉图距离向和方位向上所分布的像元数、转换到LOS向上的沉降量以及由此得到的沉降相位都不同.
(2)图2b中L波段SAR数据差分干涉图中条纹较少,在沉降梯度较大的区域(图2b的C处)的相位也保持着良好的连续性,未出现相位缺失和混叠现象,这表明L波段SAR数据更容易监测到沉降梯度和沉降量较大的矿区地面沉降,并且解缠工作容易进行.而图3b中的差分干涉条纹则非常密集,而且对应的图3b的C处的差分干涉条纹还出现了大量混叠现象,这是因为此处的高条纹率超出了相位饱和度的极大值,这种高密集以及重叠的条纹反映在差分干涉图中使相位解缠工作很难正确进行,因而也很难监测到正确的沉降结果.
图1 矿区地面沉降Fig.1 Mining land subsidence
(3)由于图2b中A和B处仅出现了轻微的地面起伏,未形成一条完整的差分干涉条纹,而图3b中相应的A和B处则形成了完整的连续差分干涉条纹,这表明由于波长较短,C波段SAR数据对微小地面沉降的反映更加敏感.
本文所选研究区位于济宁煤田西北部,矿井范围南北长4~4.5km,东西宽4~6km.矿区内地形平坦,地面标高为+37.04~+41.28m,平均为+38m.开采前期沉降幅度较小,后期较大;前1~2个月工作面沉降一般为25~66cm,一年后沉降严重的甚至可达80~120cm.结合该矿区地面沉降情况,分别选取了成像时间最相近的一对ALOS PALSAR影像和ENVISAT ASAR影像进行对比实验,影像对基本参数如表2所示.
表2 ALOS PALSAR和ENVISAT ASAR影像对的基本参数Table 2 Basic parameters of ALOS PALSAR and ENVISAT ASAR images pair
按照构成SAR影像对的方式不同,差分干涉测量分为双轨、三轨和四轨法.鉴于卫星传感器的重返周期及研究区的数据收集情况,本文采用了双轨DInSAR技术获取矿区地面沉降信息,并对影像对的配准、差分干涉图的噪声去除、相位解缠、研究区DEM的选取等关键步骤和方法进行了精化,获取矿区更为精确的地面沉降信息.如影像对配准过程中,不同的配准方法其配准精度不同,同一配准方法中窗口和相应匹配指标阈值设置不同,配准精度也不相同.为达到较好的配准效果,在四景SAR影像配准过程中,我们利用不同配准方法、在同一方法中不断改变窗口和匹配指标阈值的大小来进行配准,对配准结果进行归纳和总结,使配准尽可能精确,减少配准误差对监测结果的影响.本文所用双轨DInSAR差分干涉数据处理流程如图4所示.其中外部参考DEM选用的是研究区25m高分辨率的DEM.
5.2.1 实验结果
在进行双轨差分干涉处理过程中,首先将辅影像和外部参考DEM分别与主影像配准,然后作差分干涉处理,去除平地效应并做自适应滤波处理后得到差分干涉图和相干图.图5分别为ALOS PALSAR和ENVISAT ASAR增强后的差分干涉图和幅度图叠加的示意图和相应相干图.图6—7分别为图5a和图5c所标注的A和B附近的局部放大图及其相应的差分干涉相位.原始SAR影像的视数是1,影像如同被拉伸,与实际地物的差别较大,无法对影像进行正确的判断和解译.因此,在成像处理过程中需要进行多视处理,一般情况下,ALOS PALSAR和ENVISAT ASAR雷达干涉数据的多视系数(多视尺度)分别为1:2和1:5.多视系数的确定是由数据获取时被拉伸程度的大小决定的[21-22].多视处理会降低影像的分辨率,由公式(1)—(3)可以看出,这就不可避免地会降低InSAR监测到的最大沉降梯度和沉降量,而且不同的多视尺度(多视系数)造成SAR影像分辨率下降的程度不同,对雷达干涉数据监测到的最大沉降梯度和沉降量的能力影响也就不相同.图5—7分别是对ALOS PALSAR和ENVISAT ASAR原始数据进行了1∶2和1∶5多视处理后的SAR影像差分干涉所得到的干涉图.
图4 双轨差分干涉测量监测矿区地面沉降的数据处理流程Fig.4 Data processing flow of two-pass D-InSAR used to monitor mining subsidence in this study
对增强后的差分干涉相位利用最小费用流相位解缠算法进行相位解缠并进行相高转换和地理编码后即可得到研究区的矿区沉降图.图8和图9为ALOS PALSAR和ENVISAT ASAR经双轨差分干涉处理后得到的矿区地面沉降图及相应的局部放大图.
表3 ALOS PALSAR和ENVISAT ASAR影像监测矿区地面沉降的能力比较Table 3 Comparison of ALOS PALSAR and ENVISAT ASAR images to monitor land subsidence in mining areas
5.2.2 结果分析
(1)由于ALOS PALSAR的地面分辨率较高,图5a中研究区内地物类型、形状、纹理和分布都要比图5c清晰得多.而在相干图5b和图5d上也充分反映了这一特点.
(2)相干图是最常用的相位质量图,相干图上相干系数的大小反映了两幅雷达影像相干性的高低.由图5-9能获取研究区各像元的相干性、差分干涉条纹(形变条纹)的连续性和沉降值,表3总结了图6-9中C、D、E、F处ALOS PALSAR和ENVISAT ASAR监测矿区沉降的相干性、条纹的连续性、能够监测到的垂直向上沉降最大值.
从表3可以看出:(1)在相干性都较高的沉降区C、D两处,ALOS PALSAR和 ENVISAT ASAR差分干涉图(图6)上均出现了连续的干涉条纹,监测得到的沉降值也是可信的,但是由于波长和地面分辨率的限制,ENVISAT ASAR监测到的沉降值明显偏小;(2)在沉降区C处,ALOS PALSAR差分干涉图(图6b)未形成一条完整的干涉条纹,而ENVISAT ASAR差分干涉图(图6d)上则有完整的干涉条纹出现,其监测到的沉降值也偏小,这说明ENVISAT ASAR对微小的地面沉降具有较高的敏感度,比较适合监测幅度较小的地面沉降;(3)在沉降区E处,由于植被覆盖或沉降梯度和幅度太大,超出了ALOS PALSAR和ENVISAT ASAR的监测能力,影像对严重失相干,差分干涉图(图6b和图6d)均未出现任何干涉条纹,二者均无法监测到该处的沉降情况;(4)在沉降区F处,ALOS PALSAR和ENVISAT ASAR影像对的相干性均不太高,但ALOS PALSAR由于其波长较长的特点,其相干性要好一些.其差分干涉条纹(图6b)中仅有少量的不连续点出现,监测到的最大沉降值为766mm.而ENVISAR ASAR的差分干涉条纹(图6d)中则出现大量不连续点,条纹时断时续,监测的最大沉降值仅为83mm.由于此处失相干较严重,相位解缠误差较大,监测到的沉降结果误差较大,可信度较差;(5)综合比较这四处的沉降监测情况,可以看出,雷达影像能够监测到的最大沉降值不仅与波长、沉降影响半径、地面分辨率和入射角有关(公式(2)、(3)),还严重受两景SAR影像相干性的影响.可监测的最大沉降梯度同样也受相干性的影响,影像中相干性高的区域,InSAR监测最大沉降梯度和沉降量的能力也就越强,反之亦然.
论文通过理论推导、相位模拟和矿区地面沉降监测实验,对比分析了L和C波段雷达干涉数据矿区地面沉降的监测能力.得出:
(1)在不考虑雷达影像去相关因素的情况下,理论推导得出L和C波段雷达干涉数据监测地面沉降的最大沉降梯度分别为7.9×10-3和9.3×10-4.而当假定下沉盆地的主要影响半径为150m时,监测到的垂直向上的最大沉降量分别为756mm和76mm.但在实际应用中,由于受影像的相干性和多视尺度的影响,L和C波段雷达干涉数据可监测的地面沉降梯度和沉降量小于理论值.
(2)相位模拟结果显示在沉降梯度和沉降量都较大的区域,L波段雷达干涉数据的保相能力较强,相位连续性和清晰度都较高,并且解缠工作容易进行,沉降监测结果可靠;而C波段雷达干涉数据的保相能力相对较差,差分干涉相位会出现大量混叠现象,沉降监测结果的可信度较差.但L波段雷达干涉数据对微小地面沉降的敏感程度则不如C波段的雷达干涉数据.
(3)真实ALOS PALSAR和ENVISAR ASAR雷达干涉数据的矿区地面沉降监测实验则表明,L和C波段雷达干涉数据均能精确确定矿区地面沉降的位置,对矿区地面沉降进行定性分析,但L波段SAR数据由于其地面分辨率高,波长较长,能够更好地降低失相干和相位不连续性的影响,更适合监测沉降范围较小而沉降梯度和沉降量都较大的矿区地面沉降.这也进一步验证了“L波段可以监测范围小形变梯度大的形变,而C波段适合监测范围大形变梯度小的形变”这一事实.
(4)本文的研究结果也表明,由于受植被覆盖等其他失相关因素的影响及现有数据处理方法的限制,L和C波段雷达干涉数据监测地面沉降的实际应用能力被弱化,但二者仍表现出较强的矿区地面沉降监测能力,随着各项技术的发展和完善,相信将来雷达干涉数据矿区地面沉降的监测能力会大大提高.
(References)
[1]董玉森,Ge L L,Chang H C等.基于差分雷达干涉测量的矿区地面沉降监测研究.武汉大学学报 (信息科学版),2007,32(10):888-891.Dong Y S,Ge L L,Chang H C,et al.Mine subsidence monitoring by differential InSAR.Geomatics and Information Science of Wuhan University (in Chinese),2007,32(10):888-891.
[2]Herrera G,Tómas R,Lopez-Sanchez J M,et al.Advanced DInSAR analysis on mining areas:La Union case study(Murcia,SE Spain).Engineering Geology,2007,(90):148-159.
[3]吴立新,高均海,葛大庆等.工矿区地表沉陷D-InSAR监测试验研究.东北大学学报 (自然科学版),2005,26(8):778-782.Wu L X,Gao J H,Ge D Q,et al.Experimental study on surface subsidence monitoring with D-InSAR in mining area.Journal of Northeastern University (Natural Science)(in Chinese),2005,26(8):778-782.
[4]Ge L L,Rizos C,Han S,et al.Mining subsidence monitoring using the combined InSAR and GPS approach.The 10thFIG Int Symp.on Deformation Measurements orange,California,2001:1-10.
[5]Demirel N,Emil M K,Duzgun H S.Surface coal mine area monitoring using multi-temporal high-resolution satellite imagery.International Journal of Coal Geology,2011,86(1):3-11.
[6]Ng A H-M,Ge L L,Yan Y G,et al.Mapping accumulated mine subsidence using small stack of SAR differential interferograms in the southern coalfield of New South Wales,Australia.Engineering Geology,2010,115(1-2):1-15.
[7]Yang C S,Zhang Q,Zhao C Y,et al.Monitoring mine collapse by D-InSAR.Mining Science and Technology(China),2010,20(5):696-700.
[8]王艳,廖明生,李德仁等.利用长时间序列相干目标获取地面沉降场.地球物理学报,2007,50(2):598-604.Wang Y,Liao M S,Li D R,et al.Circulation characteristics of interannual and interdecadal anomalies of summer rainfall in north Xinjiang.Chinese J.Geophys.(in Chinese),2007,50(2):598-604.
[9]Hu Z L,Li H Q,Du P J.Case study on the extraction of land cover information from the SAR image of a coal mining area.Mining Science and Technology,2009,19(6):829-834.
[10]Wang Z Y,Zhang J Z.Use of D-InSAR technique for monitoring ground subsidence in the Yanzhou coal mining area(China).Applied Mechanics and Materials Journal,2010,34-35:756-760.
[11]Oha H J,Lee S.Integration of ground subsidence hazard maps of abandoned coal mines in Samcheok,Korea.International Journal of Coal Geology,2011,86(1):58-72.
[12]Herrera G, Tomás R, Vicente F.Mapping ground movements in open pit mining areas using differential SAR interferometry.International Journal of Rock Mechanics &Mining Sciences,2010,47(7):1114-1125.
[13]Sarychikhina O,Glowacka E, Mellors R,et al.Land subsidence in the Cerro Prieto Geothermal Field,Baja California,Mexico,from 1994to 2005:An integrated analysis of DInSAR,leveling and geological data.Journal of Volcanology and Geothermal Research,2011,204(1-4):76-90.
[14]Wang Z Y,Zhang J Z,Liu G L.Measuring land subsidence by PALSAR interferometry in Yanzhou coal mine area.Proc.of SPIE,2010,7820:1-8,doi:10.1117/12.867024.
[15]Lubis A M,Sato T,Tomiyama N,et al.Ground subsidence in Semarang-Indonesia investigated by ALOS-PALSAR satellite SAR interferometry.Journal of Asian Earth Sciences,2010,40(5):1079-1088.
[16]陶秋香,刘国林,孙翠羽.SAR影像中PS点的识别与选取.应用科学学报,2009,27(5):508-513.Tao Q X,Liu G L,Sun C Y.Identification and selection of persistent scatterer pixels from SAR images.Journal of Applied Sciences(in Chinese),2009,27(5):508-513.
[17]Massonnet D,Feigl K L.Radar interferometry and its application to changes in the Earth′s surface.Reviews of Geophysicsl,1998,36(4):441-500.
[18]赵超英.差分干涉雷达技术用于不连续性形变的监测研究[博士论文].西安:长安大学,2009.Zhao C Y.Research on the Monitoring of Discontinuous Deformation by Differential SAR Interferometry[Ph.D.thesis](in Chinese).Xi′an:Chang'an University,2009.
[19]Chen C W.Zebker H A.Phase unwrapping for Large SAR interferograms:statistical segmentation and generalized network models.IEEE Transactions on Geoscience and Remote Sensing,2002,40(8):1709-1719.
[20]阎跃观.DInSAR监测地表沉陷数据处理理论与应用技术研究[博士论文].北京:中国矿业大学,2010:23-24.Yan Y G.Data processing theory and applied technology of surface subsidence monitoring based on DInSAR[Ph.D.thesis](in Chinese).Beijing:China University of Mining &Technology,2010.
[21]Baran I,Stewart M,Claessens S.A new functional model for determining minimum and maximum detectable deformation gradient resolved by satellite radar interferometry.IEEE Transactions on Geoscence and Remote Sensing,2005,43(4):675-682.
[22]蒋弥,李志伟,丁晓利等.InSAR可检测的最大最小变形梯度的函数模型研究.地球物理学报,2009,52(7):1715-1724.Jiang M,Li Z W,Ding X L,et al.A study on the maximum and minimum detectable deformation gradients resolved by InSAR.Chinese J.Geophys.(in Chinese),2009,52(7):1715-1724.