CN102194127B - 一种多频率sar数据农作物遥感分类方法 - Google Patents

Abstract

本发明涉及农作物遥感图像分类技术领域，提出了一种多频率SAR数据农作物遥感分类方法。本发明的方法中，以多景C波段ENVISAT ASAR影像与X波段TerraSAR-X影像为数据源，利用支持向量机(SVM)分类方法，通过加入纹理特征信息，实现对农作物的精确分类。该方法为农作物精确分类提供了技术支持，同时也为农作物种植面积遥感估算及产量估测提供了重要的数据支撑和基础。

Description

一种多频率SAR数据农作物遥感分类方法

技术领域

本发明涉及农作物遥感图像分类技术领域，特别涉及一种多频率SAR数据农作物遥感分类方法。

背景技术

农作物遥感监测是遥感技术在农业领域应用的重要内容。遥感技术作为地球信息科学的前沿技术，可以在短期内连续获取大范围的地面信息，实现农业信息的快速收集和定量分析，反应迅速，客观性强，是目前最为有效的对地观测技术和信息获取手段。尤其近年来新的高空间分辨率、高光谱、雷达等遥感技术的发展，为农业现代化管理提供了新的机遇。我国幅员辽阔，农作物种类丰富，如何及时、客观、准确地收集我国农作物信息，对科学指导农业生产具有重要意义。

农作物遥感识别是遥感技术在农业领域的重要应用内容，是作物种植面积估算和产量估测的基础，具有重要的意义。其中，农作物遥感估产是应用遥感信息和遥感方法估算作物产量的过程。遥感信息是指在各种遥感平台上，使用各种传感器获取作物及其环境背景的反射、辐射信息的瞬时记录。经计算机处理、识别、分类、信息提取等遥感方法，并结合数理统计分析和地学分析，最后估测出农作物的最终产量。根据遥感资料来源的不同，农作物遥感估产可分为空间遥感作物估产和地面遥感作物估产。前者又包括以应用卫星资料为主的航天遥感作物估产和以应用飞机航测资料为主的航空遥感作物估产，估产的范围广、宏观性强。后者是根据地面遥感平台获取的农作物光谱信息进行估产，估产范围较小。农作物遥感估产包括对农作物生长过程的动态监测、种植面积测算、单位面积产量估测和总产量估测。在空间遥感估产中，农作物种植面积估算是遥感估产的重要依据之一。

但是由于中国自然环境和种植制度的复杂性，作物插花种植现象普遍存在，使农作物的遥感分类难度较大，即使采用Quickbird等高空间分辨率影像仍然存在混合像元和分类精度的问题。利用微波信号与地物的相互作用，获取地表的后向散射特征，具有一定的穿透性，在对地观测中具有独特的优势，而且微波传感器可以全天时、全天候地获取数据，在农作物遥感分类上拥有较大优势，在多云雨地区则优势更加明显。然而，微波传感器接收的信号受相干斑噪声的干扰，影响了对地物目标的提取和分析，利用单景SAR(Synthetic Aperture Radar，合成孔径雷达)数据进行农作物遥感分类往往很难达到令人满意的分类精度。微波与地物相互作用的过程会随着SAR数据的频率、极化方式和入射角等因素的影响，而多频率SAR数据能够提供多频率、多极化方式和多种入射角信息。由于微波对植被的穿透力随波长的改变而变化，因而不同频率的SAR可以获取作物从冠层到茎干不同高度的信息及其下层土壤的含水量等信息，因此基于多频率SAR数据在提高农作物遥感分类精度方面具备发展潜力。

同时，由于多频率SAR数据会为农作物的分类引入更多的数据冗余，SAR数据数量的增加，不一定会对最终分类结果的精度有很明显的提高，同时还需要更多的数据处理时间和数据购买费用。因此必须发明一种能够有效地利用多频率SAR数据提高农作物分类精度的方法，这也是一项难度较大又急需解决的任务。

发明内容

(一)要解决的技术问题

针对现有技术的缺点，本发明为了解决现有技术中农作物遥感分类精度难以提高的问题，提供了一种多频率SAR数据农作物遥感分类方法，以多频率SAR影像数据为数据源，利用支持向量机(SVM)分类方法，通过加入纹理特征信息，实现对农作物的精确分类。

(二)技术方案

为实现上述目的，本发明采用如下技术方案：

一种多频率SAR数据农作物遥感分类方法，所述方法包括步骤：

S1，根据作物物候期并结合地面观测数据，确定合适时间的ASAR及TerraSAR-X数据作为数据源，并对选定的SAR数据进行预处理；

S2，利用灰度共生矩阵GLCM方法提取SAR影像的纹理特征；

S3，利用SVM方法对SAR数据进行分类；

S4，使用地面验证数据对SVM方法的分类结果进行验证。

优选地，步骤S1中，所述数据源的确定具体为：参考历史物候数据，分别在拔节初期、拔节后期和开花期获取三景ASAR数据，并获取一景与开花期的ASAR数据同时相的TerraSAR-X数据。

优选地，所述对选定的SAR数据进行预处理具体为：对所述ASAR及TerraSAR-X数据进行辐射定标、配准、几何粗纠正和滤波处理后，再进行几何精纠正。

优选地，对ASAR数据进行辐射定标按照σ⁰＝A²sinα/K进行，其中，σ⁰是定标后的后向散射系数，A是ASAR影像上像元的数字量化输出值，α是对应像元的入射角，K是原始文件的头文件中给出的绝对定标系数；

对TerraSAR-X数据进行辐射定标按照σ⁰＝[K_s(DN)²-(NEBN)]sinθ及

进行，

其中，K_s是定标系数，DN为像元的数字量化输出值，NEBN是雷达亮度β⁰的等效噪音，θ是对应像元的入射角，deg是多项式的维数，c_i是多项式第i次项的系数，τ_ref是参考时间点，τ_min和τ_max分别为距离向上第一个像元和最后一个像元记录噪音的时间，τ为距离向上对应像元记录噪音的时间。

优选地，使用NEST 3A软件对所述ASAR数据进行所述辐射定标、配准、几何粗纠正和滤波处理。

优选地，对所述ASAR及TerraSAR-X数据进行滤波均采用5*5窗口的Gamma自适应滤波器进行。

优选地，所述几何精纠正采用二次多项式方法进行，两种数据使用相同的RapidEye参考影像进行所述几何精纠正。

优选地，步骤S2中，提取所述纹理特征包括提取同质性、对比度、熵以及能量四种纹理信息，将这些纹理信息加入SAR数据中。

优选地，所述纹理信息的提取公式为：

同质性 $\mathrm{HOM}=\underset{i=0}{\overset{N-1}{\mathrm{\Σ}}}\underset{j=0}{\overset{N-1}{\mathrm{\Σ}}}\frac{P(i,j)}{1+{(i-j)}^2},$

对比度 $\mathrm{CON}=\underset{i=0}{\overset{N-1}{\mathrm{\Σ}}}\underset{j=0}{\overset{N-1}{\mathrm{\Σ}}}P(i,j)\×{(i-j)}^2,$

熵 $\mathrm{ENT}=\underset{i=0}{\overset{N-1}{\mathrm{\Σ}}}\underset{j=0}{\overset{N-1}{\mathrm{\Σ}}}-P(i,j)\×{\log }_N\left(P(i,j)\right),$

能量 $\mathrm{ASM}=\underset{i=0}{\overset{N-1}{\mathrm{\Σ}}}\underset{j=0}{\overset{N-1}{\mathrm{\Σ}}}P{(i,j)}^2,$

其中，N为灰度级别，P为N×N归一化对称GLCM矩阵的维数，P(i，j)为归一化的自相关矩阵。

优选地，步骤S3中，基于实验区农作物类型的实际分布情况，选择每种作物类型的训练样本，对SVM分类器进行训练，进而对农作物类型进行识别分类。

优选地，步骤S4中，获取与所述数据源成像时间同步的地物的真实分布数据，使用总体分类精度和Kappa系数对分类结果的精度进行评价。

(三)有益效果

本发明的方案中，以结合纹理信息的多频率SAR数据进行农作物遥感分类，充分利用了多频率SAR数据揭示地物结构信息的优势，与已有的基于遥感的农作物分类方法相比，其分类精度更高，该方法不但为农作物精确分类提供了技术支持，同时也为农作物种植面积遥感估算及产量估测提供了重要的数据支撑和基础。

附图说明

图1为本发明中多频率SAR数据农作物遥感分类方法的流程图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明的一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。

本发明涉及一种多频率SAR数据农作物遥感分类方法，下面结合附图对本发明作进一步描述。参照图1，本发明的方法具体包括如下步骤：

步骤S1中，根据作物物候期并结合地面观测数据，确定合适时间的ENVISAT ASAR及TerraSAR-X数据作为数据源，并对选定的SAR数据进行预处理，具体包括：

参考历史物候数据，分别在拔节初期、拔节后期和开花期获取一景ASAR数据(通常为C波段)，并获取一景与开花期ASAR数据同时相的TerraSAR-X数据(通常为X波段)，为了捕捉农作物生长季内的冠层结构变化，三景ASAR数据成像时间间隔要均匀，结合拔节初期、拔节后期和开花期的时间，三景影像时间间隔为30-35天为最优，TerraSAR-X数据成像时间尽量与开花期的ASAR数据的时相相一致。

确定了数据源之后，需要对SAR数据进行数据预处理，对ASAR数据的预处理在NEST 3A软件中实现，包括辐射定标、配准、几何粗纠正以及滤波等操作。ASAR的辐射定标按下式进行：

σ⁰＝A²sinα/K

式中，σ⁰是定标后的后向散射系数，A是ASAR影像上像元的数字量化输出值，α是对应像元的入射角，K是原始文件的头文件中给出的绝对定标系数；

利用NEST 3A软件对ASAR数据进行配准、几何粗纠正和滤波处理时，配准使用自动coregistration函数对数据进行配准，并利用卫星轨道参数对ASAR数据进行几何粗纠正，分辨率重采样为12.5m，滤波器采用5*5窗口的自适应Gamma滤波器，对三景ASAR影像进行滤波处理。

获取的TerraSAR-X数据在方位向上和距离向的分辨率分别为6.6m和0.9m，为了能和ASAR数据形成对应，同样将TerraSAR-X数据重采样成12.5m，对TerraSAR-X数据的预处理包括辐射定标、配准、几何粗纠正以及滤波等操作，可采用自主开发的程序进行。TerraSAR-X数据的辐射定标按下式进行：

σ⁰＝[K_s(DN)²-(NEBN)]sinθ

$\mathrm{NEBN}=K_s\underset{i=0}{\overset{\mathrm{deg}}{\mathrm{\Σ}}}{c}_i{(\mathrm{\τ}-{\mathrm{\τ}}_{\mathrm{ref}})}^i,\mathrm{\τ}\∈[{\mathrm{\τ}}_{\min },{\mathrm{\τ}}_{\max }]$

式中，K_s是定标系数，DN(Digital Number)为像元的数字量化输出值，NEBN(Noise equivalent beta naught)是雷达亮度β⁰的等效噪音，这一参数反映了各种噪音对雷达信号的影响，θ是对应像元的入射角，deg是多项式的维数，c_i是多项式第i次项的系数，τ_ref是参考时间点，τ_min和τ_max分别为距离向上第一个像元和最后一个像元记录噪音的时间，τ为距离向上对应像元记录噪音的时间，K_s、τ_ref、τ_min以及τ_max这些参数均可以从原始文件的头文件中获取。

TerraSAR-X数据的滤波与ASAR相一致，同样采用5*5窗口的自适应Gamma滤波器完成滤波处理。

在经过辐射定标、配准、几何粗纠正及滤波等处理之后，需要进行影像的几何精纠正。为了使TerraSAR-X数据与ASAR数据能够在地理位置上实现准确对应，利用同样的RapidEye影像作为参考影像，对ASAR数据及TerraSAR-X数据进行几何精纠正，纠正在ERDASIMAGINE 8.6软件中采用二次多项式方法进行，控制误差在0.5个像元以内。

步骤S2中，利用灰度共生矩阵(GLCM)方法提取SAR影像的纹理特征，具体包括：

对经过预处理后的开花期的ASAR和TerraSAR-X两景数据进行纹理特征提取，包括同质性、对比度、熵以及能量四种纹理信息，同质性、对比度、熵以及能量四种信息的计算公式按下式进行：

同质性： $\mathrm{HOM}=\underset{i=0}{\overset{N-1}{\mathrm{\Σ}}}\underset{j=0}{\overset{N-1}{\mathrm{\Σ}}}\frac{P(i,j)}{1+{(i-j)}^2},$

对比度： $\mathrm{CON}=\underset{i=0}{\overset{N-1}{\mathrm{\Σ}}}\underset{j=0}{\overset{N-1}{\mathrm{\Σ}}}P(i,j)\×{(i-j)}^2,$

熵： $\mathrm{ENT}=\underset{i=0}{\overset{N-1}{\mathrm{\Σ}}}\underset{j=0}{\overset{N-1}{\mathrm{\Σ}}}-P(i,j)\×{\log }_N\left(P(i,j)\right),$

能量： $\mathrm{ASM}=\underset{i=0}{\overset{N-1}{\mathrm{\Σ}}}\underset{j=0}{\overset{N-1}{\mathrm{\Σ}}}P{(i,j)}^2,$

其中，N灰度级别，P为N×N归一化对称GLCM矩阵的维数，P(i，j)归一化的自相关矩阵，该自相关矩阵是通过5*5大小的窗口计算得到的，灰度级别被设定为64以获得纹理信息的均值；

具体地，灰度共生矩阵GLCM是本领域中通过研究灰度的空间相关特性来描述纹理的一种常用方法，其定义为：对于取定的方向θ和距离d，在方向为θ的直线上，一个像元灰度为i，另一个与其相距为d像元的灰度为j的点对出现的频数即为灰度共生矩阵第(i，j)阵元的值。

从SAR图像中提取纹理信息之后，将这些纹理信息加入SAR数据中，以提高农作物遥感分类的精度。

其中，步骤S3包括：

利用SVM方法对SAR数据进行分类，具体包括：

与数据获取同步，开展了地面观测及填图工作，获取了地物的实际分布知识。基于对实验区农作物类型的实际分布情况，选择每种作物类型的训练样本，对支持向量机分类器进行训练，进而对农作物类型进行识别分类。

具体地，现代的支持向量机SVM理论是由Cortes和Vapnik提出的，支持向量机基于结构风险最小化理论在特征空间中根据训练样本集构建最优分类超平面，使得学习器得到全局最优化。

其中，步骤S4包括：

利用地面验证数据对SVM方法的分类结果进行验证，具体包括：

与开花期的ASAR和TerraSAR-X数据成像时间同步的，还开展了地面观测及填图工作，获取了地物的真实分布数据。将分类结果与地面实测获取的地物类型进行比较后，导出分类结果的混淆矩阵，使用总体分类精度和Kappa系数对分类结果的精度进行评价，总体分类精度与Kappa系数通过检验样本的导出的混淆矩阵来估算得到。

举例来说，表1是根据某次观测时获取的地物真实分布数据导出的分类结果混淆矩阵：

表1 多频率SAR数据分类结果混淆矩阵

总体分类精度及Kappa系数公式如下：

总体分类精度： $p_c=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{p}_{\mathrm{ii}}/N,$

Kappa系数计算公式： $K_{\mathrm{hat}}=\frac{N\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{p}_{\mathrm{ii}}-\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{p}_{i+}{p}_{+i}}{N^2-\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{p}_{i+}{p}_{+i}},$

式中，p_ii为混淆矩阵第i行，第i列对应的像元数，n代表分类结果的类别数量，N为所有用于精度评价的像元总数，p_i+和p_+i和分别是混淆矩阵第i行和第i列的总像元数。

利用表1中的数据，计算多频率SAR数据农作物分类结果的总体精度为91.83％，Kappa系数K_hat为85.72％，以此作为评价结果对支持向量机分类器进行验证和反馈。

与现有技术相比，本发明具有以下明显的技术优势：

1.充分利用SAR数据具有一定的穿透性，并且可以全天时、全天候地获取数据的优势，在多云雨地区进行农作物分类优势明显；

2.采用了结合纹理信息的多频率SAR数据农作物遥感分类方法，充分利用了多频率SAR数据揭示地物结构信息的优势，与已有的基于遥感的农作物分类方法相比，其分类精度更高；

3.以ENVISAT ASAR数据及TerraSAR-X数据为数据源，具有较高的空间分辨率，在种植制度复杂的地区，仍然能够保持较高的农作物遥感分类精度。

以上实施方式仅用于说明本发明，而并非对本发明的限制，有关技术领域的普通技术人员，在不脱离本发明的精神和范围的情况下，还可以做出各种变化和变型，因此所有等同的技术方案也属于本发明的范畴，本发明的专利保护范围应由权利要求限定。

Claims (8)

1.一种多频率SAR数据农作物遥感分类方法，其特征在于，所述方法包括步骤：

S1，根据作物物候期并结合地面观测数据，确定合适时间的ASAR及TerraSAR-X数据作为数据源，并对选定的SAR数据进行预处理；

S2，利用灰度共生矩阵GLCM方法提取SAR影像的纹理特征；

S3，利用SVM方法对SAR数据进行分类；

S4，使用地面验证数据对SVM方法的分类结果进行验证；

在步骤S1中，所述数据源的确定具体为：参考历史物候数据，分别在拔节初期、拔节后期和开花期获取三景ASAR数据，并获取一景与开花期的ASAR数据同时相的TerraSAR-X数据；

所述对选定的SAR数据进行预处理具体为：对所述ASAR及TerraSAR-X数据进行辐射定标、配准、几何粗纠正和滤波处理后，再进行几何精纠正；

在步骤S2中，提取所述纹理特征包括提取同质性、对比度、熵以及能量四种纹理信息，将这些纹理信息加入SAR数据中。

2.根据权利要求1所述的方法，其特征在于，对ASAR数据进行辐射定标按照σ⁰＝A²sinα/K进行，其中，σ⁰是定标后的后向散射系数，A是ASAR影像上像元的数字量化输出值，α是对应像元的入射角，K是原始文件的头文件中给出的绝对定标系数；

对TerraSAR-X数据进行辐射定标按照σ⁰＝[K_s(DN)²-(NEBN)]sinθ及 $\mathrm{NEBN}=K_s\underset{i=0}{\overset{\mathrm{deg}}{\mathrm{\Σ}}}{c}_i{(\mathrm{\τ}-{\mathrm{\τ}}_{\mathrm{ref}})}^i,$ τ∈[τ_min,τ_max]进行，

其中，K_s是定标系数，DN为像元的数字量化输出值，NEBN是雷达亮度β⁰的等效噪音，θ是对应像元的入射角，deg是多项式的维数，c_i是多项式第i次项的系数，τ_ref是参考时间点，τ_min和τ_max分别为距离向上第一个像元和最后一个像元记录噪音的时间，τ为距离向上对应像元记录噪音的时间。

3.根据权利要求1所述的方法，其特征在于，使用NEST 3A软件对所述ASAR数据进行所述辐射定标、配准、几何粗纠正和滤波处理。

4.根据权利要求1所述的方法，其特征在于，对所述ASAR及TerraSAR-X数据进行滤波均采用5*5窗口的Gamma自适应滤波器进行。

5.根据权利要求1所述的方法，其特征在于，所述几何精纠正采用二次多项式方法进行，两种数据使用相同的RapidEye参考影像进行所述几何精纠正。

6.根据权利要求1所述的方法，其特征在于，所述纹理信息的提取公式为：

同质性 $\mathrm{HOM}=\underset{i=0}{\overset{N-1}{\mathrm{\Σ}}}\underset{j=0}{\overset{N-1}{\mathrm{\Σ}}}\frac{P(i,j)}{1+{(i-j)}^2},$

对比度 $\mathrm{CON}=\underset{i=0}{\overset{N-1}{\mathrm{\Σ}}}\underset{j=0}{\overset{N-1}{\mathrm{\Σ}}}P(i,j)\×{(i-j)}^2,$

熵 $\mathrm{ENT}=\underset{i=0}{\overset{N-1}{\mathrm{\Σ}}}\underset{j=0}{\overset{N-1}{\mathrm{\Σ}}}-P(i,j)\×{\log }_N\left(P(i,j)\right),$

能量 $\mathrm{ASM}=\underset{i=0}{\overset{N-1}{\mathrm{\Σ}}}\underset{j=0}{\overset{N-1}{\mathrm{\Σ}}}P{(i,j)}^2,$

其中，N为灰度级别，P为N×N归一化对称GLCM矩阵的维数，P(i，j)为归一化的自相关矩阵。

7.根据权利要求1所述的方法，其特征在于，步骤S3中，基于实验区农作物类型的实际分布情况，选择每种作物类型的训练样本，对SVM分类器进行训练，进而对农作物类型进行识别分类。

8.根据权利要求1所述的方法，其特征在于，步骤S4中，获取与所述数据源成像时间同步的地物的真实分布数据，使用总体分类精度和Kappa系数对分类结果的精度进行评价。

Images