祝传广,邓喀中,张继贤,张永红,范洪冬,张立亚
(1.中国矿业大学国土环境与灾害监测国家测绘局重点实验室,江苏徐州 221116;2.中国测绘科学研究院,北京 100036;3.煤炭资源清洁利用与环境保护湖南省重点实验室,湖南湘潭 411100)
基于多源SAR影像矿区三维形变场的监测
祝传广1,2,邓喀中1,张继贤2,张永红2,范洪冬1,张立亚3
(1.中国矿业大学国土环境与灾害监测国家测绘局重点实验室,江苏徐州 221116;2.中国测绘科学研究院,北京 100036;3.煤炭资源清洁利用与环境保护湖南省重点实验室,湖南湘潭 411100)
针对DInSAR(differential interferometric synthetic aperture radar)技术仅能获取雷达视线向(line of sight,LoS)形变的不足,研究了融合多卫星平台求解三维形变场的模型与算法。该算法基于多卫星轨道模式下具有不同成像几何的多源SAR影像联合求解矿区地表形变场。研究结果表明:采用该算法反演的下沉值与水准测量结果相互吻合,均方根误差为±4 mm,吻合程度优于单一影像源反演结果;垂直向位移场与等值线均表明下沉盆地向老采空区偏移,说明老采空区可能活化;东西向水平位移场与等值线符合开采沉陷地表移动规律,而且对于不同的成像模式,东西向水平移动的影响亦不同;由于卫星航向角的正弦值近乎为0,使得三维算法对南北向位移不敏感。
合成孔径雷达;多源影像;雷达视线向;水平移动;开采沉陷
Key words:synthetic aperture radar;multi-source imagery;line of sight;horizontal displacement;mining subsidence
开采地下煤炭资源会破坏岩层内部应力平衡力,使得地表出现下沉、倾斜等变形。为保护地面建筑物并指导地下开采工作,进行地表变形观测是必要的。目前,观测方法主要为水准测量与GPS测量,虽然精度较高(静态GPS测量可达mm级精度,水准测量可优于mm级),但是观测点的空间密度低、费用高。在Gabriel成功利用合成孔径雷达差分干涉(differential SAR interferometrical,DInSAR)技术获取地表形变[1],证明DInSAR技术进行变形监测的可行性之后,该技术凭借广空间覆盖、高空间采样率以及快速重复观测等优点,广泛应用于地震、火山、滑坡等自然灾害[2-6]的监测中。随后,为研究长时间跨度的形变信息,基于时间序列的DInSAR技术被陆续提出,如PSInSAR技术[7-8]、SBAS-InSAR技术[9-10]等,并且在监测地下水资源[11-12]开采等引起的大范围地表位移方面取得了成功。
然而,DInSAR技术仅能获取目标在雷达视线向(line of sight,LoS)的一维形变。在研究目标的垂直向形变时,通常忽略水平移动,然后根据雷达波入射角将LoS向位移投影到垂直向。这对于大范围沉降是可行的,如抽取地下水导致的地表沉降范围通常比较大,使得小范围内的水平移动可以忽略。但是在小范围内发生大的水平移动是矿区地表形变的特点,因此需要采用新的技术手段研究矿区的变形规律。
DInSAR技术获取的LoS向位移可以分解成垂直向、东西向与南北向3个位移分量,当利用至少3种具有不同成像几何的影像源获取目标在不同LoS向的位移后,就可以联合解算目标的三维形变场。基于这种思想,根据两种卫星传感器ASAR与PalSAR,选取了3种具有不同成像几何的SAR影像源,分别独立解算出LoS向位移,然后联合求解研究区域的三维形变场,并研究分析了东西向水平移动对单一影像源反演的LoS向位移的影响。
对于重复轨道干涉测量模式,如果地面目标k在两次成像期间发生了位移,则其干涉相位φ可记为
其中,ω表示缠绕算子;φlos为LoS向位移相位;φatm, φflat,φtopo与φn分别为大气变化、平地效应、地形与噪声相位。结合各相位[13-14]自身的特点,可以利用DInSAR技术[15]从式(1)中分离出φlos,也就获取了目标在LoS向的位移Dlos。在三维空间中,Dlos可以分解成垂直向分量Dv、东向分量De及北向分量Dn[16-17],如图1所示。
图1 LoS向位移三维分解模型Fig.1 Decomposition of the displacement along the LoS
图1中,航向角α指卫星飞行方向与北方向间的夹角,由北方向开始指向飞行方向。θ是研究区域中心处目标k处的雷达波入射角,下标sl表示目标k处的LoS在地面上的投影,Dn,sl表示北向位移分量Dn在sl上的投影,De,sl表示东向位移分量De在sl上的投影。Dsl则是Dn,sl与De,sl的矢量和,其在LoS向的投影为Dsl,los。Dv,los表示垂直向位移分量Dv在LoS向的投影,Dsl,los与Dv,los的矢量和即为Dlos。因此,根据图1所示的三维分解原理,目标k的LoS向位移Dlos可记为
式(2)中有3个未知量即Dv,De与Dn,并且这3个未知量是目标k在三维空间中的真实位移,与雷达波的入射角、波长等无关。因此,当获取同一目标k在n≥3个不同LoS向的位移时,式(2)可记为
2.1 区域概况
研究区域在3种SAR影像中的位置以及开采工作面、水准点的分布如图2所示,工作面采用走向长壁全部垮落法开采。其中,工作面wf1的开采时间为2009-06-10—2009-12-31,走向长度510 m,倾角5°,煤层厚2.6 m,平均采深1 163 m;工作面wf2的开采时间为2009-08-15—2010-03-20,走向长度660 m,倾角3°,煤层厚2.5 m,平均采深1 140 m,其余工作面的开采时间在2009年之前。研究区域内有水准点43个,平均间距35 m。
图2 研究区域位置以及工作面、水准点的分布Fig.2 Location of the study area and the layout of working face and leveling points
2.2 实验数据
基于2种卫星传感器ASAR与PalSAR,笔者收集了3种轨道模式下的SAR影像,其成像时间与部分几何参数见表1。由表1可以看出,ASAR数据组成的2个干涉对,其垂直基线与时间基线均较短,可以有效提高相干性以及测量精度;对于PalSAR数据组成的干涉对,其时间基线较短,虽然空间基线达580 m,但是相对于6 km的临界基线仍可认为比较短。
表1 3种SAR影像的成像参数Table 1 Imaging paremeters of three types of SAR data
首先,表1中的3个干涉对分别利用DInSAR技术提取研究区域在3个LoS向的位移。因为3个干涉对的时间跨度并不完全重叠,需要通过插值获得相同时间段(2010-01-13—2010-02-28)内的位移值。然后,根据式(4)计算研究区域内所有目标在上述时间段内的三维形变场,如图3所示。在三维空间中,以向上、向东、向北为正方向。
图3 垂直向、东西向与南北向位移场与等值线Fig.3 The displacement and contour in vertical,eastwest and north-south direction
图3(a)为垂直向位移场,图3(d)为与之对应的下沉等值线图,可以看出,在2010-01-13—2010-02-28,研究区域下沉量值在0~60 mm,下沉盆地中心位于工作面wf1的上方,但是下沉盆地并没有表现为类似工作面wf1的椭圆形状,而是向两侧偏移,分析其原因可能是受到工作面wf1的影响导致其附近的老采空区活化,从而使得下沉盆地范围向老采空区一侧偏移,其中向北方向的偏移还有可能受到工作面wf2的开采影响。
由东西向位移场图3(b)以及东西向位移等值线图3(e)可以看出:研究区域的水平移动在-15~20 mm。在下沉盆地中心,水平移动量值比较小。自盆地中心向西与向东,水平移动量值均表现为先增大后减小的趋势,并且在下沉盆地西侧表现为向东移动,在东侧则表现为向西移动,符合开采沉陷地表移动规律。
由图3(c)与3(f)代表的南北向位移场与等值线图可以看出,只在工作面wf1的西侧与南侧有5~10 mm的北向位移,其余大部分区域的南北向位移不明显且较为杂乱,究其原因,由于现在的SAR卫星均采用极轨飞行方式,飞行方向与北方向的夹角很小,式(2)中航向角α的正弦值近乎为0,因而使得本文算法对南北向位移不敏感,反演结果受噪声影响严重。
3.1 与实测下沉值的对比分析
利用研究区域内43个水准点的观测资料对本文采用的三维算法获取的垂直向位移进行验证分析,如图4(a)所示。可以看出,两者反映的下沉趋势相互一致,并且在量值上也互相吻合,最大相差10 mm,绝对值相差在5 mm以内的点占77%,基于43个水准点统计两者间的均方根误差RMSE为±4 mm。由此可以说明本文所采用的模型与算法是有效的,同时也表明本文利用DInSAR技术提取的3个LoS向位移结果也是可靠的。为进一步对比分析,根据传统算法忽略水平移动,采用式(5)计算垂直向位移。
式中,入射角θ与卫星轨道模式以及目标所在位置有关,但由于本文研究区域较小使得因位置差异引起的θ的变化可以忽略不计,因此对同一轨道模式下的SAR影像数据采用相同的θ值。
根据式(5)将3个LoS向位移投影到垂直向,结果如图4(b)所示。可以看出:传统计算结果与实际下沉趋势是一致的,但是根据PalSAR数据计算的下沉盆地向西发生了偏移,而由2个ASAR数据计算的下沉盆地形状以及下沉量值均相似,并且与实际下沉的吻合程度优于PalSAR数据结果。为定量比较2种算法,本文基于最大相差、绝对值相差在5 mm以内点的比例以及均方根误差3项标准对2种方法进行了统计分析,结果见表2。可以看出:①定量分析的结果亦表明三维算法优于传统计算方法;②由式(5)可知,随着雷达波入射角θ的增大,3组SAR数据反演结果与实际下沉的偏离应该越大。统计结果则证实了该结论的正确性。
图4 垂直向位移与水准资料对比Fig.4 Comparison between the vertical displacements and leveling data
表2 不同算法得到的垂直向位移间的统计Table 2 Statistics of vertical displacements calculated using different methods
3.2 不同LoS向位移间的对比分析
图5显示的是43个水准点在3个LoS向的位移以及利用三维算法得到的东西向水平移动。从图中可以看出:①利用2组ASAR数据提取的LoS向位移趋势互相一致,量值也近乎相同,这是由该2种SAR数据的成像几何参数相差较小所致。②ASAR与PalSAR两种数据结果相差则较大,这是因为两者的雷达波入射角与卫星航向角相差均比较大所致。并且结合实际下沉曲线可以发现,雷达视角越大,LoS向位移与实际值相差也越大。③由东西向水平移动曲线可以看出,在盆地中心,水平移动近乎为0,并且3个LoS向位移值均小于实测值。这点可以由式(2)证明,即在盆地中心,东西向与南北向位移均较小,使得LoS向位移仅是垂直向位移的一个分量。④在下沉盆地中心西侧地表向东移动,并表现为先增大后减小的趋势,在盆地中心东侧地表向西移动,也表现为先增大后减小的趋势,符合开采沉陷地表移动规律。并且由于本实验选用的传感器均采用右侧视成像,使得在下沉盆地中心西侧,东向的水平移动使得目标在ASAR影像上表现为靠近卫星,在PalSAR影像上则表现为远离卫星;而在下沉盆地东侧,西向的水平移动使得目标在ASAR影像上表现为远离卫星,在PALSAR影像上则表现为靠近卫星。同时由图5可以看出,在下沉盆地西侧利用ASAR数据获取的LoS向位移量值小于实测下沉值,利用PALSAR数据获取的LoS向位移量值则大于实测下沉值,而在下沉盆地中心东侧情况则相反。究其原因,当地表呈现下沉即Dv朝下时,由图1中Dlos的矢量合成原理可知,靠近卫星的东西向水平移动会减小Dlos的量值,而远离卫星的东西向水平移动则会使之增大。
图5 三维算法得到的东西向位移以及不同LoS向的位移Fig.5 The displacement in east-west direction calculated using 3-D method and along different LoS
(1)提出了一种基于多源SAR影像联合求解矿区地表三维变形的模型与算法。通过与实测水准资料的对比证明,三维算法获取的垂直向位移曲线与水准测量结果的吻合程度优于传统算法,为利用多源SAR影像监测矿区变形提供了理论基础。
(2)由三维算法获取的垂直向位移场与等值线可以发现,下沉盆地向老采空区偏移,说明受工作面wf1的开采影响,附近的老采空区可能发生活化。
(3)利用三维算法获取的东西向位移场与等值线符合开采沉陷地表移动规律,并且结合东西向位移对不同LoS向的位移进行对比分析,结果表明,当地表呈现为下沉时,朝向卫星的东西向水平移动会减小LoS向的位移量值,而远离卫星的东西向水平移动则会使之增大。并且,雷达波入射角越大,对东西向水平移动的影响越大。
(4)因为卫星采用极轨飞行方式,使得航向角α的正弦值近乎为0,本文提出的三维算法对南北向位移不敏感,反演精度低。欧洲空间局(ESA)与日本宇宙航空研究开发机构(JAXA)为本文提供了SAR数据(C1P.8371, C1P.10185),国家科学数据服务平台为本文提供了SRTM-3数据,在此一并致谢。
[1] Gabriel A K,Goldstein R M,Zebker H A.Mapping small elevation changes over large areas:differentialdar interferometry[J].Journal of Geophysical Research,1989,94(B7):9183-9191.
[2] Currenti Gilda,NapoliRosalba,Del Nero Ciro.Toward a realistic deformation model of the 2008 magmatic intrusion at Etna from combined DInSAR and GPS observations[J].Earth and Planetary Science Letters,2011,312(1-2):22-27.
[3] Cascini Leonardo,Fornaro Gianfranco,Peduto Dario.Advanced lowand full-resolution DInSAR map generation for slow-moving landslide analysis at different scale[J].Engineering Geology,2010,112(1-4):29-42.
[4] Shirzaei M,Buergmann R,Oncken O,et al.Response of forearc crustal to the megathrust earthquake cycle:InSAR evidence from Mejillones Peninsula,Northern Chile[J].Earth and Planetary Science letters,2012,333:157-164.
[5] Parks M M,Biggs J,Mather T A,et al.Co-eruptive subsidence at GalerasidentifiedduringanInSARsurveyofColombian Volcanoes(2006-2009)[J].Journal of Volcanology and Geothermal Research,2011,202(3-4):228-240.
[6] Hsieh C S,Shih T Y,Hu J C,et al.Using differential SAR interferometry to map land subsidence:A case study in the Pingtung Plain of SW Taiwan[J].Naural Hazards,2011,58(3):1311-1322.
[7] Ferretti A,Prati C,Rocca F.Non-linear subsidence rate estimation using permanent scatterers in differential SAR interferometry[J].IEEE Transactions on Geoscience and Remote Sensing,2000,38: 2202-2212.
[8] Ferretti A,Prati C,Rocca F.Permanent scatterers in SAR interferometry[J].IEEE Transactions on Geoscience and Remote Sensing, 2001,39:8-30.
[9] Mora O,Mallorqui J J,Duro J,et al.Longterm subsidence monitoring of urban areas using differential interferometric SAR technique[A].Preoceedings of IEEE Geoscience and Remote Sensing Symposium 2001[C].Sydney:IEEE,2001:1104-1106.
[10] Berardino P,Fornaro G,Lanari R,et al.A new algorithm for surface deformation monitoring based on small baseline differential interferograms[J].IEEE Transactions on Geoscience and Remote Sensing,2002,40:2375-2383.
[11] Samsonov S,Beavan J,Gonzalez,P J,et al.Ground deformation in the Taupo Volcanic Zone,New Zealand,observed by ALOS PALSAR interferometry[J].Geophysical Journal International,2011,187(1):147-160.
[12] 张永红,吴宏安,孙光通.时间序列InSAR技术中的形变模型研究[J].测绘学报,2012,41(6):864-869.
Zhang Yonghong,Wu Hongan,Sun Guangtong.Deformation model of time series interferometric SAR techniques[J].Acta Geodaetica et Cartographica Sinica,2012,41(6):864-869.
[13] Via G D,Crosetto M,Crippa B.Resolving vertical and east-west horizontal motion from differential interferometric synthetic aperture radar:the L’Aquila earthquake[J].Journal of Geophysical Research,2012,117(B02310):1-14.
[14] Hanssen R F.Radar interferometry data interpreation and error analysis[M].New York:Springer,2001.
[15] 范洪冬,邓喀中,祝传广,等.基于时序SAR技术的采空区上方高速公路变形监测及预测方法[J].煤炭学报,2012,37(11): 1841-1846.
Fan Hongdong,Deng Kazhong,Zhu Chuanguang,et al.Study on deformation monitoring and predicting methods for expressway above goaf based on time series SAR technique[J].Journal of China Coal Society,2012,37(11):1841-1846.
[16] Fialko Y,Simons M.The complete(3-D)surface displacement field in the epicentral area of the 1999Mw 7.1 Hector Mine earthquake, California,from space geodetic observations[J].Geophysical Research Letters,2001,28(16):3063-3066.
[17] Aly M H,Cochran E S.Spatio-temporal evolution of yellowstone deformation between 1992 and 2009 from InSAR and GPS observations[J].Bulletin of Volcanology,2009,73(9):1407-1419.
Three-dimensional deformation field detection based on multi-source SAR imagery in mining area
ZHU Chuan-guang1,2,DENG Ka-zhong1,ZHANG Ji-xian2, ZHANG Yong-hong2,FAN Hong-dong1,ZHANG Li-ya3
(1.Key Laboratory for Land Environment and Disaster Monitoring of SBSM,China University of Mining and Technology,Xuzhou 221116,China;2.Chinese Academy of Surveying and Mapping,Beijing 100036,China;3.Hunan Province Key Laboratory of Coal Resources Clean-utilization and Mine Environment, Xiangtan 411100,China)
Only one dimensional displacement in the line of sight(LoS)can be detected using DInSAR technique.In order to overcome this limitation,the three-dimensional(3-D)model and algorithm which fusing multi-source SAR images acquired with different geometric parameters in different satellite platform was proposed and validated to determine the 3-D displacement due to underground mining.Comparing with leveling observation indicates that the displacement derives from the 3-D algorithm coincide with leveling data more closely,and the root mean square error (RMSE)is±4 mm.May be due to the activation of old goafs,both the displacement and contour in vertical direction indicates that the subsidence basin extends to the old goaf.The displacement and contour in east-west direction accord with the law of surface movement caused by mining,and the effects on the displacement along LoS will be significantly different when imaging mode different.Due to the small sine value of heading angle,the 3-D algorithm is not sensitive to the displacement in the north-south direction.
P237;TD173
A
0253-9993(2014)04-0673-06
祝传广,邓喀中,张继贤,等.基于多源SAR影像矿区三维形变场的监测[J].煤炭学报,2014,39(4):673-678.
10.13225/j.cnki.jccs.2013.0424
Zhu Chuanguang,Deng Kazhong,Zhang Jixian,et al.Three-dimensional deformation field detection based on multi-source SAR imagery in mining area[J].Journal of China Coal Society,2014,39(4):673-678.doi:10.13225/j.cnki.jccs.2013.0424
2013-04-07 责任编辑:王婉洁
江苏高校优势学科建设工程资助项目(SZBF2011-6-B35);国家自然科学基金资助项目(41272389,41071273)
祝传广(1984—),男,山东成武人,博士研究生。E-mail:zhuchuanguang@163.com。通讯作者:邓喀中(1957—),男,四川资中人,教授,博士生导师。E-mail:kzdeng@cumt.edu.cn