极化SAR图像分解方法及存储介质与流程

本发明涉及微波遥感技术领域，尤其涉及一种极化sar图像分解方法及存储介质。

背景技术：

随着微波遥感技术的发展，极化合成孔径雷达(极化sar)已成为微波遥感领域的一个重要研究方向，对于极化sar的研究逐渐由单极化、单波段发展到多极化、多波段，相继出现了星载、机载等全极化sar系统。极化sar通过发射具有水平极化h或者垂直极化v的电磁波对地面进行照射，并接收该地面区域具有水平极化h和垂直极化v的散射回波，从而获得观测区域的极化散射矩阵。极化散射矩阵中包含了观测目标的详细信息，包括表面粗糙度、方位取向、对称性、几何形状结构等电磁散射特性。

极化目标分解理论最早由huynen提出，是一种对极化sar图像数据进行目标散射特性提取的主要方法。目前，极化sar图像数据分解方法主要有描述纯目标的相干目标分解和描述分布式目标的非相干目标分解两大类。相干目标分解通过对纯目标的极化散射矩阵进行加法或乘积分解，提取出一系列能够较好描述目标的极化特征，该分解方法主要包括pauli分解、cameron分解、krogager分解等方法。非相干目标分解主要分为基于mueller矩阵的非相干目标分解、基于协方差矩阵c或相干矩阵t的特征值非相干目标分解和基于散射模型的非相干目标分解，该分解方法主要包括：cloude分解、freeman分解、yamaguchi分解等。

目标分解是极化sar的一个重要部分，它是目标分类、遥感应用的前提基础，可以帮助分析散射目标的物理机制，有利于极化sar数据的地物分类、目标识别等。1998年，freeman和durden首先提出了三成分的极化目标分解方法，根据获取的极化sar数据确定观测区域中占主导地位的散射机制，将原始协方差矩阵c或相干矩阵t分解为三个散射成分(表面散射、偶次散射和体散射)的线性组合，该方法假设偶极子的旋转是均匀的以及反射对称条件的成立。2005年，yamaguchi等，在freeman分解基础上提出了四成分的极化目标分解方法，增加了螺旋体散射成分。2009年，cloude提出了一种混合freeman/eigenvalue分解方法，随后将其分解思想引入四分量模型分解中，非相干目标分解得到了迅速发展和广泛应用，但这些方法大部分是基于cloude特征值分解和freeman三成分分解，其应用都有一定的局限性。

现有的基于散射模型分解方法最大的缺陷就是负功率问题。在雷达系统中，给定发射机和接收机，那么与雷达截面积成正比的接收到的散射功率永远不可能是负的，所以在任何情况下散射功率都不可能是负值。负散射功率产生的原因是由于在分解过程中对其中一个散射分量的估计过高，由于所有散射分量的贡献之和等于接收到的相干矩阵的总功率，即总功率始终保持不变，如果某一散射分量的功率在某一像素处被高估，那么在该像素处剩余的散射分量就有可能存在负的贡献。另外，交叉极化功率的存在会降低目标分解的准确性，如二次面散射与体散射耦合产生的交叉极化功率，表面散射与体散射耦合产生的交叉极化功率；此外，传统的体散射模型不能随极化sar数据的不同而进行自适应的调整，不能满足林地中hh分量和vv分量比例的连续变化需求。

技术实现要素：

鉴于现有技术中存在的上述技术问题，本发明实施例所要解决的技术问题是提供一种极化sar图像分解方法，其能够降低体散射功率的过高估计，避免偶次散射和表面散射产生负功率的问题，以解决现有基于散射模型分解中的负功率问题。

为了解决上述技术问题，本发明实施例采用了如下技术方案：

一种极化sar图像分解方法，包括如下步骤：

s1：获取全极化sar图像数据的极化相干矩阵t；

s2：判断所述全极化sar图像数据中每个像素点处的主导散射机制，并根据所述主导散射机制对所述像素点处的所述相干矩阵t进行修正，获得修正后的相干矩阵t′；

s3：对经修正后的所述相干矩阵t′进行四分量目标分解，先分别求解螺旋体散射分量和体散射分量的功率值，然后，对分解出螺旋体散射分量和体散射分量后的剩余矩阵tr进行特征值分解，求解表面散射分量和偶次散射分量的功率值，完成四分量目标分解。

进一步地，所述相干矩阵t为3阶相干矩阵t3，每个像素点包含9个元素：

其中，t11，t12，t13，t21，t22，t23，t31，t32，t33分别表示相干矩阵t3中的9个元素。

进一步地，所述像素点处的主导散射机制的判断方法包括：

令q＝t11-t22，计算像素点处的q值，根据q值的大小判断每个所述像素点处的主导散射机制，

进一步地，步骤s2中，根据所述主导散射机制对所述相干矩阵t进行修正包括：对所述相干矩阵t进行方位角补偿和相位角旋转。

进一步地，当表面散射主导时，所述方位角补偿和所述相位角旋转包括：

利用公式(2)对所述相干矩阵t3进行方位角补偿，得到相干矩阵t3_1，

式中，[]^h表示矩阵的共轭转置，

其中，re(t13)表示取复数的实部，tan^-1表示反正切函数；

利用公式(3)对经所述方位角补偿后的所述相干矩阵t3_1进行相位角旋转，得到相干矩阵t3_2，

式中，[]^h表示矩阵的共轭转置，

其中，im(t13)表示取复数的虚部，tan^-1表示反正切函数，t11_1与t33_1是t3_1中的元素；

当偶次散射主导时，所述方位角补偿和相位角旋转包括：

利用公式(4)对所述相干矩阵t3进行方位角补偿，得到相干矩阵t3_3，

式中，[]^h表示矩阵的共轭转置，r3如下所示：

其中，re(t23)表示取实部，tan^-1表示反正切函数；

利用公式(5)对经所述方位角补偿后的所述相干矩阵t3_3进行相位角旋转，得到相干矩阵t3_4，

式中，[]^h表示矩阵的共轭转置，r4如下所示：

其中，im(t23)表示取复数的虚部，tan^-1表示反正切函数；

经所述方位角补偿和所述相位角旋转修正后的所述相干矩阵t′3的最终表达式为：

进一步地，在步骤s3中，根据公式(6)对所述相干矩阵t′3进行四分量目标分解：

式中，ts代表表面散射分量，td代表偶次散射分量，代表体散射分量，tc代表螺旋体散射分量，ms、md、mv和mc分别表示表面散射、偶次散射、体散射和螺旋体散射的贡献值，也即所占的功率值。

进一步地，在步骤s3中，根据螺旋体散射模型的计算公式(7)，计算螺旋体散射的功率值mc：

mc＝2|im(t′23)|(7)

其中，t′23为t′3中的元素，螺旋体散射分量对应的螺旋体散射模型为：

进一步地，在步骤s3中，根据体散射模型的计算公式(9)计算体散射分量

其中，γ＝<|shh|²>/<|svv|²>，表示水平极化分量与垂直极化分量之间的比值；

进一步地，根据体散射模型的功率计算公式(10)，计算体散射分量的功率值mv：

其中，<|shv|²>为经步骤s2修正后的相干矩阵t′3中t′33元素的1/4倍。

进一步地，在步骤s3中，对分解出螺旋体散射分量和体散射分量后的剩余矩阵tr进行特征值分解：

并根据所述步骤s2中确定的像素点处的主导散射机制的判断结果，求解出不同主导散射机制下的表面散射分量的功率值ms和偶次散射分量的功率值md；

其中，根据所述像素点处的主导散射机制的判断结果进行的特征值分解包括：

求解矩阵tr的特征值，求得两个特征值λ1和λ2，令λ1＞λ2，

若步骤s2中判断该像素点处q＞0，即该像素点处表面散射主导时，ms＝λ1，md＝λ2；

若步骤s2中判断该像素点处q≤0，即该像素点处偶次散射机制主导时，ms＝λ2，md＝λ1。

本发明的实施例还提供了一种计算机存储介质，其上存储有计算机可执行指令，所述计算机可执行指令由处理器执行时，实现上述的极化sar图像分解方法。

与现有技术相比，本发明实施例的有益效果在于：本发明实施例提供的极化sar图像分解方法，通过判断像素点的主导机制并进行相应的方位角补偿和相位角旋转操作，大幅度降低了体散射的过度估计现象；分解出螺旋体散射分量和体散射分量后，采用非负特征值法求解表面散射和偶次散射的功率，能够有效的解决表面散射和偶次散射产生的负功率现象，进而得到更加精确的目标散射特性分解图。另外，目标分解时，引入的广义体散射模型不再是固定不变，其能够随着极化sar图像数据的不同进行自适应的调整，满足林地中hh分量和vv分量比例的连续变化，可以在不同地物类型的分类中得到更加准确的分类类别。

附图说明

图1为本发明实施例的极化sar图像分解方法的流程图；

图2为本发明实施例的极化sar图像分解方法的对相干矩阵进行修正的流程图；

图3为本发明实施例的极化sar图像分解方法的求解四分量功率值的流程图。

具体实施方式

为使本领域技术人员更好的理解本发明的技术方案，下面结合附图和具体实施方式对本发明作详细说明。

应理解的是，可以对此处公开的实施例做出各种修改。因此，上述说明书不应该视为限制，而仅是作为实施例的范例。本领域的技术人员将想到在本申请的范围和精神内的其他修改。

包含在说明书中并构成说明书的一部分的附图示出了本申请的实施例，并且与上面给出的对本申请的大致描述以及下面给出的对实施例的详细描述一起用于解释本申请的原理。

通过下面参照附图对给定为非限制性实例的实施例的优选形式的描述，本申请的这些和其它特性将会变得显而易见。

还应当理解，尽管已经参照一些具体实例对本申请进行了描述，但本领域技术人员能够确定地实现本申请的很多其它等效形式，它们具有如权利要求所述的特征并因此都位于借此所限定的保护范围内。

当结合附图时，鉴于以下详细说明，本申请的上述和其他方面、特征和优势将变得更为显而易见。

此后参照附图描述本申请的具体实施例；然而，应当理解，所公开的实施例仅仅是本申请的实例，其可采用多种方式实施。熟知和/或重复的功能和结构并未详细描述以避免不必要或多余的细节使得本申请模糊不清。因此，本文所公开的具体的结构性和功能性细节并非意在限定，而是仅仅作为权利要求的基础和代表性基础用于教导本领域技术人员以实质上任意合适的详细结构多样地使用本申请。

图1为本发明实施例的极化sar图像分解方法的流程图；如图1所示，本申请提供了一种极化sar图像分解方法，包括步骤s1至s3：

步骤s1：获取全极化sar图像数据的极化相干矩阵t；

在步骤s1中，通过对获取的原始全极化sar图像数据进行数据处理获得相干矩阵t。以获取3阶相干矩阵t3为例对极化相干矩阵的获取进行说明，具体获取步骤包括：

步骤s11：获取原始全极化sar图像数据的极化散射矩阵。

全极化sar图像数据可以以极化散射矩阵的形式进行表示，该数据大小为m*n阶，其中，m为数据行数，n为数据的列数；全极化sar图像数据中某一位置处的像素点可以以第i行第j列位置表示，其中i∈(1，m)，i∈(1，n)。

具体地，全极化sar图像数据表示为极化后向散射矩阵s，

其中，shh和svv包含了同极化通道的回波功率，shv和svh包含了交叉极化通道的回波功率，且满足互易定理shv＝svh，h表示波的传播方向为水平(水平极化)，v表示波的传播方向位为垂直(垂直极化)。

步骤s12：从后向散射矩阵s中构建散射目标矢量在复pauli旋转矩阵基集合下得到相应的“三维目标矢量”为：

其中代表矩阵的转置；

步骤s13：求取每个像素点处的三阶相干矩阵t3，每个像素点包含9个元素，其表达式为：

式中，t11，t12，t13，t21，t22，t23，t31，t32，t33表示相干矩阵中的9个元素，[]^*表示矩阵的共轭，[]^h表示矩阵的共轭转置，<·>表示视数平均，表示目标矢量自身共轭转置矢量的外积。

优选地，在步骤s12之前，即在对极化散射矩阵s进行数据处理之前，还包括对极化散射矩阵s进行lee滤波，以去除噪声的影响。另外，由于该滤波能较好的保持边缘特性，因此它还可以减小由于滤波导致的负功率问题。

步骤s2：判断全极化sar图像数据中每个像素点处的主导散射机制，并根据该像素点处的主导散射机制对该像素点处的相干矩阵t进行修正；

图2为本发明实施例的极化sar图像分解方法的对相干矩阵进行修正的流程图，如图2所示，在步骤s2中，对第i行第j列位置处像素点的主导散射机制进行判断的方法包括：

令q＝t11-t22，计算该像素点处的q值，根据q值的大小判断每个像素位置处的主导散射机制，具体判断方法如下式所示：

进一步地，根据像素点处的主导散射机制对相干矩阵进行修正包括方位角补偿和相位角旋转。根据像素点处的主导散射机制的不同，利用不同的公式对该相素点处的相干矩阵依次进行方位角补偿和相位角旋转，对相干矩阵t3进行修正。

具体地，若判断第i行第j列位置处像素点的q＞0，其中i∈(1，m)，i∈(1，n)，则该像素点处表面散射占主导，对该像素点处的相干矩阵t3的方位角补偿和相位角旋转具体包括：

(1)方位角补偿

具体地，对相干矩阵t3进行如下操作，得到方位角补偿后的相干矩阵t3_1，具体操作如公式(2)所示：

其中，[]^h表示矩阵的共轭转置，r1如下所示：

其中，re(t13)表示取复数的实部，tan^-1表示反正切函数；

(2)相位角旋转

具体地，对经上述方位角补偿后的相干矩阵t3_1进行如下操作，得到相位角旋转后的相干矩阵t3_2，具体操作如公式(3)所示：

其中，[]^h表示矩阵的共轭转置，r2如下所示：

其中，im(t13)表示取复数的虚部，tan^-1表示反正切函数，t11_1与t33_1是t3_1中的元素。

若判断第i行第j列位置处像素点的q≤0，其中i∈(1，m)，i∈(1，n)，则该像素点处偶次散射占主导，对该像素点处的相干矩阵t3的方位角补偿和相位角旋转具体包括：

(1)方位角补偿

具体地，对相干矩阵t3进行如下操作，得到方位角补偿后的相干矩阵t3_3，具体操作如公式(4)所示：

其中，[]^h表示矩阵的共轭转置，r3如下所示：

其中，re(t23)表示取实部，tan^-1表示反正切函数；

(2)相位角旋转

对经上述方位角补偿后的相干矩阵t3_3进行如下操作，得到相位角旋转后的相干矩阵t3_4，具体操作如公式(5)所示：

其中，[]^h表示矩阵的共轭转置，r4如下所示：

其中，im(t23)表示取复数的虚部，tan^-1表示反正切函数。

经过步骤s2修正处理之后的相干矩阵t′3最终可以表示为：

其中，

具体地，相干矩阵t3_2和相干矩阵t3_4的表达式分别为：

采用判断法对不同像素点处的主导散射机制进行判断，并通过相位角旋转和方位向补偿对相干矩阵进行修正，可以更加有针对性的消除交叉极化分量hv，降低表面散射和偶次散射的负功率问题。

图3示出了本发明实施例的极化sar图像分解方法的步骤s3中四分量功率的求解流程图。如图3所示，步骤s3包括：对经修正后的极化相干矩阵t′进行四分量目标分解，先分别求解螺旋体散射分量和体散射分量的功率值，然后，对分解出螺旋体散射分量和体散射分量后的剩余矩阵tr进行特征值分解，求解表面散射分量和偶次散射分量的功率值，完成四分量目标分解。

相干矩阵t′3分解的表达式如下：

其中，ts代表表面散射分量，td代表偶次散射分量，代表体散射分量，tc代表螺旋体散射分量，ms、md、mv和mc分别表示表面散射、偶次散射、体散射和螺旋体散射的贡献值，也即所占的功率值。

在步骤s3中，对经修正后的极化相干矩阵t′3进行四分量目标分解的方法具体包括如下步骤：

步骤s31：分别求解螺旋体散射分量的功率值mc和体散射分量的功率值mv。

具体地，根据螺旋体散射模型的计算公式(7)，计算螺旋体散射的功率值mc：

mc＝2|im(t′23)|(7)

其中，t′23为t′3中的元素，螺旋体散射分量对应的螺旋体散射模型为：

根据步骤s2经修正处理之后的相干矩阵t′3，利用体散射模型的计算公式(9)，计算像素点的广义体散射模型分量

其中，γ＝<|shh|²>/<|svv|²>，表示水平极化分量与垂直极化分量之间的比值；

进一步地，根据体散射模型的功率计算公式(10)，计算体散射分量的功率值mv：

其中，<|shv|²>为经步骤s2修正后的相干矩阵t′3中t′33元素的1/4倍。

在步骤s31中，求解螺旋体散射分量的功率值mc和体散射分量的功率值mv的顺序不具体限定。采用上述公式(9)的体散射模型对极化分解中的体散射进行目标分解，在一定程度上缓解了由于体散射功率被高估而导致的表面散射和二面角散射出现负功率问题；另外，上述的体散射模型中的伽马参数γ可以随着修正的相干矩阵t3′的变化而变化，即能够随着极化sar图像数据的不同进行自适应的调整，满足了林地中hh分量和vv分量比例的连续变化，可以在不同地物类型的分类中得到更加准确的分类类别，提高分类精度，使得目标分解结果更接近实际地物目标的散射机理。

步骤s32：对分解出螺旋体散射分量和体散射分量后的剩余矩阵tr进行特征值分解，求解表面散射分量的功率值ms和偶次散射分量的功率值md，完成四分量目标分解。

经步骤s2修正后的相干矩阵t′3减去步骤s31分解出的体散射分量和螺旋体散射分量mctc即为剩余矩阵tr，剩余矩阵tr为表面散射和偶次散射分量之和：

在确定螺旋体散射分量和体散射分量的贡献之后，求解剩余矩阵tr的特征值，并根据步骤s2中确定的像素点处的主导散射机制的判断结果，求解出不同主导散射机制下的表面散射分量的功率值ms和偶次散射分量的功率值md。

具体地，求解矩阵tr的特征值，求得λ1和λ2两个特征值，令λ1＞λ2，

若步骤s2中判断该像素点处q＞0，即表面散射主导时，ms＝λ1，md＝λ2；

若步骤s2中判断该像素点处q≤0，即偶次散射机制主导时，ms＝λ2，md＝λ1。

上述矩阵tr的求解过程中，表面散射模型ts和偶次散射模型td是hybridfreeman分解中已知的模型，在此不再赘述。

步骤s31中的体散射模型公式(9)以及求解体散射分量的功率值mv的计算公式(10)，可以保证剩余矩阵tr是一个秩为2，且保持正半定的矩阵，即可以保证求解的表面散射和偶次散射均为非负数。

在步骤s3计算出极化sar图像中四种散射的功率值后，根据计算得到整个极化sar图像中表面散射、偶次散射、体散射和螺旋体散射的功率ms、md、mv和mc值，可以生成rgb伪彩色合成图。

具体地，用红色r表示二次次散射分量的功率md，绿色g表示体散射分量的功率mv和螺旋体散射分量的功率mc，蓝色b表示表面散射分量的功率ms，使用rgb合成的方法得到一幅极化sar图像分解后的伪彩色图。

由于螺旋体散射属于特殊的体散射，且其功率值较小，因此，合成图像时，可以将螺旋体散射分量的功率mc并入体散射分量的功率mv中。

下面利用德国oberpfaffenhofen地区l波段esar机载数据作为试验数据，并采用matlab软件对本发明的极化四分量分解方法进行验证。验证试验中，分别采用freeman分解方法和混合freeman/eigenvalue分解方法进行对比试验，并计算负功率的像素数量，验证结果如表1所示。

表1三种分解方法的负功率像素比例

由上述验证试验结果可知，采用本发明提供的极化sar图像分解方法，能够有效解决目标分解过程中的负功率问题。

本发明实施例提供的极化sar图像分解方法，通过判断像素点的主导机制并进行相应的方位角补偿和相位角旋转操作，大幅度降低了体散射的过度估计现象；分解出螺旋体散射分量和体散射分量后，采用非负特征值法求解表面散射和偶次散射的功率，能够有效的解决表面散射和偶次散射产生的负功率现象，进而得到更加精确的目标散射特性分解图。另外，目标分解时，引入的广义体散射模型不再是固定不变，其能够随着极化sar图像数据的不同进行自适应的调整，满足林地中hh分量和vv分量比例的连续变化，可以在不同地物类型的分类中得到更加准确的分类类别。

本发明的实施例还提供了一种计算机存储介质，其上存储有计算机可执行指令，所述计算机可执行指令由处理器执行时，实现上述根据本发明的实施例中的极化sar图像分解方法。

在一些实施例中，执行算机可执行指令处理器可以是包括一个以上通用处理设备的处理设备，诸如微处理器、中央处理单元(cpu)、图形处理单元(gpu)等。更具体地，该处理器可以是复杂指令集计算(cisc)微处理器、精简指令集计算(risc)微处理器、超长指令字(vliw)微处理器、运行其他指令集的处理器或运行指令集的组合的处理器。该处理器还可以是一个以上专用处理设备，诸如专用集成电路(asic)、现场可编程门阵列(fpga)、数字信号处理器(dsp)、片上系统(soc)等。

在一些实施例中，计算机可读存储介质可以为存储器，诸如只读存储器(rom)、随机存取存储器(ram)、相变随机存取存储器(pram)、静态随机存取存储器(sram)、动态随机存取存储器(dram)、电可擦除可编程只读存储器(eeprom)、其他类型的随机存取存储器(ram)、闪存盘或其他形式的闪存、缓存、寄存器、静态存储器、光盘只读存储器(cd-rom)、数字通用光盘(dvd)或其他光学存储器、盒式磁带或其他磁存储设备，或被用于储存能够被计算机设备访问的信息或指令的任何其他可能的非暂时性的介质等。

在一些实施例中，计算机可执行指令可以实现为多个程序模块，多个程序模块共同实现上述根据本发明的实施例中的极化sar图像分解方法。

以上实施例仅为本发明的示例性实施例，不用于限制本发明，本发明的保护范围由权利要求书限定。本领域技术人员可以在本发明的实质和保护范围内，对本发明做出各种修改或等同替换，这种修改或等同替换也应视为落在本发明的保护范围内。