基于特征值分解的极化SAR数据多元散射能量分解方法

基于特征值分解的极化sar数据多元散射能量分解方法
技术领域
1.本发明属于极化sar数据分析技术领域，涉及一种基于特征值分解的极化sar数据多元散射能量分解方法。


背景技术：

2.因为极化合成孔径雷达(polarimetric synthetic aperture radar，极化sar)具有全天时、全天候、多波段和高分辨率的成像特点，所以极化sar图像在多个领域具有广泛的应用。近年来各种空载和机载极化sar传感器的研制和发射，需要解释和处理的极化sar数据量急剧增加。由于易于实现和具有强大的物理意义，极化目标分解方法已经成为极化sar数据解释和预处理的主要途径。对于极化sar图像，极化目标分解方法可分为两类：基于散射模型的分解方法和基于特征值的分解方法。优化极化sar图像的极化目标分解算法是提高极化sar图像解译效果的重要途径之一。
3.freeman-durden分解(fd分解)是一种经典的基于散射模型的极化目标分解方法。fd分解将极化sar数据的协方差矩阵或相干矩阵分解为三个主要组成部分：表面散射、二次散射和体散射，也可以称为极化sar数据的三元散射能量分解方法。fd分解有效地描述了极化sar数据的散射特征，并且每个分解得到的散射能量具有明确的物理含义。同时，fd分解也存在一个主要缺点，即存在负的表面散射能量和二次散射能量。原因是fd分解假设极化sar数据符合反射对称理论，即共极化项与交叉极化项不相关，导致相干矩阵[t]中的t
13
和t
23
为零。因此，交叉极化项(t
33
)仅对体散射分量有贡献，导致体散射能量被过分高估，甚至可能大于总能量，最终导致出现负的表面散射能量和二次散射能量。


技术实现要素：

[0004]
本发明的目的是提供一种基于特征值分解的极化sar数据多元散射能量分解方法，该方法不要求极化sar数据符合反射对称理论假设，解决了表面散射能量和二次散射能量的负值问题，提供了多种类的散射能量。
[0005]
本发明所采用的技术方案是，基于特征值分解的极化sar数据多元散射能量分解方法，具体包括如下步骤：
[0006]
步骤1，输入极化sar数据的相干矩阵；
[0007]
步骤2，对每个像素点的极化sar数据相干矩阵进行特征值分解；
[0008]
步骤3，使用七种散射模型分解每个像素点的极化sar数据的相干矩阵；
[0009]
步骤4，计算具有较高熵和较低反熵的像素点的多元散射能量；
[0010]
步骤5，获取实数方位角补偿后的相干矩阵并分解，得到每个像素点的六种散射能量；
[0011]
步骤6，获取两次方位角补偿后的相干矩阵并分解，得到每个像素点的五种散射能量。
[0012]
本发明的特点还在于：
[0013]
步骤1中，极化sar数据的相干矩阵表示为：
[0014][0015]
其中，《t》是一个半正定的赫米特矩阵，k
p
表示单视数据的pauli向量。
[0016]
步骤2中，对每个像素点的极化sar数据相干矩阵进行特征值分解的公式如下所示：
[0017][0018]
其中，λi表示特征值，ki表示对应的特征向量，αi表示散射机制角，λ1≥λ2≥λ3≥0；
[0019]
熵h如公式(3)所示，其中pi＝λi/(λ1+λ2+λ3)，i∈{1，2，3}，反熵a如公式(4)所示：
[0020][0021][0022]
步骤3中，相干矩阵进行七元散射能量分解如下公式(5)所示：
[0023]
《t》＝msts+mdtd+mvtv+mhth+m
md
t
md
+m
cd
t
cd
+m
od
t
od (5)；
[0024]
其中，ts表示表面散射模型，td表示二次散射模型，tv表示体散射模型，th表示螺旋体散射模型，t
md
表示混合偶极子散射模型，t
cd
表示复合偶极子散射模型，t
od
表示偏振偶极子散射模型，ms表示表面散射能量，md表示二次散射能量，mv表示体散射能量，mh表示螺旋体散射能量，m
md
表示混合偶极子散射能量，m
cd
表示复合偶极子散射能量，m
od
表示偏振偶极子散射能量。
[0025]
步骤4中，当一个像素的熵h和反熵a满足关系h-a＞0.4时，ms、md和mv的计算方法如下：
[0026][0027]
步骤5中，实数方位角补偿后的相干矩阵《t(θ)》如公式(7)所示：
[0028][0029]
其中，实数的方位角补偿矩阵r(θ)如公式(8)所示：
[0030][0031]
采用《t(θ)》代替相干矩阵《t》依次执行步骤3和步骤4，由于t
23
(θ)的实部为0，此时，m
md
＝0，因此，《t(θ)》的分解为六元散射能量分解。
[0032]
步骤6中，采用实数的方位角补偿后再进行复数的方位角补偿，复数的方位角补偿矩阵如公式(9)所示，补偿后的相干矩阵如公式(10)所示：
[0033][0034][0035]
所述步骤6中，采用代替相干矩阵(t》依次执行步骤3和步骤4，因为此时，由于mh＝0和mmmd＝0；的分解为五元散射能量分解。
[0036]
本发明的有益效果是，本发明将极化sar数据的相干矩阵的特征值分解用于多余散射能量分解方法，首先，将极化sar数据的相干矩阵进行特征值分解，提取特征值空间中的散射特征，包含散射机制角、散射熵和反熵等特征；设定七种散射模型对极化sar数据的相干矩阵进行七元散射能量分解，得到所有像素点的七种散射能量；以每个像素点的散射机制角、散射熵、反熵和各项异性参数为基础，对于满足特定条件的像素点更新其表面散射能量、二次散射能量和体散射能量的值，至此得到了极化sar数据的七元散射能量分解结果；然后，利用实数的方位角补偿矩阵对相干矩阵进行方位角补偿；基于散射机制角、散射熵和反熵等特征对极化sar数据的相干矩阵进行六元散射能量分解；获取极化sar数据的六元散射能量分解结果；最后，依次利用实数方位角补偿矩阵和复数方位角补偿矩阵对相干矩阵进行方位角补偿；基于散射机制角、散射熵和反熵等特征对极化sar数据的相干矩阵进行五元散射能量分解；获取极化sar数据的五元散射能量分解结果。
附图说明
[0037]
图1是本发明基于特征值分解的极化sar数据多元散射能量分解方法中采用的极化sar实例图像；
[0038]
图2是本发明基于特征值分解的极化sar数据多元散射能量分解方法中的熵的图像；
[0039]
图3是本发明基于特征值分解的极化sar数据多元散射能量分解方法中的反熵的图像；
[0040]
图4是本发明基于特征值分解的极化sar数据多元散射能量分解方法中的七元散射能量分解后表面散射能量、二次散射能量、体散射能量的合成图像；
[0041]
图5是本发明基于特征值分解的极化sar数据多元散射能量分解方法中的七元散射能量分解后的海洋区域的平均散射机制分布图；
[0042]
图6是本发明基于特征值分解的极化sar数据多元散射能量分解方法中的七元散射能量分解后的城区区域的平均散射机制分布图；
[0043]
图7是本发明基于特征值分解的极化sar数据多元散射能量分解方法中的七元散射能量分解后的森林区域的平均散射机制分布图；
[0044]
图8是本发明基于特征值分解的极化sar数据多元散射能量分解方法中的六元散射能量分解后表面散射能量，二次散射能量，体散射能量的合成图像；
[0045]
图9是本发明基于特征值分解的极化sar数据多元散射能量分解方法中的六元散射能量分解后的海洋区域的平均散射机制分布图；
[0046]
图10是本发明基于特征值分解的极化sar数据多元散射能量分解方法中的六元散射能量分解后的城区区域的平均散射机制分布图；
[0047]
图11是本发明基于特征值分解的极化sar数据多元散射能量分解方法中的六元散射能量分解后的森林区域的平均散射机制分布图；
[0048]
图12是本发明基于特征值分解的极化sar数据多元散射能量分解方法中的五元散射能量分解后表面散射能量、二次散射能量、体散射能量的合成图像；
[0049]
图13是本发明基于特征值分解的极化sar数据多元散射能量分解方法中的五元散射能量分解后的海洋区域的平均散射机制分布图；
[0050]
图14是本发明基于特征值分解的极化sar数据多元散射能量分解方法中的五元散射能量分解后的城区区域的平均散射机制分布图；
[0051]
图15是本发明基于特征值分解的极化sar数据多元散射能量分解方法中的五元散射能量分解后的森林区域的平均散射机制分布图。
具体实施方式
[0052]
下面结合附图和具体实施方式对本发明进行详细说明。
[0053]
实施例1
[0054]
本发明基于特征值分解的极化sar数据多元散射能量分解方法，具体包括如下步骤：
[0055]
步骤1，输入极化sar数据的相干矩阵可以表示为：
[0056][0057]
在极化sar系统中，公式(1)中的k
p
表示单视数据的pauli向量，为3
×
1的列向量，表示k
p
的共轭转置向量，符号《
·
》表示一个分辨单元内的多个单视数据的平均处理，所以《t》是一个半正定的赫米特矩阵，即特征值分解《t》得到的特征值都是非负的实数。图1为美国旧金山地区的l波段airsar全极化sar数据的合成图像，用相干矩阵的对角元素t
11
、t
22
和t
33
合成的图像，图像的大小为700
×
600，含有海洋、城区、森林等地物，适合检验极化sar数据分解算法的特征表达能力。在图1中使用矩形框由上至下分别选取海洋、城区和森林三
种地物类型，用于散射能量的展示。
[0058]
步骤2，对每个像素点的极化sar数据相干矩阵进行特征值分解，并计算散射熵和反熵特征，步骤2的具体过程为：
[0059]
对每个像素点的极化sar数据相干矩阵进行特征值分解的公式如下所示：
[0060][0061]
公式(2)中，λi表示特征值，ki表示对应的特征向量，αi表示散射机制角，βi表示交叉极化的相关幅度，φi、δi和ρi表示相位角，上标t表示转置操作，上标*表示共轭操作。因为《t》是一个半正定的赫米特矩阵，所以有λ1≥λ2≥λ3≥0。
[0062]
通过特征值，计算出熵h如公式(3)所示，其中pi＝λi/(λ1+λ2+λ3)，i∈{1，2，3}，反熵a如公式(4)所示。熵表示目标的随机程度，取值范围为0≤h≤1。h＝0表示目标随机程度最低，为确定性目标；h＝1表示目标的随机程度最高，为完全随机目标。图2为熵的图像，像素取值范围为0至1。反熵a作为熵的一个补偿参数，当熵h的值较高时(比如h≥0.7)才发挥作用，a≥0.5表示目标的异质性较高，通常发生于地物目标为人工建筑上；a＜0.5表示目标的异质性较低，通常发生于地物目标为植被区域。图3为反熵a的图像，像素取值范围为0至1。
[0063][0064][0065]
步骤3，使用七种散射模型分解每个像素点的极化sar数据的相干矩阵，结合特征值分解的特征，得到每个像素点的七种散射能量；
[0066]
步骤4，使用步骤2中的熵和反熵确定具有高熵和低反熵的像素点；重新计算该像素点的多元散射能量；
[0067]
步骤5，使用实数的方位角补偿矩阵对步骤1中的极化sar数据的相干矩阵进行方位角补偿，获取实数方位角补偿后的相干矩阵；
[0068]
步骤6，使用六种散射模型分解步骤5中得到的实数方位角补偿后相干矩阵，得到每个像素点的六种散射能量；
[0069]
步骤7，依次使用实数的方位角补偿矩阵和复数的方位角补偿矩阵对步骤1中的极化sar数据的相干矩阵进行方位角补偿，获取两次方位角补偿后的相干矩阵；
[0070]
步骤8，使用五种散射模型分解步骤7中得到的方位角补偿后相干矩阵，得到每个像素点的五种散射能量。
[0071]
实施例2
[0072]
在实施例1的基础上，步骤3的具体过程为：
[0073]
对步骤1中输入的相干矩阵进行七元散射能量分解，表示如下：
[0074]
《t》＝msts+mdtd+mvtv+mhth+m
md
t
md
+m
cd
t
cd
+m
od
t
od (5)；
[0075]
公式(5)中，ts表示表面散射模型，td表示二次散射模型，tv表示体散射模型，th表示螺旋体散射模型，t
md
表示混合偶极子散射模型，t
cd
表示复合偶极子散射模型，t
od
表示偏振偶极子散射模型。相应的，ms表示表面散射能量，md表示二次散射能量，mv表示体散射能
量，mh表示螺旋体散射能量，m
md
表示混合偶极子散射能量，m
cd
表示复合偶极子散射能量，m
od
表示偏振偶极子散射能量。
[0076]
七种散射模型分别表示如公式(6)-公式(12)所示，其中β为实数，α为复数，分别表示表面散射和二次散射模型中散射系数之间的关系。公式(8a)、公式(8b)，公式(8c)和公式(8d)表示不同情况下的四种体散射模型：
[0077][0078][0079][0080][0081][0082][0083][0084][0085][0086][0087]
将上述散射模型带入到公式(1)，如果体散射矩阵为公式(8a)所示，可以得到相干矩阵中的各个元素如下所示：
[0088][0089][0090][0091][0092][0093][0094]
通过公式(17)和公式(18)，解得mh、m
md
、m
cd
、m
od
分别为：
[0095]
mh＝2imag(t
23
)
ꢀꢀ
(19)；
[0096]mmd
＝2real(t
23
)
ꢀꢀ
(20)；
[0097]mcd
＝2imag(t
13
)
ꢀꢀ
(21)；
[0098]mod
＝2real(t
13
)
ꢀꢀ
(22)；
[0099]
在公式(19)-公式(22)中，imag表示取复数的虚部值，real表示取复数的实部值。
[0100]
去掉螺旋体散射分量，混合偶极子散射分量，复合偶极子散射分量，偏振偶极子散射分量后，剩余三个分量的极化sar的相干矩阵《t
re
》表示为公式(23)，因为和则《t
re
》如公式(24)所示。
[0101][0102][0103]
《t
re
)中的每个元素表示如下：
[0104][0105]
为保证ms、md和mv为非负值，首先必须确保和均为非负值。如果存在一些像素的和小于零的情况，必须适当限制mh、m
md
、m
cd
和m
od
以确保
和均为非负值。若均为非负值。若则需同时逐渐降低mh、m
md
、m
cd
和m
od
，直至满足如公式(26)所示，其中r表示四种散射能量减小的系数。然后，如果此时则令m
cd
＝0，m
od
＝0。类似地，如果则令mh＝0，m
md
＝0。
[0106][0107]
根据上述方法确定mh、m
md
、m
cd
和m
od
的值后，只与mv有关，设条件参数如果c1＜0，则体散射模型选择公式(8d)中的t
v4
；如果c1≥0，体散射模型选择公式(8a)、(8b)和(8c)中的一个，根据同极化比r来确定所用的体散射模型，同极化比r的表示如公式(27)所示。如果|r|≤2，体散射模型选择公式(8a)的t
v1
；如果r＞2，体散射模型选择公式(8b)的t
v2
；如果r＜-2，体散射模型选择公式(8c)的t
v3
。确定体散射模型后，体散射能量mv的解如公式(28)所示。
[0108][0109][0110]
《t
re
》中含有表面散射分量、二次散射分量和体散射分类，去掉体散射分类后的相干矩阵表示为《t
sd
》，如公式(29)所示。根据不同的体散射模型，《t
sd
》的元素可以表示为公式(30)。
[0111]
[0112][0113]
对《t
sd
》进行特征值分解，获得两个特征值和特征向量，如公式(31)和公式(32)所示。λ1和λ2为特征值，令λ1≥λ2，k1和k2为特征向量，利用散射机制角α1确定表面散射能量和二次散射能量，具体方法如公式(33)所示。
[0114]
《t
sd
》＝λ1k1+λ2k2ꢀꢀ
(31)；
[0115][0116][0117]
其中，φi和δi表示相位角。如果λ2小于零，则逐步降低mv的值，mv＝mv*r，其中r∈{0.9，0.8，0.7，...，0.1，0}的取值由大到小，直至满足λ2≥0，然后再利用公式(31)，公式(32)和公式(33)求取ms和md。
[0118]
实施例3
[0119]
在实施例2的基础上，步骤4的具体过程为：
[0120]
为了获得非负特征值λ1和λ2，在步骤3中降低了体散射能量mv，这会导致体散射能量mv被低估，尤其是在高熵的植被区域。为了解决该问题，重新定义了高熵植被区域的ms、md和mv的计算方法。
[0121]
当一个像素的熵h和反熵a满足关系h-a＞0.4时(熵h和反熵a由步骤2中得到)，表示当前目标具有较高的熵h和较低的反熵a，则认为该像素属于植被区域，ms、md和mv的计算方法如下：
[0122][0123]
至此步为止，极化sar图像的七元散射能量分解完成。图4为七元散射能量分解后的表面散射能量、二次散射能量和体散射能量的合成图像。图5为图1选取的海洋区域的七种散射能量的平均值，图6为图1选取的城区区域的七种散射能量的平均值，图7为图1选取的森林区域的七种散射能量的平均值。
[0124]
步骤5的具体过程为：
[0125]
实数的方位角补偿矩阵如公式(35)所示：
[0126][0127]
使用r(θ)对步骤1中的极化sar数据的相干矩阵《t》进行方位角补偿，获取实数方位角补偿后的相干矩阵《t(θ)》，如公式(36)所示，求解θ的值如公式(37)所示，《t(θ)》中的每个元素如公式(38)所示。
[0128][0129][0130][0131]
步骤6的具体过程为：
[0132]
用步骤5中的《t(θ)》代替相干矩阵《t》依次执行步骤3和步骤4。根据公式(36)可以得到t
23
(θ)的实部为(θ)的实部为将公式(37)的θ代入，得到real(t
23
(θ))＝0。根据公式(20)，可知m
md
＝0；因此《t(θ)》的分解为六元散射能量分解。图8为六元散射能量分解后的表面散射能量、二次散射能量和体散射能量的合成图像。图9为图1选取的海洋区域的六种散射能量的平均值，图10为图1选取的城区区域的六种散射能量的平均值，图11为图1选取的森林区域的六种散射能量的平均值。
[0133]
步骤7的具体过程为：
[0134]
实数的方位角补偿矩阵如公式(35)所示，补偿后的相干矩阵《t(θ)》如公式(36)和公式(38)所示。
[0135]
实数的方位角补偿后再进行复数的方位角补偿，复数的方位角补偿矩阵如公式(39)所示，补偿后的相干矩阵如公式(40)所示，求解得到的值如公式(41)所示，中的每个元素如公式(42)所示。
[0136][0137][0138][0139][0140]
步骤8的具体过程为：
[0141]
用步骤7中的代替相干矩阵《t》依次执行步骤3和步骤4。
[0142]
因为将公式(41)的代入，可以得到根据公式(19)和公式(20)，可知mh＝0和mmd＝0；因此的分解为五元散射能量分解。图12为五元散射能量分解后的表面散射能量、二次散射能量和体散射能量的合成图像。图13为图1选取的海洋区域的五种散射能量的平均值，图14为图1选取的城区区域的五种散射能量的平均值，图15为图1选取的森林区域的五种散射能量的平均值。

技术特征：
1.基于特征值分解的极化sar数据多元散射能量分解方法，其特征在于：具体包括如下步骤：步骤1，输入极化sar数据的相干矩阵；步骤2，对每个像素点极化sar数据相干矩阵进行特征值分解；步骤3，分解每个像素点的极化sar数据的相干矩阵；步骤4，计算具有较高熵和较低反熵的像素点的多元散射能量；步骤5，获取实数方位角补偿后的相干矩阵并分解，得到每个像素点的六种散射能量；步骤6，获取两次方位角补偿后的相干矩阵并分解，得到每个像素点的五种散射能量。2.根据权利要求1所述的基于特征值分解的极化sar数据多元散射能量分解方法，其特征在于：所述步骤1中，极化sar数据的相干矩阵表示为：其中，<t>是一个半正定的赫米特矩阵，k
p
表示单视数据的pauli向量。3.根据权利要求2所述的基于特征值分解的极化sar数据多元散射能量分解方法，其特征在于：所述步骤2中，对每个像素点的极化sar数据相干矩阵进行特征值分解的公式如下所示：其中，λ
i
表示特征值，k
i
表示对应的特征向量，α
i
表示散射机制角，λ1≥λ2≥λ3≥0；熵h如公式(3)所示，其中p
i
＝λ
i
/(λ1+λ2+λ3)，i∈{1，2，3}，反熵a如公式(4)所示：)，i∈{1，2，3}，反熵a如公式(4)所示：4.根据权利要求3所述的基于特征值分解的极化sar数据多元散射能量分解方法，其特征在于：所述步骤3中，相干矩阵进行七元散射能量分解如下公式(5)所示：<t>＝m
s
t
s
+m
d
t
d
+m
v
t
v
+m
h
t
h
+m
md
t
md
+m
cd
t
cd
+m
od
t
od
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(5)；其中，t
s
表示表面散射模型，t
d
表示二次散射模型，t
v
表示体散射模型，t
h
表示螺旋体散射模型，t
md
表示混合偶极子散射模型，t
cd
表示复合偶极子散射模型，t
od
表示偏振偶极子散射模型，m
s
表示表面散射能量，m
d
表示二次散射能量，m
v
表示体散射能量，m
h
表示螺旋体散射能量，m
md
表示混合偶极子散射能量，m
cd
表示复合偶极子散射能量，m
od
表示偏振偶极子散射能量。5.根据权利要求4所述的基于特征值分解的极化sar数据多元散射能量分解方法，其特征在于：所述步骤4中，当一个像素的熵h和反熵a满足关系h-a＞0.4时，m
s
、m
d
和m
v
的计算方法如下：
6.根据权利要求5所述的基于特征值分解的极化sar数据多元散射能量分解方法，其特征在于：所述步骤5中，实数方位角补偿后的相干矩阵<t(θ)>如公式(7)所示：其中，实数的方位角补偿矩阵r(θ)如公式(8)所示：采用<t(θ)>代替相干矩阵<t>依次执行步骤3和步骤4，由于t
23
(θ)的实部为0，此时，m
md
＝0，因此，<t(θ)>的分解为六元散射能量分解。7.根据权利要求6所述的基于特征值分解的极化sar数据多元散射能量分解方法，其特征在于：所述步骤6中，采用实数的方位角补偿后再进行复数的方位角补偿，复数的方位角补偿矩阵如公式(9)所示，补偿后的相干矩阵如公式(10)所示：如公式(10)所示：所述步骤6中，采用代替相干矩阵<t>依次执行步骤3和步骤4，因为此时，由于m
h
＝0和m
md
＝0；的分解为五元散射能量分解。

技术总结
本发明公开了一种基于特征值分解的极化SAR数据多元散射能量分解方法，具体包括如下步骤：步骤1，输入极化SAR数据的相干矩阵；步骤2，对每个像素点的极化SAR数据相干矩阵进行特征值分解，并计算散射熵和反熵等特征；步骤3，使用七种散射模型分解每个像素点的极化SAR数据的相干矩阵，结合特征值分解的特征，得到每个像素点的七种散射能量；步骤4，进行极化SAR数据的七元散射能量分解；步骤5，获取实数方位角补偿后的相干矩阵；步骤6，进行极化SAR数据的六元散射能量分解；步骤7，获取两次方位角补偿后的相干矩阵；步骤8，进行极化SAR数据的五元散射能量分解。元散射能量分解。元散射能量分解。


技术研发人员：张爽 王璐 崔真 罗静蕊
受保护的技术使用者：西安理工大学
技术研发日：2023.07.10
技术公布日：2023/10/15