CN104122554B - 高分辨率sar图像目标的属性散射中心特征提取方法 - Google Patents

Abstract

本发明公开一种高分辨率SAR图像目标属性散射中心特征提取方法，主要解决现属性散射中心提取方法的模型失配、计算复杂度高、参数估计不准确和无法满足目标识别自动化要求的问题。其实现过程是：1.提取雷达图像中目标的支撑连通区域，并自动初始化参数；2.根据属性散射中心模型，利用字典拆分实现字典降维，构建位置参数字典和属性参数字典；3.通过求解0范数优化问题得到目标属性散射中心参数集合。本发明能有效提取目标属性散射中心，可用于精确估计目标及其重要部件几何尺寸，满足目标识别自动化要求。

Description

高分辨率SAR图像目标的属性散射中心特征提取方法

技术领域

本发明属于雷达技术领域，涉及一种雷达目标属性散射中心特征提取方法，可用于估计目标及其重要部件的几何尺寸，为目标分类识别提供重要的特征信息。

背景技术

合成孔径雷达SAR技术是从20世纪50年代发展起来的，在近七十年里蓬勃发展。提高SAR的空间分辨率并从中提取有效的目标信息一直是SAR领域的研究热点。高分辨率SAR图像是目标的二维散射图像，可近似看成由少量强散射中心的后向散射回波之和产生的。目前用于提取目标散射中心的模型主要分为三种：理想点散射中心模型，衰减指数模型，以及属性散射中心模型。其中属性散射中心模型最符合SAR目标的散射特性，它描述的属性散射中心具有丰富的几何意义和物理意义。

属性散射中心模型是1999年Michael J.Gerry和Lee C.Potter基于几何绕射理论和物理光学理论提出的一个适用于高分辨SAR图像数据的参数化模型。见[M.J.Gerry,L.C.Potter,I.J.Gupta,and A.van der Merwe,A parametric model for syntheticaperture radar measurements[J].IEEE Transactions on Antennas and Propagation,1999,Vol.47,NO.7,pp.1179-1188]。该模型不仅描述了散射中心的位置信息和强度信息，还描述了散射中心的几何尺寸，以及散射中心对频率和方位的依赖性，并可以由此反演出目标的几何结构。属性散射中心模型相比于点散射模型和衰减指数模型，包含更丰富的可用于目标分类识别的特征。

属性散射中心特征提取实质上是一个从目标回波数据中估计各个散射中心参数的过程。由于属性散射中心模型包含的散射中心信息较多，导致模型结构复杂以及参数维数较高，增加了模型参数估计的复杂性。现有方法主要分为两种，分别是基于图像域的散射中心提取和基于频域数据的散射中心提取。基于图像域数据的散射中心提取通过对SAR图像进行图像分割，得到阶数较低的目标散射区或者是孤立的散射中心，利用近似最大似然方法估计目标的属性散射中心参数。这种方法对图像质量要求较高，对于某些散射中心由于倾斜而导致的局部强度较弱的情况难以估计准确。此外由于该类方法需要求解非凸优化问题，并且存在很多局部最小解，所以存在参数初始化、模型阶数选择和散射中心结构类别判别的问题，使得最终参数估计精度较低。基于频域数据的散射中心提取避免了模型失配和散射中心类别判别的问题，但是由于属性散射中心的参数维度较大，计算复杂度高，现有的基于频域数据的散射中心提取方法在估计参数时都采取分步估计的方法以降低计算复杂度，这样难免会降低参数的估计精度，此外这些方法在初始化参数时不能自动选取参数支撑区，无法满足目标识别的自动化要求。

发明内容

本发明的目的在于提出一种高分辨率SAR图像目标的属性散射中心特征提取方法，以解决现有方法中复杂度高，参数估计不准确以及无法自动化处理的问题。

本发明是这样实现的：

一.技术思路

在高频区的雷达回波中，目标散射回波可近似看成由少量强散射中心的后向散射回波叠加产生，说明雷达回波在属性散射中心的参数空间具有很强的稀疏性。考虑雷达接收回波中掺杂一些噪声信号，利用图像处理的技术降低噪声的影响；通过在图像域提取目标支撑区实现参数的自动初始化；考虑到属性散射中心参数空间维数较高，导致构造的冗余字典维数远远大于实际可处理的维数，基于属性散射中心模型，通过拆分冗余字典降低字典维度，并通过结合属性散射中心模型特点实现对字典的再次降维，再通过求解0范数优化问题提取属性散射中心特征。

二.技术方案

根据上述思路本发明的实现步骤包括如下：

(1)提取原始图像I₀中的局部极大值区域，得到目标图像I₁；对目标图像I₁进行分割，得到目标的轮廓区域R，用轮廓区域R对原始图像I₀进行掩膜处理，得到目标连通支撑域图像I₂；

(2)确定散射中心的位置集合Θ₁和散射中心的属性集合Θ₂。

2a)根据轮廓区域R中的非零元素位置得到散射中心距离维坐标参数x的取值范围X和散射中心方位维坐标参数y的取值范围Y，确定散射中心的位置集合Θ₁，Θ₁＝{(x,y)|x∈X,y∈Y}；

2b)由散射中心方位维坐标参数y的取值范围确定散射中心长度参数L的取值范围L；由雷达回波数据录取的方位角域确定散射中心方位角取值范围Φ；由雷达回波数据录取的中心频率确定散射中心方位依赖因子γ取值范围Γ；设定散射中心的频率依赖因子α取值范围为Λ＝{-1,-0.5,0,0.5,1}；最终确定散射中心的属性集合Θ₂：

${\mathrm{\Θ}}_2=\left\{(L,\stackrel{\‾}{\mathrm{\φ}},\mathrm{\α},\mathrm{\γ})\right|L\∈L,\stackrel{\‾}{\mathrm{\φ}}\∈\mathrm{\Φ},\mathrm{\α}\∈\mathrm{\Λ},\mathrm{\γ}\∈\mathrm{\Γ}\};$

(3)根据所述的位置集合Θ₁和属性集合Θ₂，利用属性散射中心模型分别构建位置信息字典D₁(x,y|Θ₁)和属性信息字典

(4)利用正交匹配追踪OMP法得到目标属性散射中心参数的初步估计集合Θ₀；

(5)对初步估计集合Θ₀进行优化，得到目标属性散射中心的特征集合Θ'：

5a)令迭代集合Θ等于初步估计集合Θ₀，记Θ＝{θ₁,...,θ_i,...,θ_p}，其中θ_i是Θ中的第i组参数，1≤i≤p，p为迭代集合Θ中参数总组数，令缓存集合Θ₃等于初步估计集合Θ₀，子迭代次数q＝1；

5b)从迭代集合Θ中去掉第q组参数，得到子迭代集合Θ'_q为：

Θ'_q＝{θ₁,...,θ_q-1,θ_q+1,...,θ_p}；

5c)根据子迭代集合Θ'_q，利用属性散射中心模型构建子迭代字典D'_q；

5d)由子迭代字典D'_q得到子迭代字典的空间投影算子

其中(·)^-1表示对矩阵求逆，(·)^H表示对矩阵进行共轭转置；

5e)利用以上参数计算余量相关矩阵C'：

$C^{\&prime;}=D_1{(x,y|{\mathrm{\&Theta;}}_1)}^H\&CenterDot;\mathrm{diag}\left((I-P_{D_q^{\&prime;}})s\right)\&CenterDot;D_2(L,\stackrel{\&OverBar;}{\mathrm{\&phi;}},\mathrm{\&alpha;},\mathrm{\&gamma;}|{\mathrm{\&Theta;}}_2)---<1>$

其中s是目标连通支撑域图像I₂的频域观测信号，diag(·)表示对角化操作，I是单位矩阵；

5f)取出余量相关矩阵C'中模值最大元素对应的行数n'和列数m'，得到一组子参数其中Θ₁(n')表示位置集合Θ₁中第n'组参数，Θ₂(m')表示属性集合Θ₂中第m'组参数；

5g)用子参数更新迭代集合Θ中第q组参数，即：

${\mathrm{\&theta;}}_q={\mathrm{\&theta;}}_q^{\&prime;}(x,y,L,\stackrel{\&OverBar;}{\mathrm{\&phi;}},\mathrm{\&alpha;},\mathrm{\&gamma;});$

5h)判断q＜p是否成立，若成立，则令q＝q+1，返回执行5b)；若不成立，则设q＝1，执行步骤5i)；

5i)判断Θ₃＝Θ是否成立，若不成立，则令Θ₃＝Θ，返回步骤5b)；若成立，则此时的迭代集合Θ就是目标属性散射中心的特征集合Θ'，即Θ'＝Θ。

本发明与现有技术相比具有以下优点：

1.解决了属性散射中心特征提取方法的模型失配问题：

现有的属性散射中心特征提取方法存在模型失配问题，并且由于图像分割会导致目标的一些重要散射中心易于丢失。本发明提出的高分辨率SAR图像目标的属性散射中心特征提取方法从频域观测数据出发，通过求解0范数优化问题得到频域观测数据的稀疏表示，可以有效提取目标散射中心；

2.实现了对属性散射中心参数的自动初始化：

现有的属性散射中心特征提取方法存在参数初始化范围不能自动选取的问题，无法满足目标识别的自动化要求。本发明方法通过结合目标SAR图像的图像域信息，分割出目标的支撑连通区域，实现了对参数的自动初始化，满足了雷达目标识别的自动化要求，并且通过图像分割减小了噪声对属性散射中心特征提取的影响；

3.实现了属性散射中心参数联合估计：

现有的属性散射中心特征提取方法由于属性散射中心参数空间维度过高，采用参数分布估计的方法提取属性散射中心参数集合，降低了属性散射中心参数的估计准确度。本发明方法基于属性散射中心模型，通过拆分属性散射中心冗余字典，降低了字典维度，实现了属性散射中心参数联合估计，提高了属性散射中心参数的估计准确度；

4.能得到目标及其重要部件的精确几何尺寸特征：

由于本发明方法能够有效提取目标属性散射中心参数，所以本发明根据提取的属性散射中心参数集合，可以得到目标及其重要部件的精确几何尺寸特征。

附图说明

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

图2是本发明使用的MSTAR数据库中T72坦克的原始雷达图像；

图3是图2的局部极大值区域图像；

图4是用本发明方法得到的目标连通支撑域图像。

具体实施方式

以下结合附图对本发明的技术方案和效果做进一步详细描述。

一、技术原理

高分辨率SAR图像是目标的二维散射图像，可近似看成由少量强散射中心的后向散射回波之和产生的。目前提取目标散射中心的模型主要分为三种：理想点散射中心模型，衰减指数模型，以及属性散射中心模型。其中属性散射中心模型最符合SAR目标的散射特性，不仅描述了散射中心的位置信息和强度信息，还描述了散射中心的几何尺寸，以及散射中心对频率和方位的依赖性，并可以由此反演出目标的几何结构。属性散射中心模型相比于点散射模型和衰减指数模型，包含更丰富的可用于目标分类识别的特征。它描述的属性散射中心具有丰富的几何意义和物理意义。

根据属性散射中心模型可知，目标中第i个散射中心频率-方位二维回波信号为：

$\begin{array}{c}{E}_i(f,\mathrm{\&phi;};{\mathrm{\&theta;}}_i)=A_i\&CenterDot;{\left(j\frac{f}{f_c}\right)}^{{\mathrm{\&alpha;}}_i}\exp \left(\frac{-j4\mathrm{\&pi;f}}{c}(x_i\cos \mathrm{\&phi;}+y_i\sin \mathrm{\&phi;})\right)\&CenterDot;\\ \sin c\left(\frac{2\mathrm{\&pi;f}}{c}{L}_i\sin (\mathrm{\&phi;}-\stackrel{\&OverBar;}{{\mathrm{\&phi;}}_i})\right)\exp (-2\mathrm{\&pi;f}{\mathrm{\&gamma;}}_i\sin \mathrm{\&phi;})\end{array}---<2>$

其中，i表示散射中心序号，f为雷达发射信号频率，φ为雷达方位角，exp(·)为自然指数函数，sinc(·)为辛克函数，c为光速，θ_i表示第i个散射中心的参数向量，包含A_i为散射中心的散射强度，x_i为距离维坐标，y_i为方位维坐标，L_i为分布式散射中心方位维的长度，为分布式散射中心的方位角，α_i为频率依赖因子，一般α_i∈{-1,-0.5,0,0.5,1}，γ_i为局部式散射中心的方位依赖因子。

由各个散射中心的回波信号之和，即可构成目标频率-方位二维回波信号：

$E(f,\mathrm{\&phi;};\mathrm{\&Theta;})=\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}{E}_i(f,\mathrm{\&phi;};{\mathrm{\&theta;}}_i),{\mathrm{\&Theta;}}^T=[{\mathrm{\&theta;}}_i^T,...,{\mathrm{\&theta;}}_i^T,...,{\mathrm{\&theta;}}_M^T]---<3>$

其中i表示散射中心序号，E_i(f,φ；θ_i)为第i个散射中心的回波信号，M为散射中心个数，Θ表示M个散射中心参数矩阵，(·)^T表示转置操作；

将目标回波信号用矩阵形式表达，其表达式为：

s＝D(Θ)σ+n <4>

其中s为目标回波信号E(f,φ；Θ)的列向量化，D(Θ)为散射中心参数矩阵Θ对应的字典，σ为散射系数向量，n表示噪声。

在雷达回波中，由于目标散射场的绝大部分能量仅由少量强散射中心贡献，故说明雷达回波在属性散射中心的参数空间具有很强的稀疏性。根据雷达回波在属性散射中心的参数空间的稀疏性，通过求解0范数优化问题即可得到观测数据s的稀疏表示，以及散射中心参数估计集合，通过散射中心参数集合可以估计目标及其重要部件的几何尺寸特征。

由于属性散射中心参数空间维数较高，导致构造的冗余字典维数远远大于实际可处理的维数，因此可通过分析属性散射中心模型对字典进行简化，降低参数空间维数。

二、实现步骤

参照图1，本发明的实现步骤如下：

步骤1，目标连通支撑域提取：

从MSTAR数据库中提取T72坦克原始图像I₀，如图2所示；提取原始图像I₀的局部极大值区域，得到目标图像I₁，如图3所示；对目标图像I₁进行分割，得到目标的轮廓区域R，用轮廓区域R对原始图像I₀进行掩膜处理，得到目标连通支撑域图像I₂，如图4所示。

步骤2，确定目标散射中心的位置集合Θ₁和散射中心的属性集合Θ₂。

2a)根据目标的轮廓区域R中的非零元素位置得到散射中心距离维坐标参数x的取值范围X和散射中心方位维坐标参数y的取值范围Y，确定散射中心的位置集合Θ₁，Θ₁＝{(x,y)|x∈X,y∈Y}；

例如，当雷达图像数据为公开MSTAR实测数据时，根据轮廓区域R中的非零元素位置可以得到散射中心距离维坐标参数x的取值范围X＝{x|-2.0700≤x≤1.2937}，散射中心方位维坐标参数y的取值范围Y＝{y|-6.50≤y≤3.12}，所以散射中心的位置集合Θ₁＝{(x,y)|x∈X,y∈Y}；

2b)由散射中心方位维坐标参数y的取值范围确定散射中心长度参数L的取值范围L；由雷达回波数据录取的方位角域确定散射中心方位角取值范围Φ；由雷达回波数据录取的中心频率确定散射中心方位依赖因子γ取值范围Γ；设定散射中心的频率依赖因子α取值范围为Λ＝{-1,-0.5,0,0.5,1}；最终确定散射中心的属性集合Θ₂：

${\mathrm{\&Theta;}}_2=\left\{(L,\stackrel{\&OverBar;}{\mathrm{\&phi;}},\mathrm{\&alpha;},\mathrm{\&gamma;})\right|L\&Element;L,\stackrel{\&OverBar;}{\mathrm{\&phi;}}\&Element;\mathrm{\&Phi;},\mathrm{\&alpha;}\&Element;\mathrm{\&Lambda;},\mathrm{\&gamma;}\&Element;\mathrm{\&Gamma;}\};$

例如，当雷达图像数据为公开MSTAR实测数据时，散射中心长度参数L的取值范围L＝{L|0≤L≤9.8800}，散射中心方位角取值范围散射中心方位依赖因子γ取值范围Γ＝{γ|0≤γ≤1.0417e-9}，频率依赖因子α取值范围为Λ＝{-1,-0.5,0,0.5,1}，所以散射中心的属性集合 ${\mathrm{\&Theta;}}_2=\left\{(L,\stackrel{\&OverBar;}{\mathrm{\&phi;}},\mathrm{\&alpha;},\mathrm{\&gamma;})\right|L\&Element;L,\stackrel{\&OverBar;}{\mathrm{\&phi;}}\&Element;\mathrm{\&Phi;},\mathrm{\&alpha;}\&Element;\mathrm{\&Lambda;},\mathrm{\&gamma;}\&Element;\mathrm{\&Gamma;}\};$

步骤3，根据所述的位置集合Θ₁和属性集合Θ₂，利用属性散射中心模型分别构建位置信息字典D₁(x,y|Θ₁)和属性信息字典

3a)输入位置集合Θ₁和属性集合Θ₂；

3b)由位置集合Θ₁和属性集合Θ₂，分别产生位置原子d_u(f,φ)和属性原子d'_l(f,φ)；

$d_u(f,\mathrm{\&phi;})=\mathrm{vec}\left(\exp \left(\frac{-j4\mathrm{\&pi;f}}{c}(x_u\cos \mathrm{\&phi;}+y_u\sin \mathrm{\&phi;})\right)\right),u=1,...,N_1$

$d_l^{\&prime;}(f,\mathrm{\&phi;})=\mathrm{vec}({\left(j\frac{f}{f_c}\right)}^{{\mathrm{\&alpha;}}_l}\&CenterDot;\sin c\left(\frac{2\mathrm{\&pi;f}}{c}{L}_l\sin (\mathrm{\&phi;}-\stackrel{\&OverBar;}{{\mathrm{\&phi;}}_l})\right)\&CenterDot;\exp (-2\mathrm{\&pi;f}{\mathrm{\&gamma;}}_l\sin \mathrm{\&phi;}))---<5>$

l＝1,...,N₂

其中，vec(·)表示列向量化操作，N₁表示位置集合Θ₁的参数总组数，N₂表示属性集合Θ₂的参数总组数，exp(·)为自然指数函数，sinc(·)为辛克函数，(x_u,y_u)为位置集合Θ₁的第u组参数，为属性集合Θ₂的第l组参数，f为雷达发射信号频率，f_c为雷达发射信号中心频率，φ为雷达波束方位角，c为光速；

3c)由位置原子d_u(f,φ)和属性原子d'_l(f,φ)，得到归一化的位置原子和归一化的属性原子为：其中，||·||₂为2范数算子；

3d)由归一化的位置原子和归一化的属性原子分别构建位置信息字典D₁(x,y|Θ₁)和属性信息字典为：

$\begin{array}{c}{D}_1(x,y|{\mathrm{\&Theta;}}_1)=[{\hat{d}}_1,...,{\hat{d}}_u,...,{\hat{d}}_{N_1}]\\ D_2(L,\stackrel{\&OverBar;}{\mathrm{\&phi;}},\mathrm{\&alpha;},\mathrm{\&gamma;}|{\mathrm{\&Theta;}}_2)=[{\hat{d}}_1^{\&prime;},...,{\hat{d}}_l^{\&prime;},...,{\hat{d}}_{N_2}^{\&prime;}]\end{array}.---<6>$

步骤4，利用正交匹配追踪OMP法，得到目标属性散射中心参数的初步估计集合Θ₀。

4a)输入位置信息字典D₁(x,y|Θ₁)，属性信息字典以及重构能量误差约束因子ε，ε由目标连通支撑域图像I₂的信噪比确定，本实例取ε＝0.001；

4b)设初始迭代次数h＝1，将信号余量r初始化为频域观测信号s，s是目标连通支撑域图像I₂的频域观测信号，将散射系数向量σ初始化为0，将临时集合Θ'₀初始化为空；

4c)利用以上参数计算相关系数矩阵C： $C=D_1{(x,y|{\mathrm{\&Theta;}}_1)}^H\&CenterDot;\mathrm{diag}\left(r\right)\&CenterDot;D_2(L,\stackrel{\&OverBar;}{\mathrm{\&phi;}},\mathrm{\&alpha;},\mathrm{\&gamma;}|{\mathrm{\&Theta;}}_2),$ 其中，diag(·)表示对角化操作；找出相关系数矩阵C中模值最大元素对应行数n和列数m，得到一组估计参数Θ₁(n)表示位置集合Θ₁中第n组参数，Θ₂(m)表示属性集合Θ₂中第m组参数；

4d)将估计参数添加到临时集合Θ'₀；Θ'₀＝Θ'₀∪{θ_h}；

4e)根据临时集合Θ'₀，利用属性散射中心模型构建初步估计字典D'(Θ'₀)；

4f)由初步估计字典D'(Θ'₀)计算系数向量σ；σ＝pinv(D'(Θ'₀))·s，其中pinv(·)表示伪逆；

4g)由系数向量σ估计重构信号

4h)由重构信号更新信号余量r；

4i)判断||r||₂≤ε是否成立，若不成立，则令h＝h+1，返回步骤4c)；若成立，则此时的临时集合Θ'₀就是初步估计集合Θ₀，即Θ₀＝Θ'₀。

步骤5，对初步估计集合Θ₀进行优化，得到目标属性散射中心的特征集合Θ'。

5a)令迭代集合Θ等于初步估计集合Θ₀，记Θ＝{θ₁,...,θ_i,...,θ_p}，其中θ_i是Θ中的第i组参数，1≤i≤p，p为迭代集合Θ中参数总组数，令缓存集合Θ₃等于初步估计集合Θ₀，子迭代次数q＝1；

5b)从迭代集合Θ中去掉第q组参数，得到子迭代集合Θ'_q为：

Θ'_q＝{θ₁,...,θ_q-1,θ_q+1,...,θ_p}；

5c)根据子迭代集合Θ'_q，利用属性散射中心模型构建子迭代字典D'_q；

5c1)输入子迭代集合Θ'_q；

5c2)由子迭代集合Θ'_q，产生子迭代原子

$\begin{array}{c}{d}_t^0(f,\mathrm{\&phi;})=\mathrm{vec}({\left(j\frac{f}{f_c}\right)}^{{\mathrm{\&alpha;}}_t}\&CenterDot;\exp \left(\frac{-j4\mathrm{\&pi;f}}{c}(x_t\cos \mathrm{\&phi;}+y_t\sin \mathrm{\&phi;})\right)\&CenterDot;\\ \sin c\left(\frac{2\mathrm{\&pi;f}}{c}{L}_t\sin (\mathrm{\&phi;}-\stackrel{\&OverBar;}{{\mathrm{\&phi;}}_t})\right)\&CenterDot;\exp (-2\mathrm{\&pi;f}{\mathrm{\&gamma;}}_t\sin \mathrm{\&phi;}))\end{array}---<7>$

1≤t≤N₃

其中N₃是子迭代集合Θ'_q的参数总组数，为子迭代集合Θ'_q的第t组参数；

5c3)由子迭代原子得到归一化的子迭代原子为：

5c4)由归一化的子迭代原子构建子迭代字典D'_q为：

5d)由子迭代字典D'_q得到子迭代字典的空间投影算子

其中(·)^-1表示对矩阵求逆，(·)^H表示对矩阵进行共轭转置；

5e)利用以上参数计算余量相关矩阵C'：

$C^{\&prime;}=D_1{(x,y|{\mathrm{\&Theta;}}_1)}^H\&CenterDot;\mathrm{diag}\left((I-P_{D_q^{\&prime;}})s\right)\&CenterDot;D_2(L,\stackrel{\&OverBar;}{\mathrm{\&phi;}},\mathrm{\&alpha;},\mathrm{\&gamma;}|{\mathrm{\&Theta;}}_2)---<8>$

其中s是目标连通支撑域图像I₂的频域观测信号，diag(·)表示对角化操作，I是单位矩阵；

5f)取出余量相关矩阵C'中模值最大元素对应的行数n'和列数m'，得到一组子参数其中Θ₁(n')表示位置集合Θ₁中第n'组参数，Θ₂(m')表示属性集合Θ₂中第m'组参数；

5g)用子参数更新迭代集合Θ中第q组参数，即；

${\mathrm{\&theta;}}_q={\mathrm{\&theta;}}_q^{\&prime;}(x,y,L,\stackrel{\&OverBar;}{\mathrm{\&phi;}},\mathrm{\&alpha;},\mathrm{\&gamma;});$

5h)判断q＜p是否成立，若成立，则令q＝q+1，返回执行5b)；若不成立，则设q＝1，执行步骤5i)；

5i)判断Θ₃＝Θ是否成立，若不成立，则令Θ₃＝Θ，返回步骤5b)；若成立，则此时的迭代集合Θ就是目标属性散射中心的特征集合Θ'，即Θ'＝Θ。

本发明的效果通过以下实测数据的实验进一步说明：

1)实验场景：

实验所用的实测数据为公开的MSTAR数据中方位角为80.774185°俯仰角为15°的T72坦克的合成孔径雷达SAR数据，雷达的中心频率f_c＝9.599GHz，带宽B＝591MHz。T72坦克真实几何尺寸为：车体长6.41m，车宽3.52m，炮筒长6.155米，炮筒向前时车身总长9.445米。

2)实验内容：

对于方位角为80.774185°俯仰角为15°的T72坦克MSTAR数据，利用本发明提取目标属性散射中心的特征集合，如表1和表2所示。

表1本发明方法处理T72的炮筒散射中心参数

表2本发明方法处理T72的坦克体散射中心参数

3)实验结果分析：

计算表1中散射中心1和散射中心2端点的空间距离，可得到本发明估计出T72的炮筒总长度为6.24米，

计算表2中散射中心3和散射中心4端点的空间距离，可得到本发明估计出T72的坦克体长度为6.23米，

计算表1和表2中散射中心2和散射中心3端点的空间距离，可得到本发明估计出T72的坦克体宽度为3.27米；计算表1和表2中，散射中心1和散射中心4，计算两散射中心端点的空间距离，可得到本发明估计出T72坦克的总长为9.66米，

结合T72坦克的真实几何尺寸可知，本发明方法估计的目标及其重要部件的几何尺寸相对误差在7.1％以内。

实验结果说明，本发明能够有效地准确提取目标属性散射中心特征，本发明提取的目标属性散射中心特征可以用来精确估计目标及其重要部件的几何尺寸。

Claims (3)

1.一种高分辨率SAR图像目标的属性散射中心特征提取方法，包括如下步骤：

(1)提取原始图像I₀中的局部极大值区域，得到目标图像I₁；对目标图像I₁进行分割，得到目标的轮廓区域R，用轮廓区域R对原始图像I₀进行掩膜处理，得到目标连通支撑域图像I₂；

(2)确定散射中心的位置集合Θ₁和散射中心的属性集合Θ₂；

2a)根据轮廓区域R中的非零元素位置得到散射中心距离维坐标参数x的取值范围X和散射中心方位维坐标参数y的取值范围Y，确定散射中心的位置集合Θ₁，Θ₁＝{(x,y)|x∈X,y∈Y}；

2b)由散射中心方位维坐标参数y的取值范围确定散射中心长度参数L的取值范围L；由雷达回波数据录取的方位角域确定散射中心方位角取值范围Φ；由雷达回波数据录取的中心频率确定散射中心方位依赖因子γ取值范围Γ；设定散射中心的频率依赖因子α取值范围为Λ＝{-1,-0.5,0,0.5,1}；最终确定散射中心的属性集合Θ₂：

(3)根据所述的位置集合Θ₁和属性集合Θ₂，利用属性散射中心模型分别构建位置信息字典D₁(x,y|Θ₁)和属性信息字典

(4)利用正交匹配追踪OMP法得到目标属性散射中心参数的初步估计集合Θ₀：

4a)输入位置信息字典D₁(x,y|Θ₁)，属性信息字典以及重构能量误差约束因子ε，ε由目标连通支撑域图像I₂的信噪比确定；

4b)设初始迭代次数h＝1，将信号余量r初始化为频域观测信号s，s是目标连通支撑域图像I₂的频域观测信号，将散射系数向量σ初始化为0，将临时集合Θ'₀初始化为空；

4c)利用以上参数计算相关系数矩阵C：

$C=D_1{(x,y|{\mathrm{\&Theta;}}_1)}^H\&CenterDot;diag\left(r\right)\&CenterDot;D_2(L,\stackrel{\&OverBar;}{\mathrm{\&phi;}},\mathrm{\&alpha;},\mathrm{\&gamma;}|{\mathrm{\&Theta;}}_2)---<4>$

找出相关系数矩阵C中模值最大元素对应行数n和列数m，得到一组估计参数Θ₁(n)表示位置集合Θ₁中第n组参数，Θ₂(m)表示属性集合Θ₂中第m组参数；

4d)将估计参数添加到临时集合Θ'₀；Θ'₀＝Θ'₀∪{θ_h}；

4e)根据临时集合Θ'₀，利用属性散射中心模型构建初步估计字典D'(Θ'₀)；

4f)由初步估计字典D'(Θ'₀)计算系数向量σ；σ＝pinv(D'(Θ'₀))·s，其中pinv(·)表示伪逆；

4g)由系数向量σ估计重构信号

4h)由重构信号更新信号余量r；

4i)判断||r||₂≤ε是否成立，若不成立，则令h＝h+1，返回步骤3c)；若成立，则此时的临时集合Θ'₀就是初步估计集合Θ₀，即Θ₀＝Θ'₀；

(5)对初步估计集合Θ₀进行优化，得到目标属性散射中心的特征集合Θ'：

5a)令迭代集合Θ等于初步估计集合Θ₀，记Θ＝{θ₁,...,θ_i,...,θ_p}，其中θ_i是Θ中的第i组参数，1≤i≤p，p为迭代集合Θ中参数总组数，令缓存集合Θ₃等于初步估计集合Θ₀，子迭代次数q＝1；

5b)从迭代集合Θ中去掉第q组参数，得到子迭代集合Θ'_q为：

Θ'_q＝{θ₁,...,θ_q-1,θ_q+1,...,θ_p}；

5c)根据子迭代集合Θ'_q，利用属性散射中心模型构建子迭代字典D'_q；

5d)由子迭代字典D'_q得到子迭代字典的空间投影算子

其中(·)^-1表示对矩阵求逆，(·)^H表示对矩阵进行共轭转置；

5e)利用以上参数计算余量相关矩阵C'：

$C^{\&prime;}=D_1{(x,y|{\mathrm{\&Theta;}}_1)}^H\&CenterDot;diag\left(\right(I-P_{D_q^{\&prime;}}\left)s\right)\&CenterDot;D_2(L,\stackrel{\&OverBar;}{\mathrm{\&phi;}},\mathrm{\&alpha;},\mathrm{\&gamma;}|{\mathrm{\&Theta;}}_2)---<1>$

其中s是目标连通支撑域图像I₂的频域观测信号，diag(·)表示对角化操作，I是单位矩阵；

5f)取出余量相关矩阵C'中模值最大元素对应的行数n'和列数m'，得到一组子参数其中Θ₁(n')表示位置集合Θ₁中第n'组参数，Θ₂(m')表示属性集合Θ₂中第m'组参数；

5g)用子参数更新迭代集合Θ中第q组参数，即：

${\mathrm{\&theta;}}_q={\mathrm{\&theta;}}_q^{\&prime;}(x,y,L,\stackrel{\&OverBar;}{\mathrm{\&phi;}},\mathrm{\&alpha;},\mathrm{\&gamma;});$

5h)判断q<p是否成立，若成立，则令q＝q+1，返回执行5b)；若不成立，则设q＝1，执行步骤5i)；

5i)判断Θ₃＝Θ是否成立，若不成立，则令Θ₃＝Θ，返回步骤5b)；若成立，则此时的迭代集合Θ就是目标属性散射中心的特征集合Θ'，即Θ'＝Θ。

2.根据权利要求1所述的方法，其中步骤(3)所述的根据位置集合Θ₁和属性集合Θ₂，利用属性散射中心模型分别构建位置信息字典D₁(x,y|Θ₁)和属性信息字典按如下步骤进行：

2a)输入位置集合Θ₁和属性集合Θ₂；

2b)由位置集合Θ₁和属性集合Θ₂，分别产生位置原子d_u(f,φ)和属性原子d'_l(f,φ)；

$d_u(f,\mathrm{\&phi;})=vec\left(\exp \right(\frac{-j4\mathrm{\&pi;}f}{c}\left(x_{u}cos\mathrm{\&phi;}+y_u\sin \mathrm{\&phi;}\right)\left)\right),u=1,...,N_1$

$d_l^{\&prime;}(f,\mathrm{\&phi;})=vec({\left(j\frac{f}{f_c}\right)}^{{\mathrm{\&alpha;}}_l}\&CenterDot;\sin c(\frac{2\mathrm{\&pi;}f}{c}{L}_{l}sin\left(\mathrm{\&phi;}-\stackrel{\&OverBar;}{{\mathrm{\&phi;}}_l}\right))\&CenterDot;\exp (-2{\mathrm{\&pi;f\&gamma;}}_{l}sin\mathrm{\&phi;}\left)\right)---<2>$

l＝1,...,N₂

其中，vec(·)表示列向量化操作，N₁表示位置集合Θ₁的参数总组数，N₂表示属性集合Θ₂的参数总组数，exp(·)为自然指数函数，sinc(·)为辛克函数，(x_u,y_u)为位置集合Θ₁的第u组参数，为属性集合Θ₂的第l组参数，f为雷达发射信号频率，f_c为雷达发射信号中心频率，φ为雷达波束方位角，c为光速；

2c)由位置原子d_u(f,φ)和属性原子d'_l(f,φ)，得到归一化的位置原子和归一化的属性原子为：其中，||·||₂为2范数算子；

2d)由归一化的位置原子和归一化的属性原子分别构建位置信息字典D₁(x,y|Θ₁)和属性信息字典为：

$\begin{array}{c}{D}_1(x,y|{\mathrm{\&Theta;}}_1)=\&lsqb;{\hat{d}}_1,...,{\hat{d}}_u,...,{\hat{d}}_{N_1}\&rsqb;\\ D_2(L,\stackrel{\&OverBar;}{\mathrm{\&phi;}},\mathrm{\&alpha;},\mathrm{\&gamma;}|{\mathrm{\&Theta;}}_2)=\&lsqb;{\hat{d}}_1^{\&prime;},...,{\hat{d}}_l^{\&prime;},...,{\hat{d}}_{N_2}^{\&prime;}\&rsqb;\end{array}---<3>.$

3.根据权利要求1所述的方法，其中步骤5c)所述的根据子迭代集合Θ'_q，利用属性散射中心模型构建子迭代字典D'_q，按如下步骤进行：

5c1)输入子迭代集合Θ'_q；

5c2)由子迭代集合Θ'_q，产生子迭代原子

$\begin{array}{c}{d}_t^0(f,\mathrm{\&phi;})=vec({\left(j\frac{f}{f_c}\right)}^{{\mathrm{\&alpha;}}_t}\&CenterDot;\exp (\frac{-j4\mathrm{\&pi;}f}{c}(x_t\cos \mathrm{\&phi;}+y_t\sin \mathrm{\&phi;}))\&CenterDot;\\ \sin c\left(\frac{2\mathrm{\&pi;}f}{c}{L}_t\sin \right(\mathrm{\&phi;}-{\stackrel{\&OverBar;}{\mathrm{\&phi;}}}_t))\&CenterDot;\exp (-2{\mathrm{\&pi;f\&gamma;}}_t\sin \mathrm{\&phi;}))\\ 1\&le;t\&le;N_3\end{array}---<5>$

其中N₃是子迭代集合Θ'_q的参数总组数，为子迭代集合Θ'_q的第t组参数；

5c3)由子迭代原子得到归一化的子迭代原子为：

${\hat{d}}_t^0=\frac{d_t^0(f,\mathrm{\&phi;})}{\left|\right|d_t^0(f,\mathrm{\&phi;})|{|}_2};$

5c4)由归一化的子迭代原子构建子迭代字典D'_q为：

$D_q^{\&prime;}=\&lsqb;{\hat{d}}_1^0,...,{\hat{d}}_t^0,...,{\hat{d}}_{N_3}^0\&rsqb;.$