利用Sentinel-1A数据的纳木错湖面月际变化监测

常亚茹，张继贤，韩文立，丁夏萌，张越

(1.山东科技大学 测绘与空间信息学院，山东 青岛 266590；2.国家测绘产品质量检验测试中心，北京 100036；3.兰州交通大学 测绘与地理信息学院，兰州 730070)

0 引言

气候变化已成为目前全球面临的一大关键性问题，而湖泊作为参与水文循环和气候变化过程的重要组成部分，是全球环境变化的敏感地带和典型区域，它的形成与消失、扩张与收缩均受季节气候模式和长期环境变化影响[1]。青藏高原被誉为“全球气候变化的驱动机与放大器”“亚洲水塔”，分布着地球上海拔最高、数量最多、面积最大的高原湖泊群[2]，因此受到学者们的高度关注。2019年7月在亚洲水塔国际研讨会上，中国科学家发起以青藏高原冰川和湖泊为核心的第三极水塔计划，第三极变化情况已成为国际水资源与生态环境研究热点之一。青藏高原地区海拔较高、人口相对稀少，湖泊的自然条件保持较好，湖面的面积变化受人类活动影响较小[3]，能够较为真实地反映区域气候变化情况，对监测冰川融化[4]影响和反映气候模式变化[5]有着重要作用。纳木错湖是青藏高原第二大湖泊，是一个封闭的湖盆区，空气稀薄，温差较大，生态环境基本为天然状态，人为扰动因素较小，其湖泊的变化主要受自然环境和气候变化的影响[6]。因此，分析和研究纳木错湖面积变化趋势与气候、地形等因素的关系，对反馈高海拔区域生态环境和气候变化有着深远意义[7]。

近年来，许多学者逐渐开展对青藏高原湖泊变化的研究。相较于传统人工测量手段成本高、数据缺乏等问题，利用遥感技术能更全面、更快速高效地获取湖泊动态变化情况[8]，尤其在自然条件恶劣、缺乏地面观测数据和水文观测不连续的青藏高原高海拔地区。然而光学遥感观测受天气和时间的影响较大[9]，因此目前研究主要是利用有限的几景光学遥感影像对高原湖泊提取[10]，缺乏时间上的连续性，在一定程度上难以反映高原湖泊逐年、逐季度甚至是逐月的动态变化情况。合成孔径雷达(synthetic aperture radar，SAR)的优势在于不受云、雨、雾的影响，具有全天候和全天时的探测能力，且能获取不同于光学传感器的遥感地物信息，可以实现对湖泊连续性变化监测研究[11]。目前，基于SAR影像提取水体的方法已有广泛应用，Gong等[12]利用灰度和纹理特征的变化对洪水检测进行研究，取得了较高的精度；庞科臣等[13]利用改进的Otsu法分割DEM，去除山体阴影对水体提取，能有效降低水体提取的虚警率；邓滢等[14]基于结合纹理与极化分解的面向对象方法对极化SAR图像识别水体区域，能够提取完整水陆边界和提高小型水体检测率；王庆等[15]引入纹理特征和不同极化通道之间的极化差、极化比等参数与第一主成分阈值相结合，研究鄱阳湖水体的水域面积动态变化情况；湛南渝等[16]提出将简化的SLIC超像素分割算法运用到SAR影像水体提取，实现基于SAR数据2019年株洲洪水监测。由于高原地区受云雾和海拔高的影响较大，现有研究多为长时间序列的高原湖泊面积变化，缺少针对一年间连续12月湖泊面积变化的研究，而星载SAR能够获取高原多云地区的短周期数据。基于此，本文利用2018年12期Sentinel-1A雷达影像数据，以青藏高原纳木错湖为研究区域，采用面向对象分割方法，通过最大类间方差法分割阈值，并结合灰度共生矩阵提取的纹理特征，提取湖泊边界信息，深入研究2018年纳木错湖泊逐月面积变化情况及驱动因素，为相似气候条件的高海拔湖泊群水体提取和变化监测研究提供依据。

1 研究区域及实验数据

1.1 研究区域

纳木错湖位于西藏自治区中部(30°30′N～30°55′N，90°16′E～91°03′E)，面积约1 920 km2，湖面平均海拔达4 718 m，最大水深超过120 m，是世界上海拔最高的大型湖泊。纳木错东部和南部是终年积雪的念青唐古拉山脉和冈底斯山脉，北部和西部为地形起伏缓和的藏北高原丘陵和湖群，湖区地势略微低于周边区域(图1)。纳木错湖处于水流的聚集中心处，沿湖水系发达，有罗萨、查哈苏太河、打尔古藏布等河流。该地区属于亚寒带季风半干旱气候区，环境寒冷干燥多风，光、热、水资源充足，年日照时数达3 000小时左右，日照强烈，雨、旱季节分明，每年6月至9月为雨季，多年平均降水量为410 mm，11月至翌年5月为旱季，对气候波动较为敏感[17]。

图1 纳木错湖高程分布

1.2 数据及其预处理

1)研究数据源。本文主要利用了Sentinel-1A雷达数据、DEM数据和气象资料数据。Sentinel-1A数据来源于美国地质调查局网站(http：//glovis.usgs.gov/)，选取2018年1月至12月的12期Sentinel-1干涉宽幅模式(IW)的地距(GRD)影像为研究数据源，极化方式为VV/VH，其中GRD数据为经过多视处理、采用WGS-84椭球投影到地距的影像。图2为该影像示例。Sentinel-1A数据及参数见表1。DEM 数据使用SRTMDEM 30 m的高程数据(数据来源于地理空间数据云)。气象数据包括纳木错湖附近的西藏当雄、班戈、那曲和申扎四个气象站2018年的日平均气温、日降水量、日平均风速、日日照时数、日平均蒸发量气象数据，用四个气象站点的气象数据取平均数值表示纳木错湖地区的气象值。

图2 纳木错湖泊2018年Sentinel-1A影像(以1、4、7、11月为例)

表1 Sentinel-1A数据及参数表

2)数据预处理。本文利用SARscape工具对12景Sentinel-1A影像进行预处理。原始的SAR影像数据为幅度数据，需对其进行辐射定标、地形校正、裁剪、滤波和地理编码等预处理，将实验区域的雷达强度图像转换为后向散射系数图像。

2 研究方法

2.1 最大类间方差法

表面近乎光滑的水体在雷达影像中后向散射值小，表现为暗区，灰度阈值分割法依据此特征可以解算得到图像灰度直方图的极值，从而确定水体的分割阈值，将图像分为背景和目标两部分，小于阈值部分的图像归为水体，大于阈值部分的归为背景，形成二值图。最大类间方差法是由日本学者 Otsu提出来的一种全局最优阈值确定方法，也称为Otsu法[18]，是应用较广泛的灰度阈值分割法之一。Otsu法的基本思想是利用图像的灰度值，将原图像分为背景和目标两部分，以两者之间的最大方差确定图像的分割阈值。由于方差是图像像元灰度分布均匀性的一种度量，因此，方差值越大说明图像中目标和背景两部分区域的差别也越大，所对应的分割效果就越好[19]。

2.2 纹理特征提取

在已有的纹理特征提取方法中，灰度共生矩阵(gray-level co-occurrence matrix，GLCM)应用最为广泛，且其提取精度高，因此采用灰度共生矩阵提取水体纹理信息并结合水体几何信息特征，在一定程度上能够减少其他地物对提取水体信息的干扰。灰度共生矩阵是利用窗口内特征值作为该窗口中心像素特征值的统计方法，其通过研究灰度的空间相关特性来描述纹理特征[20]，反映了图像中任意两点灰度的相关性和纹理特征的统计性质。在SAR影像上，灰度共生矩阵反映图像空间结构特征，表现区域内一致性和区域间相对性[21]，其意义是，细纹理灰度空间变化很快，粗纹理随距离增大仅有细微变化。GLCM纹理特征值较多，其中，均值反映图像灰度均匀性和纹理规则程度，水域、道路、建筑物区域的灰度较为均匀，在均值图像上表现为较高亮度，适用于水体提取；均值越小表示相关地物的纹理杂乱无章、难以描述；反之值越大代表纹理规律性强、易于描述。综上，本研究选取均值特征描述纹理信息。

2.3 结合纹理特征与阈值分割法提取湖泊

目前基于SAR数据提取水体的方法主要有灰度阈值分割法、机器学习法、滤波法和结合辅助信息的提取方法等[22]，本文结合灰度共生矩阵提取纹理特征和最大类间方差法阈值分割的面向对象分割方法提取水体。首先，对Sentinel-1A影像数据进行预处理(图3)；其次，通过最大类间方差法确定每个影像的最佳阈值，得到水体和其他背景物二值图；最后，将初水体提取二值图和由灰度共生矩阵提取的纹理特征结合，补充添加的纹理信息，弥补了只采用最大类间方差法产生的相干斑噪声和图像各处灰度信息不均匀的缺陷。另外，青藏高原地区由于地形的原因会在SAR图像上产生地形阴影干扰信息，结合DEM高程辅助数据，以消除由于山体阴影造成的错分混淆现象。

图3 研究技术路线

3 结果与分析

3.1 纳木错湖面积提取结果与评价

1)湖泊面积提取结果。本实验中对预处理后的12景SAR影像通过最大类间方差法确定的最佳阈值分别为：-14.5、-15、-14.7、-16.5、-18.5、-18、-15.8、-16、-16.5、-15.8、-14.8、-14，根据这些阈值分别得到湖泊提取的二值图。由图4(a)中Otsu法提取湖泊结果发现，水体和陆地未能区分，湖泊的轮廓不清晰，存在虚警率较大的问题。值得注意的是，阈值分割方法简单、运算速度较快，能进行水体的快速提取，但此方法缺乏空间特征考量，在区分与水体散射特征相似的地物时较容易发生误分的情况，精度和鲁棒性较差，同时由于受相干斑噪声和图像各处灰度信息不均匀现象的影响，难以获取准确的阈值，尤其是高海拔地区湖泊表面有湖冰存在，只通过简单的阈值法无法完全正确地提取湖泊的整体边界，提取的水体边缘线不够光滑连续，水体和河岸边缘线也不能明确分离，存在边缘误差特征。

随着SAR影像的空间分辨率不断提高，图像亮度范围随之增大，纹理结构信息也更加丰富。本文在提取水体时引入纹理特征[23]，增强图像的空间信息，抑制SAR图像的斑点噪声。将纹理灰度共生矩阵与图像中的灰度信息相结合，使提取结果的“椒盐现象”明显减少。从图4(b)湖泊提取结果发现，基于GLCM-均值方法水体提取结果边界清晰但仍存在误提现象，尤其在影像的右上侧部分明显，未能提取完整湖体边界。通过最大类间方差法确定水体分割初阈值，并结合灰度共生矩阵提取纹理信息，弥补只通过阈值法在处理大面积影像时获取的结果主观性强、可追踪性差的缺点，湖泊提取结果如图4(c)所示，该方法能更好地区分出水域和其他地形，提取边界清晰明了。

图4 不同方法提取结果

2)精度分析与评价。

(1)不同水体提取方法精度分析。为了验证方法的有效性和精度，采用地理国情数据人工标注的解译水体范围作为真实水体值，将本文方法、Otsu阈值法和基于纹理特征提取的水体范围与真实水体提取值对比分析。引用误差统计方法进行精度评估，准确率F较大，说明提取方法效果较好；准确率F较小代表提取方法效果较差。

如表2所示，以纳木错湖为研究区域，Otsu方法提取湖泊面积准确率最低，仅有54.73%，这是由于Otsu法未能将水体和陆地区分开，存在大量误提现象；基于纹理特征的提取精度88.26%，基本能提取湖泊的完整边界，但噪声影响依旧存在；结合纹理特征与Otsu方法提取精度明显提高，达到95.12%，能较好地获取湖泊提取结果。

表2 不同水体提取方法精度对比

(2)不同季节的水体提取精度分析。根据SAR影像的辐射特性，水体的后向散射系数由波长、斜视角、复介电常数、粗糙度等雷达参数决定。在水体提取分析中，Sentinel-1A SAR影像主要考虑地表粗糙度对后向散射系数的影响，表面越粗糙，后射回波信号越强，后向散射系数越大；平静的水面属于平滑面，主要表现为镜面反射，其后向散射强度低。图5为纳木错不同月份的水体后向散射系数变化图，各阶段后向散射强度层次明显，夏季风速较小且湖面处于完全消融状态，近似为镜面反射，其后向散射强度值明显低于其他季节；冬季后向散射系数值最大，由于温度降低湖面开始结冰，且冬季风速较大，导致纳木错冬季时湖面呈粗糙形态，后向散射回波信号增强。基于此，以1、4、7、10月为例，分析不同季节、不同湖面形态的水体提取精度，准确率分别为：80.15%、83.62%、89.37%、85.26%。夏季7月份时水体提取精度最高，湖面状态在一年中最为平静，提取边界清晰；秋季10月份精度次之，温度逐渐降低至零下，湖面开始结冰，表面形态发生变化，影响水体提取精度；春季4月份时精度由于受风速的影响不是很高；冬季1月份水体提取精度最低，湖面既是结冰状态又受大风影响，导致水体边界有所误提。

图5 纳木错水体后向散射系数时变特征

3.2 纳木错湖湖泊面积变化监测分析

基于本文提取水体方法得到2018年12个月份纳木错湖的空间分布范围，分析了纳木错湖在不同月份湖面面积的变化情况，具体如下。

由图6发现，2018年12个月内纳木错湖的面积变化呈现起伏状态，主要经历三个阶段：从3—6月迅速增加，到6—10月平稳变化，最终10月至次年2月收缩减小。从数值上来看，9月份纳木错湖面积达到峰值，为2 020.639 km2，3月份时湖泊面积最小，只有1 980.179 km2，二者相差40.46 km2。由于纳木错地处藏北高原，每年冰封期长达5个月(完全封冻时间近3个月)，结合影像数据可直观发现1—3月份湖泊大部分处于冰冻状态，直到4月份开始消融，此时湖水面积开始明显增加；夏季由于温度上升导致冰川融水和湖冰自身消融，加之降水量增加，使得9—10月份湖泊的面积达到峰值；从11月进入冬季，温度降低降水减少，湖泊面积发生收缩变化。

图6 2018年纳木错湖面积统计图

通过图7纳木错湖12个月份提取的湖面矢量叠加图可以观察得到，纳木错湖泊周围都存在边界扩张的现象，其中有四处明显变化区域，分别在湖区的东北部、西北部和西南部，主要是由南部的冈底斯山脉和东部的念青唐古拉山脉冰雪消融补给，以及周围的罗萨、查哈苏太河、打尔古藏布等河流注入引起变化。

图7 2018年纳木错湖空间分布及面积变化

3.3 湖面面积变化驱动因素分析

1)气象资料分析。结合2018年纳木错湖周围的当雄、班戈、那曲和申扎四个气象站点的日均气温、日均降水量、日均风速、日均蒸发量、日均日照时数等气象数据(图8)，分析了纳木错湖区2018年的气候变化。从图8气象数据和其变化趋势线能够看出，2018年期间气温的变化起伏最为明显，1—4月气温均处于0 ℃以下，5—10月温度均处于0 ℃以上，6—9月日均温度都可达到10 ℃，11月起气温再次降到0 ℃以下，且温度与湖泊的面积的变化走向大致相似。2018年月降水的变化趋势也较为明显，1—3月几乎无降水，5—9月的降水明显增加进入雨季，其中7—8月降水最为丰富达到日均5 mm，10—12月再次进入枯水期几乎无降水。纳木错地区风资源丰富，2018年期间大风日主要集中在12月至翌年5月，大风日占全年的一半时间以上，最大日均风速可达到7 m/s。夏季蒸发量最大，冬季蒸发量最小，蒸发量整体波动不大。湖区及附近区域日照强烈，2018年日照时数接近3 000小时，年均日照率大于65%。

图8 2018年纳木错区域气象数据

2)湖面面积变化与气象资料相关性分析。通过对两个或多个具有相关性的变量进行分析，能够得到两个变量之间的线性关系程度，即Pearson相关性分析。纳木错湖泊水位的变化直接反映为湖水表面积的变化，而影响纳木错湖水位变化的因素有各支流、雪山、气温和降雨等原因。通过将统计得到的纳木错湖湖面面积统计值与各气象数据之间进行相关性分析，由面积-气象数据相关性分析表3的结果分析得出：纳木错湖面积变化与气温(相关系数r=0.680，p<0.05)和降水(相关系数r=0.638，p<0.05)存在较强的正相关性，与平均风速存在强的负相关(相关系数r=-0.796，p<0.01)；湖面面积变化与蒸发量的相关系数为r=-0.300，p=0.344(p>0.05)，与日照时数的相关系数为r=-0.171，p=0.596(p>0.05)，因此湖面面积变化与蒸发量和日照时数变化的相关性不显著。

表3 面积-气象数据相关性分析表

由气象数据和相关性分析结果可得到：纳木错湖区雨、旱季节分明，热季和雨季均为每年的6—9月，湖水水位在9—10月达到年内稳定最大值；多风是纳木错流域及其附近地区气候的显著特点，风速增大加快湖泊水面蒸散发。随着气温的不断升高，可能会使湖泊补给流域的上游冰川加速消融和湖冰自身融化，温度升高是纳木错湖泊面积增加的主要原因，其中，纳木错湖东部和南部的念青唐古拉山脉、冈底斯山脉雪冰融化为湖泊增加了水域补给；雨季降水量增加也是影响湖面积变化的重要因素，风速对湖泊面积变化起负相关作用，蒸发量和日照时数两个变量对纳木错湖面积变化影响较小。

3)地形和相对水深分析。对纳木错湖及周围区域地形分布进行分析发现：湖区西北侧是人类居住活动区域，地势相对湖泊较高；东南侧是念青唐古拉山脉，四周多雪山，湖泊位于地势最为低洼的区域，是水流的汇集中心，且不易流散。叠加DEM高程图中可以观察到，纳木错湖地处区域地势相对较低，因此能够推断湖水水位增高部分会掩盖地势相对较低地区，水位继续上涨，纳木错湖仍会向地势相对较低的区域扩张。

利用2018年3月Landsat 8 OLI遥感影像数据反演得到整个纳木错湖泊的相对水深，从相对水深反演结果(图9)发现：相对水深最小的位置是陆地，值为0为图中白色区域，值越大表示相对水深越大，即在图中呈现的颜色越深；对纳木错湖水深反演结果进行不同的颜色分级显示，可以清晰明了地看出整个纳木错湖的面积变化过程，即从深红色变到橘色再到淡黄色最后到绿色，边界较为清晰；纳木错湖泊向水深相对较浅的区域扩张消融，主要的扩张边界在西侧和东北侧方向；相对水深反演结果与2018年纳木错湖面积变化过程(图7)基本一致。

图9 纳木错湖相对水深反演

4)湖泊景观形状指数变化分析。景观形状指数用来描述自然景观受外界条件包含人为因素干扰的难易程度，也能反映景观的复杂程度和易损程度，此值越小代表湖泊景观越容易收到外界干扰，反之亦然。其中，湖泊景观形状指数(lake landscape shape index，LLSI)的计算如式(1)所示。

(1)

式中：LLSI为湖泊景观形状指数；L为湖泊的周长；A为湖泊的面积。

由图10得到，1、2、3、4、12月LLSI值较小，由于期间温度降低湖面结冰和受大风影响，湖泊易受到外界干扰，LLSI值均低于2.6；5—11月LLSI值较大，是由于期间湖泊的气候环境和湖面状态相对稳定，受外界干扰影响降低。

图10 2018年纳木错区域LLSI变化

4 结束语

本文基于Sentinel-1A合成孔径雷达数据，以纳木错湖为研究区域，在Otsu阈值法对影像粗分割提取的基础上，利用灰度共生矩阵分析影像纹理特征，并结合DEM辅助数据消除山体阴影造成的错分混淆现象，对粗提取结果进一步优化整合，实现纳木错湖泊的准确提取，最终获得湖泊2018年逐月面积变化情况。结合周围四个气象站点的气象数据，分析湖区月际变化趋势及原因，结论如下。

1)通过实验对比分析得到，Otsu法提取精度为54.73%，处理速度快但精度较低；基于纹理特征方法提取精度为88.26%，可以明显获取湖泊边界，但存在边界不清晰且误提的现象；结合纹理特征与阈值分割的方法提取精度达到95.12%，效果得到明显提升。实验结果证明，SAR影像结合本文方法可用于高海拔地区湖泊水体信息的提取，并能实现逐年、逐月的常态化动态变化监测。

2)2018年内，纳木错湖面积呈现起伏变化，3至6月面积增加，6至10月平稳变化，10月至次年2月收缩减小；在9至10月份时面积达到了峰值2 020.639 km2，3月份时湖泊面积最小，只有1 980.179 km2。通过湖泊面积变化与气象数据的相关性分析，气温升高和降水增加等因素是纳木错湖泊面积变化的主要正相关因素，风速是影响湖面面积变化的主要负相关因素，蒸发量和日照时数变量对湖泊逐月面积变化结果影响较小。另外，对整个湖区进行相对水深反演和对湖区周围地形分析，反演结果与面积变化过程很好地对应，在西北、西南和东北方向有明显扩张。

综上所述，利用合成孔径雷达SAR数据对纳木错湖观测，可以有效地获取其水体短期内的连续变化信息，掌握其变化规律，分析影响湖面面积变化的因素。本文仍存在不足之处，目前只针对纳木错湖2018年逐月数据进行了水体提取，今后需扩大时间跨度，研究年度、季度甚至月度的湖泊面积变化情况，并进一步加强对相似气候条件下湖泊的研究及对比分析。