[发明专利]一种图像熵最优的压缩传感SAR稀疏自聚焦成像方法

摘要

本发明公开了一种图像熵最优的压缩传感SAR稀疏自聚焦成像方法，它是针对SAR回波信号测量模型中方位向相位误差对压缩传感SAR成像的影响，以及压缩传感SAR成像模型中的未知相位误差估计与补偿问题，利用成像模型中方位向相位误差特性和稀疏目标特征，采用压缩传感SAR图像熵作为评价准则，在每一次迭代处理过程中利用SAR方位向回波与观测目标的关系估计方位向相位误差，然后对压缩传感成像模型进行相位误差补偿，接着再进行压缩传感SAR成像，并采用逐次迭代方法使得压缩传感SAR成像的图像熵最优，从而提高了压缩传感SAR成像质量。

主权项

一种图像熵最优的SAR稀疏自聚焦成像方法，其特征是它包括如下步骤：步骤1、初始化SAR系统参数：初始化SAR系统参数包括：平台速度矢量，记做V；天线初始位置矢量，即方位向0慢时刻天线位置，记做P(0)；雷达工作中心频率，记做f_c；雷达载频波长，记做λ；雷达发射基带信号的信号带宽，记做B_r；雷达发射信号脉冲宽度，记做T_P；雷达发射信号的调频斜率，记做f_dr；雷达脉冲重复频率，记做PRF；方位向等效天线长度，记做D_a；雷达接收系统的采样频率，记做f_s；光在空气中的传播速度，记做C；距离向快时刻序列，记做t，t＝1,2,…,N_R，N_R为距离向快时刻总数；方位向慢时刻序列，记做l，l＝1,2,…,N_A，N_A为方位向慢时刻总数；上述参数均为SAR系统标准参数，在SAR系统设计和观测过程中已经确定；根据SAR成像系统方案和观测方案，SAR成像方法需要的初始化系统参数均为已知；步骤2、初始化SAR成像空间参数以及获取原始回波信号：初始化SAR的成像空间参数，包括：以雷达波束照射场区域地平面所构成的空间坐标作为SAR的成像空间，该成像空间记为Ω；平台到SAR成像空间中心的参考斜距，记做R_ref；将成像空间Ω均匀划分成大小相等的平面单元格，也称为分辨单元，单元网格在水平横向、水平纵向边长分别记为d_x和d_y，单元网格大小选择为SAR系统传统理论成像分辨率或SAR系统传统理论成像分辨率的二分之一；观测场景目标空间Ω中第m个单元格的坐标矢量，记做P_m，m表示SAR成像空间Ω中第m个单元格，m＝1,2,…,M，M为成像空间Ω中的单元格总数；SAR成像空间Ω中所有单元格的散射系数按位置顺序排列组成向量，记做α，向量α由M行1列组成；散射系数向量α中第m个元素，记做α_m，m＝1,2,…,M；根据步骤1中初始化的平台速度矢量V，天线初始位置矢量P(0)和雷达系统的脉冲重复频率PRF，采用公式P(l)＝P(0)+V·l/PRF，l＝1,2,…,N_A，计算得到天线在第l个方位向慢时刻的位置矢量，记为P(l)，N_A为步骤1的方位向慢时刻总数；采用公式R(P(l),P_m)＝||P(l)‑P_m||₂，l＝1,2,…,N_A，m＝1,2,…,M，计算得到在第l个方位向慢时刻SAR成像空间Ω中第m个单元格到天线的距离，记为R(P(l),P_m)，其中||·||₂表示定义2中的向量L2范数，P_m为初始化得到成像空间Ω中第m个单元格的坐标矢量，M为初始化的成像空间Ω中单元格总数；采用公式τ_m(l)＝2·R(P(l),P_m)/C，l＝1,2,…,N_A，m＝1,2,…,M，计算得到在第l个方位向慢时刻SAR成像空间Ω中第m个单元格到天线的时间延时，记为τ_m(l)，其中C为步骤1中初始化得到的光在空气中的传播速度；第l个方位向慢时刻和第t个距离向快时刻中SAR天线的原始回波数据记为s(t,l)，t＝1,2,…,N_R，l＝1,2,…,N_A，其中N_R为步骤1中初始化的距离向快时刻总数；在SAR实际成像中，回波数据s(t,l)，t＝1,2,…,N_R，l＝1,2,…,N_A，可由SAR系统中雷达数据接收机提供；步骤3、建立SAR回波信号的线性测量模型：将步骤2中获取所有SAR原始回波信号s(t,l)按顺序排列组成向量，记为回波信号向量S，回波信号向量S由D行1列组成，其中D＝N_A·N_R，N_R为步骤1中初始化的距离向快时刻总数，N_A为步骤1中初始化的方位向慢时刻总数；采用公式t＝1,2,…,N_R，l＝1,2,…,N_A，m＝1,2,…,M，d＝1,2,…,D，计算得到成像空间Ω中第m个单元格在回波信号向量S第d个元素信号对应的时延函数，记为φ_d(m)，其中exp(·)表示e指数运算符号，f_c为步骤1初始化得到的雷达工作中心频率，f_dr为步骤1初始化得到的发射信号调频斜率，τ_m(l)为步骤2得到的在第l个方位向慢时刻SAR成像空间Ω中第m个单元格到天线的时间延时，t为距离向的第t个快时刻，j为虚数单位(即‑1的开根值)，π为圆周率；SAR原始回波信号向量S与成像空间Ω中所有单元格散射系数向量α之间的测量矩阵，记为A；测量矩阵A由SAR成像空间Ω中所有单元格对应的时延函数构成，A为D行M列的二维矩阵，具体表达式为其中，φ₁(1)为成像空间Ω中第1个单元格在回波信号向量S第1个元素信号对应的时延函数，φ₁(2)为成像空间Ω中第2个单元格在回波信号向量S第1个元素信号对应的时延函数，φ₁(M)为成像空间Ω中第M个单元格在回波信号向量S第1个元素信号对应的时延函数，φ₂(1)为成像空间Ω中第1个单元格在回波信号向量S第2个元素信号对应的时延函数，φ₂(2)为成像空间Ω中第2个单元格在回波信号向量S第2个元素信号对应的时延函数，φ₂(M)为成像空间Ω中第M个单元格在回波信号向量S第2个元素信号对应的时延函数，φ_D(1)为成像空间Ω中第1个单元格在回波信号向量S第D个元素信号对应的时延函数，φ_D(2)为成像空间Ω中第2个单元格在回波信号向量S第D个元素信号对应的时延函数，φ_D(M)为成像空间Ω中第M个单元格在回波信号向量S第D个元素信号对应的时延函数，φ₁(1),φ₁(2),…,φ₁(M)分别为成像空间Ω中第1,2,…,M个单元格在回波信号向量S第1个元素信号对应的时延函数向量，φ₂(1),φ₂(2),…,φ₂(M)分别为成像空间Ω中第1,2,…,M个单元格在回波信号向量S第2个元素信号对应的时延函数向量，φ_D(1),φ_D(2),…,φ_D(M)分别为成像空间Ω中第1,2,…,M个单元格在回波信号向量S第D个元素信号对应的时延函数向量；步骤4、对SAR成像空间的散射系数向量进行初步估计：采用公式计算得到SAR成像空间Ω中散射系数向量的初始估计值，记做其中A^H为矩阵A的共轭转置，A为步骤3中得到的SAR测量矩阵，上标H为共轭转置运算符号，S为步骤3中得到的SAR原始回波信号向量；步骤5、初始化稀疏自聚焦成像算法所需的参数：初始化稀疏自聚焦成像算法所需参数，包括：成像处理迭代的最大迭代次数，记做Maxiter；成像处理第i次迭代过程中SAR成像模型中的相位误差向量，记为相位误差向量Φ⁽ⁱ⁾为D行1列的向量，i为自然数，表示为成像处理的第i次迭代，i＝1,2,…,Maxiter，元素为向量Φ⁽ⁱ⁾中的第1个元素，元素为向量Φ⁽ⁱ⁾中的第2个元素，元素为向量Φ⁽ⁱ⁾中的第D个元素，右上角正体符号T为转置运算符号；成像处理第i次迭代过程中第l个方位向慢时刻的相位误差值，记为l＝1,2,…,N_A；成像处理迭代的阈值，记做τ；采用公式W⁽ⁱ⁾＝diag(exp(j·Φ⁽ⁱ⁾))，i＝1,2,…,Maxiter，计算得到成像处理第i次迭代过程中SAR测量模型的初始相位误差矩阵值，记做W⁽ⁱ⁾，i＝1,2,…,Maxiter，其中相位误差矩阵W⁽ⁱ⁾为D行D列的对角矩阵，diag(·)为将向量元素作为对角矩阵中对角元素的运算符号；成像处理第i次迭代过程中SAR成像图像的图像熵，记为L⁽ⁱ⁾，i＝1,2,…,Maxiter；在成像迭代处理之前，成像处理第0次迭代中SAR测量模型相位误差矩阵的初值，记为W⁽⁰⁾，成像处理第0次迭代中SAR成像图像熵的初值，记为L⁽⁰⁾；步骤6、估计SAR方位向的相位误差向量：包括以下步骤：步骤6.1、获取目标回波的响应向量：采用合成孔径雷达标准距离压缩方法对SAR原始回波数据s(t,l)，t＝1,2,…,N_R，l＝1,2,…,N_A，进行距离压缩处理，得到距离压缩后的SAR回波数据，记做s_r(t,l)，t＝1,2,…,N_R，l＝1,2,…,N_A，其中s(t,l)为步骤2得到的第l个方位向慢时刻和第t个距离向快时刻中SAR天线的原始回波数据；采用公式$\mathrm{\β}(l,m)=s_r(\mathrm{ceil}\left(\frac{R(P\left(l\right),P_m)-R_{\mathrm{ref}}}{2\·C\·f_s}\right),l)\·\exp (j\·\frac{4\·\mathrm{\π}\·R(P\left(l\right),P_m)}{\mathrm{\λ}}),$l＝1,2,…,N_A，m＝1,2,…,M，和公式s_l＝[β(l,1),β(l,2),…,β(l,M)]，l＝1,2,…,N_A，计算得到第l个方位向目标回波的响应向量，记为s_l，其中β(l,1)为第l个方位向慢时刻SAR成像空间Ω中第1个单元格对应的β(l,m)，β(l,2)为第l个方位向慢时刻SAR成像空间Ω中第2个单元格对应的β(l,m)，β(l,M)为第l个方位向慢时刻SAR成像空间Ω中第M个单元格对应的β(l,m)，ceil(·)表示上取整运算符号，R_ref为步骤2初始化得到的平台到SAR成像空间中心参考斜距，C为步骤1初始化得到的光在空气中的传播速度，f_s为步骤1初始化得到的雷达接收系统的采样频率，R(P(l),P_m)为步骤2中得到的第l个方位向慢时刻SAR成像空间Ω中第m个单元格到天线的距离，λ为步骤1中初始化得到的波长；令q为自然数，q取值范围为q＝1,2,…,N_A；当q＝1时，采用公式计算得到第q个方位向慢时刻的数据向量，记为X；当q＝N_A时，采用公式计算得到第q个方位向慢时刻的数据向量，记为X；当1＜q＜N_A时，采用公式计算得到第q个方位向慢时刻的数据向量，记为X，向量X的第m个元素记为x_m，其中为成像处理第i次迭代过程中第l个方位向慢时刻的相位误差值，l＝1,2,…,N_A，s_l为第l个方位向目标回波的响应向量；采用公式Y＝s_q计算得到第q个方位向慢时刻的数据向量，记为Y，向量Y的第m个元素记为y_m，其中s_q为l等于q时第l个方位向目标回波的响应向量，即s_q＝s_l；采用公式m＝1,2,…,M，计算得到元素f_m；再采用公式F＝[f₁,f₂,…,f_M]计算得到的数据向量，记为F，其中|·|²为绝对值平方运算符号，f₁为m＝1时对应的元素f_m，f₂为m＝2时对应的元素f_m，f_M为m＝M时对应的元素f_m，x_m为向量X的第m个元素，y_m为向量Y的第m个元素；步骤6.2、利用施密特正交化方法获取单位正交向量：采用公式和m＝1,2,…,M，计算得到元素a_m和b_m，其中x_m为步骤6.1得到的向量X中第m个元素，y_m为步骤6.1得到的向量Y中第m个元素；再采用a＝[a₁,a₂,…,a_M]、b＝[b₁,b₂,…，b_M]，m＝1,2,…,M，计算得到向量a和向量b，a₁为m＝1时对应的元素a_m，a₂为m＝2时对应的元素a_m，a_M为m＝M时对应的元素a_m，b₁为m＝1时对应的元素b_m，b₂为m＝2时对应的元素b_m，b_M为m＝M时对应的元素b_m；利用施密特正交化方法对向量a和向量b进行标准正交化，得到向量a和向量b构成的平面对应的单位正交向量其中为向量的第1个元素，为向量的第2个元素，为向量的第1个元素，为向量的第2个元素；步骤6.3、计算二次型矩阵及参数估计：采用公式$r_1=\frac{{{\tilde{a}}_2}^2+{{\tilde{b}}_2}^2}{{\tilde{a}}_2\·{\tilde{b}}_1-{\tilde{a}}_1\·{\tilde{b}}_2},r_2=\frac{{{\tilde{a}}_1}^2+{{\tilde{b}}_1}^2}{{\tilde{a}}_2\·{\tilde{b}}_1-{\tilde{a}}_1\·{\tilde{b}}_2},r_3=-\frac{{\tilde{a}}_1\·{\tilde{a}}_2+{\tilde{b}}_1\·{\tilde{b}}_2}{{\tilde{a}}_2\·{\tilde{b}}_1-{\tilde{a}}_1\·{\tilde{b}}_2}$和$R=\left[\begin{array}{cc}{r}_1& r_3\\ r_3& r_2\end{array}\right],$计算得到二次型矩阵R，其中为向量的第1个元素，为向量的第2个元素，为向量的第1个元素，为向量的第2个元素，和分别为步骤6.2得到的单位正交向量；采用公式Γ·Λ·Γ^T＝R对二次型矩阵R进行特征值分解，得到二次型矩阵R的特征向量矩阵记为Γ，Γ^T为矩阵Γ的转置，得到二次型矩阵R的特征值矩阵记为Λ；采用公式η＝(Γ·x₀‑I)·R^‑1计算得到参数的估计值，记为η，其中x₀为向量F在向量a和向量b所张成的二维平面的垂足点，I为单位矩阵，R^‑1为二次型矩阵R的逆矩阵；步骤6.4、估计方位向相位误差：采用公式σ＝[a,b]^‑1·(η·R+I)^‑1·x₀计算得到的向量，记为σ，其中向量σ维数为2×1，向量σ中第1个元素记为σ₁，向量σ中第2个元素记为σ₂，a和b为步骤6.2得到的向量，η为步骤6.3得到的估计值向量，R为步骤6.3得到的二次型矩阵，x₀为向量f在向量a和b所张成的二维平面Ω内的垂足点，I为单位矩阵，上标^‑1表示矩阵求逆运算符号；采用公式计算更新成像处理第i次迭代过程中第q个方位向慢时刻的相位误差向量，其中atan(·)为求解正切函数反函数运算符号；将相位误差向量Φ⁽ⁱ⁾中第(q‑1)·N_R+1至q·N_R个元素的值全部赋值为其中N_R为距离向快时刻总数，为相位误差向量Φ⁽ⁱ⁾中的第(q‑1)·N_R+1个元素，为相位误差向量Φ⁽ⁱ⁾中的第q·N_R个元素；步骤6.5、方位向相位误差值逐个估计：对于所有方位向序号q，q＝1,2,…,N_A，采用步骤6.1到步骤6.4逐个方位向直到估计所有方位向相位误差，最终得到方位向的相位误差向量Φ⁽ⁱ⁾；步骤7、压缩传感稀疏成像：采用公式W⁽ⁱ⁾＝diag(exp(j·Φ⁽ⁱ⁾))计算得到成像处理第i次迭代过程中的相位误差矩阵，记做W⁽ⁱ⁾，其中Φ⁽ⁱ⁾为步骤6.5得到的成像处理第i次迭代过程中的相位误差向量，diag(·)为将向量元素作为对角矩阵对角元素的运算符号；采用公式和标准压缩传感稀疏重构方法计算得到成像处理第i次迭代过程中SAR目标成像空间的散射系数向量，记做其中表示求取满足括号中最小值时对应自变量α的最优值，向量S为步骤3得到的SAR数据回波信号向量，表示向量L2范数的平方运算符号，||·||₁表示向量L1范数运算符号；步骤8、计算SAR图像熵：采用公式m＝1,2,…,M，和公式计算得到成像处理第i次迭代过程中SAR成像图像熵，记做L⁽ⁱ⁾，i＝1,2,…,Maxiter，其中为步骤8中得到的散射系数向量的第m个元素，i为成像处理的第i次迭代，i＝1,2,…,Maxiter，M为步骤2初始化得到的SAR成像空间Ω中划分的单元格总数，为向量L2范数的平方运算符号，Σ(·)表示向量元素求和运算符号，log₂(·)表示底数为2的对数运算符号；步骤9、成像算法迭代条件判定：若满足条件||L⁽ⁱ⁾‑L^(i‑1)||₂≥τ且i≤Maxiter，则令成像处理迭代次数i的值加1，然后重复执行步骤6到步骤9，其中L⁽ⁱ⁾为成像处理第i次迭代过程中SAR成像图像熵，L^(i‑1)为成像处理第i‑1次迭代过程中SAR成像图像熵，τ为步骤5中初始化得到的算法迭代阈值，Maxiter为步骤5中初始化得到的成像处理最大迭代次数，||·||₂为向量L2范数运算符号；若满足条件||L⁽ⁱ⁾‑L^(i‑1)||₂＜τ或者i＞Maxiter，则成像处理第i次迭代过程得到的散射系数向量和相位误差向量Φ⁽ⁱ⁾为最终的SAR成像散射系数向量和相位误差估计结果。