张永红,刘 冰,2,吴宏安,程 霞,3,康永辉
(1.中国测绘科学研究院摄影测量与遥感研究所,北京 100830; 2.山东科技大学 测绘科学与工程学院,山东 青岛 266590; 3.中南大学 地球科学与信息物理学院,湖南 长沙 410083)
2017年4月1日,中共中央、国务院印发通知,决定设立河北雄安新区。雄安新区是继深圳经济特区和上海浦东新区之后又一具有全国意义的新区,是千年大计、国家大事。雄安新区被明确要建设成为绿色生态宜居新城区、创新驱动发展引领区、协调发展示范区、开放发展先行区[1]。雄安新区的环境承载容量是关乎其历史使命能否实现的重要因素。
雄安新区地处京津保腹地,位于华北平原中心。华北平原因长期地下水超采是世界上地面沉降漏斗最大、沉降覆盖面积最大的地区[2]。雄安新区的地面沉降状况是评估其环境承载容量的重要因子。已发表的地面沉降研究文献中完整覆盖雄安新区的主要有葛大庆等[3]及张永红等[4]。前者展示了中国中东部地区2012年1月至2015年10月RADARSAT-2影像的地面沉降速率图,该图比例尺较小,显示了华北平原大的沉降中心分布,但缺乏相关地面沉降状况的定量描述[3]。后者利用ERS、ENVISAT、RADARSAT-2多颗卫星时序SAR数据研究了京津冀地区1992~2014年的地面沉降,给出了京津冀平原区(包括雄安新区)2012年2月至2014年7月的沉降速率图[4],但对雄安新区地面沉降没有专门的论述。Zhu等对更早时段也有研究[5-7],但研究区都在雄安新区外围。中国地质调查局在最新公布的雄安新区地质调查结论中指出“区内场地稳定性和工程建设适宜性总体较好,…,但应关注地面沉降问题”[8]。由此可见,地面沉降对雄安新区的建设有直接影响,但从公开发布的文献资料来看,目前针对雄安新区地面沉降的专门研究还是空白。
本文以2012年1月至2016年11月获取的28期RADARSAT-2影像为遥感数据源,采用自主研发的多主影像相干目标小基线InSAR(MCTSB-InSAR)技术[4,9],提取了雄安新区在该时段内的平均沉降速率和时间序列累计沉降量,分析了该区域的地面沉降原因,对该区域地面沉降危险性进行了评价,以期揭示雄安新区较详细的地面沉降状况。
本研究使用了28期时间序列RADARSAT-2宽模式影像,数据格式为单视复数据(SLC数据),像元大小分别为11.83 m(距离向)和5.21 m(方位向),雷达波中心入射角约为26.2°。影像的获取时间跨度从2012年1月28日至2016年11月14日,RADARSAT-2影像时空基线参数见表1。影像覆盖区面积约为170 km×155 km,包括天津南部,河北沧州、衡水、保定、廊坊部分区域(图1)。图1中方框为RADARSAT-2影像覆盖区,以Landsat TM影像为背景。为了降低斑点噪声,增加相干性,对单视复数据做了5×1的多视处理,多视处理后像元大小为25 m(距离向)和25 m(方位向)。本文选用了30 m格网间距的SRTM 数字高程数据去除地形相位。另外,从天津市测绘院获取了研究区内天津市域2012~2015年154个水准点的观测数据,用于精度评价。水准点分布位置见图1。
表1 RADARSAT-2影像时空基线参数Tab.1 Spatial and Temporal Baselines of RADARSAT-2 Images
图1 雄安新区地理位置Fig.1 Location of Xiong’an New Area
MCTSB-InSAR技术主要包括小基线干涉像对选择、高相干点提取、高相干点网络连接和形变参数反演4个主要步骤。在此仅介绍本文数据的处理过程,更多细节请参阅文献[4]和[9]。
图2 干涉像对连接图Fig.2 Connections of Interferograms
在干涉像对选择方面,依据影像集的时空基线分布状况,以干涉像对连接图必须连通[4]和干涉相干性最优为原则,确定了干涉像对的时间基线不超过180 d、垂直基线小于250 m的时空基线阈值,共形成了66个干涉像对(图2)。本文采用平均相干系数、幅度离差[9]、平均幅度等参数阈值原则,共提取了约6.78×105个高相干点,这些高相干目标主要分布在建筑物、道路、桥梁、铁路及公共设施等地物上。对于提取的高相干点,采用局部Delaunay三角网连接[10],分块大小为2 074×1 978,重叠度为50个像元,共得到4 509 200边。
将地表形变分解为线性形变(主形变)和非线性残余形变,对第i幅干涉图上两相邻高相干点(设为m、n)构成的边上干涉相位差可表示为
Δφm,n(Ti)=Δδm,n(Ti)+Δβm,n(Ti)+
Δαm,n(Ti)+ΔNm,n(Ti)
(1)
(2)
为求解未知数Δhm,n与Δvm,n,建立模型相干系数方程[10-11]。其表达式为
(3)
式中:γm,n为相邻高相干点的模型相干系数;M为干涉图数目;j为中间参数。
式(3)中模型相干系数γm,n取得最大值的高程误差之差与线性形变速率之差为正确解。求解这两个参数的方法有二维周期图法[12]或空间搜索法[13-14]。对Delaunay三角网上所有边完成最大化求解后,将γm,n≥0.7的边作为可靠的连接关系,并将模型相干系数γm,n作为权值,利用某一参考点的线性形变速率和高程误差,采用加权最小二乘方法解算得到各高相干点上线性形变速率和高程误差的绝对量。从差分干涉相位中去除高相干点的线性形变相位和高程误差相位,得到高相干点的残余相位,包括大气影响相位、非线性形变相位以及噪声相位。根据残余相位3个分量的不同时空频谱特征,利用时空域的滤波方法可以将三者分离出来[12],并得到每个SAR成像时刻对应的非线性形变(以第一期影像为时间参考)。将线性形变和非线性形变累加到一起,可以得到每个时刻的雷达视线方向(LOS)形变量。为了方便和水准数据进行比较,假设沉降是地面形变的主要分量,将雷达视线方向形变量转化为垂向沉降量。转化公式为
(4)
式中:dV为垂向沉降量;dR为雷达视线方向形变量;θ为雷达视线入射角。
通过MCTSB-InSAR技术提取了整景SAR RADARSAT-2影像覆盖区2012年1月至2016年11月平均沉降速率信息(图3),图3中正值代表地面抬升,负值代表地面沉降。整个区域内最大沉降速率为每年184 mm,位于河北省廊坊市胜芳镇;与廊坊市沉降中心相邻的是天津市武清区王庆坨镇的沉降漏斗,最大沉降速率为每年181 mm;河北平原的保定东南部、衡水北部及沧州西部有连接成片的大面积沉降带,沉降速率大都在每年45 mm以上。相对而言,雄安新区大部分区域的沉降速率在每年40 mm以内。
将本文RADARSAT-2影像覆盖区监测结果与相关文献进行了比较。张玲等研究了河北沧州地区2012年地面沉降速率,指出沧州西部存在任丘—肃宁沉降带,肃宁最大沉降速率超过每年100 mm[15];从本文监测结果来看,任丘、肃宁是沧州西部的主要沉降中心,肃宁2012~2016年最大沉降速率为每年89 mm,这可能暗示2012年后肃宁地面沉降有所减缓。Zhu等利用2009年3月至2013年3月的TerraSAR-X影像获取了天津西部沉降信息,其中最大沉降速率出现在王庆坨镇,为每年137 mm[5];唐嘉等利用2007年1月至2010年10月的PALSAR影像和2007年2月至2009年6月的ASAR影像分别提取了胜芳镇和王庆坨镇的沉降速率,胜芳镇最大沉降速率为每年206.9 mm(PALSAR影像)、211.6 mm(ASAR影像),王庆坨镇最大沉降速率为每年120.2 mm(PALSAR影像)、125.5 mm(ASAR影像)[7];本文结果显示,2012~2016年胜芳镇和王庆坨镇的最大沉降速率分别为每年184、181 mm。比较三者,虽然时间范围不一致,但都显示胜芳镇和王庆坨镇是该地区沉降最严重的区域,进而可以推断胜芳镇沉降2010年以后呈现一定程度的减缓,王庆坨镇沉降在2013年以后加剧。总体来说,这些用不同时段不同卫星SAR影像提取的沉降信息显示了较好的一致性,这也从侧面验证了本文研究结果的可靠性。
图3 SAR RADARSAT-2影像覆盖区2012年1月至2016年11月平均沉降速率Fig.3 Average Subsidence Rates Derived from SAR RADARSAT-2 Image from January 2012 to November 2016
图4 85个水准点与最邻近高相干点的沉降速率差值直方图Fig.4 Histogram of 85 Difference Values of Subsidence Rates Derived from InSAR and Levelling Measurement
利用从天津市测绘院获取的154个水准数据对MCTSB-InSAR技术获取的地面沉降结果进行精度检验。这些水准点的分布见图1,天津市测绘院每年第四季度(10月至12月)施测一次。需要说明的是,本文使用的154个水准数据为2012~2015年的,与SAR数据的获取时段并不完全重合,对精度评估会有一定影响。将水准点逐年的沉降量累加然后除以时间,得到2011~2015年的平均沉降速率,然后选择那些半径80 m范围内有至少一个高相干点的水准点,将水准点的平均沉降速率与最邻近高相干点的线性沉降速率进行对比。这样的水准点共85个,二者沉降速率差值直方图分布见图4,其中绝大部分的差值在每年[-5 mm,5 mm]区间内,差值绝对值超过每年8 mm的点共有9个。分析发现,这些差值较大点上的沉降速率均在每年75 mm以上。造成这些点上InSAR监测与水准测量结果相差较大的原因除前述的观测时段不一致外,水准测量本身的时点误差也是重要因素。据了解,天津测绘院的水准测量每年大约从10月开始,12月底结束,全部完成后要经过一次平差处理,对某个水准点而言,它的测量时刻有约2个月的时差,对于沉降速率较大的情况,2个月的时差将对水准测量的沉降速率估算产生较大影响。差值的平均值已用于图3所示的沉降速率的标定,差值的标准差为每年5.2 mm,表明MCTSB-InSAR技术反演的沉降速率精度较高。
图5 2012年1月至2016年11月平均沉降速率Fig.5 Average Subsidence Rates from January 2012 to November 2016
图6 2012年1月至2013年8月累计沉降量序列Fig.6 Sequence of Accumulative Subsidence from January 2012 to August 2013
图8 2012年1月至2016年11月累计沉降量序列Fig.8 Sequence of Accumulative Subsidence from January 2012 to November 2016
将RADARSAT-2影像覆盖区内的雄安新区(包括雄县、容城县、安新县)沉降监测结果抽取出来,分别显示于图5(平均沉降速率)和图6~8(累计沉降量)。从沉降区域的空间分布来看,雄县沉降区域较大,主要分布在大营镇、雄州镇、北沙口乡、米家务镇、朱各庄乡及昝岗镇等地区,最大沉降速率为每年75 mm,位于雄县大营镇西昝村;容城县容城镇沉降较大,最大沉降速率达每年28 mm,其周边区域沉降速率较小,不超过每年15 mm;安新县绝大部分区域没有发生地面沉降,只在南部的芦庄乡和刘李庄镇形成了相对较大的地面沉降,沉降速率最大值分别为每年35 mm和37 mm。从地面沉降的演化进程来看,雄县大营镇至北沙口乡一带的地面沉降从2013年4月开始逐步加剧,此后随时间推移,沉降速率和沉降范围不断扩展,至2016年11月,最大累计沉降量达到411 mm;而安新县南部的地面沉降发展较缓慢,直到2013年底才有明显的沉降,到2016年11月,安新县刘李庄镇和芦庄乡最大沉降量分别为235 mm和218 mm。由此可见,雄安新区地面沉降随时间推移呈现不断加剧的趋势。
经分析有关资料,推断地热资源的开发是引发雄县出现较大沉降的最主要原因。雄县位于华北平原牛驼镇地热田范围内,全县约六成面积下蕴藏着地热资源,地热水储量达821×108m3[16]。雄县从2003年就开始使用地热资源,随着经济的发展,地热井的开凿量也随之增加,至2017年2月,雄县共有地热井68眼[16-17]。其中,大营镇、雄州镇地热开采条件较好,地热井数量较多[18],沉降也较大;北沙口乡、米家务镇、朱各庄乡及昝岗镇西部也较适宜开发地热资源,地热井数量次之;而位于雄县东部的双堂乡、张岗乡及龙湾乡地热可持续开发条件较差,地热井数目也较少,地面几乎无沉降。图9(a)为截止至2009年底的雄县地热井位置分布示意图[18],图9(b)为2012~2016年雄县累计沉降量,可见地热井的位置分布和累计沉降量分布之间存在着很强的空间相关性。
造成雄安新区地面沉降的另一个因素可能是塑料包装企业及纺织企业消耗了大量的地下水。雄县是华北地区最大的塑料包装产业基地[19],全县的塑料包装企业有3 000多家,主要集中在包装印刷、塑料制品和颗粒销售。包装印刷及塑料制品都是高耗水行业。鉴于华北平原长期处于水资源短缺状态[2],推测雄县规模巨大的塑料包装产业很可能导致地下水开采量增多,进一步加剧了雄县的地面沉降。安新县刘李庄镇位于白洋淀南部,也拥有多达47家塑料加工企业;芦庄乡是安新县纺织厂的聚集地,乡镇内分布着大量的纺织厂,而纺织业也是高耗水行业之一。因此,有理由推断安新县刘李庄镇及芦庄乡两地的地面沉降主要是塑料包装业及纺织业大量消耗地下水所致。
图9 雄县地热井位置分布及累计沉降量Fig.9 Geothermal Well Location and Acuumulative Subsidence of Xiongxian
在雄安新区高相干点沉降监测结果的基础上,采用克里金(Kriging)插值方法将高相干点上的沉降速率值内插至每个像元,可以得到整个雄安新区完整的沉降速率图,然后根据杨艳等对北京市郊区地面沉降危险性分级评价标准[20],将雄安新区地面沉降危险性评价分为严重区(沉降速率大于每年50 mm)、较严重区(沉降速率为每年30~50 mm)、一般区(沉降速率为每年10~30 mm)和轻微区(沉降速率小于每年10 mm)(图10)。从图10可以看出,雄安新区约78%以上的区域为地面沉降轻微区,地面沉降较严重区和严重区都集中在雄县,面积分别为55.76 km2和12.14 km2,二者合计仅占雄安新区面积的4.38%(表2)。
图10 雄安新区地面沉降危险性分区Fig.10 Ground Subsidence Risk Map of Xiong’an New Area
危险性分级面积/km2占比/%严重区 12.140.78较严重区55.763.60一般区263.1917.01轻微区1 216.0378.61
(1)利用2012年1月至2016年11月28期RADARSAT-2影像,基于多主影像相干目标小基线InSAR技术提取了雄安新区及周边的平均沉降速率和累计沉降量信息。利用85个水准观测数据对InSAR获取的平均沉降速率进行了精度验证,标准差为每年5.2 mm。
(2)雄安新区约78%以上的区域为地面沉降轻微区(沉降速率小于每年10 mm)。雄县部分区域发生了较大的地面沉降,容城县容城镇及其周边区域发生了较小的地面沉降,安新县南部的刘李庄镇及芦庄乡小部分地区发生了相对较大的地面沉降。雄安新区地面沉降危险性评价为较严重和严重的区域都集中在雄县,共67.9 km2,占雄安新区面积的4.38%,最大沉降速率为每年75 mm。整体而言,雄安新区地面沉降呈现不断加大的趋势。
(3)雄安新区地面沉降的主要原因是地热资源开采和包装印刷业及纺织业等导致的地下水开采。建议未来对雄县的地热资源开采和地下水开采加以控制,以遏制雄安新区地面沉降的进一步发展,为雄安新区发展战略的顺利实施提供良好的生态环境。
本文所用水准数据为天津市测绘院提供,特此致谢!