不同地表覆盖对InSAR技术的相干性影响研究

周鑫城,左小清,李勇发,杨 栩,邓云龙,周家厚

(昆明理工大学 国土资源工程学院,云南 昆明 650093)

为了解地表的沉降情况,众多学者采用持续性的形变监测技术对地表进行观测,从而揭示地表的活动状态,预测与防范地质灾害的发生[1]。当前,常规的地表形变监测技术主要有:水准测量[2]、全球卫星导航定位技术(global positioning system,GPS)[3]等。近些年来,随着地表监测技术水平和科学技术的发展,诸如全球导航卫星系统(global navigation satellite system,GNSS)[4]、三维激光扫描[5]、地理信息系统(geographic information system,GIS)[6]等高新技术的引入,使得地表形变监测的工作效率有显著的提高。这些技术的应用,为区域的沉降情况提供了可靠的数据基础,但存在一定的缺陷,如空间分辨率低、观测周期长和监测成本高等问题。随着雷达卫星(synthetic aperture radar,SAR)逐步发育成熟,具有高时空分辨率、不受云雨条件限制、覆盖范围广、数据获取周期较短、安全性高等优势的干涉合成孔径雷达(interferometric synthetic aperture radar,InSAR)成为地表形变监测的主要手段之一[7-9]。相干性[10-11]是干涉测量的基础,不仅直接影响形变监测数据的处理难度,而且失相干也是地表形变误差的主要来源。InSAR失相干[12]包括时间失相干和空间失相干,而失相干主要是由时间基线或空间基线过长引起的。目前,多数学者从定性的角度进行研究,发现InSAR相干性与时空基线的长度呈负相关[13-14]。欧洲航天局的哨兵双星系统的出现极大地克服了大部分空间失相干所带来的误差,而引起时间失相干的主要原因是研究区分辨单元内目标的物理和化学特性随时间的变化而变化。

当前,研究不同地表对相干性的影响较少,余祥伟等[15]分析植被及坡度对InSAR相干性的影响,发现相较于坡度,植被是导致高覆盖山区InSAR影像时间失相干的主要因素,相干性和植被覆盖度之间存在一定的相关性。佀文娜等[16]研究典型地物的时间失相干影响,发现风速对相干性没有影响,而降雨量是导致地物失相干的重要因素。王天祥等[17]使用D-InSAR技术对地表进行形变监测,在获取大量数据后进行仿真试验,发现空间基线和时间基线与InSAR相干性之间存在指数关系,而且当空间基线的长度不超过临界基线的三分之一时,可以保证得到较好的InSAR相干性。

现有研究主要对植被、坡度、降雨量等多方面进行综合分析,尚未针对不同地表覆盖对InSAR的相干性影响进行深入研究。因此,本文选取昆明市滇池及周边环境作为研究区,采用小基线集(small baselines subset InSAR,SBAS-InSAR)技术对欧洲航天局哥白尼计划(global monitoring for environment and security,GMES)中的地球观测卫星哨兵1号(Sentinel-1)卫星数据进行处理。从定量角度分析不同地表覆盖对相干性的影响,深入分析在不同月份、不同季节和不同降雨量下,地表覆盖类型对相干性的影响,为提高InSAR形变精度提供有力的理论依据和数据支撑。

1 研究方法

1.1 SBAS-InSAR技术

SBAS-InSAR技术是在InSAR技术基础上发展起来的一种多主影像的时序InSAR方法。该技术的测量原理是利用两幅具有相同观测几何的影像构成干涉像对,结合卫星和观测目标之间的几何关系,对干涉像对所产生的相位差进行相位滤波、相位解缠、地理配准等步骤,获取观测目标的三维信息和形变信息[18]。

在评价干涉质量的好坏时,分析相位相干性常被作为一个重要的标准。相干性被定义为一个干涉像对中主、副影像是否存在相似,而两者之间的相似程度目前常使用相干系数进行描述。相干系数(γ)数值分布于[0,1]区间内,相干系数为0表示主、副影像之间无相似程度,相干系数为1表示主、副影像之间完全相同。相干系数(γ)计算模型如下:

式中,N和M为最终获取相干性像元的尺寸大小;n和m为数据内像元的行列号;μ*为共轭复数;μ1(n,m),μ2(n,m)为主、副影像数据内像元再坐标(n,m)处的复数值;|·|为复数的绝对值。

由于雷达两次成像期间观测区域内目标的物理和化学特性随时间变化而变化,使得相干性在目标条件不稳定的区域通常较低,而在目标条件稳定的区域(建筑区域及工矿用地)相干性通常较高。

1.2 归一化水指数

归一化水指数是用遥感影像中的特定波段进行归一化差值处理,以突显影像中水体的信息[19-20]。该指数的主要优点之一是能够很好的区分水体区与非水体区,其利用水体吸收的绿光(492～577 nm)和反射近红外(780～2 526 nm)的原理建立表达式,计算公式如下:

(1)

式中,B3为遥感影像中的绿波段;B5为遥感影像中的近红外波段。INDWI的取值范围为[-1,1]。

1.3 垂直不透水面指数

垂直不透水面指数是利用所选样本的光谱特性建立较高的关联表达式,以凸显影像中不透水面(建筑区域及工矿用地)的信息[21-23]。该指数相较于归一化建筑指数能够较好的区分不透水面和透水面,计算公式如下:

(2)

式中,a为选取样本构造出的关联表达式的斜率;b为选取样本构造出的关联表达式的截距;B2为遥感影像中的蓝波段。

1.4 归一化植被指数

归一化植被指数是反应植被生长趋势的重要参数之一,被广泛应用于干旱季的监测、农作物生产情况的监测和预测等[24-25]。该指数的主要优点之一是能够较好的区分植被区域与非植被区域,其利用植被吸收红光(625～740 nm)和反射近红外(780～2 526 nm)的原理建立表达式,计算公式如下:

(3)

式中,B4为遥感影像中的红波段。INDVI的取值范围为[-1,1]。

2 研究区与数据集

2.1 研究区

云南省昆明市总体的地貌趋势为北部高,南部低,呈阶梯式降低,高程主要分布在1 500～2 800 m之间;以亚热带高原季风气候为主,因此,植被以常绿阔叶林为主,年平均气温为15℃左右,年日照平均时长为2 200 h左右,年平均降水量为1 035 mm。经全国第三次全国国体调查,昆明市耕地面积约为3 865.73 km2;林地面积约为11 317.67 km2;草地面积约为1 244.40 km2;建筑区域及工矿用地面积约为1 480.53 km2。本次研究区包括云南省昆明市滇池及周边区域,包括安宁区的东部地区,地处中国西南地域、云贵高原中部地区,其覆盖研究区覆盖面积约为2 287 km2,位于东经:102°28′05″E—102°57′40″E,北纬:24°36′34″N—25°13′52″N,如图1所示。

图1 研究区位置Fig.1 Location of study area

2.2 数据集

本文使用的数据集包括从欧空局(https://search.asf.alaska.edu/)获取2020年1月11日—2021年2月22日的35景降轨Sentinel 1A数据,该数据的波段为C波段,波长为5.6 cm,分辨率为5 m×20 m,观测模式为干涉宽幅 (interferometric wide swath,IW),极化方式为VV(单极化),单视复数(single look comple,SLC)。土地利用分类数据来源于中国科学院资源环境科学与数据中心数据平台 (http://www.resdc.cn)。从地理空间数据云(http://www.gscloud.cn/)网站下载得到的2020年2月5日Landsat 8影像。2020年全年MODIS系列NDVI产品的原始数据是从Nasa(https://www.nasa.gov/)网站下载得到的2020年全年MODIS数据。降雨数据是从中国地面气候资料日值数据集中获取。

2.3 数据预处理

本文通过使用SBAS-InSAR技术对35景Sentinel 1A降轨数据进行处理,将SBAS-InSAR的时间基线阈值设置为120 d,经过数据导入、研究区裁剪、地理配准、反演等步骤,最终得到204幅相干性图,并将获取相干性的两个日期对相干性图进行命名,每幅影像的分辨率为17 m×17 m,地理坐标为WGS_1984。

本文对地理空间数据云中2020年全年的影像进行筛选,最终使用2020年2月5日的Landsat 8影像(云量小于5%)作为研究影像,利用光谱特性与三种指数(INDWI、IPII和INDVI)的计算公式,计算出研究区内各像元在2020年2月5日的三种指数值,其分辨率为30 m*30 m,投影坐标为WGS_1984_48N。

为研究不同类型的地表覆盖与InSAR相干性之间的关系,以三种指数的像元分辨率和坐标为准,对相干性数据进行匹配、投影转换、重采样、掩膜提取等步骤使4种数据的行列号保持一致。对各类型的地表覆盖随机选择1 200个点,分析所选点的指数数据与相干性之间的关系,以确定各类型的地表覆盖是否存在研究的意义,并进一步了解不同地表覆盖对相干性的影响。

MODIS系列的归一化植被指数产品在获取时间序列上的间隔为16 d,因此,需对该数据进行预处理,获取每月的月度归一化植被指数数据。降雨量是日值数据,同样也要对此进行月份划分,获取月度降雨量数据。相干性数据是由每两幅不同时间段的影像计算得到的,而本文所使用的时间基线为120 d,因此,会出现跨月时间影像生成的相干性图,针对该现象,本文将命名中较早的日期作为分类依据,进行归纳各月的相干性图,以方便计算各月的相干性平均值(如:2020年3月11日的影像和2020年4月4日的影像所计算出来的相干性图归纳为三月份的相干性图,2020年3月11日的影像和2020年6月3日的影像所计算出来的相干性图也归纳为三月份的相干性图)。最终结合降雨量数据与季节数据,研究分析在降雨量与季节变化的条件下,不同地表覆盖对InSAR相干性的影响。

3 结果与分析

3.1 三种指数结果与分析

图2(a)为求取平均值后的相干性图,从该图可知,研究区内的相干性取值范围为0～1,其中,建筑区域的相干性最高,其次是植被覆盖区域的相干性,水体的相干性最低。本文的土地利用一级类别有五类,分别为:水域、林地、草地、耕地和建筑区域及工矿用地,其二级分类如图2(b)所示。统计一级类别的像元,可以发现建筑区域及工矿用地的像元最多,其次为林地。

(a)求取平均值后相干性图 (b)二级分类图2 相干性与土地利用分类图Fig.2 Coherence and land use classification map

图3(a)为研究区的INDWI指数图,由图可知,其INDWI取值范围为-0.61～1。图3(b)为研究区的IPII指数图,IPII指数的获取与INDWI指数不同。在影像中选取暗色不透水面(居民建筑区等)、亮色不透水面(工厂等)、蓝色不透水面(厂房等)、透水面(林地、草地、耕地等)5种类别作为样本,根据样本生成的特征空间,计算出(2)式中的a和b分别为1.514和2 118.066,因此,(2)式可以表示为

图3 INDWI、IPII、INDVI指数图Fig.3 INDWI, IPII and INDVI index diagram

由图3(b)可知,IPII数值的取值范围为-7 299.6～8 065.1。图3(c)为研究区的INDVI影像图,其取值范围为-1～0.82。

图3中,三种指数可以明显的区分出水体区域、建筑区域和植被覆盖区域。为定量分析不同地表覆盖与相干性之间的关系,利用随机点绘制散点图,并通过最小二乘法进行拟合,最终得到三种指数与相干性之间的关联程度图,如图4所示。由于每种地表覆盖类型的光谱特性不同,用不同的指数表达不同地表覆盖可以更好的体现地表覆盖的差异性,因此,本文对不同地表覆盖的随机点赋予不同的指数。INDWI指数能够很好的区分水体区域与非水体区域,所以对水体类型的随机点赋予INDWI数据;IPII指数主要用于区分透水面和不透水面,而透水面主要是植被覆盖区域,不透水面主要是建筑区域,因此,对两者区域赋予IPII数值能够进一步了解两者与相干性之间的关系,由于建筑区域及工矿用地的建设需符合国家绿化面积要求,其周边附着一定的植被及绿化带,而在农村居民点及城镇周边,多为植被覆盖区与建筑区的交汇,影像中存在混淆像元,因此,为减少此方面的误差,对建筑区域及工矿用地与植被覆盖区的随机点赋予INDVI数值。

图4 INDWI、IPII、INDVI指数与相干性的关联程度Fig.4 Correlation degree of INDWI, IPII, INDVI index and coherence

从整体上看,图4中存在一些离散点,出现离散点的主要原因是三种指数由影像通过波段差计算得到的,而相干性的计算方式是由两幅Sentinel 1A影像通过干涉处理得到的,在影像源的使用和计算方法上存在差异;各类型地表覆盖的随机点是在土地利用的分类基础上进行选择的,而土地利用在分类上存在误差,因此,在数值上会出现异常值,在图像上出现无规律分布的离散点。如图4(a)所示,结合光学遥感影像和土地利用类型可以发现,该点的利用类型为水体,而在遥感影像中,该点位于昆明南绕城高速上,且从相干性和INDWI的数值大小可以断定,该点属于异常点。因此,研究三种指数与相干性之间的关系前,需对离散点进行剔除。

图4(a)为INDWI指数相关性图,图4(b)为IPII指数相关性图,图4(c)为INDVI指数相关性图。由图4(a)可知,INDWI数值集中分布在0.75～0.85之间,其对应的相干性分布在0.05～0.58之间,两者的相关系数R2为0.180 6,呈负相关。且从图中的点分布可知,INDWI数值与相干性之间无关联程度,故研究水体对相干性的影响没有意义,后续将不对水体进行研究分析。

由图4(b)可知,不透水面(建筑区域及工矿用地)的IPII数值主要分布在0～2 000之间,其对应的相干性主要分布在0.4～0.9之间;透水面(植被覆盖区)的IPII数值主要分布在-2 000～0之间,其对应的相干性数值主要分布在0.4～0.8之间,两者的相关系数R2为0.180 3,呈正相关。导致IPII指数与InSAR相干性之间的关联程度低有两个原因,一是本文所选取的随机点多为植被覆盖类型,所以在图中点分布较为密集,呈现出上部分稀疏下部分密集的分布情况;二是IPII数值取值范围远大于InSAR相干性。虽然在图中呈现出随着IPII数值的增加,相干性增加,但是两者之间的关联程度低,无法显著的体现出不同地表对相干性的影响。

从图4(c)可知,植被覆盖区域的INDVI数值主要分布在0～0.55之间,相干性主要分布在0.4～0.8之间;建筑区域及工矿用地的INDVI数值主要分布在-0.2～0之间,相干性主要分布在0.6～0.98之间,两者的相关系数R2为0.399 8,呈负相关。结合INDVI指数特性,不同地表覆盖的INDVI取值范围不同,即可以通过INDVI指数区分不同地物。故本文将通过使用INDVI作为主要的研究指数,研究不同地表覆盖对相干性的影响。

3.2 不同地表覆盖类型对相干性影响分析

为研究不同地表覆盖对相干性的影响,本文对4种地表覆盖类型的1 200个随机点赋予2020年月度平均INDVI数值和月度平均相干性,将随机点的INDVI数值以0.2的步长进行划分,绘制各类型的地表覆盖在不同INDVI数值下相干性的分布图,如图5所示。

图5 不同地表覆盖类型的相干性分布图Fig.5 Coherence distribution map of different surface Types

对比不同地表覆盖的整体相干性可以发现,建筑区域及工矿用地的整体相干性较高,植被覆盖区的整体相干性较低。由图5可知,建筑区域及工矿用地的相干性在INDVI数值为0.34～0.57的区间内下降程度最大,而在0.18～0.34区间内,相干性的衰减程度平稳,由此可见,在建筑区域及工矿用地中,相干性受INDVI数值变化影响显著,特别是在INDVI数值超过0.34之后,随着INDVI数值增加,相干性以较大的趋势下降。林地的相干性在INDVI数值为0.35～0.66的区间内快速下降,而在0.27～0.35和0.66～0.81的区间内下降程度比较平稳,可以确定林地的相干性受INDVI数值的变化影响较大,特别是在INDVI数值为0.35～0.66之间,之后,随着INDVI数值增加,相干性将以较平稳的趋势下降。耕地的相干性在INDVI数值为0.36～0.45区间下降趋势较小,在0.27～0.36和0.45～0.81的区间内下降趋势程度最大,由此可知,耕地的相干性受INDVI数值变化影响明显;而草地的相干性受INDVI数值的增加而降低,且以较相同的趋势下降,由此可见,草地的相干性变化趋势受INDVI数值的影响较小。

进一步对比4种地表覆盖类型的相干性变化趋势,可以发现,在INDVI数值为0.18～0.55之间,不同地表覆盖的相干性变化趋势不相同,建筑区域及工矿用地的变化趋势最大,其次是林地,而耕地与草地的变化趋势相似;在INDVI数值为0.55～0.65区间内,草地的变化趋势低于林地与耕地;在INDVI数值为0.65～0.75区间内,草地的变化趋势高于林地与耕地;INDVI数值大于0.75,三种植被类型的相干性变化趋势一致,由此可知,当INDVI数值大于0.75之后,相干性不再受地表覆盖类型的影响,而只受植被覆盖率影响,且随着INDVI数值增加而降低。

3.3 基于各月不同地表覆盖对相干性的影响分析

为研究不同月份下,地表覆盖类型对相干性的影响,本文提取2020年月度平均INDVI数据和平均相干性数据绘制图6,图中的数值刻度为相干性。

图6 不同地表覆盖的月季变化Fig.6 The change of Rose in different surface covers

从图6中可以发现,在各月份中,建筑区域及工矿用地的相干性最高,且不低于0.7,其次是耕地和草地,林地的相干性最低,其主要原因是建筑区域及工矿用地的INDVI最低,林地的INDVI最高。从整体上看,该图的整体趋势靠近1月与12月,偏离6月与7月。从整体的相干性变化趋势上可以发现,在11月至2月,建筑区域及工矿用地的相干性变化趋势较为平缓,其次是林地,草地与耕地的变化趋势最大;而在5月至8月,建筑区域及工矿用地、林地和草地的变化稳定,耕地的变化较大。结合不同地表覆盖类型的植被覆盖度可以发现,在11月至2月建筑区域及工矿用地的整体植被覆盖度低,林地的整体覆盖度高,因此,整体变化趋势较小,而耕地与草地的植被覆盖度高于建筑区,低于林地,因此,在该阶段内,植被覆盖率增加较快,相干性变化大;在5月至8月,建筑区、林地与草地的植被生长情况稳定,耕地在该阶段属于丰收期,因此,相干性变化较大。由此可知,植被覆盖度对相干性影响较大,且对于不同地表覆盖,对相干性的影响不同。

为进一步分析不同地表覆盖类型的相干性在各季节中情况,按照节气的划分规则,将12个月份划分成春、夏、秋、冬四个季节,提取对应的相干性求取平均值,并加入研究区内的降雨量数据绘制出图7。

图7 不同地表覆盖的四季相干性变化Fig.7 Seasonal coherence change of different surface covers

图7中(a)图为耕地,(b)图为林地,(c)图为草地,(d)图为建筑区域及工矿用地。由图7可知,建筑区域及工矿用地的相干性变化趋势较小,耕地、林地和草地的相干性变化趋势较大。从INDVI数值的季节分布情况,可以发现秋季的INDVI数值最高,其主要原因是研究区以亚热带高原季风气候为主,植被以常绿阔叶林为主[26]。对比不同地表覆盖的四季相干性分布情况,可以发现不同地表覆盖的相干性在冬季最高,其次是春季,秋季和夏季最低。其主要原因是冬季的植被处于落叶期,植被覆盖率低,而且降雨量少,因此,冬季相干性高;春季的植被处于生长期,植被覆盖率和降雨量均大于冬季的情况,因此,冬季的相干性大于春季的相干性。由此可以发现,相干性受降雨量和植被覆盖率共同影响。对比夏季和秋季的情况可以发现,夏季的植被覆盖率低于秋季的植被覆盖率,而夏季的降雨量高于秋季的降雨量,从图中可以发现夏季的相干性略低于秋季的相干性。为进一步了解植被覆盖率和降雨量对相干性的影响,本文使用灰色关联度进行对比分析,如表1所示。

表1 不同地表覆盖的INDVI与降雨量之间的灰色关联度分析Tab.1 Grey correlation analysis between INDVI and rainfall under different surface covers

由表1可以发现,不同地表覆盖的相干性受INDVI影响较大,降雨量较小,对比4种不同地表覆盖的降雨量关联度可知,4种不同地表覆盖的相干性受降雨量的影响一致,数值在0.4～0.412之间;INDVI对相干性的影响相差较大,林地受INDVI影响较大,其余类型的地表覆盖类型受INDVI的影响大致一致,数值在0.87左右。

4 结论

相干性是InSAR技术广泛应用于形变监测与三维重建的基础,相干性的好坏决定着InSAR技术最终结果的精度。本研究以覆盖昆明市滇池周边的Sentinel 1A数据为例,基于不同指数研究地表覆盖与相干性之间的关联程度,将效果最好的指数为依据,进一步分析不同地表覆盖对相干性的影响及在不同月份、季节和降雨量下对相干性的影响,得出以下结论:

1)INDVI指数能够显著体现出不同地表覆盖与相干性之间的关系。随着INDVI增加,不同地表覆盖的相干性以不同的趋势下降。

2)在不同月份下,所有地表覆盖类型的相干性在1月与12月最高,6月、7月和8月的相干性最低。且在一年的时间内,建筑区域及工矿用地的相干性最高,其次是耕地与草地,林地的相干性最低。

3)在四季的变化中,建筑区域及工矿用地的相干性变化趋势较平缓,植被覆盖区域的相干性变化趋势较大,且春、冬两季的相干性较高,夏、秋两季的相干性低。

4)通过灰色关联度的分析可知,INDVI对不同地表覆盖的相干性影响较大,降雨量对相干性的影响较小。相比于其他地表覆盖,INDVI对林地的相干性影响最大。