CN104008552A - 基于动态时间弯曲的时序sar影像耕地提取方法 - Google Patents

Abstract

本发明涉及一种基于动态时间弯曲的时序SAR影像耕地提取方法，步骤包括：构建时序SAR影像；提取耕地参考时间序列；计算待分类像元时间序列与耕地参考时间序列之间的动态时间弯曲距离；计算结果阈值分割，待分类像元归类为耕地与非耕地；分割结果空域滤波，滤除孤立的耕地像元，填补连片耕地之间的缝隙，得到耕地的最终提取结果。本发明考虑到耕地时间序列特有的“时间轴弯曲”现象，使用动态时间弯曲距离（DTW）作为相似性度量标准，从而实现耕地像元与非耕地像元的划分，解决了传统方法无法适应时间轴畸变的时间序列相似性度量这一问题，提高了耕地的提取精度。本发明方法适应性强，提取精度可达82%以上，能够满足实际生产的需要。

Description

基于动态时间弯曲的时序SAR影像耕地提取方法

技术领域

本发明涉及一种遥感影像信息提取方法，特别涉及一种基于动态时间弯曲的时序SAR影像耕地提取方法。

背景技术

耕地是粮食生产的基本条件，是确保粮食安全的物质基础，掌握实时的耕地信息对国家农业政策的制定具有重要意义。遥感技术经过几十年的发展，已经形成了“星——空——地”一体的对地观测网络，并产生了海量堆积数据，因此使用遥感技术进行大范围耕地信息的快速提取领域具有独到的优势。

耕地表面附着的作物拥有随季节变化十分明显的物候特征，不同生长周期(如播种期、营养生长期、生殖生长期和成熟期等)在遥感影像上有着不同的表现，因此使用单一时间节点的遥感影像往往不能有效提取出耕地信息。时间序列遥感影像(以下简称时序遥感影像)是卫星传感器对同一地区不同时间进行观测所获得的数值序列，通过构建时序遥感影像，可以识别观测数值的时间变化特征，并利用该特征进行耕地信息的有效提取。然而受观测条件多变、农作物物候差异较大等因素的影响，现阶段利用时序遥感影像反演耕地的工作仍存在不少问题，大致可归纳为两个方面：

1)所使用遥感数据的不足。当前用于耕地反演的时序遥感数据大多数为SPOT/VEGETATION、NOAA/AVHRR、MODIS等影像数据构建的NDVI、EVI等植被指数产品(Sakamoto,2005；Wardlow,2008；Dheeravath,2010)，它们拥有较高的时间分辨率，然而空间分辨率很低(最高仅250m)，复合像元问题十分严重；同时，卫星传感器的光学观测过程中常受到云层遮挡，得到的时间序列存在大量噪声；此外，NDVI本身会受到土壤背景和含水量等因素较大的影响，且在高植被覆盖区域容易达到饱和(Baush,1993；唐怡，2006)。合成孔径雷达(SAR)为主动式传感器，发射的微波可以穿透云层，能够全天时不受天气影响完成对地观测任务。随着技术的发展，SAR传感器的空间分辨率以及时间分辨率都有了极大的提升，欧空局的ENVISAT ASAR数据空间分辨率最高可达10m，中国HJ-1C星SAR数据空间分辨率最高可达4m，德国TerraSAR-X、加拿大RADARSAT-2的SAR数据空间分辨率最高可达1m；同时，SAR可侧视成像的特性保证了影像的高时间分辨率，许多卫星传感器能够提供时间间隔小于10天的SAR影像数据(TerraSAR-X对全球95％的地区可实现2天重访)。然而尽管SAR影像应用范围越来越广，却少有将其构建时间序列进行耕地信息提取的研究。

2)所采用方法的不足。利用时序遥感影像提取耕地的关键技术是时间序列相似性度量，即定量评估待分类像元时间序列与耕地参考时间序列的相似性程度。对大区域耕地而言，其种植的农作物种类多样，不同种类农作物的物候期开始时间不一，持续长短不同，致使耕地时间序列在时间轴上存在一定的平移和弯曲(Calendar Shift)。而当前进行遥感时间序列相似性度量的方法大体上可分为三类：距离法(欧氏距离、马氏距离等)、代数法(作差法、作商法等)和相关法(Pearson相关系数、Spearman相关系数等)，它们均无法适用于“时间轴移位”的时间序列。动态时间弯曲(Dynamic time warping,DTW)是一个在计算机科学中常用的时间序列相似度算法，它可以在一定程度上自动寻找发生平移和弯曲的时间节点，并对序列的时间轴进行动态弯曲以实现自适应匹配，从而计算动态弯曲距离，作为两条时间序列的相似度，从机理上看十分适于耕地时间序列的相似性分析，然而相应的研究或应用却并不多见。

发明内容

本发明要解决的技术问题是：克服现有技术缺点，提出一种基于动态时间弯曲的时序SAR影像耕地提取方法。

本发明基于动态时间弯曲的时序SAR影像耕地提取方法，包括以下步骤：

第一步、时序SAR影像构建——对研究区研究时段内的原始SAR影像数据进行预处理和配准，并按照时间先后顺序构建时序SAR影像；

第二步、耕地参考时间序列提取——选择研究区内的若干耕地像元作为训练样本，针对所述时序SAR影像中的每一幅影像，分别计算训练样本中耕地像元的后向散射系数的均值和标准差，若所述标准差之和小于设定的可靠性阈值，则利用所述均值依据对应影像的成像时间构建耕地参考时间序列；否则，重新选择训练样本并重复本步骤，直至标准差小于可靠性阈值；

第三步、时间序列相似度计算——针对每个待分类像元，构建其后向散射系数时间序列，并计算所述后向散射系数时间序列与耕地参考时间序列之间的动态时间弯曲距离，作为两时间序列之间的相似度；

第四步、计算结果阈值分割——针对相似度计算结果，使用预先设定的分割阈值进行二值分割，将动态时间弯曲距离大于分割阈值的对应待分类像元归类为非耕地像元，其余归为耕地像元，获得阈值分割结果；

第五步、分割结果空域滤波——使用空间操作算子对所述阈值分割结果进行空域滤波，滤除孤立的耕地像元，填补连片耕地之间的缝隙，得到耕地的最终提取结果。

本发明基于动态时间弯曲的时序SAR影像耕地提取方法，还具有如下改进：

1、所述第一步中，预处理至少包括轨道参数更新、辐射校正、地形校正、几何校正和去噪滤波。

2、所述第一步中，影像的配准误差不大于0.5像素。

3、所述第二步中，构成训练样本的耕地像元通过人工识别提取，训练样本中的耕地像元数量不低于9个。

4、第二步中的可靠性阈值Th1通过公式Th1＝N×c×R计算获得，其中，N为时序SAR影像中的影像数量，c为比率系数，取值范围为[0.3,0.5]，R为所有影像的训练样本中耕地像元的后向散射系数均值的最大值与最小值之差。

5、所述第三步中，利用待分类像元在每一幅影像中的后向散射系数按照对应影像的成像时间构建待分类像元的后向散射系数时间序列。

6、所述第四步中，分割阈值Th2通过公式Th2＝c₁×Σ计算获得，其中，c₁为调节系数，其取值范围为0.3～1，Σ为第二步中所有标准差的求和。

7、所述第五步中，空间操作算子为众数滤波算子、数学形态学开闭算子、中值滤波算子中的一种，空间操作算子的尺寸size根据下式获得：

其中，re为SAR影像的空间分辨率，ps为SAR影像的像元大小，int′()为四舍五入取整函数。

本发明的有益效果如下：

1)、本发明使用较高分辨率合成孔径雷达(SAR)后向散射系数构建时间序列，耕地提取结果的空间分辨率可以达到30米(如ENVISAT ASAR、ERS-1/2SAR等卫星传感器)甚至更高，对空间的辨析程度较传统的MODIS/NDVI时间序列等产品(最高空间分辨率250米)提高了8倍以上。

2)、考虑到耕地时间序列特有的“时间轴弯曲”现象，本发明使用动态时间弯曲距离(DTW)作为相似性度量标准，从而实现耕地像元与非耕地像元的划分，解决了传统方法无法适应时间轴畸变的时间序列相似性度量这一问题，从理论上提高了耕地的提取精度。

3)、本发明针对训练样本提出了可靠性阈值，确保训练样本的可靠性，并且针对动态距离分割提出了分割阈值的选取方法，确保最终耕地提取的精度。

4)、本发明方法适应性强，经实践证明，使用本发明方法提取耕地，提取精度可达82％以上，能够满足实际生产的需要，相比传统的耕地提取方法具有更强实用性。

附图说明

下面结合附图对本发明作进一步的说明。

图1为本发明实施例中实验区区位图。

图2为本发明实施例中耕地参考时间序列示意图。

图3-a为“时间轴弯曲”时间序列示意图。

图3-b为本发明实施例中动态弯曲路径示意图。

图4为本发明实施例中研究区耕地提取结果图(深色表示耕地)。

具体实施方式

下面根据附图详细阐述本发明，使本发明的技术路线和操作步骤更加清晰。

本实施例选择江苏省徐州市贾汪区西部地区为研究区，经度范围为东经117.28°～117.56°，纬度范围为北纬34.30°～34.54°，面积614.57km²(见图1)。该研究区属暖温带半湿润季风气候，四季分明，年平均气温14.2℃，平均降雨量802.4mm，研究区内耕作方式为一年两熟，5月至10月种植一季稻或夏玉米，10月至次年5月种植冬小麦。

实验数据为覆盖该区域的2009年12月至2010年12月共11景ENVISAT ASAR图像模式数据，L1级产品，经过了系统几何校正，空间分辨率30米，像元大小20.25米，VV极化，每景影像获取时间见表1。

表1  实验所用ENVISAT ASAR影像获取时间

本发明实施例基于动态时间弯曲的时序SAR影像耕地提取方法，包括以下步骤：

第一步、时序SAR影像构建——对研究区研究时段内的原始SAR影像数据进行预处理和配准，并按照时间先后顺序构建时序SAR影像。

首先使用欧空局提供的NEST软件依次对原始SAR影像数据进行轨道参数更新、辐射校正、地形校正、几何校正和去噪滤波等预处理操作，其中地形校正使用Range-Doppler模型和ASTER GDEM数据，去噪滤波采用3*3大小的EnhancedLee滤波算子。然后使用局部自适应配准方法(Cheng,L.et al.,Generationof Pixel-Level SAR Image Time Series Using a Locally Adaptive MatchingTechnique,2014)对所有影像进行配准，影像的配准误差不大于0.5像素，最后，按照时间影像获取时间的先后顺序构建时序SAR影像。

第二步、耕地参考时间序列提取——选择研究区内的若干耕地像元作为训练样本，针对所述时序SAR影像中的每一幅影像，分别计算训练样本中耕地像元的后向散射系数的均值和标准差，若所述标准差之和小于设定的可靠性阈值，则利用所述均值依据对应影像的成像时间构建耕地参考时间序列；否则，重新选择训练样本并重复本步骤，直至标准差小于可靠性阈值。

本步骤中，在ENVI中打开构建好的时序SAR影像，对照Google Earth历史影像在时序SAR影像中选取耕地训练样本，根据设定的可靠性阈值Th1来对训练样本进行可靠性评判。

可靠性阈值Th1通过公式Th1＝N×c×R计算获得，其中，N为时序SAR影像中的影像数量，本例中N＝11，c为比率系数，取值范围为[0.3,0.5]，本例中取0.3，R为所有影像的训练样本中耕地像元的后向散射系数均值的最大值与最小值之差。

对所有影像的训练样本耕地像元后向散射系数标准差求和，若标准差之和大于可靠性阈值Th1时，需重新选择样本并重复上述可靠性评价过程，直至标准差之和不大于Th1，将各时间节点均值按先后顺序排列，得到耕地参考时间序列。

本实施例共选取了研究区1000个耕地像元(为了确保提取精度，最少不能低于9个)，得到的耕地时间序列S＝[-9.762，-13.898，-11.319，-12.304，-15.811，-14.235，-11.285，-8.761，-8.753，-9.585，-11.945](参见图2)。每幅影像训练样本的后向散射系数标准差按照时间先后排序，STD＝[3.498，2.144，3.169，2.640，2.939，2.813，2.592，2.716，3.021，3.134，3.165]。

第三步、时间序列相似度计算——针对每个待分类像元，构建其后向散射系数时间序列Q，并计算所述后向散射系数时间序列Q与耕地参考时间序列S之间的动态时间弯曲距离，作为两时间序列之间的相似度。

本实施例使用Microsoft Visual Studio2010和C#编程，计算待分类像元后向散射系数时间序列Q与耕地参考时间序列S的动态弯曲距离。

时间轴弯曲的同类时间序列(见图3-a)之间的动态时间弯曲距离计算，属于现有成熟算法，本例中给出一种算法，对动态时间弯曲距离进行说明。但本发明中所涉动态弯曲距离不限于这种算法获得。

计算时间序列Q(假设Q＝q₁,q₂,...,q_n)与时间序列S(假设S＝s₁,s₂,...,s_m)的动态弯曲距离(本例中下标n＝m)，具体步骤如下：

a1)、计算Q与S的距离矩阵D，距离矩阵D的元素d(i,j)＝|q_i-s_j|，即：

$D=\left[\begin{array}{cccc}d(q_1,s_m)& d(q_2,s_m)& ...& d(q_n,s_m)\\ d(q_1,s_{m-1})& d(q_2,s_{m-1})& ...& d(q_n,s_{m-1})\\ ...& ...& ...& ...\\ d(q_1,s_1)& d(q_2,s_1)& ...& d(q_n,s_1)\end{array}\right]$

a2)、定义与距离矩阵D维度相同的最小代价矩阵C，矩阵C中元素的值c(i,j)表示在距离矩阵D中从原点(1,1)处经过一条弯曲路径到达位置(i,j)处的最小路径代价；

a3)、根据距离矩阵D计算矩阵C中各元素的值，其中最后一个元素c(n,m)即为时间序列Q与时间序列S的动态弯曲距离。

上述步骤a2)中定义的弯曲路径为距离矩阵D中一组连续元素的有序集合P＝p₁,p₂,p₃,...,p_k(见图3-b)，且需满足以下条件：

Ⅰ)边界条件：p₁＝d(1,1)，p_k＝d(n,m)，即该路径的起止元素为距离矩阵的斜对角线的两端元素；

Ⅱ)连续性：若有p_k＝d(i,j)和p_k-1＝d(i',j')，则必有|i-i'|≤1且|j-j'|≤1。即路径中相邻元素必须连续；

Ⅲ)单调性：若有p_k＝d(i,j)和p_k-1＝d(i',j')，则必有i-i'≥0且j-j'≥0。即若路径P通过位置(i,j)，则其必须至少通过(i-1,j)、(i,j-1)和(i-1,j-1)三个位置中的一个，保证路径的后一位置在前一位置的基础上保持在行方向和列方向不减；

Ⅳ)有界性：max(m,n)≤k≤m+n-1。即该路径所经过的矩阵元素个数存在上限和下限。

上述步骤a3)中计算最小代价矩阵C中各个元素值的方法如下：

b1)、将距离矩阵D中原点位置(1,1)处元素值赋给最小代价矩阵C中原点位置的元素，即c(1,1)＝d(1,1)；

b2)、计算第一行元素值，c(1,j)＝c(1,j-1)+d(1,j)，其中j＝2,3,...,m；

b3)、计算第一列元素值，c(i,1)＝c(i-1,1)+d(i,1)，其中i＝2,3,...,n；

b4)、计算矩阵其余元素值，计算公式如下：

c(i,j)＝min(c(i-1,j),c(i,j-1),c(i-1,j-1))+d(i,j)

式中，c(i,j)代表位置(i,j)的元素，min()为最小值取值函数，i＝2,3,...,n，j＝2,3,...,m。

第四步、计算结果阈值分割——针对相似度计算结果，使用预先设定的分割阈值进行二值分割，将动态时间弯曲距离大于分割阈值的对应待分类像元归类为非耕地像元，其余归为耕地像元，获得阈值分割结果。

本实施例中，分割阈值Th2通过公式Th2＝c₁×Σ计算获得，其中，c₁为调节系数，其取值范围为0.3～1，本例中c₁取0.3，Σ为第二步中所有标准差的求和，本例中Σ＝31.832。于是，计算得到分割阈值Th2为9.55。最终获得的分割结果可以使用二值图像来表示，例如，1代表耕地，0代表非耕地。

第五步、分割结果空域滤波——使用空间操作算子对所述阈值分割结果进行空域滤波，滤除孤立的耕地像元，填补连片耕地之间的缝隙，得到耕地的最终提取结果，见图4。图4中深色的像元为提取的耕地。

第4步中使用众数滤波算子(Marjority)进行分割结果的空域滤波(此外还可以使用数学形态学开闭算子、中值滤波算子等)，空间操作算子的尺寸size根据下式获得：

其中，re为SAR影像的空间分辨率，此处为30，ps为SAR影像的像元大小，此处为20.25，int′()为四舍五入取整函数。最终计算得到空间滤波算子的尺寸为3*3。

利用地面真实数据对所提取耕地的精度验证结果为：生产者精度84.37％，用户精度82.19％，表明本发明方法有效进行较高空间分辨率的耕地信息提取，且精度较高，满足生产需要。

除上述实施例外，本发明还可以有其他实施方式。凡采用等同替换或等效变换形成的技术方案，均落在本发明要求的保护范围。

Claims (8)

1.一种基于动态时间弯曲的时序SAR影像耕地提取方法，包括以下步骤：

第一步、时序SAR影像构建——对研究区研究时段内的原始SAR影像数据进行预处理和配准，并按照时间先后顺序构建时序SAR影像；

第二步、耕地参考时间序列提取——选择研究区内的若干耕地像元作为训练样本，针对所述时序SAR影像中的每一幅影像，分别计算训练样本中耕地像元的后向散射系数的均值和标准差，若所述标准差之和小于设定的可靠性阈值，则利用所述均值依据对应影像的成像时间构建耕地参考时间序列；否则，重新选择训练样本并重复本步骤，直至标准差小于可靠性阈值；

第三步、时间序列相似度计算——针对每个待分类像元，构建其后向散射系数时间序列，并计算所述后向散射系数时间序列与耕地参考时间序列之间的动态时间弯曲距离，作为两时间序列之间的相似度；

第四步、计算结果阈值分割——针对相似度计算结果，使用预设的分割阈值进行二值分割，将动态时间弯曲距离大于分割阈值的对应待分类像元归类为非耕地像元，其余归为耕地像元，获得阈值分割结果；

第五步、分割结果空域滤波——使用空间操作算子对所述阈值分割结果进行空域滤波，滤除孤立的耕地像元，填补连片耕地之间的缝隙，得到耕地的最终提取结果。

2.根据权利要求1所述的基于动态时间弯曲的时序SAR影像耕地提取方法，其特征在于：所述第一步中，预处理至少包括轨道参数更新、辐射校正、地形校正、几何校正和去噪滤波。

3.根据权利要求1所述的基于动态时间弯曲的时序SAR影像耕地提取方法，其特征在于：所述第一步中，影像的配准误差不大于0.5像素。

4.根据权利要求1所述的基于动态时间弯曲的时序SAR影像耕地提取方法，其特征在于：所述第二步中，构成训练样本的耕地像元通过人工识别提取，训练样本中的耕地像元数量不低于9个。

5.根据权利要求1所述的基于动态时间弯曲的时序SAR影像耕地提取方法，其特征在于：第二步中的可靠性阈值Th1通过公式Th1＝N×c×R计算获得，其中，N为时序SAR影像中的影像数量，c为比率系数，取值范围为[0.3,0.5]，R为所有影像的训练样本中耕地像元的后向散射系数均值的最大值与最小值之差。

6.根据权利要求1所述的基于动态时间弯曲的时序SAR影像耕地提取方法，其特征在于：所述第三步中，利用待分类像元在每一幅影像中的后向散射系数按照对应影像的成像时间构建待分类像元的后向散射系数时间序列。

7.根据权利要求1所述的基于动态时间弯曲的时序SAR影像耕地提取方法，其特征在于：所述第四步中，分割阈值Th2通过公式Th2＝c₁×Σ计算获得，其中，c₁为调节系数，其取值范围为0.3～1，Σ为第二步中所有标准差的求和。

8.根据权利要求1所述的基于动态时间弯曲的时序SAR影像耕地提取方法，其特征在于：所述第五步中，空间操作算子为众数滤波算子、数学形态学开闭算子、中值滤波算子中的一种，空间操作算子的尺寸size根据下式获得：

其中，re为SAR影像的空间分辨率，ps为SAR影像的像元大小，int′()为四舍五入取整函数。