本发明涉及基于纹理特征辅助的SAR影像冰川识别方法，利用Sentinel‑1A数据的干涉测量方法，结合纹理信息识别冰川。其步骤包括利用Sentinel‑1A数据进行干涉处理得到相干系数图。提取相干系数图的均值、方差、同质性、反差、相异性、熵、相关性7个纹理特征。对不同纹理特征组合之间的冰川提取效果进行比较。选择分类精度最高的VV极化方式下均值、方差、同质性、相异性四个纹理特征，利用SVM分类器进行训练，得到冰川、岩石及积雪等类别。将分类后的数据转换成shapefile格式，通过山脊线进行分割，获得冰川边界数据集。与传统的基于相干系数图阈值分割法相比，本发明引入纹理特征，冰川识别更准确。
1.基于纹理特征辅助的SAR影像冰川识别方法，包括以下步骤:第一步、准备遥感影像数据和数字高程模型数据；
步骤1a、下载某地区两景的SAR影像，两者具有相同的轨道号和图幅号；
步骤1b、下载与SAR影像同地区同时期的光学卫星影像，且冰川上方无明显云层覆盖；
第二步、将两景SAR影像进行干涉处理，包括：步骤2a、影像配准，使两景SAR影像坐标达到一致；
步骤2b、干涉图计算，通过像元的复数共轭，计算两幅SAR影像像元的相干系数；
步骤2c、去除平地效应和相位解缠，使干涉相位在2π的范围内，生成解缠后的相干系数图；
步骤2d、对解缠后的相干系数图进行地理编码，将相干系数图的坐标转化为地理坐标，以对应地面真实位置；
第三步、对地理编码后的相干系数图生成纹理特征图像，计算水平方向上相隔1个像元M距离的一对像素灰度值分别为i和j出现的频次，进而生成灰度共生矩阵，i和j是0到2之间的整数，M取值为3,4,5或6，基于灰度共生矩阵提取均值、方差、同质性、反差、相异性、熵、相关性7个纹理特征图像；
第四步、基于纹理特征进行分类器的训练，包括以下步骤：
3a、对光学卫星影像目视解译，选取若干冰川和非冰川；
3b、选取的冰川和非冰川按比例划分为训练组和验证组；
3c、分别使用训练组、验证组中的冰川多边形和非冰川多边形裁剪所述纹理特征图像，得到基于纹理特征的训练样本和验证样本；
3d、选择不同的纹理特征组合所对应的训练样本进行支持向量机分类器训练，使用训练后的模型对该纹理特征组合所对应的验证样本进行分类，并记录每个支持向量机分类器的分类精度；
3e、分类精度最高的支持向量机分类器为最优分类器，对应的纹理特征组合为最佳纹理特征组合；
第五步、使用所述最优分类器对最佳纹理特征组合所对应的纹理特征图像进行分类，得到分类结果，完成冰川识别。
2.根据权利要求1所述基于纹理特征辅助的SAR影像冰川识别方法，其特征在于：所述两景SAR影像的时间间隔不超过36天。
3.根据权利要求1所述基于纹理特征辅助的SAR影像冰川识别方法，其特征在于：所述SAR影像为Sentinel‑1A影像，两景Sentinel‑1A影像的时间间隔为12天，所述光学卫星影像为Landsat影像，所述DEM影像为SRTM DEM影像。
4.根据权利要求1所述基于纹理特征辅助的SAR影像冰川识别方法，其特征在于：第三步中，均值、方差、同质性、反差、相异性、熵、相关性的计算公式如下：(1)均值mean：代表图像亮度
(2)方差variance：代表图像像元值与均值之间的偏差2
variance＝∑∑(f(i,j)‑mean)(3)同质性homogeneity：代表图像灰度的均匀性(4)反差contrast和相异性dissimilarity：代表图像灰度的对比性
contrast＝∑∑(i‑j) f(i,j)dissimilarity＝∑∑|i‑j|f(i,j)(5)熵entropy：代表图像平均信息量entropy＝‑∑∑f(i,j)log[f(i,j)](6)相关性correlation：代表不同相对位置上相似性程度和相关匹配程度2
其中，i，j为像素灰度值，f(i,j)为像素灰度值分别为i和j出现的概率，n为窗口大小，σ2
5.根据权利要求1所述基于纹理特征辅助的SAR影像冰川识别方法，其特征在于：第三步完成后，对7个纹理特征图像合成一景多特征影像，然后进行第四步的步骤3c中，使用多边形裁剪多特征影像完成对所有纹理特征图像的裁剪。
6.根据权利要求1所述基于纹理特征辅助的SAR影像冰川识别方法，其特征在于：所述第四步的步骤3b中，选取的冰川和非冰川按照7:3的比例随机划分为训练组和验证组。
7.根据权利要求1所述基于纹理特征辅助的SAR影像冰川识别方法，其特征在于：第五步分类完成后，进行分类精度评价，若不满足分析要求，则转至第四步。
8.根据权利要求1所述基于纹理特征辅助的SAR影像冰川识别方法，其特征在于：还具有第六步、对于分类结果填补空洞并且手工修正明显的错分像元，改进分类结果。
9.根据权利要求8所述基于纹理特征辅助的SAR影像冰川识别方法，其特征在于：第六步完成后，基于坡向变率法从DEM数据中提取山脊线，使用山脊线对识别结果进行分割，得到冰川边界数据集。
基于纹理特征辅助的SAR影像冰川识别方法
技术领域
[0001] 本发明涉及基于纹理特征辅助的SAR影像冰川识别方法，属于遥感地学应用技术领域。
背景技术
[0002] 冰川资源是地球上重要的淡水资源，作为全球气候变化的重要指示器，利用遥感方法进行冰川识别对于监测冰川面积及其储量，并研究全球气候变化具有重要意义。但是由于表碛覆盖型冰川表面有一层岩屑混合物，与冰川周围的岩石在光学影像上具有相似的光谱信息，波段比值法及归一化雪盖指数等光学遥感识别方法只能将识别较干净的纯净冰川，且光学遥感受到天气条件的影响较大，难以对同一地区进行同时相的观测。SAR能够全天时、全天候的监测地物，且由于冰川的运动造成其在干涉测量过程中相对其他地物具有较低相干性。同时，由于冰川表面失相干现象的存在，在相干系数图上冰川区域没有形成周期性变化的条纹，从而具有不同于其他地物的纹理特征，可以利用这一特征识别冰川。
发明内容
[0003] 本发明要解决的技术问题是：提出基于纹理特征辅助的SAR影像冰川识别方法。考虑到光学遥感数据只能将识别较干净的纯净冰川，而SAR影像的相干系数阈值分割法会在冰川末端产生很多大小不一的空洞，需要较多的人工修补。
[0004] 为了解决上述技术问题，本发明提出的技术方案是：基于纹理特征辅助的SAR影像冰川识别方法，包括以下步骤:
[0005] 第一步、准备遥感影像数据和数字高程模型数据；
[0006] 步骤1a、下载某地区两景的SAR影像，两者具有相同的轨道号和图幅号；
[0007] 步骤1b、下载与SAR影像同地区同时期的光学卫星影像，且冰川上方无明显云层覆盖；
[0008] 步骤1c、下载对应地区DEM数据；
[0009] 第二步、将两景SAR影像进行干涉处理，包括：
[0010] 步骤2a、影像配准，使两景SAR影像坐标达到一致；
[0011] 步骤2b、干涉图计算，通过像元的复数共轭，计算两幅SAR影像像元的相干系数；
[0012] 步骤2c、去除平地效应和相位解缠，使干涉相位在2π的范围内，生成解缠后的相干系数图；
[0013] 步骤2d、对解缠后的相干系数图进行地理编码，将相干系数图的坐标转化为地理坐标，以对应地面真实位置；
[0014] 第三步、对地理编码后的相干系数图生成纹理特征图像，计算水平方向上相隔1个M像元距离的一对像素灰度值分别为i和j出现的频次，进而生成灰度共生矩阵，i和j是0到2之间的整数，M取值为3,4,5或6，基于灰度共生矩阵提取均值、方差、同质性、反差、相异性、熵、相关性7个纹理特征图像；
[0015] 第四步、基于纹理特征进行分类器的训练，包括以下步骤：
[0016] 3a、对光学卫星影像目视解译，选取若干冰川和非冰川；
[0017] 3b、选取的冰川和非冰川按比例划分为训练组和验证组；
[0018] 3c、分别使用训练组、验证组中的冰川多边形和非冰川多边形裁剪所述纹理特征图像，得到基于纹理特征的训练样本和验证样本；
[0019] 3d、选择不同的纹理特征组合所对应的训练样本进行支持向量机分类器训练，使用训练后的模型对该纹理特征组合所对应的验证样本进行分类，并记录每个支持向量机分类器的分类精度；
[0020] 3e、分类精度最高的支持向量机分类器为最优分类器，对应的纹理特征组合为最佳纹理特征组合；
[0021] 第五步、使用所述最优分类器对最佳纹理特征组合所对应的纹理特征图像进行分类，得到分类结果，完成冰川识别。
[0022] 本方法利雷达数据，通过相干系数图上的纹理特征对冰川进行识别，首先生成两景影像的相干系数图，以衡量像元的相位差，利用DEM进行地理编码。在相干系数图上生成灰度共生矩阵，提取7个纹理特征，并对比了不同纹理特征组合之间的提取精度，最后利用精度最高的特征组合对影像进行SVM分类，得到冰川识别结果。通过坡向变率提取山脊线的方法生成山脊线，分割冰川识别结果，可得到冰川边界数据集。与传统的基于相干系数图阈值分割法方法相比，本发明引入纹理信息，在冰川内部精度有了提升，冰川识别更为准确。
[0023] 本发明采用Sentinel‑1A、Landsat以及SRTM DEM数据，数据获取相对简单。Sentinel‑1A通过干涉测量得到相干系数图，结合灰度共生矩阵提取图像的纹理信息，使分类的特征空间中加入了纹理信息，实现相干系数信息在不同维度上的扩展。通过和Landsat影像的比对，提供冰川样本，构建SVM分类器进行分类，提取冰川边界。
[0024] 综上，本发明的执行步骤简单易行，分类效果较好。目前基于传统的光学影像的分类方法受制于青藏高原地区复杂的天气条件及难以提取表碛覆盖型冰川。相干系数阈值分割法需要大量人工参与，分类效率较差。本发明综合利用SAR相干系数图的均值、方差、同质性、相异性特征，提出识别冰川的方法。该方法只使用Sentinel‑1A、Landsat数据和SRTM DEM数据，以较高精度完成冰川的识别，为研究冰川面积及储量的变化提供了新的思路，对研究全球气候的影响有深远的科学意义。
附图说明
[0025] 图1是本发明基于纹理特征辅助的SAR影像冰川识别方法流程图。
[0026] 图2是相干系数图，(a)为VH极化下的相干系数图，(b)为VV极化下的相干系数图。
[0027] 图3是VV极化方式提取的7个纹理特征示意图，(a)为均值的纹理特征图像，(b)为方差的纹理特征图像，(c)为同质性的纹理特征图像，(d)为对比性的纹理特征图像，(e)为反差的纹理特征图像，(f)为熵的纹理特征图像、(g)为相关性的纹理特征图像。
[0028] 图4是初步分类结果图。
[0029] 图5是经手工修改后的最终分类结果图。
具体实施方式
[0030] 下面结合附图对本发明作进一步的说明。使本发明的技术路线和操作步骤更加清晰。本发明实例研究区选取喀喇昆仑山的克勒青河上游(35.1°N‑36.2°N，76.3°E‑77.9°E)区域，采用的Sentinel‑1A数据是2018年9月11日和9月23日的两景SLC数据，时间基线12天。Landsat8 OLI数据获取时间是2018年8月29日。DEM数据是SRTM Version4数据。
[0031] 如图1所示，为本实施例基于纹理特征辅助的SAR影像冰川识别方法流程图，具体步骤如下：
[0032] 第一步、准备遥感影像数据和数字高程模型(DEM)数据。
[0033] 步骤1a、下载某地区两景时间间隔为12天的Sentinel‑1A SLC数据，两者具有相同的轨道号和图幅号；
[0034] 步骤1b、下载与Sentinel‑1A数据同地区同时期的Landsat影像，且冰川上方无明显云层覆盖；
[0035] 步骤1c、下载对应地区SRTM DEM数据。
[0036] 第二步、将两景Sentinel‑1A影像SLC格式数据进行干涉处理，包括：
[0037] 步骤2a、影像配准，使两景图像坐标达到一致；
[0038] 步骤2b、干涉图计算，通过像元的复数共轭，计算两幅SAR影像像元的相干系数；
[0039] 步骤2c、去除平地效应和相位解缠，使干涉相位在2π的范围内，生成解缠后的相干系数图；
[0040] 步骤2d、对解缠后的相干系数图进行地理编码，将SAR图像坐标转化为地理坐标，以对应地面的真实位置。
[0041] Sentinel‑1A数据有两种极化方式成像模式，分别是VV极化和VH极化，本步骤中，将两种极化方式成像的图像分别进行干涉处理，以比较不同极化方式对冰川提取效果的影响，结果如图2所示。由于同极化(VV)比交叉极化(VH)渗透能力强，VV极化的相干性比VH极化的相干性高，在图像上呈现更高的亮度。本步骤选择SRTM DEM数据进行地理编码。原因有以下几点，第一、该数据覆盖范围广；第二、该数据获取全部免费，使用方便，具有一定的权威性；第三、该数据质量较好，错误较少，使用可靠。
[0042] 第三步、对地理编码后的相干系数图生成纹理特征图像。计算水平方向上相隔1个M像元距离的一对像素灰度值分别为i和j出现的频次，进而生成灰度共生矩阵，i和j是0到2之间的整数，本例中M取值为3，也可以取4、5或6，一般不超过6。基于灰度共生矩阵提取均值、方差、同质性、反差、相异性、熵、相关性7个纹理特征图像，如图3所示。并对7个纹理特征图像合成一景多特征影像(类似假彩色合成)，作为用于分类的影像；均值(a)值越小代表亮度越大，冰川的均值较小，主要分布于0‑0.36之间。方差(b)值越大代表着偏差越大，冰川的方差相对于其它地物较大，主要分布于0.0049‑0.1380。同质性(c)值大代表着图像纹理之间缺少变化，冰川的同质性集中分布于0.98‑1。对比性(d)和反差(e)值越大代表着图像差异性越大，冰川与其它地物越具有区分性。熵(f)值越大代表着信息量越丰富，冰川熵集中分布在1.887‑2.197。相关性(g)值均匀相等时，相关值就大，冰川的相关性在0.051‑0.624之间均匀分布。
[0043] 本步骤中需要提取灰度共生矩阵和纹理特征，本发明采用R.Haralick等在Journal of Signal and Information Processing发表的《Texture Filters and Fractal Dimension on Image Segmentation》一文中提出的灰度共生矩阵和特征，各特征计算式如下：
[0044] (1)均值mean：代表图像亮度
[0045]
[0046] (2)方差variance：代表图像像元值与均值之间的偏差
[0047] variance＝∑∑(f(i,j)‑mean)2
[0048] (3)同质性homogeneity：代表图像灰度的均匀性
[0049]
[0050] (4)反差contrast和相异性dissimilarity：代表图像灰度的对比性
[0051] contrast＝∑∑(i‑j)2f(i,j)2
[0052] dissimilarity＝∑∑|i‑j|f(i,j)
[0053] (5)熵entropy：代表图像平均信息量
[0054] entropy＝‑∑∑f(i,j)log[f(i,j)]
[0055] (6)相关性correlation：代表不同相对位置上相似性程度和相关匹配程度[0056]
[0057] 其中，i，j为像素灰度值，f(i,j)为像素灰度值分别为i和j出现的概率，n为窗口大2 2
小，σ＝∑∑(i‑j) f(i,j)。
[0058] 第四步、基于纹理特征进行分类器的训练，包括以下步骤：
[0059] 3a、对Landsat影像目视解译，选取若干冰川和非冰川；
[0060] 3b、选取的冰川和非冰川按照7:3的比例随机划分为训练组和验证组；
[0061] 3c、分别使用训练组、验证组中的冰川多边形和非冰川多边形裁剪所述纹理特征图像，得到基于纹理特征的训练样本和验证样本。为了方便可对多特征影像裁剪，即可实现一次完成对所有7个纹理特征图像的裁剪；
[0062] 3d、选择不同的纹理特征组合所对应的训练样本进行支持向量机(SVM)分类器训练，使用训练后的模型对该纹理特征组合所对应的验证样本进行分类，以检验样本中被正确识别的比例作为该分类器的精度；并记录每个支持向量机分类器的分类精度；
[0063] 3e、比较不同特征组合的分类精度，分类精度最高的支持向量机分类器为最优分类器，对应的纹理特征组合为最佳纹理特征组合，本实施例中，最佳纹理特征组合为：均值、方差、同质性和相异性。
[0064] 本步骤采用支持向量机(SVM)方法进行图像的分类。SVM建立在统计学习理论上，基本思想是通过一个非线性映射，将数据x映射到高维特征空间F，并在这个空间进行分类，获得最优分类面，即分类面不仅两类正确分开(训练错误率为0)，且分类间隔最大。
[0065] 第五步、使用所述最优分类器对最佳纹理特征组合所对应的纹理特征图像进行分类，得到分类结果，如图4所示。
[0066] 第六步、对于分类结果填补空洞并且手工修正明显的错分像元，改进分类结果，完成冰川识别，如图5所示。
[0067] 本步骤中考虑到即使分类精度很高，还是会存在部分错分像元，因此需要手动消除，手工消除的错分像元仅仅是可以判断并非常明显的、严重影响分析结果的像元，因此手工的工作量不大。
[0068] 第七步、通过坡向变率(SOA)法提取山脊线，使用山脊线对识别结果进行分割，获得冰川边界数据集。
[0069] 除上述实施例外，本发明还可以有其他实施方式。凡采用等同替换或等效变换形成的技术方案，均落在本发明要求的保护范围。
