[发明专利]一种合成孔径雷达高效自聚焦后向投影BP方法

权利要求书

1.一种合成孔径雷达高效自聚焦后向投影BP方法，其特征是它包括以下几个步骤： 

步骤1、用于低内存的高效自聚焦BP算法相关参数的初始化 

初始化的参数均为己知，且初始化的参数如下：光速为C；雷达发射线性调频信号，雷达发射脉冲波长为λ；雷达发射脉冲的带宽Ｂ；雷达发射脉冲的时宽为T_p；雷达脉冲重复周期为T；雷达回波距离向采样频率为f_s；雷达回波数据矩阵为Ｓ_Ｋ×L；雷达回波数据的方位向点数和距离向点数分别为K和L(K和L均为正整数)，K也称为慢时刻数；慢时刻向量为t_s=[-K／2，1-K／2，...，K／2-1]×T；几何坐标系为三维笛卡尔坐标系O-XYZ;雷达平台速度矢量为V，雷达平台速度矢量V的大小为1行3列；雷达平台在零时刻的位置向量为P_t0，P_t0的大小为1行3列；O-XY平面内的矩形场景为Θ；将Θ离散化为像素点网格，记为Ω_M×N；像素点网格Ω_M×N中X方向和Y方向的像素点点数分别为M和N；像素点网格Ω_M×N中X方向和Y方向的像素点间隔分别为dx和dy；像素点网格Ω_M×N中第m行、第n列像素点在K个慢时刻的粗聚焦后向投影结果向量为I_mn，I_mn的大小为K行1列，m=1，2，...，M，n=1,2，...，N；像素点网格Ω_M×N中第m行、第n列像素点在K个慢时刻的精聚焦后向投影结果向量为I_1mn，I_1mn的大小为K行1列，m=1,2，...，M，n=1,2，，...N；场景中心位置向量为P_c，P_c的大小为1行3列；雷达回波接收波门延迟距离为R₀；升采样数据矩阵为SS_K×P，ss_K×P的大小为K行、8×L列，K为雷达回波数据的方位向点数，L为雷达回波数据的距离向点数；粗聚焦图像矩阵为Ic_M×N，IcMf×N的大小为M行、N列；精聚焦图像矩阵为Image_M×N，Image_Ｍ×N的大小为M行、N列；天线相位中心矩阵为APC_３×K； 

步骤2、对雷达回波数据矩阵的每一行进行脉冲压缩 

取出所有雷达回波数据s_Ｋ×L，利用传统的脉冲压缩方法对雷达回波数据s_Ｋ×L的每一行进行脉冲压缩，得到脉冲压缩后的雷达回波数据矩阵PS_K×L； 

步骤3、对脉冲压缩后的数据矩阵的每一行进行频域升采样 

对步骤2得到的脉冲压缩后的雷达回波数据矩阵PS_K×Ｌ的每一行做如下8倍频域升采样处理： 

步骤3.1、取出步骤2中脉冲压缩后的雷达回波数据矩阵PS_K×Ｌ的第k行向量，记为s_k，k=1，2，...，Ｋ，K为步骤1定义的雷达回波数据的方位向点数； 

步骤3.2、对向量s_k作传统的快速傅里叶变换(FFT)，得到向量f_k； 

步骤3.3、在向量f_k的L／2+1位置插入7×L个零，得到向量z_k， 

z_k=[f_k(1)，f_k(2），...，f_k(L／2),0_1×7L，f_k(L／2+1），...f_k(L)]，f_k(1)为向量f_k中的第1个元素，f_k(2)为向量f_k中的第2个元素，f_k(L／2)为向量g_k中的第L／2个元素，0_1×7L为1行、7×L列的零向量，f_k(L／2+1)为向量f_k中的第L／2+1个元素，f_k(L)为向量f_k中的第L个元素，L为步骤1提供的雷达回波数据的距离向点数，L均为正整数； 

步骤3.4、对向量z_k作传统的快速傅里叶逆变换(IFFT)，得到向量ss_k，将向量ss_k存放到矩阵SS_K×p的第K行，SS_K×P为步骤1提供的升采样数据矩阵； 

步骤4、对整个像素点网格Ω_M×N进行粗聚焦BP成像 

对步骤1定义的像素点网格Ω_M×N中的每个像素点做如下处理： 

步骤4.1、采用公式P1=P_c+[(m-M／2)×dx，(n-N／2)×dy,0]，计算像素点网格Ω_M×N中第m行、第n列像素点的位置向量，记为P1，m=1,2，...，M，n=1,2，...，N，M和N分别为步骤1定义的像素点网格Ω_M×N中X方向和Y方向的像素点点数，P_c为步骤1定义的场景中心位置向量，dx和dy分别为步骤1定义的像素点网格Ω_M×N中X方向和Y方向的像素点间隔； 

步骤4.2、采用公式P_k1=P_t0+V×t_S(k1)计算在第k1个慢时刻的平台位置向量，记为Pk1，k1=1,2，...，K，K为步骤1定义的雷达回波数据的方位向点数，P_t0为步骤1定义的雷达平台在零时刻的位置向量，V为步骤1定义的雷达平台速度矢量，t_s(k1)为向量t_s中第k1个元素，t_s为步骤1定义的慢时刻向量； 

步骤4.3、采用公式R_k1=norm(P_k1-P1)计算步骤4.2计算出的平台位置P_k1与步骤4.1计算出的像素点位置P1的距离，记为R_k1，norm(·)表示向量取模运算； 

步骤4.4、采用公式τ₁=2R_k1／C-2R₀／c计算在第K1个慢时刻从像素点P1反射的回波延时，记为τ₁，k1=1,2，...,K，K为步骤1定义的雷达回波数据的方位向点数，R_K1为步骤4.3计算出的平台位置P_k1与像素点位置P1的距离，C为步骤1定义的光速，R₀为步骤1定义的雷达回波接收波门延迟距离； 

步骤4.5、采用公式id₁=round((τ₁×f_s+L／2)×8)计算在第k1个慢时刻从像素点P1反射的回波在距离向数据中的位置，记为id₁，id₁为正整数，τ₁为步骤4.4计算出的在第k1个慢时刻，从像素点P1反射的回波延时，f_s为步骤1定义的雷达回波距离向采样频率，L为步骤1定义的雷达回波数据的距离向点数，round(·)表示四舍五入运算； 

步骤4.6、取出步骤3计算出的升采样数据矩阵SS_K×P中第k1行、第id₁列的元素，记为s_id1，k1=1,2，...，K，K为步骤1定义的雷达回波数据的方位向点数，id₁为步骤4.5计算出的在第k1个慢时刻，从像素点P1反射的回波在距离向数据中的位置； 

采用公式I_k1=s_id1×exp(j4rrR_k1／λ)，计算在第k1个慢时刻像素点P1的粗聚焦后向投影结果，记为I_k1，k1=1,2，...，k，K为步骤1定义的雷达回波数据的方位向点数，j为虚数单位，R_k1为步骤4.3计算出的平台位置P_k1与像素点位置P1的距离，λ为步骤1定义的雷达发射脉冲的波长，exp(·)表示指数运算；将I_k1存放到向量I_mn的第k1个位置，I_mn为步骤1定义的像素点网格Ω_M×N中第m行、第n列像素点在K个慢时刻的粗聚焦后向投影结果向量； 

步骤4.7、重复步骤4.1～4.6，直到计算出像素点P1在所有K个慢时刻的粗聚焦后向投影结果向量; 

步骤4.8、采用公式计算像素点P1总的粗聚焦后向投影结果，记为II_mn，I_mn(k1)为向量I_mn中第k1个元素，I_mn为步骤1定义的像素点网格Ω_M×N中第m行、第n列像素点在K个慢时刻的粗聚焦后向投影结果向量，m=1,2，...，M，n=1,2，...,N；将II_mn存放到矩阵Ic_M×N中第m行、第n列位置，Ic_M×N为步骤1定义的粗聚焦图像矩阵； 

步骤5、在粗聚焦BP成像结果中选择用于自聚焦的小块场景区域 

首先，利用Matlab软件中的绘图函数imagesc(·）绘制步骤4中获得的粗聚焦图像矩阵Ic_Ｍ×N的幅度图像；然后，采用遍历法搜索出幅度图像中强度大的Q个像素点，并通过人眼观察，从这Q个像素点中选出最为孤立的1个像素点，最为孤立的像素点周围区域的强度很小，记为P₀，以像素点P₀为中心，选择包含像素点P₀散布能量的一个矩形像素点区域，记为Ω_Ｅ×H，矩形像素点区域Ω_Ｅ×H在X坐标方向的起始编号和终止编号分别为IDX_stａｒt和IDX_ｅｎd，矩形像素点区域Ω_E×Ｈ在Y坐标方向的起始编号和终止编号分别为IＤＹ_{ｓｔａｒt}和IDY_ｅｎd，矩形像素点区域Ω_Ｅ×H在X坐标方向和Y坐标方向的像素点点数分别为E和H，Ｅ=IDX_eｎa-IDＸ_{ｓｔａｒt}+1，H=IＤＹ_ｅnd-IDＹ_{ｓｔａｒt}+1，矩形像素点区域Ω_Ｅ×H内总的像素点点数为W，W=E×H；上述所述的“强度大”和“强度很小”属于图像处理领域普通技术人员均能领会的术语； 

步骤6、对矩形像素点区域Ω_Ｅ×H内所有像素点做后向投影 

步骤6.1、将矩形像素点区域Ω_Ｅ×H中的所有像素点以Y坐标方向为优先顺序进行编号，编号为i,2，...，W，W为步骤5定义的矩形像素点区域Ω_E×H内总的像素点点数；定义：区域后向投影结果矩阵为Ia_W×K，距离历史矩阵为R_W×K，矩形像素点区域Ω_E×H内所有像素点的位置矩阵为Q_W×3，矩形像素点区域Ω_E×H中第q个像素点在K个慢时刻的后向投影结果向量为I_q，q=1,2，...，W； 

对矩形像素点区域Ω_E×H中的每个像素点做如下处理： 

步骤6.2、采用公式P_q=P_c+[(e-M／2)×dx，(h-N／2)×dy，0]计算矩形像素点区域Ω_E×H中第q个像素点的位置向量，记为P_q，q=1,2，...，W，W为步骤5定义的矩形像素点区域Ω_E×H内总的像素点点数，P_c为步骤1定义的场景中心位置向量， 

e=IDX_start+mod((q-1)／W)，IDX_srart为步骤5定义的矩形像素点区域Ω_E×H在X坐标方向的起始编号，W为步骤5定义的矩形像素点区域Ω_E×H内总的像素点点数，mod(·)表示取余运算，M为步骤1定义的像素点网格Ω_M×N中X方向的像素点点数，dx为步骤1定义的像素点网格Ω_M×N中X方向的像素点间隔，h=IDY_start+[(q-1)／W]，IDY_start为步骤5定义的矩形像素点区域Ω_E×H在Y坐标方向的起始编号，W为步骤5定义的矩形像素点区域Ω_E×H内 总的像素点点数，表示向下取整运算，N为步骤1定义的像素点网格Ω_M×N中Y方向的像素点点数，dy为步骤1定义的像素点网格Ω_M×N中Y方向的像素点间隔；将P_q存放到矩阵Q_W×3的第q行，q=1,2，...，W，W为步骤5定义的矩形像素点区域Ω_E×H内总的像素点点数，Q_W×3为步骤6.1定义的矩形像素点区域Ω_E×H内所有像素点的位置矩阵； 

步骤6.3、采用公式P_k2=P_t0+V×t_s(k2)计算在第k2个慢时刻的平台位置向量，记为P_k2，k2=1，2，...,K，K为步骤1定义的雷达回波数据的方位向点数，P_t0为步骤1定义的雷达平台在零时刻的位置向量，V为步骤1定义的雷达平台速度矢量，t_s(k2)为向量t_s中第k2个元素，t_s为步骤1定义的慢时刻向量； 

步骤6.4、采用公式R_k2=norrn(P_k2-P_q)计算步骤6.3计算出的平台位置P_k2与步骤6.2计算出的像素点位置P_q的距离，记为R_k2，nonm(·)表示向量取模运算；将R_k2存放到矩阵R_W×K的第q行、第k2列位置，R_W×K为步骤6.1定义的距离历史矩阵，q=1,2，...，W，W为步骤5定义的矩形像素点区域Ω_E×H内总的像素点点数，k2=1，2，...,K，K为步骤1定义的雷达回波数据的方位向点数； 

步骤6.5、采用公式τ₂=2R_k2／C-2R₀／c计算在第k2个慢时刻，从像素点P_q反射的回波延时，记为τ₂，k2=1，2，...，K，K为步骤1定义的雷达回波数据的方位向点数，R_k2为步骤6.4计算出的平台位置P_k2与像素点位置P_q的距离，C为步骤1定义的光速，R₀为步骤1定义的雷达回波接收波门延迟距离； 

步骤6.6、采用公式id₂=round((τ₂×f_s+L／2)×8)计算在第k2个慢时刻从像素点P_q反射的回波在距离向数据中的位置，记为id₂，id₂为正整数，τ₂为步骤6.5计算出的在第k2个慢时刻，从像素点P_q反射的回波延时，f_s为步骤1定义的雷达回波距离向采样频率，L为步骤1定义的雷达回波数据的距离向点数，round(·)表示四舍五入运算； 

步骤6.7、取出步骤3计算出的升采样数据矩阵SS_K×P中第k2行、第id₂列的元素，记为s_id2，k2=1,2，...，k，K为步骤1定义的雷达回波数据的方位向点数，id₂为步骤6.6计算出的在第 k2个慢时刻从像素点P_q反射的回波在距离向数据中的位置；采用公式I_k2=s_id2×exp（j4πR_k2／λ)计算在第k2个慢时刻像素点P_q的后向投影结果，记为Ｉ_k2，k2=12，...,Ｋ，Ｋ为步骤1定义的雷达回波数据的方位向点数，j为虚数单位，R_k2为步骤6.4计算出的平台位置P_k2与像素点位置P_q的距离，λ为步骤1定义的雷达发射脉冲的波长，exp(·)表示指数运算；将Ｉ_k2存放到向量I_q的第k2个位置，I_q为步骤6.1定义的矩形像素点区域Ω_E×H中第q个像素点在K个慢时刻的后向投影结果向量，k2=1,2，...，K，K为步骤1定义的雷达回波数据的方位向点数； 

步骤6.8、重复步骤6.2～6.7，直到像素点P_q在所有K个慢时刻的后向投影结果向量都计算完毕；将像素点P_q在所有K个慢时刻的后向投影结果向量I_q存放到矩阵Ia_W×K的第q行，Ia_W×K为步骤6.1定义的后向投影结果矩阵，q=1,2，...,W，W为步骤5定义的矩形像素点区域Ω_E×H内总的像素点点数； 

步骤7、对矩形像素点区域Ω_E×H内的所有像素点做基于图像强度的BP自聚焦处理 

采用传统的基于图像强度的自聚焦BP算法，对步骤6中计算出的矩形像素点区域Ω_E×H内所有像素点的后向投影结果矩阵Ia_W×K做BP自聚焦处理，得到1行、K列的相位误差向量Error_1×K； 

步骤8、计算矩形像素点区域Ω_E×H内的所有像素点在所有慢时刻的距离 

定义W行、1列且元素全为1的向量为e_W×1，采用公式 

R1_W×K=(R_W×K-e_W×1×Error_1×K×λ／π／2)／2计算矩形像素点区域Ω_E×H内的所有像素点在所有慢时刻的距离矩阵，记为R1_W×K，R_W×K为步骤6计算出的距离历史矩阵，Error_1×K为步骤7计算出的相位误差向量，λ为步骤1定义的雷达发射脉冲波长； 

步骤9、计算天线相位中心 

对步骤8中得到的距离矩阵R1_W×K的每一列做如下处理： 

步骤9.1、采用公式P_k3=P_t0+V×t_s(k3)计算在第k3个慢时刻的平台位置向量，记为P_k3，k3=1,2，...，K，K为步骤1定义的雷达回波数据的方位向点数，P_t0为步骤1定义的雷 达平台在零时刻的位置向量，V为步骤1定义的雷达平台速度矢量，t_s(k3)为向量t_s中第k3个元素，t_s为步骤1定义的慢时刻向量； 

步骤9.2、取出距离矩阵R1_W×K中的第k3列向量，记为R_new； 

步骤9.3、以作为目标函数，P_k3为步骤9.1计算出的在第k3个慢时刻的平台位置向量，k3=1,2，...，K，K为步骤1定义的雷达回波数据的方位向点数，Q_W×3(p，：)为矩阵Q_W×3的第p行向量，p=1，...，W，W为步骤5定义的矩形像素点区域Ω_E×H内总的像素点点数，Q_W×3为步骤6计算出的矩形像素点区域Ω_E×H内所有像素点的位置矩阵，norm(·)为向量求模运算，R_new(p)为R_new的第p个元素，R_new为步骤9.2取出的向量，采用坐标下降法，求解出第k3个慢时刻的天线相位中心向量，记为apc_3×1，并将向量apc_3×1存放到矩阵APC_3×K的第k3列，APC_3×K为步骤1定义的天线相位中心矩阵； 

步骤10、对整个像素点网格Ω_M×NN进行精聚焦BP成像 

对步骤1定义的像素点网格Ω_M×N中的每个像素点做如下处理： 

步骤10.1、采用公式P4=P_c+[(m-M／2)×dx，（n-N／2)×dy，0]计算像素点网格Ω_M×N中第m行、第n列像素点的位置向量，记为P4，m=1,2，...，M，n=1,2，...，N，M和N分别为步骤1定义的像素点网格Ω_Ｍ×N中X方向和Y方向的像素点点数，P_c为步骤1定义的场景中心位置向量，dx和dy分别为步骤1定义的Ω_M×N中X方向和Y方向的像素点间隔； 

步骤10.2、采用公式R_k4=norm(P_k4-P4)计算在慢时刻k4的平台位置P_k4与步骤10.1计算出的像素点位置P4的距离，记为R_k4，P_k4为步骤9中得到的天线相位中心矩阵APC_3×K的第k4列向量，K4=1,2，...,K，K为步骤1定义的雷达回波数据的方位向点数，norm(·)表示向量求模运算； 

步骤10.3、采用公式τ₄=2R_k4／c-2R₀／C计算在第k4个慢时刻从像素点P4反射的回波延时，记为τ₄，k4=1,2，...，k，K为步骤1定义的雷达回波数据的方位向点数，R_k4为步骤 10.2计算出的平台位置P_k4与像素点位置P4的距离，C为步骤1定义的光速，R₀为步骤1定义的雷达回波接收波门延迟距离； 

步骤10.4、采用公式id₄=round((τ₄×f_s+L／2)×8)计算在第k4个慢时刻从像素点P4反射的回波在距离向数据中的位置，记为id₄，id₄为正整数，τ₄为步骤10.3计算出的在第k4个慢时刻，从像素点P4反射的回波延时，f_s为步骤1定义的雷达回波距离向采样频率，L为步骤1定义的雷达回波数据的距离向点数，round(·)表示四舍五入运算； 

步骤10.5、取出步骤3计算出的升采样数据矩阵ss_K×P中第k4行、第id₄列的元素，记为s_id4，k4=1,2，...，K，K为步骤1定义的雷达回波数据的方位向点数，id₄为步骤10.5计算出的在第k4个慢时刻，从像素点P4反射的回波在距离向数据中的位置；采用公式I_k4=s_id4×exp（j4πR_k4／λ)计算在第k4个慢时刻，像素点P4的精聚焦后向投影结果，记为I_k4，k4=1,2，...，K，K为步骤1定义的雷达回波数据的方位向点数，j为虚数单位，R_k4为步骤10.2计算出的平台位置P_k4与像素点位置P4的距离，λ为步骤1定义的雷达发射脉冲的波长，exp(·)表示指数运算；将I_k4存放到向量I_1mn的第k4个位置，I_1mn为步骤1定义的像素点网格Ω_M×N中第m行、第n列像素点在K个慢时刻的精聚焦后向投影结果向量； 

步骤10.6、重复步骤10.1～10.5，直到像素点P4在所有K个慢时刻的精聚焦后向投影结果都计算完毕； 

步骤10.7、采用公式计算像素点P4总的后向投影结果，记为II_1mn，I_1mn(k4)为向量I_１mn中第k4个元素，I_1mn为步骤1定义的像素点网格Ω_M×N中第m行、第n列像素点在K个慢时刻的精聚焦后向投影结果向量，m=1,2，...，M，n=1,2，...，Nr；将II_1mn存放到矩阵Image_M×Ｎ中第m行、第n列位置，Image_M×Ｎ为步骤1定义的精聚焦图像矩阵。