一种基于全极化SAR图像的植被覆盖信息提取方法

摘要

本发明公开了一种基于全极化SAR图像的植被覆盖信息提取方法，通过对获取的全极化SAR图像预处理，得到基准SAR图像，再进行相干矩阵的计算，求取其特征值和特征向量，从而进行Cloude-Pottier分解，从而获得极化熵和散射角，再依次通过H/alpha分类和Wishart H/alpha分类，识别和提取出植被覆盖的散射特点以及覆盖信息。本发明通过H/alpha分类算法和Wishart H/alpha分类算法，能够准确识别、提取植被覆盖信息，大大减少了识别误差率。

说明书

技术领域

本发明属于图像处理技术领域，更为具体地讲，涉及一种基于全极化SAR 图像的植被覆盖信息提取方法。

背景技术

植被是生态系统组成的重要部分，具有改善地方气候、防止水土流失、调 节河流流量、减轻环境污染等作用。现如今，乱砍滥伐、过度放牧以及环境污 染使得植被资源遭受到掠夺性破坏，导致生物的多样性受到破坏，水土流失严 重以及土壤荒漠化的加剧。因此，提取植被的信息，时刻关注植被的变化，对 于我国乃至世界保护生物多样性、减轻环境污染，维持社会稳定以及促进经济 发展具有重大的意义。以往，对于植被信息的提取主要是地表实测比较传统的 方法，需动用大量的人力物力收集整理资料并进行分析，还要花费大量的时间 和资金，而且结果的准确性很低，这就使得地表实测只能在很小尺度的区域提 取植被的信息；现如今，SAR技术的发展为大范围提取植被信息提供了一个新 的发展方向，速度快且准确度高。

与单极化SAR数据相比，全极化SAR数据可以提取出与数据本身无关的散 射机制的信息，该信息对于所有的全极化SAR数据均是稳定的，因此，它可以 实现非监督的分类方法，而不需要地表真实数据，或者其它地图或地理信息系 统的相关数据。同时，其含有丰富的极化特征信息，对地物几何结构特征敏感 的特点为植被信息的提取提供了新的角度和方法。如何充分利用全极化SAR数 据提取植被覆盖信息是现有技术需要解决的难题。

发明内容

本发明的目的在于克服现有技术的不足，提供一种基于全极化SAR图像的 植被覆盖信息提取方法，以实现基于SAR图像的植被覆盖信息的识别与提取。

为实现上述发明目的，本发明一种基于全极化SAR图像的植被覆盖信息提 取方法，其特征在于，包括以下步骤：

(1)、通过星载SAR获取全极化SAR图像；

(2)、对全极化SAR图像依次进行辐射定标和地理编码，得到基准SAR 图像；

(3)、计算相干矩阵K

(3.1)、用散射矩阵表示基准SAR图像的全极化特性，得到极化散射矩阵 S；

(3.2)、对极化散射矩阵S进行矢量化，得相干矩阵K，即：

$\left(\begin{array}{c}S=\left(\begin{array}{cc}{\sigma }_{\mathrm{HH}}& {\sigma }_{\mathrm{HV}}\\ {\sigma }_{\mathrm{VH}}& {\sigma }_{\mathrm{VV}}\end{array}\right)\\ K=k\times k^{*T}\\ k={[{\sigma }_{\mathrm{HH}}+{\sigma }_{\mathrm{VV}},{\sigma }_{\mathrm{HH}}-{\sigma }_{\mathrm{VV}},2\sigma ]}^T\\ \sigma =\frac{1}{2}({\sigma }_{\mathrm{HV}}+{\sigma }_{\mathrm{VH}})\end{array}\right)$

其中，σ_HH、σ_HV、σ_VH和σ_VV分别对应四个极化图像的后向散射系数值，T 表示转置，*T表示共轭转置；

(4)、对相干矩阵K进行多视处理，得到相干矩阵K^*

(5)、对相干矩阵K^*进行Cloude-Pottier分解

(5.1)、对相干矩阵K^*求取特征值λ和特征向量u；

(5.2)、根据λ和u进行Cloude-Pottier分解，得到极化熵H和平均散射角；

$\left(\begin{array}{cc}H=\underset{i=1}{\overset{3}{\Sigma }}-P_i{\log }_3{P}_i& P_i=\frac{{\lambda }_i}{\underset{i=1}{\overset{3}{\Sigma }}{\lambda }_i}\end{array}\right)$

$\left(\begin{array}{cc}\bar{\alpha }=\underset{i=1}{\overset{3}{\Sigma }}{P}_i{\alpha }_i& {\alpha }_i=\arccos \end{array},\left(\right|u_i\left(1\right)\left|\sqrt{u_i^T*u_i}\right)\right)$

其中，λ_i为相干矩阵K^*的特征值，分别为λ₁、λ₂和λ₃；u_i为对应λ_i的特征向 量，u_i(1)为特征向量u_i的第一个元素值；极化熵H是表征散射媒质的随机性， 即各散射机制在总散射过程中所占的比重，H∈[0,1]；平均散射角表征从表面 散射到二面角散射的平均散射机制，是一个从0°到90° 连续变化的参量，即

(6)、利用分类算法对基准SAR图像中的像元点初始分类

根据极化熵H和平均散射角将基准SAR图像中的像元点分为8个类别， 即：

类别1：散射特性高熵环境下的二面角散射；

类别2：散射特性高熵环境下的体散射；

类别3：散射特性中熵环境下的偶次或奇次散射；

类别4：散射特性中熵环境下的体散射；

类别5：散射特性中熵环境下的表面散射；

类别6：散射特性低熵环境下的偶次或奇次散射；

类别7：散射特性低熵环境下的体散射；

类别8：散射特性低熵环境下的表面散射；

(7)、利用分类算法对基准SAR图像中的像元点精确分类

(7.1)、计算类别m中像元点的聚类中心V_m

$V_m=\frac{1}{N_m}\underset{p\in {\omega }_m}{\Sigma }{K}_p^{*}$

其中，N_m是类别m中像元点的总数，m＝1,2,…8；ω_m是类别m中像元点的集 合，p是集合ω_m中的像元点；为像元点p的相干矩阵；

(7.2)、计算基准SAR图像中的每个像元点与聚类中心V_m的距离

$d(K^{*},V_m)=\ln \left|V_m\right|+\mathrm{Tr}\left(V_m^{-1}{K}^{*}\right)$

其中，||是求行列式，Tr是求矩阵主对角线上元素的和；

(7.3)、确定像元点的类别

选出每个像元点离聚类中心V_m的最小距离d＝min(d(K^*,V_m))，再将该像元点 归入对应的类别中；

(8)、提取植被覆盖信息

根据植被的散射类型，将分类算法所分的8个类别中，类别1、 2、4和7中的像元点视为植被，类别6中的像元点视为建筑物，类别8中的像 元点视为湖泊和裸土，类别3和5中的像元点视为建筑物或植被；

根据两类别间的距离计算公式，并结合类别3和5的聚类中心，分别计算 类别3和5的聚类中心与已经视为建筑物或植被类别的聚类中心的距离；

$D_{\hat{i}\hat{j}}=\frac{1}{2}\left[\ln \right|V_{\hat{i}}|+\ln |V_{\hat{j}}|+\mathrm{Tr}(V_{\hat{i}}^{-1}{V}_{\hat{j}}+V_{\hat{j}}^{-1}{V}_{\hat{i}})];\hat{i},\hat{j}=\mathrm{1,2}...8,\hat{i}\ne \hat{j}$

如果某个类别的聚类中心与视为植被类别的聚类中心的距离小于视为建筑 物类别的聚类中心的距离，则该类别中的像元点视为植被，反之视为建筑物。

本发明的发明目的是这样实现的：

本发明基于全极化SAR图像的植被覆盖信息提取方法，通过对获取的全极 化SAR图像预处理，得到基准SAR图像，再进行相干矩阵的计算，求取其特征 值和特征向量，从而进行Cloude-Pottier分解，从而获得极化熵和散射角，再依 次通过H/alpha分类和Wishart H/alpha分类，识别和提取出植被覆盖的散射特点 以及覆盖信息。本发明通过H/alpha分类算法和Wishart H/alpha分类算法，能够 准确识别、提取植被覆盖信息，大大减少了识别误差率。

附图说明

图1是基于全极化SAR图像的植被覆盖信息提取方法流程图；

图2是全极化SAR图像；

图3是对VV极化图像进行辐射定标后的结果图；

图4是对VV极化图像进行地理编码后的结果图；

图5是相干矩阵的计算结果图；

图6是多视处理的原理框图；

图7是多视处理前、后对比图；

图8是Cloude-Pottier分解后各像元点的极化熵H和散射角

图9是分类图；

图10是植被覆盖信息提取图。

具体实施方式

下面结合附图对本发明的具体实施方式进行描述，以便本领域的技术人员更 好地理解本发明。需要特别提醒注意的是，在以下的描述中，当已知功能和设 计的详细描述也许会淡化本发明的主要内容时，这些描述在这里将被忽略。

实施例

图1是基于全极化SAR图像的植被覆盖信息提取方法流程图。

在本实施例中，如图1所示，本发明一种基于全极化SAR图像的植被覆盖 信息提取方法，包括以下步骤：

S1、通过星载SAR获取全极化SAR图像；

在本实施例中，可以通过Radarsat-2卫星获取C波段的高分辨率全极化SAR 图像，在全极化SAR图像中，包括HH、HV、VH、VV四个极化图像，即图2 中的(a)～(d)，四幅图像为内蒙古赤峰境内的达里湖地区的四个极化图像。

S2、图像预处理

在对全极化SAR图像进行辐射定标和地理编码时，即是分别对四个极化的 图像进行辐射定标和地理编码。

其中，辐射定标是将图像的灰度值转换为后向散射系数，从而将像元点与 地物散射特性相联系；在本实施例中，以四个极化图像中的VV极化图像为例， 辐射定标后的结果如图3所示；然后，对图像进行地理编码，校正图像中的几 何畸变，得到基准SAR1图像，如图4(a)所示，同理可以得到另外三个基准 图像。为了便于处理，将5431×5712大小的基准SAR1图像裁减为901×829，如 图4(b)所示。

S3、计算相干矩阵K

S3.1)、用散射矩阵表示基准SAR图像的全极化特性，得到极化散射矩阵 S；

S3.2)、对极化散射矩阵S进行矢量化，得相干矩阵K，即：

$\left(\begin{array}{c}S=\left(\begin{array}{cc}{\sigma }_{\mathrm{HH}}& {\sigma }_{\mathrm{HV}}\\ {\sigma }_{\mathrm{VH}}& {\sigma }_{\mathrm{VV}}\end{array}\right)\\ K=k\times k^{*T}\\ k={[{\sigma }_{\mathrm{HH}}+{\sigma }_{\mathrm{VV}},{\sigma }_{\mathrm{HH}}-{\sigma }_{\mathrm{VV}},2\sigma ]}^T\\ \sigma =\frac{1}{2}({\sigma }_{\mathrm{HV}}+{\sigma }_{\mathrm{VH}})\end{array}\right)---\left(1\right)$

其中，σ_HH、σ_HV、σ_VH和σ_VV分别对应四个极化图像的后向散射系数值，T 表示转置，*T表示共轭转置；由于k为3×1的矩阵，因此K为3×3矩阵，且K有 3个特征值和3个特征向量，即λ₁、λ₂和λ₃，u₁、u₂和u₃；

将矩阵k代入K后，可表示为：

$\left(\begin{array}{c}K=k\times k^{*T}\\ \left(\begin{array}{ccc}{|{\sigma }_{\mathrm{HH}}+{\sigma }_{\mathrm{VV}}|}^2& ({\sigma }_{\mathrm{HH}}+{\sigma }_{\mathrm{VV}}){({\sigma }_{\mathrm{HH}}+{\sigma }_{\mathrm{VV}})}^{*}& ({\sigma }_{\mathrm{HH}}+{\sigma }_{\mathrm{VV}}){\left(2\sigma \right)}^{*}\\ {({\sigma }_{\mathrm{HH}}+{\sigma }_{\mathrm{VV}})}^{*}({\sigma }_{\mathrm{HH}}-{\sigma }_{\mathrm{VV}})& {|{\sigma }_{\mathrm{HH}}+{\sigma }_{\mathrm{VV}}|}^2& ({\sigma }_{\mathrm{HH}}-{\sigma }_{\mathrm{VV}}){\left(2\sigma \right)}^{*}\\ {({\sigma }_{\mathrm{HH}}+{\sigma }_{\mathrm{VV}})}^{*}\left(2\sigma \right)& {({\sigma }_{\mathrm{HH}}-{\sigma }_{\mathrm{VV}})}^{*}\left(2\sigma \right)& {\left|2\sigma \right|}^2\end{array}\right)\end{array}\right)---\left(2\right)$

相干矩阵K的得来进行详细描述，如下：

散射矩阵和相干矩阵均是以像元点为单位，即每个像元点都对应着其自身的 散射矩阵和相干矩阵，从公式(1)可知，每个相干矩阵拥有3×3个元素，因此 又可以组成一幅图像。

在本实施例中，设像元点p的相干矩阵为$K_p^{*}=\left(\begin{array}{ccc}{K}_{p}11& K_{p}12& K_{p}13\\ K_{p}21& K_{p}22& K_{p}23\\ K_{p}31& K_{p}32& K_{p}33\end{array}\right)$则元素 K_p11与基准SAR图像中其他像元点的相干矩阵的第一行第一列的元素，按照基 准SAR图像中的坐标，组成一幅基准SAR图像同样大小的图像K11，图像K11中 的各元素值为基准SAR图像中对应位置的像元点的相干矩阵的第一行第一列元 素值。同理，可得图像K12，K13，…,K33，再将图像K11，K12，…,K33作为矩阵元素， 得相干矩阵K。

在本实施例中，图5即为基准SAR图像中各像元点的相干矩阵的对应位置 的元素组成的图像，其中图5(1)～图5(9)分别对应图像K11，K12，…,K33。

S4、对相干矩阵K进行多视处理，得到相干矩阵K^*。

多视处理的原理是用邻域窗口内的均值代替矩阵中的元素值；在本实施例 中，如图6所示，以像元点p为中心(图中的深色方格)，邻域窗口选取为5×5， 计算邻域窗口中所有像元点的相干矩阵均值，再把该相干矩阵均值赋给像元点 p。在本实施中，分别对图像K11，K12，…,K33进行多视处理，再组成多视处理后 的相干矩阵，记为K^*。

在本实施例中，如图7所示，以图像K11为例，对图像K11进行多视处理， 其多视处理前、后对比如图7(a)和图7(b)所示。

S5、对相干矩阵K^*进行Cloude-Pottier分解

S5.1)、对相干矩阵K^*求取特征值λ和特征向量u；

S5.2)、根据λ和u进行Cloude-Pottier分解，得到极化熵H和平均散射角；

$\left(\begin{array}{cc}H=\underset{i=1}{\overset{3}{\Sigma }}-P_i{\log }_3{P}_i& P_i=\frac{{\lambda }_i}{\underset{i=1}{\overset{3}{\Sigma }}{\lambda }_i}\end{array}\right)$

$\left(\begin{array}{cc}\bar{\alpha }=\underset{i=1}{\overset{3}{\Sigma }}{P}_i{\alpha }_i& {\alpha }_i=\arccos \end{array},\left(\right|u_i\left(1\right)\left|\sqrt{u_i^T*u_i}\right)\right)$

其中，λ_i为相干矩阵K^*的特征值，分别为λ₁、λ₂和λ₃；u_i为对应λ_i的特征向 量，u_i(1)为特征向量u_i的第一个元素值；

极化熵H是表征散射媒质的随机性，即各散射机制在总散射过程中所占的 比重，H∈[0,1]；若H＝0，相关矩阵K^*只有一个特征值不为零，此时散射波处于 完全极化状态完全，只有一种确定性的散射机制起作用；若H＝1，有三个相等 的特征值，目标的散射退化为随机的噪声，散射波处于完全非极化状态。

平均散射角表征从表面散射到二面角散射的平均散 射机制，是一个从0°到90°连续变化的参量，即当时，散射机 制类型为各向同性的表面散射；当时，散射机制类型为偶极子散射；当 时，散射机制类型为各向同性的二面角散射。

在本实例中，图8(a)是各像元点的极化熵值，图8(b)是各像元点的平 均散射角值。

S6、利用分类算法对基准SAR图像中的像元点初始分类

根据极化熵H和平均散射角对基准SAR图像进行分类，其主要思想是根 据不同的H和值，将基准SAR图像中的像元点分为8个类别，即：

类别1：散射特性高熵环境下的二面角散射；

类别2：散射特性高熵环境下的体散射；

类别3：散射特性中熵环境下的偶次或奇次散射；

类别4：散射特性中熵环境下的体散射；

类别5：散射特性中熵环境下的表面散射；

类别6：散射特性低熵环境下的偶次或奇次散射；

类别7：散射特性低熵环境下的体散射；

类别8：散射特性低熵环境下的表面散射；

在在本实施例中，对基准SAR图像中的像元点进行分类的结果如图9 所示。

S7、利用Wishart分类算法对基准SAR图像中的像元点精确分类

Wishart分类算法需要借助分类算法的结果，从而能够更清楚的区 分开地物的类别，更符合地物的散射特性分布。下面对具体的分类过程进行描 述。

S7.1)、计算类别m中像元点的聚类中心V_m

$V_m=\frac{1}{N_m}\underset{p\in {\omega }_m}{\Sigma }{K}_p^{*}$

其中，N_m是类别m中像元点的总数，m＝1,2,…8；ω_m是类别m中像元点的集 合，p是集合ω_m中的像元点；为像元点p的相干矩阵；

此步需要计算出步骤S6中8个类别的聚类中心，得到每个类别的聚类中心 点V_m；

S7.2)、计算基准SAR图像中的像元点与聚类中心V_m的距离

$d(K^{*},V_m)=\ln \left|V_m\right|+\mathrm{Tr}\left(V_m^{-1}{K}^{*}\right)$

其中，||是求行列式，Tr是求矩阵主对角线上元素的和；

此步的运算量较大，需要计算出每个像元点到聚类中心的距离，即对于每 个像元点来说，将会计算出8个距离值；

S7.3)、确定像元点的类别

选出每个像元点离聚类中心V_m的最小距离d＝min(d(K^*,V_m))，即在8个距离 值中找出最小距离值，再结合该距离值找到聚类中心，确定出类别，最后将该 像元点归入该类别中；

通过此方法就可以将基准SAR图像中的所有像元点进行精确分类。

S8、提取植被覆盖信息

根据植被的散射类型，将Wishart分类算法所分的8个类别中，类别1、 2、4和7中的像元点视为植被，类别6中的像元点视为建筑物，类别8中的像 元点视为湖泊和裸土，类别3和5中的像元点视为建筑物或植被；

根据两类别间的距离计算公式，并结合类别3和5的聚类中心，分别计算 类别3和5的聚类中心与已经视为建筑物或植被类别的聚类中心的距离；

$D_{\hat{i}\hat{j}}=\frac{1}{2}\left[\ln \right|V_{\hat{i}}|+\ln |V_{\hat{j}}|+\mathrm{Tr}(V_{\hat{i}}^{-1}{V}_{\hat{j}}+V_{\hat{j}}^{-1}{V}_{\hat{i}})];\hat{i},\hat{j}=\mathrm{1,2}...8,\hat{i}\ne \hat{j}$

如果某个类别的聚类中心与视为植被类别的聚类中心的距离小于视为建筑 物类别的聚类中心的距离，则该类别中的像元点视为植被，反之视为建筑物。

在本实施例中，以类别3为例，计算类别3的聚类中心与类别1(植被)的 聚类中心的距离L₁，再计算类别3的聚类中心与类别6(建筑物)的聚类中心的 距离L₂，如果L₁＜L₂，则类别3中的像元点为植被，反之类别3为建筑物。

通过上述计算后，可以最终确定类别3和5中的像元点的类别，在本实施 例中，类别3和5的像元点均视为植被，最后得到如图10所示的分类图，其中 灰色区域为植被，黑色区域为建筑物，白色区域为湖泊和裸土。

尽管上面对本发明说明性的具体实施方式进行了描述，以便于本技术领域 的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对 本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定 的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发 明创造均在保护之列。