一种基于最小剩余功率的极化SAR图像目标分解方法

摘要

一种基于最小剩余功率的极化SAR图像目标分解方法，包括对极化SAR图像中的像素进行散射矢量构造并计算极化协方差矩阵的步骤；计算总功率、体散射加权系数、面散射加权系数和偶次散射加权系数的步骤；对面散射加权系数、偶次加权系数、体散射加权系数进行扩展操作的步骤；通过计算最小剩余功率对面散射加权系数、偶次散射加权系数和体散射加权系数进行修正的步骤；利用修正后的加权系数进行面散射功率、偶次散射功率和体散射功率计算的步骤；本发明计算简便高效，可有效获得极化SAR目标分解后的多种散射机制功率和物理参数，能显著地提取各种地物的散射特性和目标特征，可用于雷达目标识别技术领域。

权利要求

1.一种基于最小剩余功率的极化SAR图像目标分解方法，其特征在于：包括如下步骤：A、对极化SAR图像中的像素进行散射矢量构造的步骤；

B、由散射矢量计算像素的极化协方差矩阵的步骤；

C、由极化协方差矩阵计算像素的总功率、体散射加权系数、面散射加权系数和偶次散射加权系数的步骤；

D、对面散射加权系数和偶次加权系数进行扩展操作的步骤；

E、对体散射加权系数进行扩展操作的步骤；

F、通过计算最小剩余功率对面散射加权系数、偶次散射加权系数和体散射加权系数进行修正的步骤；

G、利用修正后的面散射加权系数、偶次散射加权系数和体散射加权系数进行面散射功率、偶次散射功率和体散射功率计算的步骤；

其中，所述步骤D的具体操作方法为：

(1) 、设置取值在(0 ,1]之间的扩展系数eco和正整数m ，依据公式计算序号i∈{1,2,…,m}对应的面散射加权系数fs,i；

(2)、设置正整数n，依据公式 计算序号j∈{1,2,…,n}对应的偶次散射加权系数fd,j；

其中，fs为由极化协方差矩阵计算得到的面散射的加权系数，fd为由极化协方差矩阵计算得到的偶次散射的加权系数；

所述步骤E的具体操作方法为：

(1)、利用公式Ps,i＝(1+|α|2)fs,i，对所有的i∈{1,2,…,m}，计算其对应的面散射功率Ps,i，利用公式Pd,j＝(1+|β|2)fd,j，对所有的j∈{1,2,…,n}，计算其对应的偶次散射功率Pd,j；

(2)、对所有的i∈{1,2,…,m}和j∈{1,2,…,n}的组合，利用公式fv,k＝8(Span-Ps,i-Pd,j)/3，计算k∈{1,2,…,m×n}集合内其对应的体散射加权系数fv,k；

其中，α为偶次散射特征参数，β为面散射特征参数，Span为由极化协方差矩阵计算得到的像素的总功率；

所述步骤F的具体操作方法为：

(1)、利用公式Crem,k＝C-fs,i×Cs-fd,j×Cd-fv,k×Cv，对所有的k∈{1,2,…,m×n}，计算其对应的剩余协方差矩阵Crem,k，其中，C为由散射矢量计算得到的像素的极化协方差矩阵，(2)、利用公式Prem,k＝trace(Crem,k)，计算其对应的剩余功率Prem,k；

(3)、利用公式 得到最小剩余功率 修正后的面散射加权系数fs'、偶次散射加权系数fd'、体散射加权系数fv'分别为取最小剩余功率 时所对应的fs,i、fd,j和fv,k。

2.根据权利要求1所述的一种基于最小剩余功率的极化SAR图像目标分解方法，其特征在于，所述步骤A的具体操作方法为：设极化SAR图像中某个像素对应的散射矩阵为在单基站电磁波互易原理下，SHV＝SVH，则该像素对应的散射矢量为

3.根据权利要求1所述的一种基于最小剩余功率的极化SAR图像目标分解方法，其特征在于：所述步骤B中，利用公式 由散射矢量k计算得到该像素对应的极化协方差矩阵C；

其中，H表示矢量转置共轭，*表示复共轭，<·>为集合平均。

4.根据权利要求3所述的一种基于最小剩余功率的极化SAR图像目标分解方法，其特征在于，所述步骤C的具体操作方法为：(1)、利用公式Span＝C11+C22+C33和fv＝3C22/2，计算该像素对应的总功率Span和体散射加权系数fv；

(2)、当ReC13>0时，α＝-1，利用公式

来计算该像素对应的

面散射的加权系数fs和偶次散射的加权系数fd，

当ReC13≤0时，β＝1，利用公式

来计算该像素对应的面

散射的加权系数fs和偶次散射的加权系数fd；

其中，ReC13为协方差矩阵元素C13的实部，ImC13协方差矩阵元素C13的虚部，real(β)为β的实部，imag(β)为β的虚部，real(α)为α的实部，imag(α)为α的虚部。

5.根据权利要求1所述的一种基于最小剩余功率的极化SAR图像目标分解方法，其特征在于：所述步骤G中，利用公式Ps＝fs'(1+|α|2)、Pd＝fd'(1+|β|2)、 分别计算面散射功率Ps、偶次散射功率Pd和体散射功率Pv。

6.根据权利要求1所述的一种基于最小剩余功率的极化SAR图像目标分解方法，其特征在于：在步骤D中，所述的扩展系数eco取0.3。

说明书

一种基于最小剩余功率的极化SAR图像目标分解方法

技术领域

[0001] 本发明属于雷达自动目标识别技术领域，具体的说是一种基于最小剩余功率的极化SAR图像目标分解方法，用于估计目标中几种不同的散射机制所占的分量，为目标分类识别提供重要的特征信息。

背景技术

[0002] 机载和星载极化合成孔径雷达发射电磁波照射地面，并接收该地面区域的散射回波，从而获得每个分辨单元对应的散射矩阵。该散射矩阵包含方位指向、物理构成、形状结构、极化特征等有效信息，能够较为完整的描述雷达照射目标的电磁散射特性。对由散射矩阵构成的极化SAR图像进行准确的目标特征识别提取，始终是极化SAR图像解译迫切需要解决的一个核心问题。

[0003] 极化目标分解是极化SAR图像特征提取的主要实现方法。极化SAR目标分解可分为描述纯目标的相干类和描述分布式目标的非相干两类方法。非相干类目标分解方法的基本思想是将极化相干矩阵或协方差矩阵拆分为若干个典型散射机制的加权和，典型散射机制包括面散射、偶次散射、体散射和螺旋体散射等。例如当前广泛使用的Freeman-Durden的3成分分解方法即是将极化SAR的协方差矩阵分解为面散射、偶次散射和体散射的叠加。这些传统的非相干类目标分解方法首先得到体散射功率，然后再计算面散射功率和偶次散射功率。由于总功率固定且优先计算获得体散射功率，高估计的体散射功率会导致偶次散射功率和面散射功率被过低估计，甚至出现负功率。

[0004] 当前解决方法大致可分为三类，第一种方法是构造通用或自适应的体散射模型准确的描述各类自然实物地貌，从而得到各典型散射机制的更合理功率值；第二类是在分解过程中添加非负功率约束，通过消除负功率进而解决体散射过估计的问题。第三类是引入极化方位角，通过最小化交叉极化项来降低体散射功率，补偿增加偶次散射或面散射功率。这三类方法中，第一类方法在构造体散射模型时，假设森林散射方位向对称或反射对称等条件成立，这种状况下部分实物地貌与此假设前提不符，因此构造的体散射模型并不是针对所有各类地貌都十分有效。第二类方法在分解过程中，需要迭代计算矩阵行列式特征值或求解方程组，计算复杂度较高。第三类方法虽然能够有效降低森林区域的体散射功率，但是对城市区域的偶次散射功率补偿增加较小，依然可能出现城市建筑物区域体散射被过高估计的问题。

[0005] 发明实用新型内容

[0006] 为了克服现有技术中体散射功率易被过高估计、偶次散射功率和面散射功率易被过低估计带来的目标分类识别误差，本发明提出了一种基于最小剩余功率的极化SAR图像目标分解方法，该方法通过最小剩余功率的约束，修正了不同散射机制的加权系数，从而解决了现有方法中体散射功率高估和偶次散射功率的低估问题。

[0007] 本发明为解决上述技术问题，所提供的技术方案是：一种基于最小剩余功率的极化SAR图像目标分解方法，包括如下步骤：

[0008] A、对极化SAR图像中的像素进行散射矢量构造的步骤；

[0009] B、由散射矢量计算像素的极化协方差矩阵的步骤；

[0010] C、由极化协方差矩阵计算像素的总功率、体散射加权系数、面散射加权系数和偶次散射加权系数的步骤；

[0011] D、对面散射加权系数和偶次加权系数进行扩展操作的步骤；

[0012] E、对体散射加权系数进行扩展操作的步骤；

[0013] F、通过计算最小剩余功率对面散射加权系数、偶次散射加权系数和体散射加权系数进行修正的步骤；

[0014] G、利用修正后的面散射加权系数、偶次散射加权系数和体散射加权系数进行面散射功率、偶次散射功率和体散射功率计算的步骤；

[0015] 其中，所述步骤D的具体操作方法为：

[0016] (1)、设置取值在(0 ,1]之间的扩展系数eco和正整数m，依据公式计算序号i∈{1,2,…,m}对应的面散射加权系数fs,i；

[0017] (2)、设置正整数n，依据公式 计算序号j∈{1,2,…,n}对应的偶次散射加权系数fd,j；

[0018] 其中，fs为由极化协方差矩阵计算得到的面散射的加权系数，fd为由极化协方差矩阵计算得到的偶次散射的加权系数；

[0019] 所述步骤E的具体操作方法为：

[0020] (1)、利用公式Ps,i＝(1+|α|2)fs,i，对所有的i∈{1,2,…,m}，计算其对应的面散射功率Ps,i，利用公式Pd,j＝(1+|β|2)fd,j，对所有的j∈{1,2,…,n}，计算其对应的偶次散射功率Pd,j；

[0021] (2)、对所有的i∈{1,2,…,m}和j∈{1,2,…,n}的组合，利用公式fv,k＝8(Span-Ps,i-Pd,j)/3，计算k∈{1,2,…,m×n}集合内其对应的体散射加权系数fv,k；

[0022] 其中，α为偶次散射特征参数，β为面散射特征参数，Span为由极化协方差矩阵计算得到的像素的总功率；

[0023] 所述步骤F的具体操作方法为：

[0024] (1)、利用公式Crem,k＝C-fs,i×Cs-fd,j×Cd-fv,k×Cv，对所有的k∈{1,2,…,m×n}，计算其对应的剩余协方差矩阵Crem,k，其中，C为由散射矢量计算得到的像素的极化协方差矩阵，

[0025]

[0026] (2)、利用公式Prem,k＝trace(Crem,k)，计算其对应的剩余功率Prem,k；

[0027] (3)、利用公式 得到最小剩余功率 修正后的面散射加权系数fs'、偶次散射加权系数f'd、体散射加权系数f'v分别为取最小剩余功率 时所对应的fs,i、fd,j和fv,k。

[0028] 所述步骤A的具体操作方法为：设极化SAR图像中某个像素对应的散射矩阵为在单基站电磁波互易原理下，SHV＝SVH，则该像素对应的散射矢量为

[0029] 所述步骤B中，利用公式 由散射矢量k计算得到该像素对应的极化协方差矩阵C；

[0030] 其中，H表示矢量转置共轭，*表示复共轭，<·>为集合平均。

[0031] 所述步骤C的具体操作方法为：

[0032] (1)、利用公式Span＝C11+C22+C33和fv＝3C22/2，计算该像素对应的总功率Span和体散射加权系数fv；

[0033] (2)、当ReC13>0时，α＝-1，利用公式

[0034]

[0035]

[0036]

[0037] 来计算该像素对应的面散射的加权系数fs和偶次散射的加权系数fd，

[0038] 当ReC13≤0时，β＝1，利用公式

[0039]

[0040]

[0041]

[0042] 来计算该像素对应的面散射的加权系数fs和偶次散射的加权系数fd；

[0043] 其中，ReC13为协方差矩阵元素C13的实部，ImC13协方差矩阵元素C13的虚部，real(β)为β的实部，imag(β)为β的虚部，real(α)为α的实部，imag(α)为α的虚部。

[0044] 所述步骤G中，利用公式Ps＝fs'(1+|α|2)、Pd＝fd'(1+|β|2)、 分别计算面散射功率Ps、偶次散射功率Pd和体散射功率Pv。

[0045] 在步骤D中，所述的扩展系数eco取0.3。

[0046] 有益效果：

[0047] 本发明的一种基于最小剩余功率的极化SAR图像目标分解方法，先计算得到各种典型散射机制的基本加权系数，然后再以最小剩余功率为约束修正各种典型散射机制的加权系数，可有效改善传统的目标分解方法造成的体散射功率高估和偶次散射功率的低估问题，使目标分类识别的结果更合理、有效和精确。

附图说明

[0048] 图1是本发明的实现流程图；

[0049] 图2是San Francisco地区的全极化数据对应的后向散射功率图像及7个典型地物选择区域；

[0050] 图3是美国Haywrd地区的极化数据对应的后向散射功率图像及5个典型地物选择区域。

具体实施方式

[0051] 极化SAR作为一种主动遥感系统，具有全天时全天候的工作能力，近年来在军事、城市规划、测绘、病虫害防治、森林防火和地质勘探等领域得到了广泛应用。基于模型的非相干极化分解能够获取多个散射机制的功率和参数，从而能有效提取极化SAR数据中地物目标的物理特性和极化特征，是近年来极化SAR领域内最重要研究方向之一。非相干极化分解过程是先对地面反射电磁波对应的散射矩阵进行统计平均构造为二阶描述式，如Mueller矩阵、协方差矩阵或相干矩阵等，再将这些二阶描述式分解为几种典型散射机制的的线性组合。

[0052] Freeman等人所提的目标分解方法假设所有地物由三种典型散射机制构成，将极化协方差矩阵C分解为体散射、面散射和偶次散射模型组成的线性组合，表达式为：

[0053] C＝fsCs+fdCd+fvCv+Crem

[0054] 其中，C为极化SAR中最小分辨单元对应的协方差矩阵，fs，fd，fv依次为面散射加权系数、偶次散射加权系数和体散射的加权系数，Cs，Cd，Cv依次为面散射、偶次散射和体散射的极化协方差矩阵，Crem为分解后剩余项。

[0055] 参照图1，本发明的一种基于最小剩余功率的极化SAR图像目标分解方法，具体实施步骤如下：

[0056] 步骤一、设极化SAR图像中某个像素对应的散射矩阵为 其中，下标H和V分别表示水平和垂直极化状态，且每个散射矩阵元素中的第1个下标代表接收信号的极化状态，第2个下标代表发射信号的极化状态，例如SHV中H代表接收信号的极化状态，V代表发射信号的极化状态，在单基站电磁波互易原理下，散射矩阵S可认为是对称的，即SHV＝SVH，用散射矩阵构造散射矢量

[0057] 步骤二、由散射矢量k计算得到该像素的极化协方差矩阵C，计算公式为其中，H表示矢量转置共轭，*表示复共轭，<·>为集合平均；

[0058] 步骤三、计算该像素的总功率Span＝C11+C22+C33和体散射加权系数fv＝3C22/2；

[0059] 步骤四、计算得到面散射的加权系数fs和偶次散射的加权系数fd，以及对应的参数β和α，α为偶次散射特征参数，β为面散射特征参数；

[0060] 具体步骤为：

[0061] 当ReC13>0时

[0062] α＝-1

[0063]

[0064]

[0065]

[0066]

[0067] 其中，ReC13为协方差矩阵元素C13的实部，ImC13协方差矩阵元素C13的虚部，real(β)表示β的实部，imag(β)表示β的虚部；

[0068] 当ReC13≤0时

[0069] β＝1

[0070]

[0071]

[0072]

[0073]

[0074] 其中，ReC13为协方差矩阵元素C13的实部，ImC13协方差矩阵元素C13的虚部，real(α)表示α的实部，imag(α)表示α的虚部；

[0075] 步骤五、设置扩展系数eco的初始值，且eco＞0，扩展系数eco的范围一般取(0,1]，实际经验来说取0.3发明效果较好，一般根据实际经验来取，或者取几个数值，如0.3，0.5，0.6，来选择一个实验效果最好的，一般取0.3即可；对偶次加权系数fd和面散射加权系数fs进行扩展操作；

[0076] 对偶次加权系数fd和面散射加权系数fs进行扩展操作的过程为：

[0077] (1)、设置取值在(0 ,1]之间的扩展系数eco和正整数m，依据公式计算序号i∈{1,2,…,m}对应的面散射加权系数fs,i；

[0078] (2)、设置正整数n，依据公式 计算序号j∈{1,2,…,n}对应的偶次散射加权系数fd,j；

[0079] 其中，i和j均表示区间中某一份的序号数值；

[0080] 步骤六、对所有的i∈{1,2,…,m}和j∈{1,2,…,n}，由总功率Span、面散射加权系数fs,i和偶次散射加权系数fd,j计算得到所有k∈{1,2,…,m×n}的体散射加权系数fv,k，其中，k表示集合k∈{1,2,…,m×n}中的一个序号数值；

[0081] 具体过程为：

[0082] 6a)对所有的i∈{1,2,…,m}，计算对应的面散射功率值Ps,i＝(1+|α|2)fs,i，对所有的j∈{1,2,…,n}，计算对应的偶次散射功率值Pd,j＝(1+|β|2)fd,j，

[0083] 6b)对所有的i∈{1,2,…,m}和j∈{1,2,…,n}的组合，计算k∈{1,2,…,m×n}集合内所有的体散射加权系数fv,k＝8(Span-Ps,i-Pd,j)/3；

[0084] 步骤七、对所有下标k∈{1,2,…,m×n}，计算剩余协方差矩阵Crem,k＝C-fs,i×Cs-fd,j×Cd-fv,k×Cv，其中， 然后再计算剩余功率为Prem,k＝trace(Crem,k)，trace(·)表示求矩阵的迹；

[0085] 步骤八、确定修正后的面散射加权系数fs'，修正后的偶次散射加权系数fd'和修正2
后的体散射加权系数fv'，并计算面散射功率Ps＝fs'(1+|α|)、偶次散射功率Pd＝fd'(1+|β|2)和体散射功率

[0086] 其中，确定fs'，fd'，fv'的具体过程为：

[0087] 对所有的k∈{1,2,…,m×n}，最小剩余功率为 同时得到取最小剩余功率 时对应的面散射加权系数为fs'＝fs,i，偶次散射加权系数fd'＝fd,i和体散射加权系数fv'＝fv,k。

[0088] 本发明中，采用最小剩余功率进行约束的原理是：剩余功率最小是从物理角度推导出的要求，传统的目标分解方法没有考虑这个约束，本发明专利添加此要求后，得出的结果必然更加合理、有效和精确，符合物理意义。就如同原来的汽车没有刹车，现在有了刹车系统，更加合理。

[0089] 下面通过实验对本发明的效果做进一步的说明：

[0090] 1、实验内容

[0091] 实验使用两组数据，第1组为由星载系统Radarsat-2获取的San Francisco(旧金山)地区C波段数据，方位向分辨率为4.82米，距离向为4.73米，包含多种地物类型诸如森林、城市建筑物、海洋等；第2组机载UAVSAR系统获取美国Haywrd地区的L波段全极化数据，方位向分辨率为7.2米，距离向分辨率为5米。

[0092] 实验过程如下：分别采用Freeman-Durden的文章[Athree-component scattering model for polarimetric SAR data]中的基于模型的三成分方法与本发明基于最小剩余功率的极化SAR图像目标分解方法对上述两组极化数据进行分解。极化协方差矩阵对极化数据采用5*5的滑窗计算得到。在两组数据中分别选择7个和5个不同的典型区域，计算每个区域中各个散射机制所占功率的百分比。参数设定为eco＝0.3，m＝5和n＝5。

[0093] 2、实验数据分析

[0094] 在表1和表2中，Sur代面散射机制对应的散射功率占总功率的比，Vol代表体散射机制对应的散射功率占总功率的比，Dbl代表偶次散射机制对应的散射功率占总功率的比。假设选取的某一个区域中共包含有N个像素点则
其中，Pd,i，Ps,i，Pv,i，Spani分别表示第i个像素点中偶次散射机制对应的功率值、面散射机制对应的功率值，体散射机制对应的功率值和总功率值。

[0095] 表1 San Francisco数据不同分解分解方法下各散射机制所占百分比[0096]

[0097]

[0098] 由San Francisco地区对应的Google Earth光学图像可知，区域1～区域3为城市建筑物。对城市建筑区，偶次散射机制对应的散射功率所占比例越高，越利于后续的检测识别。

[0099] 由表1可知，本发明方法在区域1中的偶次散射所占功率比值为51.4％，而Freeman-Durden方法偶次散射所占功率百分比为39.6％，较本发明方法功率值减少11.8％。对于区域2和区域3，本发明方法相比Freeman-Durden方法偶次散射机制功率分别高出9％和12％。这是由于本发明方法有效的抑制了城市建筑物区域体散射过高估计，增加了偶次散射机制的功率值，说明本发明方法在城市建筑物区获取散射机制功率更为准确有效。在区域4和区域5中，由光学图像可知对应于森林植被区域，本发明与Freeman-Durden方法分解结果中体散射机制所占功率值相差不大，表明在这两个区域体散射机制占主导，与实际地貌一致。对于区域6和区域7，由光学图像可知相应地貌为海洋，本发明方法分解得到的面散射功率相对于Freeman-Durden方法由小幅度提升，约为3％左右。这是由于本发明增加了最小剩余功率为约束，降低了目标分解方法的剩余功率值，增加了海洋区域作为主导散射机制的面散射功率值。

[0100] 表2 Haywrd数据不同分解分解方法下各散射机制所占百分比

[0101]

[0102] 由美国Haywrd地区对应的光学图像可知，区域1～区域3对应的是城市等人造目标地区。

[0103] 由表2可知，Freeman-Durden方法分解后在区域1中的偶次散射机制所占功率比例为42.9％，采用本发明方法结果中偶次散射功率所占比例为51.1％；同样，本发明方法对区域2和区域3的目标分解结果偶次散射功率分别比Freeman-Durden方法高3.9％和7.6％；本发明使得城市区域的偶次散射机制对应的功率比例优于Freeman-Durden方法。区域4和区域5对应的是海洋，通常面散射机制占主导，本发明使得该区域对应的主要散射功率所占比例高于采用Freeman-Durden方法求出的主导散射功率所占比例，最大增加值为8％。

[0104] 综上所述，本发明方法对极化数据进行目标分解后得到的功率值和特征参数能更有效反映出各种实际地物的散射特性和物理结构，优于传统使用的Freeman-Durden方法。

[0105] 以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的包含范围内。