CN103593669A - 一种极化合成孔径雷达影像四分量分解的方法 - Google Patents

Abstract

本发明公开了一种极化合成孔径雷达影像四分量分解的方法，它涉及极化SAR影像地物散射特性的分析，尤其涉及对现有极化SAR影像分解方法中难以处理的具有定向角不为0度的建筑物散射类型的计算方法，如图1所示，所述方法包括：一、针对影像中定向角不为0度的建筑物所对应的像素，提出筛选算法；二、针对定向角不为0度的建筑物提出能够描述其散射特征的广义偶次散射模型；三、提出非负特征值分解技术解算散射分量；通过本发明提出的方法可以有效避免原有散射模型带来的不足，对于高分辨率全极化SAR影像能够获得满意的影像分解结果，并具有新颖性、创造性和简单实用的特点。

Description

一种极化合成孔径雷达影像四分量分解的方法

技术领域

本发明涉及到极化合成孔径雷达（PolarimetricSyntheticApertureRadar，POLSAR）影像的散射信息提取，以及散射类型的建模，散射分量的估计与获取，具体为一种极化合成孔径雷达影像四分量分解的方法。

背景技术

传统的合成孔径雷达（SyntheticApertureRadar，SAR）采用单一极化通道发射和接收电磁波，只能获取地物对某一种极化方式的电磁波的散射特性，获取的地物散射信息有限。POLSAR作为SAR新的发展方向，通过两个通道发射和接收不同极化方式的电磁波，获取四组测量数据，构成一个完备的极化基，以此来探测地面目标对电磁波的散射特性。与传统的单极化SAR仅能获取目标的雷达散射截面积（RadarCrossSection，RCS）相比，散射目标对不同极化方式的电磁波具有不同的调制作用，极化SAR可以获得反映分辨单元散射特性的极化散射矩阵、Mueller矩阵或者Stokes矩阵。极化散射矩阵将目标散射的能量特性、相位特性以及极化特性统一起来，完整地描述雷达目标的电磁散射过程。极化散射矩阵蕴含着更加丰富的地物信息，使人们可以更为深入地提取与分析地物的物理特性，极大地促进了参数反演的研究。目标的极化散射特性与其形状结构有着本质的联系，可反映目标表面粗糙度、对称性和取向等其他雷达参数不能提供的信息，是完整刻画目标特性不可或缺的。

从极化SAR影像数据中，可以提取目标的极化散射特性，包括目标定向角、散射类型、螺旋性、对称性等，用以细致刻画目标的散射特征，以及进一步实现全极化数据的目标检测和目标识别等应用。其中提取目标的散射特性的理论核心是目标的极化分解技术，即对目标散射矩阵进行分解，使人们更加充分地认识与分析目标的散射机理，进而提取出目标的各种散射特征参数，反演出地面目标的信息。因此目标极化分解与散射特征提取是极化SAR应用研究的一个必要组成部分，它在分类、目标检测、目标识别、参数反演、影像增强融合这些具体研究方向中都有许多实际应用。

Yamaguchi分解是目前国际上常用的极化SAR四分量分解方法，该方法虽然能够将实际地物的散射机制与地物的物理散射模型联系起来，但受限于模型的描述能力，计算结果中存在负功率散射分量，同时该方法还有过高估计体散射分量功率的问题，特别是对于高分辨率的极化SAR影像，上述局限性更为明显。具体表现为，对于高分辨率全极化SAR影像，城市地区人工建筑物通常与雷达的方位向形成一个的夹角，造成这些人工建筑物定向角不为0度，进而导致该类地物在做Yamaguchi四分量分解时，体散射分量过高，与实际地物散射类型不符。国际上常使用定向角补偿的方法，把目标的散射矩阵的定向角旋转到0度后，再做极化分解；但随着定向角的增大，定向角补偿的效果降低，建筑物的定向角达到45度时，定向角补偿效果最差，因而分解结果中，地物体散射分量功率仍然占主导，进而无法获取该类地物的正确散射机制，给地物类型鉴别带来困难。

发明内容

有鉴于上述分析中原有Yamaguchi分解方法的局限性，本发明的目的在于避免上述背景技术中的不足而提供一种面向极化SAR散射特性的地物偶次散射模型，并基于该模型构造了一种极化SAR四分量分解方法。本发明方法实现简单，能够有效避免原有方法带来的不足，对于高分辨率全极化SAR影像能够获得满意的影像分解结果。

本发明中所述的方法是这样实现的，一种极化合成孔径雷达影像四分量分解的方法，用于提取极化合成孔径雷达影像中地物的散射机制，其特征在于包括步骤：

①对极化合成孔径雷达影像进行筛选，对于极化合成孔径雷达影像中定向角为0度的建筑物所对应的像素，转入步骤⑤；对于极化合成孔径雷达影像中定向角不为0度的建筑物所对应的像素，转入步骤②；

②根据定向角不为0度的建筑物的散射特点，建立能够描述其散射特性的广义偶次散射模型；

③基于广义偶次散射模型，针对定向角不为0度的建筑物所对应的像素，构造广义偶次散射模型的四分量分解方程；

④提出求解广义偶次散射模型的四分量分解方程的非负特征值分解方法，求取散射分量；所述的散射分量包括地物偶次散射机制的功率、体散射机制的功率、表面散射机制的功率和螺旋体散射机制的功率；

⑤使用带有定向角补偿的Yamaguchi分解方法求取其散射分量；所述的散射分量包括地物偶次散射机制的功率、体散射机制的功率、表面散射机制的功率和螺旋体散射机制的功率；

完成极化合成孔径雷达影像的四分量分解。

其中，步骤①中对于极化合成孔径雷达影像中定向角不为0度的建筑物所对应的像素的筛选具体包括以下步骤：

（101）使用H/α/A极化合成孔径雷达影像分类方法，并选取分类结果中的第1，第2，第5，和第6区域生成像素集合A，这样集合A中的像素对应着两类地物，即以体散射机制为主的茂密植被区域和交叉极化通道功率较大的定向角不为0度的建筑物；

（102）如果极化合成孔径雷达影像为全极化合成孔径雷达数据，利用极化合成技术，在电磁波的定向角

和椭圆率角χ张成的二维空间中，将共极化通道的极化响应特征最大化，在像素集合A中选取电磁波的定向角

位于5度至85度之间的像素作为定向角不为0度的建筑物所对应的像素；

如果极化合成孔径雷达影像为极化干涉数据，利用植被时间去相干的特性，在像素集合A中提取出交叉极化通道的相干性ρ大于0.8的像素，就得到定向角不为0度建筑物所对应的像素。

其中，步骤②中的广义偶次散射模型T_d(ψ)用下式表示：

$T_d\left(\mathrm{\ψ}\right)=R\left(\mathrm{\ψ}\right)T_d{R}^H\left(\mathrm{\ψ}\right)=\frac{1}{1+{\left|\mathrm{\α}\right|}^2}\left[\begin{array}{ccc}{\left|\mathrm{\α}\right|}^2& \mathrm{\α}\cos \left(2\mathrm{\ψ}\right)& \mathrm{\α}\sin \left(2\mathrm{\ψ}\right)\\ {\mathrm{\α}}^{*}\cos \left(2\mathrm{\ψ}\right)& {\cos }^2\left(2\mathrm{\ψ}\right)& \sin \left(4\mathrm{\ψ}\right)/2\\ {\mathrm{\α}}^{*}\sin \left(2\mathrm{\ψ}\right)& \sin \left(4\mathrm{\ψ}\right)/2& {\sin }^2\left(2\mathrm{\ψ}\right)\end{array}\right]$

$T_d=\frac{1}{1+{\left|\mathrm{\α}\right|}^2}\left[\begin{array}{ccc}{\left|\mathrm{\α}\right|}^2& \mathrm{\α}& 0\\ {\mathrm{\α}}^{*}& 1& 0\\ 0& 0& 0\end{array}\right]$

$R\left(\mathrm{\ψ}\right)=\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos \left(2\mathrm{\ψ}\right)& -\sin \left(2\mathrm{\ψ}\right)\\ 0& \sin \left(2\mathrm{\ψ}\right)& \cos \left(2\mathrm{\ψ}\right)\end{array}\right]$

其中，T_d是Yamaguchi分解中使用相干矩阵的形式表示的偶次散射模型，R(ψ)是相干矩阵T_d所对应的定向角旋转矩阵，旋转角度为ψ；H表示矩阵的共轭转置；复数α是表示与材质的反射率与介电性质有关的参数；α的限定条件为|α|＜1。

其中，步骤③中的基于广义偶次散射模型的四分量分解方程为：

T＝P_d·T_d(ψ)+P_v·T_v+P_s·T_s+P_c·T_h

其中，T是用相干矩阵形式表示的地物极化散射矩阵，T_d是表示广义偶次散射模型的相干矩阵，P_d表示地物偶次散射机制的功率；T_v是用相干矩阵形式表示的体散射模型，P_v表示体散射机制的功率；T_s是用相干矩阵形式表示的表面散射模型，P_s表示表面散射机制的功率；T_h是用相干矩阵形式表示的螺旋体散射模型，P_c表示螺旋体散射机制的功率。

其中，步骤④中的散射分量的求解方法包括以下步骤：

⑴极化散射矩阵归一化；

⑵根据T′₂₃的虚部确定螺旋体散射功率P′_c的值；

⑶基于NNED的二分迭代搜索方法求取偶次散射功率P′_d的值；

⑷根据剩余矩阵T′_rm2求取体散射功率P′_v的值为：

P′_v＝3min(T′_rm2(1,1),T′_rm2(3,3))

⑸根据极化分解的功率恒定准则，计算表面散射功率P′_s为：

P′_s＝span(T′_rm2)-P′_v

⑹根据极化总功率，求取的散射分量为：

$\left\{\begin{array}{c}{P}_d=P_d^{\′}\·\mathrm{span}\left(T\right)\\ P_v=P_v^{\′}\·\mathrm{span}\left(T\right)\\ P_s=P_s^{\′}\·\mathrm{span}\left(T\right)\\ P_c=P_c^{\′}\·\mathrm{span}\left(T\right)\end{array}\right.$

完成散射分量的求解。

本发明相比背景技术具有如下优点：

（i）分析现有的Yamaguchi分解方法无法正确获取定向角不为0度的地物的散射信息的缺陷，提出一种适合于描述该类地物的广义偶次散射模型。

（ii）基于本发明提出的广义偶次散射模型，构造一种极化SAR四分量分解方法，首次在国际上解决了上述存在的问题。

（iii）针对广义偶次散射模型，提出了基于非负特征值分解的定向角不为0度的地物所对应像素的筛选方法，便于有针对性地应用广义偶次散射模型。

（iv）本发明提出的方法计算简单，求解结果具有唯一性，使用了二分搜索法搜寻解空间，计算效率高。

（v）本发明所提出的方法，在国际上首次实现了在极化分解方法中，完全利用散射矩阵各个元素的信息，避免了地物散射信息的损失。

附图说明

图1是本发明的数据处理流程图；

图2是像素筛选方法流程图。

具体实施方式

下面结合附图1-2对本发明做进一步的描述：

首先，按照惯例，需要对数据做一定的预处理，即对以相干矩阵形式表示的原始影像进行滤波，以降低SAR相干斑噪声对分解结果的影响，提高SAR数据的信噪比，这里选用降噪和边缘保持效果均理想的LeeSigma滤波方法。

接下来，实施一种极化合成孔径雷达影像四分量分解的方法，如图1所示，用于提取极化合成孔径雷达影像中地物的散射机制，其特征在于包括步骤：

①对极化合成孔径雷达影像进行筛选，对于极化合成孔径雷达影像中定向角为0度的建筑物所对应的像素，转入步骤⑤；对于极化合成孔径雷达影像中定向角不为0度的建筑物所对应的像素，转入步骤②。

对于定向角不为0度的建筑物，使用现有的分解方法，其体散射功率与偶次散射功率不能达到合理比例，该方法的流程如图2所示。根据数据源的不同，可以分成两种情况考虑。首先，使用处理极化SAR影像常用的H/α/A分类方法，并选取分类结果中的第1，第2，第5，和第6区域生成像素集合A，这样集合A中的像素对应着两类地物，即以体散射机制为主的茂密植被区域，和交叉极化通道功率较大的定向角不为0度的建筑物。基于这个结果，如果使用的是全极化数据，则使用如下方法进行处理。利用极化合成技术，即在电磁波的定向角

和椭圆率角χ张成的二维空间中最大化共极化通道的极化响应特征，当共极化通道的极化响应特征最大化时，所对应的电磁波的定向角

的值，在植被区域通常接近于0度或者90度，而在其他区域没有这种现象，因此可以在像素集合A中选取那些电磁波的定向角

位于5度至85度之间的像素点作为目标像素点。该方法具体的计算过程如下：

需要在电磁波的定向角

和椭圆率角χ张成的二维空间中进行逐点搜索，以获取共极化通道能量的最大值。其中，归一化的共极化通道的能量可以使用Stokes矩阵表示成如下式所示的形式。

σ＝s^TM_s

其中，T表示向量的转制，σ表示归一化后的共极化通道的后向散射强度，s为Stokes矢量，如下式所示。

M表示Stokes矩阵，如下式所示。

$M=\left[\begin{array}{cccc}{m}_{11}& m_{12}& m_{13}& m_{14}\\ m_{12}& m_{22}& m_{23}& m_{24}\\ m_{13}& m_{23}& m_{33}& m_{34}\\ m_{14}& m_{24}& m_{34}& m_{44}\end{array}\right]$

如果记b如下式所示，

$b=\left[\begin{array}{c}{m}_{12}\\ m_{13}\\ m_{14}\end{array}\right]$

记C如下式所示，

$C=\left[\begin{array}{ccc}{m}_{22}& m_{23}& m_{24}\\ m_{23}& m_{33}& m_{34}\\ m_{24}& m_{34}& m_{44}\end{array}\right]$

记a为如下式所示的形式

则σ可以表示为下式所示的形式，并且a满足a^Ta＝1。

σ＝m₁₁+2b^Ta+a^TCa

上式可以认为是共极化响应通道的极化响应特征的目标函数，它是未知参数

和χ的函数，现在需要将上式给出的目标函数最大化，以求取对应的未知参数

和χ对应的值。现在使用最速上升法进行求解，首先对上式求关于a的梯度，如下式所示。

$\frac{\∂\mathrm{\σ}}{\∂a}=2(b+\mathrm{Ca})$

接下来可以按照如下步骤进行求解。

步骤201：设置向量a的初始值，如下式所示。

a=[  -1  0  0  ]^T

步骤202：计算梯度值

步骤203：计算a'的值，如下式所示。

$a^{\′}=a+\mathrm{\Δ}\·\frac{\∂\mathrm{\σ}}{\∂a}$

步骤204：将a′进行归一化处理，如下式所示。

$a^{\left(\mathrm{new}\right)}=\frac{a^{\′}}{a^{\′T}{a}^{\′}}$

步骤205：使用下式把向量a更新为a^(new)，并回到步骤202。

a＝a^(new)

步骤206：当满足下式时，停止上述迭代过程，其中ε为一个给定的小正数，通常可以取为0.0001。

Max(|a^(new)-a|)＜ε

从而，使用上述的步骤201至步骤206，可以解算当共极化通道能量最大化时，所对应的电磁波的定向角

和椭圆率角χ的值。

如果使用的是极化干涉数据，问题就变得简单，仅需要在像素集合A中找到那些交叉极化通道的相干性ρ大于0.8的像素即可，这是因为像素集合A中由于植被区域的去相干，因此相干性很低，而建筑物属于稳定目标，因而相干性很高。

②根据定向角不为0度的建筑物的散射特点，建立能够描述其散射特性的广义偶次散射模型；

$T_d\left(\mathrm{\ψ}\right)=R\left(\mathrm{\ψ}\right)T_d{R}^H\left(\mathrm{\ψ}\right)=\frac{1}{1+{\left|\mathrm{\α}\right|}^2}\left[\begin{array}{ccc}{\left|\mathrm{\α}\right|}^2& \mathrm{\α}\cos \left(2\mathrm{\ψ}\right)& \mathrm{\α}\sin \left(2\mathrm{\ψ}\right)\\ {\mathrm{\α}}^{*}\cos \left(2\mathrm{\ψ}\right)& {\cos }^2\left(2\mathrm{\ψ}\right)& \sin \left(4\mathrm{\ψ}\right)/2\\ {\mathrm{\α}}^{*}\sin \left(2\mathrm{\ψ}\right)& \sin \left(4\mathrm{\ψ}\right)/2& {\sin }^2\left(2\mathrm{\ψ}\right)\end{array}\right]$

$T_d=\frac{1}{1+{\left|\mathrm{\α}\right|}^2}\left[\begin{array}{ccc}{\left|\mathrm{\α}\right|}^2& \mathrm{\α}& 0\\ {\mathrm{\α}}^{*}& 1& 0\\ 0& 0& 0\end{array}\right]$

$R\left(\mathrm{\ψ}\right)=\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos \left(2\mathrm{\ψ}\right)& -\sin \left(2\mathrm{\ψ}\right)\\ 0& \sin \left(2\mathrm{\ψ}\right)& \cos \left(2\mathrm{\ψ}\right)\end{array}\right]$

其中，T_d是Yamaguchi分解中使用相干矩阵的形式表示的偶次散射模型，R(ψ)是相干矩阵T_d所对应的定向角旋转矩阵，旋转角度为ψ；H表示矩阵的共轭转置；复数α是表示与材质的反射率与介电性质有关的参数；α的限定条件为|α|＜1。

③基于广义偶次散射模型，针对定向角不为0度的建筑物所对应的像素，构造广义偶次散射模型的四分量分解方程；

本发明构造的四分量分解方法可以表示如下：

$\begin{array}{c}T=P_d\·T_d\left(\mathrm{\ψ}\right)+P_v\·T_v+P_s\·T_s+P_c\·T_h\\ =P_d\·T_d\left(\mathrm{\ψ}\right)+\frac{P_v}{3}\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right]+\frac{P_s}{1+{\left|\mathrm{\β}\right|}^2}\left[\begin{array}{ccc}1& \mathrm{\β}& 0\\ {\mathrm{\β}}^{*}& {\left|\mathrm{\β}\right|}^2& 0\\ 0& 0& 0\end{array}\right]+\frac{P_c}{2}\left[\begin{array}{ccc}0& 0& 0\\ 0& 1& \stackrel{-}{+j}\\ 0& \±j& 1\end{array}\right]\end{array}$

如上式所示，将目标的相干矩阵T分解为四种基本散射机制的加权和，权值即为相应的散射机制的功率。四种散射机制分别为偶次散射，体散射，表面散射，以及螺旋体散射。其中，表面散射模型使用Bragg表面散射模型T_s，参数β是与目标的材质与介电性质有关的复数，并且满足|β|＜1；螺旋体散射模型T_h与Yamaguchi分解中使用的体散射模型一致；由于该类特殊地物仍然属于城市区域的散射目标，因此体散射模型选用熵最大的体散射模型T_v；最后，本发明中提出的广义偶次散射模型。从而，可以获得如下所示的方程组。

$T_{11}=\frac{P_d\·{\left|\mathrm{\α}\right|}^2}{1+{\left|\mathrm{\α}\right|}^2}+\frac{P_s}{1+{\left|\mathrm{\β}\right|}^2}+\frac{P_v}{3}$

$T_{22}=\frac{P_s\·{\left|\mathrm{\β}\right|}^2}{1+{\left|\mathrm{\β}\right|}^2}+\frac{P_d\·{\cos }^2\left(2\mathrm{\ψ}\right)}{1+{\left|\mathrm{\α}\right|}^2}+\frac{P_v}{3}+\frac{P_c}{2}$

$T_{33}=\frac{P_d\·{\sin }^2\left(2\mathrm{\ψ}\right)}{1+{\left|\mathrm{\α}\right|}^2}+\frac{P_v}{3}+\frac{P_c}{2}$

$T_{12}=\frac{P_s\·\mathrm{\β}}{1+{\left|\mathrm{\β}\right|}^2}+\frac{P_d\·\mathrm{\α}}{1+{\left|\mathrm{\α}\right|}^2}\cos \left(2\mathrm{\ψ}\right)$

$T_{13}=\frac{P_d\·\mathrm{\α}}{1+{\left|\mathrm{\α}\right|}^2}\sin \left(2\mathrm{\ψ}\right)$

$T_{23}=\frac{P_d}{1+{\left|\mathrm{\α}\right|}^2}\sin \left(2\mathrm{\ψ}\right)\cos \left(2\mathrm{\ψ}\right)+\frac{P_c}{2}\left(\stackrel{-}{+}j\right)$

从上述方程组可知，本发明的四分量分解方法利用了目标的相干矩阵T的全部信息，特别是元素T₁₃和T₂₃的信息能够被利用起来，这两个元素的信息在国际上现有的分解方法中均被忽略掉。

④提出求解广义偶次散射模型的四分量分解方程的非负特征值分解方法，求取散射分量；所述的散射分量包括地物偶次散射机制的功率、体散射机制的功率、表面散射机制的功率和螺旋体散射机制的功率；

由于引入了定向角参数ψ，因此无法直接计算出各种散射机制的功率值。基于上述方程组，在图1中给出了每种散射机制功率的求解方法。由于求解过程中需要迭代搜索，因此在分解最初阶段，对相干矩阵T进行功率归一化处理，获得归一化后的相干矩阵T'，后续的求解过程均基于T'进行。

首先，由最后一个方程可知，根据T′₂₃的虚部确定螺旋体散射功率P′_c的大小，如图1所示，其中，Im表示取复数的虚部，|·|是取绝对值运算。从T'中减去螺旋体散射成分后获得剩余矩阵T′_rm1，接下来主要问题是要确定偶次散射功率P′_d的大小。

在这里提出一种基于非负特征值分解（Nonnegative EigenvalueDecomposition，NNED）的迭代搜索方法，用于求解偶次散射功率P′_d，并避免求解过程中出现负功率值。如图1所示，首先对目标的相干矩阵T进行功率归一化，从而可以确定偶次散射功率P′_d的搜索范围为[0,1]，并将定向角参数ψ的搜索范围设定为[ψ₁,ψ₂]。从而，偶次散射功率P′_d具体的迭代求解过程如下。

步骤301：在[0,1]内，由小到大给定P′_d一个固定的值，基于这个固定值，在区间[ψ₁,ψ₂]内对定向角参数ψ进行循环，即在[ψ₁,ψ₂]内给定一个ψ的值，计算对应的剩余矩阵T′_rm2。

步骤302：对T′_rm2进行特征值分解并找到三个特征值中最小的特征值λ₃。

步骤303：由于T′_rm2是Hermitian矩阵，因此其所有特征值均非负。如果λ₃大于0，则返回步骤301，尝试下一个P′_d的值。

步骤304：当P′_d上升到某个值时，总能够使得无论ψ在[ψ₁,ψ₂]内取何值，剩余矩阵T′_rm2的最小的那个特征值λ₃均小于等于0，此时终止迭代。

步骤305：记下终止迭代时，P′_d所对应的值，即为需要求解的最优值，此时定向角参数ψ的值也与地物的实际定向角最接近，与此同时，剩余矩阵T′_rm2的功率也达到了唯一的全局最小值。

进一步分析发现，可以对步骤301中P′_d的取值方法进行优化，即在搜索过程中，随着偶次散射功率P′_d的逐渐增大，在[ψ₁,ψ₂]内所能够获得的T′_rm2的最小特征值λ₃的最大值将逐渐减小，因此在搜索P′_d的最优值过程中可以使用二分搜索法来降低计算量。例如，如果将P′_d的搜索步长设定为0.001，则使用二分搜索法后，在最坏情况下经过10次循环就能搜索到P′_d的最优值，如果按照穷举搜索，在最坏情况下，需要1000次循环。

如图1所示，在确定了偶次散射功率P′_d的最优值之后，需要确定其他类型的散射分量的功率，如步骤307～步骤309所述。

步骤306：从剩余矩阵T′_rm1中减去偶次散射分量，从而获得剩余矩阵T′_rm2。

步骤307：现在剩余矩阵T′_rm2的功率包含了其余两种基本散射成分，即P′_v和P′_s。其中，体散射功率P′_v的计算方法如下式所示。

P′_v＝3min(T′_rm2(1,1),T′_rm2(3,3))

步骤308：根据极化分解的功率恒定准则，认为T′_rm2减去体散射功率P′_v后的剩余功率是表面散射功率P′_s，如下式所示。

P′_s＝span(T′_rm2)-P′_v

步骤309：将所求出的各个散射分量的功率乘以总功率才能获得散射功率的绝对值，如下式所示。

$\left\{\begin{array}{c}{P}_d=P_d^{\′}\·\mathrm{span}\left(T\right)\\ P_v=P_v^{\′}\·\mathrm{span}\left(T\right)\\ P_s=P_s^{\′}\·\mathrm{span}\left(T\right)\\ P_c=P_c^{\′}\·\mathrm{span}\left(T\right)\end{array}\right.$

接下来给出定向角参数ψ的搜索区间[ψ₁,ψ₂]的计算方法。基于上述分析，可以获得如下的表达式。

$\frac{P_d^{\′}\sin \left(4\mathrm{\ψ}\right)}{2(1+{\left|\mathrm{\α}\right|}^2)}=\mathrm{Re}\left(T_{23}^{\′}\right)$

$\frac{\mathrm{\α}{P}_d^{\′}\sin \left(2\mathrm{\ψ}\right)}{1+{\left|\mathrm{\α}\right|}^2}=T_{13}^{\′}$

$\mathrm{\α}=\frac{T_{13}^{\′}\cos \left(2\mathrm{\ψ}\right)}{\mathrm{Re}\left(T_{23}^{\′}\right)}$

考虑到|α|＜1，可以进一步获得：

-τ＜cos(2ψ)＜τ

$\mathrm{\τ}=\frac{\left|\mathrm{Re}\left(T_{23}^{\′}\right)\right|}{\left|T_{13}^{\′}\right|}$

如果τ≥1，则区间[ψ₁,ψ₂]取为固定区间[0，π]；如果τ＜1，区间[ψ₁,ψ₂]可以由下面两个公式给出。

$\mathrm{\ψ}\∈[b,\frac{\mathrm{\π}}{2}-b]\∪[\frac{\mathrm{\π}}{2}+b,\mathrm{\π}-b]$

$b=\frac{{\cos }^{-1}\left(\mathrm{\τ}\right)}{2}$

进一步可以观察到，如果ψ在上式中的其中任何一个区间取值，解算出的广义偶次散射模型的表达式均能与另一个区间内对应，因此，可以将定向角参数ψ的搜索范围进一步缩减，以降低计算量，如下式所示。

$[{\mathrm{\&psi;}}_1,{\mathrm{\&psi;}}_2]=\left\{\begin{array}{cc}[0,\frac{\mathrm{\&pi;}}{2}]& \mathrm{\&tau;}\&GreaterEqual;1\\ [b,\frac{\mathrm{\&pi;}}{2}-b]& \mathrm{\&tau;}<1\end{array}\right.$

⑤使用带有定向角补偿的Yamaguchi分解方法求取其散射分量；所述的散射分量包括地物偶次散射机制的功率、体散射机制的功率、表面散射机制的功率和螺旋体散射机制的功率；

⑥完成极化合成孔径雷达影像的四分量分解。

Claims (5)

1.一种极化合成孔径雷达影像四分量分解的方法，用于提取极化合成孔径雷达影像中地物的散射机制，其特征在于包括步骤：

①对极化合成孔径雷达影像进行筛选，对于极化合成孔径雷达影像中定向角为0度的建筑物所对应的像素，转入步骤⑤；对于极化合成孔径雷达影像中定向角不为0度的建筑物所对应的像素，转入步骤②；

②根据定向角不为0度的建筑物的散射特点，建立能够描述其散射特性的广义偶次散射模型；

③基于广义偶次散射模型，针对定向角不为0度的建筑物所对应的像素，构造广义偶次散射模型的四分量分解方程；

④提出求解广义偶次散射模型的四分量分解方程的非负特征值分解方法，求取散射分量；所述的散射分量包括地物偶次散射机制的功率、体散射机制的功率、表面散射机制的功率和螺旋体散射机制的功率；

⑤使用带有定向角补偿的Yamaguchi分解方法求取其散射分量；所述的散射分量包括地物偶次散射机制的功率、体散射机制的功率、表面散射机制的功率和螺旋体散射机制的功率；

完成极化合成孔径雷达影像的四分量分解。

2.根据权利要求1所述的一种极化合成孔径雷达影像四分量分解的方法，其特征在于：步骤①中对于极化合成孔径雷达影像中定向角不为0度的建筑物所对应的像素的筛选具体包括以下步骤：

（101）使用H/α/A极化合成孔径雷达影像分类方法，并选取分类结果中的第1，第2，第5，和第6区域生成像素集合A，这样集合A中的像素对应着两类地物，即以体散射机制为主的茂密植被区域和交叉极化通道功率较大的定向角不为0度的建筑物；

（102）如果极化合成孔径雷达影像为全极化合成孔径雷达数据，利用极化合成技术，在电磁波的定向角

和椭圆率角χ张成的二维空间中，将共极化通道的极化响应特征最大化，在像素集合A中选取电磁波的定向角

位于5度至85度之间的像素作为定向角不为0度的建筑物所对应的像素；

如果极化合成孔径雷达影像为极化干涉数据，利用植被时间去相干的特性，在像素集合A中提取出交叉极化通道的相干性ρ大于0.8的像素，就得到定向角不为0度建筑物所对应的像素。

3.根据权利要求1所述的一种极化合成孔径雷达影像四分量分解的方法，其特征在于步骤②中的广义偶次散射模型T_d(ψ)用下式表示：

$T_d\left(\mathrm{\&psi;}\right)=R\left(\mathrm{\&psi;}\right)T_d{R}^H\left(\mathrm{\&psi;}\right)=\frac{1}{1+{\left|\mathrm{\&alpha;}\right|}^2}\left[\begin{array}{ccc}{\left|\mathrm{\&alpha;}\right|}^2& \mathrm{\&alpha;}\cos \left(2\mathrm{\&psi;}\right)& \mathrm{\&alpha;}\sin \left(2\mathrm{\&psi;}\right)\\ {\mathrm{\&alpha;}}^{*}\cos \left(2\mathrm{\&psi;}\right)& {\cos }^2\left(2\mathrm{\&psi;}\right)& \sin \left(4\mathrm{\&psi;}\right)/2\\ {\mathrm{\&alpha;}}^{*}\sin \left(2\mathrm{\&psi;}\right)& \sin \left(4\mathrm{\&psi;}\right)/2& {\sin }^2\left(2\mathrm{\&psi;}\right)\end{array}\right]$

$T_d=\frac{1}{1+{\left|\mathrm{\&alpha;}\right|}^2}\left[\begin{array}{ccc}{\left|\mathrm{\&alpha;}\right|}^2& \mathrm{\&alpha;}& 0\\ {\mathrm{\&alpha;}}^{*}& 1& 0\\ 0& 0& 0\end{array}\right]$

$R\left(\mathrm{\&psi;}\right)=\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos \left(2\mathrm{\&psi;}\right)& -\sin \left(2\mathrm{\&psi;}\right)\\ 0& \sin \left(2\mathrm{\&psi;}\right)& \cos \left(2\mathrm{\&psi;}\right)\end{array}\right]$

其中，T_d是Yamaguchi分解中使用相干矩阵的形式表示的偶次散射模型，R(ψ)是相干矩阵T_d所对应的定向角旋转矩阵，旋转角度为ψ；H表示矩阵的共轭转置；复数α是表示与材质的反射率与介电性质有关的参数；α的限定条件为|α|＜1。

4.根据权利要求3所述的一种极化合成孔径雷达影像四分量分解的方法，其特征在于：步骤③中的基于广义偶次散射模型的四分量分解方程为：

T＝P_d·T_d(ψ)+P_v·T_v+P_s·T_s+P_c·T_h

其中，T是用相干矩阵形式表示的地物极化散射矩阵，T_d是表示广义偶次散射模型的相干矩阵，P_d表示地物偶次散射机制的功率；T_v是用相干矩阵形式表示的体散射模型，P_v表示体散射机制的功率；T_s是用相干矩阵形式表示的表面散射模型，P_s表示表面散射机制的功率；T_h是用相干矩阵形式表示的螺旋体散射模型，P_c表示螺旋体散射机制的功率。

5.根据权利要求4所述的一种极化合成孔径雷达影像四分量分解的方法，其特征在于：步骤④中的散射分量的求解方法包括以下步骤：

⑴极化散射矩阵归一化；

⑵根据T′₂₃的虚部确定螺旋体散射功率P′_c的值；

⑶基于NNED的二分迭代搜索方法求取偶次散射功率P′_d的值；

⑷根据剩余矩阵T′_rm2求取体散射功率P′_v的值为：

P′_v＝3min(T′_rm2(1,1),T′_rm2(3,3))

⑸根据极化分解的功率恒定准则，计算表面散射功率P′_s为：

P′_s＝span(T′_vm2)-P′_v

⑹根据极化总功率，求取的散射分量为：

$\left\{\begin{array}{c}{P}_d=P_d^{\&prime;}\&CenterDot;\mathrm{span}\left(T\right)\\ P_v=P_v^{\&prime;}\&CenterDot;\mathrm{span}\left(T\right)\\ P_s=P_s^{\&prime;}\&CenterDot;\mathrm{span}\left(T\right)\\ P_c=P_c^{\&prime;}\&CenterDot;\mathrm{span}\left(T\right)\end{array}\right.$

完成散射分量的求解。

Images