CN103440489A - 基于像素级sar影像时间序列相似性分析的水体提取方法 - Google Patents

Abstract

本发明涉及一种像素级SAR影像时间序列的水体提取方法，该方法首先，对SAR影像数据集进行预处理，经过高精度匹配，构建像素级SAR影像时间序列，生成时间序列文本数据；其次，采样选取纯净水体像元和混合水体像元的时间序列，选取DTW作为时间序列的相似性度量，计算其DTW值作为最大阈值；然后计算所有像元的像素级SAR影像时间序列与纯净水体像元时间序列的DTW值，采用最大阈值方法分割SAR影像，获取二值图像；最后，采用8邻域搜索方法对二值图像进行操作以提高水体识别精度。该方法能够准确提取稳定的水资源分布范围，提取结果不受山体阴影、雨季积水及部分植被的影响，能够满足水体制图的要求。

Description

基于像素级SAR影像时间序列相似性分析的水体提取方法

技术领域

本发明涉及一种SAR影像的水体提取方法，特别是涉及一种利用像素级SAR影像时间序列相似性分析的水体提取方法。

背景技术

水资源是自然生态环境资源的重要组成部分，水资源的分布和供需状况已经成为评价一个国家或者地区综合竞争力的重要指标之一。遥感技术在水体分布提取方面得到了广泛的应用，如MODIS、TM、SPOT等。在微波遥感领域，合成孔径雷达（SAR）能够不受日照和天气条件的限制完成对地观测任务，并对某些地物具有一定的穿透能力，具有多极化、多入射角、多模式、短重访周期等特点，在水资源提取与监测中具有独特的优势。然而，由于SAR斜距成像的特点，山体阴影、雨季积水及部分植被严重影响了水体的提取精度。针对SAR影像水体目标，学者们提出了很多提取方法，比如2006年，朱俊杰等在《水科学进展》17卷第4期中发表“基于纹理与成像知识的高分辨率SAR图像水体检测”一文，综合利用阈值法、小波纹理法和建筑物的SAR成像知识对水体进行了检测，实现了高分辨率SAR图像的水体检测；2008年，胡德勇等在《中国图象图形学报》13卷第2期中发表“单波段单极化SAR图像水体和居民地信息提取方法研究”一文，针对单波段单极化SAR影像，基于纹理分析方法提取水体及居民地目标，并利用支持向量机进行分类；2009年，程明跃等在《计算机工程》35卷第2期中发表“基于模糊加权SVM的SAR图像水体自动检测”一文，使用树型小波方法对SAR图像进行纹理分析，有效地利用SAR图像中的纹理信息，结合图像中的灰度信息获得判别水体的决策函数，通过该决策函数能够检测出图像的水体区域；2010年，李景刚等在《自然灾害学报》19卷第3期中发表“ENVISAT卫星先进合成孔径雷达数据水体提取研究—改进的最大类间方差阈值法”一文，提出一种综合考虑类间和类内方差两个因素的改进的最大类间方差法，对ENVISAT ASAR数据进行水体提取；2011年，王栋等在《计算机工程与应用》47卷第2期中发表“盲信号分离和序贯滤波的SAR影像水体提取”一文，针对SAR影像水体目标的亮度及形状分布特征，进一步采用序列非线性滤波处理，从分离出的图像分量中提取出水体目标。

然而这些方法都存在以下缺点：（1）这些方法大多是针对单幅影像提取水体，无法去除雨季积水及部分植被（如水稻）的影响；（2）多时相的水体提取方法虽能够在一定程度上克服以上影响，但还不够彻底；（3）制图要求中明确规定，湖泊、水库等水体制图应以水涯线为准，即一年中大部分时间平稳的常水位，而以往的方法只能提取瞬时影像中的水体分布，无法提取常年稳定的水体面积。因此这些方法无法准确提取稳定的水资源分布范围，同时受山体阴影、雨季积水及部分植被的影响，提取结果尚不能满足水体制图的要求，无法应用于大范围区域水体分布的提取与制图中。

发明内容

本发明解决的技术问题是：提出一种精度高的基于像素级SAR影像时间序列相似性分析的水体提取方法，该方法能够准确提取稳定的水资源分布范围，提取结果不受山体阴影、雨季积水及部分植被的影响，能够满足水体制图的要求，从而可以应用于大范围区域水体分布的提取与制图中。

为了解决上述技术问题，本发明提出的技术方案是：一种基于像素级SAR影像时间序列相似性分析的水体提取方法，包括以下步骤:

第一步、对SAR影像数据集进行预处理，所述SAR影像数据集包含有N幅不同时相的SAR影像；

第二步、对SAR影像数据集的N幅SAR影像进行几何配准，得到配准后的SAR影像数据集；

第三步、针对坐标值相同的SAR影像像元的后向散射系数，根据时间顺序建立与坐标值相关的像元后向散射系数时间序列；

第四步、在SAR影像数据集的图幅范围内目视选取M个纯净水体像元，相应的获得M个纯净水体坐标，对每幅SAR图像中M个纯净水体坐标位置的像元的后向散射系数求均值，得到SAR图像集的N个纯净水体后向散射系数均值，根据时间顺序建立纯净水体后向散射系数均值时间序列；

第五步、在SAR影像数据集的图幅范围内目视选取K个水体边缘像元，相应的获得K个水体边缘坐标，计算所述纯净水体后向散射系数均值时间序列与每个水体边缘坐标对应的所述像元后向散射系数时间序列之间的DTW距离值，并将计算得到的K个DTW距离值作平均，得到用于划分水体像元的最大阈值；

第六步、构建与SAR影像数据集的图幅范围相同的空白影像，第三步得到的像元后向散射系数时间序列与纯净水体后向散射系数均值时间序列之间的DTW距离值小于或等于所述最大阈值，则空白影像中相应坐标值的像元为水体像元；否则为非水体像元。

为了更精确的提取水体像元，在第六步完成之后可以去除水体内的噪点。

作为优选的方案，本发明提供了一种去除水体内噪点的方法如下：

1）将第六步得到的影像中的水体像元设为1，非水体像元设为0，从而将SAR影像赋值为0-1二值图像，假设二值图像为m行*n列，将二值影像四周扩充至(m+1)行*（n+1）列，四周像元赋值为0，得到新的二值图像；

2）选取3*3滑动窗口，从新的二值图像左上角像元（1,1）开始，依次判断滑动窗口中心点（i,j）的值，若该中心点（i，j）的八邻域像素（i-1，j-1）、(i-1,j)、(i-1,j)、(i,j-1)、(i,j+1)、（i+1，j-1）、（i+1，j）、(i+1,j+1)均为1，则将中心点像元（i,j）赋值为1；

3）将所述滑动窗口按从左到右、从上到下顺序移动到下一个像元，依次重复执行步骤2），直到全幅影像扫描一遍结束。

另外，本发明第一步中对SAR影像数据集的预处理包括影像的裁剪和拼接、元数据更新、辐射校正、几何校正和地形校正。

作为优选的方案，第二步中进行几何配准的方法为：选取SAR影像数据集中任一幅SAR影像作为主影像，其他影像为从影像进行几何配准。

作为优选的方案，本发明中SAR影像中第i行第j列像元的后向散射系数

${\mathrm{\σ}}_{\mathrm{ij}}^0=\frac{{{\mathrm{DN}}_{\mathrm{ij}}}^2}{K}\sin \left({\mathrm{\α}}_{\mathrm{ij}}\right)$

式中，DN_ij为第i行第j列像元的辐射强度值；α_ij为第i行第j列像元的的入射角；K为绝对定标因子。

本发明带来的有益效果是：本发明利用SAR影像数据，提出了一种基于像素级SAR影像时间序列的水体提取方法。本方法的核心在于充分利用像素级SAR影像时间序列的相似性，即采用动态时间弯曲距离（DTW）作为相似性的度量，利用最大阈值法提取水体像元。该方法能够准确提取稳定的水资源分布范围，提取结果不受山体阴影、雨季积水及部分植被的影响，能够应用于大范围区域水体分布的提取与制图中。该方法提取的水体分布精确度很高，能够满足水体制图的要求。

附图说明

下面结合附图对本发明的基于像素级SAR影像时间序列相似性分析的水体提取方法作进一步说明。

图1是本发明实施例的研究区的原始SAR影像数据集。

图2是本发明实施例的2008年7月8号的原始SAR影像。

图3是图2预处理前的示意图。

图4是图2预处理后的示意图。

图5是本发明实施例的纯净水体后向散射系数均值时间序列的曲线示意图。

图6是本发明实施例的水体像元提取结果示意图。

图7是本发明实施例的水体像元提取结果去除水体内噪点的示意图。

具体实施方式

实施例

本实施例采用的实验数据为2008年1～12月不同时期的25景宽幅模式（WSM）的ENVISAT ASAR影像，如图1所示。该数据是标准格式的1B级产品，空间分辨率为150m，像素间距为75m，极化方式有HH和VV两种。原始影像包含振幅（Amplitude）和密度（Intensity）两种信息，以2008年7月8号影像为例，原始ENVISAT ASRA影像如图2所示。

研究区选取深圳市宝安区西乡河中游段，区域范围为东经113.857度～113.991度，北纬22.584～22.719度，研究区覆盖面积为15km*15km。研究区南临南海，处于亚热带季风气候区，夏季高温多雨，其余季节较为干燥。主要地表覆盖类型包括水体、建筑物、植被、裸地等。

本实施例采用欧空局的NEST（Next ESA SAR Toolbox）软件完成像素级SAR影像时间序列构建的预处理工作，形成高精度匹配的SAR影像数据集。采用MATLAB读写SAR影像数据集，计算DTW相似度，利用最大阈值分割提取水体像元，并利用8邻域搜索方法提高精度。本实施例选用的MATLAB版本为7.0。

本实施例的基于像素级SAR影像时间序列相似性分析的水体提取方法，包括以下步骤:

第一步、对SAR影像数据集进行预处理，所述SAR影像数据集包含有N幅不同时相的SAR影像。

本实施例的SAR影像数据集含有25幅SAR影像，利用NEST软件对25幅SAR影像进行预处理，预处理包括影像的裁剪和拼接、元数据更新、辐射校正、几何校正和地形校正等。2008年7月8号的SAR影像预处理前后对比如图3和图4所示。本实施例中地形校正采用距离—多普勒地形校正方法，根据各影像元数据中的轨道状参数和下载的90m SRTM DEM数据对SAR影像进行地形纠正。

第二步、对SAR影像数据集的N幅SAR影像进行几何配准，得到配准后的SAR影像数据集。

本实施例采用NEST软件进行SAR影像的几何配准，以2008年7月8号的SAR影像为主影像，其他SAR影像为从影像，对所有从影像分别做几何配准，其中采用窗口设置为32*32，采用二阶多项式，总体RMS阈值设置为0.4，保证几何配准的误差精度保持在0.5个像素以内。

作为进一步的改进，本实施例还可以采用3×3窗口的Enhanced Lee滤波方法对SAR影像剔除光斑噪声，去除异常值。

第三步、针对坐标值相同的SAR影像像元的后向散射系数，根据时间顺序建立与坐标值相关的像元后向散射系数时间序列。

本实施例利用MATLAB按时序建立SAR影像文件列表，读取影像数据。按时序建立影像文件列表，以行为单位，依次读取每幅影像，根据像元坐标标识提取该像元在所有影像中的后向散射系数值。写出时间序列数据，每个像元的时间序列形成一条记录，将所有像素级SAR影像时间序列数据写入excel文件中,其中每一行代表一个坐标值相关的像元后向散射系数时间序列。

本实施例中SAR影像中第i行第j列像元的后向散射系数

${\mathrm{\σ}}_{\mathrm{ij}}^0=\frac{{{\mathrm{DN}}_{\mathrm{ij}}}^2}{K}\sin \left({\mathrm{\α}}_{\mathrm{ij}}\right)$

式中，DN_ij为第i行第j列像元的辐射强度值；α_ij为第i行第j列像元的的入射角；K为绝对定标因子。

第四步、在SAR影像数据集的图幅范围内目视选取M个纯净水体像元，相应的获得M个纯净水体坐标，对每幅SAR图像中M个纯净水体坐标位置的像元的后向散射系数求均值，得到SAR图像集的N个纯净水体后向散射系数均值，根据时间顺序建立纯净水体后向散射系数均值时间序列。

本实施例中M取50，将50个纯净水体坐标位置的像元的后向散射系数求均值，如图4所示，得到纯净水体后向散射系数均值时间序列

X=[-26.15343,-21.79253,-23.23705,-26.23199,-27.34405,-26.20714,-25.36589,-24.02563,-26.0951,-25.69172,-26.29467,-27.59024,-25.18612,-25.21841,-27.38671,-23.54009,-26.38829,-25.83653,-26.83594,-24.67267,-27.69025,-28.12293,-26.05999,-25.4612,-26.57747]。

本实施例的纯净水体后向散射系数均值时间序列的曲线示意图如5所示。

第五步、在SAR影像数据集的图幅范围内目视选取K个水体边缘像元，相应的获得K个水体边缘坐标，计算所述纯净水体后向散射系数均值时间序列与每个水体边缘坐标对应的所述像元后向散射系数时间序列之间的DTW距离值，并将计算得到的K个DTW距离值作平均，得到用于划分水体像元的最大阈值。

本实施例中K取50，纯净水体后向散射系数均值时间序列与每个水体边缘坐标对应的像元后向散射系数时间序列之间的DTW距离值的计算方法为现有技术范畴，本例的具体方法如下：

假设Y是水体边缘像元j坐标对应的像元后向散射系数时间序列，其长度为m，而X的长度为n：

X=a₁,a₂,…,a_n;Y=b₁,b₂,…,b_m

1）为了利用DTW将时间序列X和Y对准，需要定义一个距离矩阵D，其元素为d(i,j)=(a_i-b_j)²，即

$D=\left[\begin{array}{cccc}d(a_1,b_m)& d(a_2,b_m)& ...& d(a_n,b_m)\\ d(a_1,b_{m-1})& d(a_2,b_{m-1})& ...& d(a_n,b_{m-1})\\ ...& ...& ...& ...\\ d(a_1,b_1)& d(a_2,b_1)& ...& d(a_n,b_1)\end{array}\right]$

令Cost为最低成本矩阵，即距离矩阵中某个元素沿着弯曲路径到达初始元素的最低成本，Cost[m,n]就代表时间序列X和Y的最短距离。

2）设置Cost计算的约束条件。

在距离矩阵D中，定义一组连续矩阵元素的集合P=p₁,p₂,…,p_k，称之为弯曲路径，该路径需满足以下条件：

边界条件：p₁=D[1,1]，p_k=D[n,m]，即该路径必须从距离矩阵的第一个元素开始，并在最后一个元素截止。

连续性：若p_k=D[r,c]，p_k-1=D[r’,c’]，则必有|r-r’|≤1且|c-c’|≤1，即路径中相邻元素必须是连续的。

单调性：若p_k=D[r,c]，p_k-1=D[r’,c’]，则必有|r-r’|≥0且|c-c’|≥0，即路径的后一位置必须在前一位置的基础上保持在行方向和列方向不减，也即是说若路径P通过点[i,j]，则其必须至少通过[i-1,j]、[i,j-1]和[i-1,j-1]三个点中的一个。

有界性：max(m,n)≤k≤m+n-1，由前三条性质即可推导出此性质，即该路径所经过的矩阵元素个数存在上限和下限。

3）计算d(1,1)=(q₁–s₁)²，并将该值赋值给Cost[1,1]。

4）令i=2，3，…，m，计算

Cost[i,1]=Cost[i-1,1]+d[i,1]；式中，d[i,1]=(q_i–s₁)²。

令j=2，3，…，n，计算

Cost[1,j]=Cost[1,j-1]+d[1,j]；式中，d[1,j]=(q₁–s_j)²。

5）令i=2，3，…，m，j=2，3，…，n，循环计算Cost[i-1,j],Cost[i,j-1],Cost[i-1,j-1]的值，并将三者中的最小值作为d[i,j]的累加值，即Cost[i,j]=d[i,j]+min(Cost[i-1,j],Cost[i,j-1],Cost[i-1,j-1]；

式中d[i,j]=(q_i–s_j)²;最终计算所得的Cost[m,n]即为所求的DTW距离值。

本实施例中m,n均为25，并将50个DTW距离值作平均，得到用于划分水体像元的最大阈值为650。

第六步、构建与SAR影像数据集的图幅范围相同的空白影像，第三步得到的像元后向散射系数时间序列与纯净水体后向散射系数均值时间序列之间的DTW距离值小于或等于所述最大阈值，则空白影像中相应坐标值的像元为水体像元；否则为非水体像元。最终结果如图6所示。

本实施例还可以作以下改进，在第六步完成之后去除水体内的噪点。去除水体内噪点的判断规则为若一个像元的周围8个像元均为水体像元，即将该像元划分水体。具体方法如下：

1）将第六步得到的影像中的水体像元设为1，非水体像元设为0，从而将SAR影像赋值为0-1二值图像，假设二值图像为m行*n列，将二值影像四周扩充至(m+1)行*（n+1）列，四周像元赋值为0，得到新的二值图像；

2）选取3*3滑动窗口，从新的二值图像左上角像元（1,1）开始，依次判断滑动窗口中心点（i,j）的值，若该中心点（i，j）的八邻域像素（i-1，j-1）、(i-1,j)、(i-1,j)、(i,j-1)、(i,j+1)、（i+1，j-1）、（i+1，j）、(i+1,j+1)均为1，则将中心点像元（i,j）赋值为1；

3）将所述滑动窗口按从左到右、从上到下顺序移动到下一个像元，依次重复执行步骤2），直到全幅影像扫描一遍结束。

经过以上操作判断出的所有水体像元进行标记显示，如图7所示。

对本实施例提取的水体像元结果进行精度评价，选取完整率和正确率两个指标，其中完整率为正确提取像元数与真实水体像元总数之比；正确率为正确提取像元数与提取水体像元总数之比。利用2008年深圳市的土地利用类型变更矢量数据为参考，将其经过类型归并，采用ArcGIS软件最大面积法进行栅格化，同时归并孤立像元，作为本文典型地物提取精度评价的参考图像。通过统计，真实水体像元为4146个，本方法提取水体像元3865个，误判像元376个，经过计算，完整率为84.15%；正确率为90.27%。

本发明的基于像素级SAR影像时间序列相似性分析的水体提取方法不局限于上述实施例所述的具体技术方案，凡采用等同替换形成的技术方案均为本发明要求的保护范围。

Claims (6)

1.一种基于像素级SAR影像时间序列相似性分析的水体提取方法，包括以下步骤:

第一步、对SAR影像数据集进行预处理，所述SAR影像数据集包含有N幅不同时相的SAR影像；

第二步、对SAR影像数据集的N幅SAR影像进行几何配准，得到配准后的SAR影像数据集；

第三步、针对坐标值相同的SAR影像像元的后向散射系数，根据时间顺序建立与坐标值相关的像元后向散射系数时间序列；

第四步、在SAR影像数据集的图幅范围内目视选取M个纯净水体像元，相应的获得M个纯净水体坐标，对每幅SAR图像中M个纯净水体坐标位置的像元的后向散射系数求均值，得到SAR图像集的N个纯净水体后向散射系数均值，根据时间顺序建立纯净水体后向散射系数均值时间序列；

第五步、在SAR影像数据集的图幅范围内目视选取K个水体边缘像元，相应的获得K个水体边缘坐标，计算所述纯净水体后向散射系数均值时间序列与每个水体边缘坐标对应的所述像元后向散射系数时间序列之间的DTW距离值，并将计算得到的K个DTW距离值作平均，得到用于划分水体像元的最大阈值；

第六步、构建与SAR影像数据集的图幅范围相同的空白影像，第三步得到的像元后向散射系数时间序列与纯净水体后向散射系数均值时间序列之间的DTW距离值小于或等于所述最大阈值，则空白影像中相应坐标值的像元为水体像元；否则为非水体像元。

2.根据权利要求1所述的基于像素级SAR影像时间序列相似性分析的水体提取方法，其特征在于：第六步完成之后去除水体内的噪点。

3.根据权利要求2所述的基于像素级SAR影像时间序列相似性分析的水体提取方法，其特征在于：去除水体内噪点的方法如下：

1）将第六步得到的影像中的水体像元设为1，非水体像元设为0，从而将SAR影像赋值为0-1二值图像，假设二值图像为m行*n列，将二值影像四周扩充至(m+1)行*（n+1）列，四周像元赋值为0，得到新的二值图像；

2）选取3*3滑动窗口，从新的二值图像左上角像元（1,1）开始，依次判断滑动窗口中心点（i,j）的值，若该中心点（i，j）的八邻域像素（i-1，j-1）、(i-1,j)、(i-1,j)、(i,j-1)、(i,j+1)、（i+1，j-1）、（i+1，j）、(i+1,j+1)均为1，则将中心点像元（i,j）赋值为1；

3）将所述滑动窗口按从左到右、从上到下顺序移动到下一个像元，依次重复执行步骤2），直到全幅影像扫描一遍结束。

4.根据权利要求1所述的基于像素级SAR影像时间序列相似性分析的水体提取方法，其特征在于：所述第一步中，对SAR影像数据集的预处理包括影像的裁剪和拼接、元数据更新、辐射校正、几何校正和地形校正。

5.根据权利要求1所述的基于像素级SAR影像时间序列相似性分析的水体提取方法，其特征在于，第二步中进行几何配准的方法为：选取SAR影像数据集中任一幅SAR影像作为主影像，其他影像为从影像进行几何配准。

6.根据权利要求1所述的基于像素级SAR影像时间序列相似性分析的水体提取方法，其特征在于：SAR影像中第i行第j列像元的后向散射系数

${\mathrm{\σ}}_{\mathrm{ij}}^0=\frac{{{\mathrm{DN}}_{\mathrm{ij}}}^2}{K}\sin \left({\mathrm{\α}}_{\mathrm{ij}}\right)$

式中，DN_ij为第i行第j列像元的辐射强度值；α_ij为第i行第j列像元的的入射角；K为绝对定标因子。

Images