基于雷达长时间观测的天体目标转动估计与三维重建方法
技术领域
本发明涉及基于雷达长时间观测的天体目标转动估计与三维重建方法，属于逆合成孔径雷达三维成像领域。
背景技术
逆合成孔径雷达(ISAR)可以对非合作目标进行全天时、全天候的成像。其中，三维ISAR可以对目标的三维结构进行重构，获取目标丰富的特征信息，具有广泛的应用前景。
三维ISAR成像方法可以分为两类，即干涉逆合成孔径雷达(InISAR)方法和单天线ISAR三维重构方法。其中，InISAR方法依赖于多基地或多通道系统，此方法可以在高度维上形成合成孔径，从而获得高度维的分辨率。但是，受到不同雷达或通道之间同步精度与相干性的影响，此方法需投入较高的硬件成本。单天线ISAR三维重构方法通过目标自身的运动获得高度维的合成孔径，此方法无需考虑天线阵列的布局，硬件复杂度较低。基于单天线的ISAR图像序列，通过分析目标三维结构在成像平面的投影方程，可以实现三维重建，但这仅限于运动已知的合作目标。此外，基于光学图像处理的三维重建理论，可以依据ISAR图像序列中特显点的航迹矩阵实现三维重建，但此方法需要预先对ISAR图像进行准确地定标。
近年来，深空探测成为了人们的研究热点。受到地基雷达数量的限制，基于InISAR的天体三维重建方法实用性较低。此外，由于天体距离远、双程衰减严重，回波的信噪比(SNR)极低，ISAR图像序列将受到噪声的严重干扰，难以实现精确的定标与三维重建。因此，目前现有的三维ISAR成像方法在天体观测中的实现具有一定的局限性。
发明内容
有鉴于此，针对低信噪比下天体目标三维ISAR成像难的问题，本发明提出一种基于雷达长时间观测的天体目标转动估计与三维重建方法，可以在低信噪比下实现天体目标的转轴、转速联合估计，同时实现天体目标的三维重建。
1、天体目标转动估计与三维重建方法，其特征在于，包括如下步骤：
步骤一、天体目标运动建模：
在以天体中心为原点的天体惯性坐标系O-XYZ下，天体中心坐标为[0,0,0]T，设雷达的坐标为[x0(t),y0(t),z0(t)]T；
设天体在惯性坐标系下的转轴为单位向量的方向：
其中，a、b和c分别为单位向量在天体惯性坐标系的三个坐标轴的投影；γ为的方向与Z轴正方向的夹角；为惯性坐标系X轴与OM之间的夹角，OM为向量在惯性坐标系的XOY平面的投影；天体上的第p个散射点Qp在以天体为中心的惯性坐标系下的初始直角坐标为[xp0,yp0,zp0]T，p＝1,2,…,P，P表示散射点个数；将上述以天体中心为原点的惯性直角坐标系转换为极坐标系，则第p个散射点Qp对应的初始极坐标可以表示为[αp,βp,Hp]T，其中，αp为通过天体惯性坐标系Z轴和Qp的半平面与天体惯性坐标系XOZ平面所构成的角；βp为OQp与Z轴正方向的夹角；Hp为第p个散射点到原点的距离，即为待搜索的物理量；则第p个散射点对应的天体惯性坐标系下的直角坐标为：
第p个散射点绕转轴以转速ω转动后的坐标为：
其中，r＝[a,b,ω]T为天体转动参数；t为慢时间；
则平动补偿后第p个散射点回波的斜距历程为：
其中，RT(t)为补偿的平动；式(5)中，待搜索的参数为天体转动参数r和高度值Hp；
步骤二：设置参数搜索空间：
设置天体转动参数r的搜索空间和散射点高度的搜索空间；
步骤三：首先，将雷达回波分成若干子孔径，在每个子孔径内利用GRFT获得相参积累结果，然后将子孔径间的结果进行非相干积累，从而获得最终的积累输出值；其中，积累过程表示为：
其中，sr表示脉冲压缩后的回波；Nc表示子孔径个数；Tc表示每个子孔径的时间；λ为波长，c为光速，k为子孔径序号；
然后，在三维参数搜索过程中，在固定每组待搜索转动参数r＝[a,b,ω]T的情况下，对每个散射点的高度值Hp进行P次一维遍历搜索，每次遍历搜索根据公式(6)得到一系列的积累输出值，保留P个积累输出值的最大值与对应的P个高度值；接着，将上述P个输出最大值叠加，得到目标函数s1(r)，表示为：
最后，按照式(7)对参数搜索空间中的转动参数进行启发式寻优，当式(7)的目标函数s1(r)取最大值时，搜索到的转动参数和高度值为目标的真实值；
步骤四：利用获得的P个散射点的高度值对天体的三维形状进行重建。
较佳的，所述步骤四中，三维重建的公式为：
其中，为第p个散射点三维坐标的估计值。
较佳的，所述步骤三中，非相干积累采用基于混合广义瑞登傅里叶变换(HGRFT)方法实现。
较佳的，所述启发式寻优方法采用粒子群优化算法实现。
较佳的，所述转动参数搜索空间为U＝[1h,4h]×[-0.5,0.5]×[-0.5,0.5]，高度搜索空间为UH＝[550m,800m]。
较佳的，高度搜索间隔为0.5m。
本发明具有如下有益效果：
(1)本发明提出了一种基于混合广义瑞登傅里叶变换(HGRFT)的目标转动参数估计与三维重建方法，基于长时间积累，首先对天体绕转轴的自旋运动进行了建模，然后利用HGRFT对回波的能量实现相参与非相参联合的积累；
(2)当运动模型与目标真实运动吻合时，可以获得积累的峰值，并根据峰值在参数搜索空间中的位置快速并准确的估计目标的转轴、转速与散射点的高度；
(3)本发明方法还可以与启发式寻优方法相结合，进一步提升参数的搜索效率。
附图说明
图1表示本发明建立的天体目标三维ISAR几何模型；
图2表示1999KW4主星三维形状；
图3表示网格剖分后的小行星形状；
图4表示信噪比为-10dB下的脉冲压缩结果；
图5表示重建后的小行星形状；
图6表示高度值估计偏差。
具体实施方式
下面结合附图和实施例对本发明方法的实施方式做详细说明。
若雷达发射线性调频信号，则接收到回波的可以表示为：
其中，c为光速；fc为载频；λ＝c/fc为波长；τ为快时间；t为慢时间；Tp为脉冲宽度；kr为线性调频信号的调频率；σp为第p(p＝1,2,…,)个散射点回波的幅度；Rp(t)为第p个散射点的斜距历程。
脉冲压缩后的回波可以表示为
其中，B为信号带宽；Ap＝σp为脉冲压缩后的幅度。
本发明的实施步骤如下：
步骤一、目标运动建模：
目标与雷达间的距离可以根据接收到回波的时间延迟获得，目标相对于雷达的方位角和俯仰角可以根据天线的波束指向确定。则在以天体中心为原点的惯性坐标系下，天体中心坐标为[0,0,0]T，雷达的坐标为[x0(t),0(t),0(t)]T，如图1所示。
设天体在惯性坐标系下的转轴为单位向量的方向， 其中，a、b和c分别为单位向量在三个坐标轴的投影；γ为的方向与Z轴正方向的夹角；为X轴与OM之间的夹角，OM为向量在XOY平面的投影。天体上的第p(p＝1,2,…,P)个散射点在以天体为中心的惯性坐标系下的初始直角坐标为[xp0,yp0,zp0]T，则其对应的初始极坐标可以表示为[αp,βp,Hp]T，其中，αp和βp为第p个散射点在极坐标系下的坐标，需要预先设定；Hp为待搜索的高度。则第p个散射点对应的天体惯性坐标系下的直角坐标为[xp0,yp0,zp0]T
第p个散射点绕转轴以转速ω转动后的坐标为
其中，r＝[a,b,ω]T为天体转动参数，
则平动补偿后第p个散射点回波的斜距历程为
其中，RT(t)为补偿的平动，可以通过现有的参数化平动补偿方法获得。式(6)中，待搜索的参数为天体转动参数r和三维高度Hp(1,2,…,P)。
步骤二：设置参数搜索空间
设置三维参数搜索空间U＝[ωmin,ωmax]×[amin,amax]×[bmin,bmax]用于搜索转动参数r，其中，(amin,bmin)和(amax,bmax)为天体转轴的搜索边界，ωmin和ωmax为天体转速的搜索边界；此外，设置P个相同的一维参数搜索空间UH＝[Hmin,Hmax]用于搜索P个散射点的高度，其中，Hmin和Hmax为高度的搜索边界。
步骤三：基于混合广义瑞登傅里叶变换(HGRFT)，联合遍历与启发式寻优方法搜索参数估计值，具体为：
HGRFT方法首先将回波划分为Nc个子孔径，使得每个子孔径内有Nd个脉冲，每个子孔径的时间Tc＝dTr。其中，混合广义瑞登傅里叶变换HGRFT方法采用申请号为“202011196999.3”，名称为“基于混合广义瑞登傅里叶变换的雷达目标检测方法”专利的方法。在每个子孔径内利用GRFT获得相参积累结果，然后将子孔径间的结果进行非相干积累，从而获得最终的积累输出值。其中，在本发明中，基于HGRFT的积累过程可以表示为：
其中，k为子孔径序号。
在本发明中，参数搜索为三维搜索嵌套一维搜索的过程。由于启发式寻优方法可以加速多维参数搜索过程，因此，三维搜索采用启发式寻优方法，一维搜索采用遍历搜索方法。搜索过程的目标函数为HGRFT的积累输出值，具体的步骤如下。
首先，在三维参数搜索过程中，在固定每组待搜索转动参数r＝[a,b,ω]T的情况下，对每个散射点的高度值Hp进行P次一维遍历搜索，保留P个HGRFT的输出最大值与对应的P个高度值。接着，将上述P个输出最大值叠加，可以表示为
按照式(8)对三维参数搜索空间U中的转动参数进行启发式寻优。当搜索的转动参数和高度值搜索至目标的真实值时，式(8)的输出达到最大值s1(r)。最终估计的转动参数值与P个散射点的高度值Hp可以表示为
步骤四：重建目标的三维形状
获得了P个散射点的高度估计值后天体的三维形状可以重建为
其中，为第p个散射点三维坐标的估计值。
实施例1
为了进一步验证本发明的可行性与有效性，开展了对近地小行星1999KW4主星的转动参数估计与三维重建仿真。
1999KW4主星的旋转周期为2.8h，直径为1.5km，其三维形状如图2所示。以200米为间距对小行星的表面进行剖分，得到的由390个散射点组成的三维形状如图3所示。仿真参数如表1所示
表1：仿真参数表
参数
参数值
载频
2380MHz
带宽
30MHz
采样率
30MHz
脉宽
2us
脉冲重复频率
1Hz
子孔径数
100
积累时间
2000s
雷达坐标
[1000m,1000m,500m]
转轴
(0,0,1)
在参数搜索过程中，采用粒子群优化算法(PSO)作为启发式寻优方法，其中，PSO粒子数设置为100，最大迭代次数设置为500，迭代门限设置为10。三维参数搜索空间为U＝[ωmin,ωmax]×[amin,amax]×[bmin,vmax]＝[1h,4h]×[-0.5,0.5]×[-0.5,0.5]，高度搜索空间为UH＝[550m,800m]，高度搜索间隔为0.5m。
向脉冲压缩后的回波中加入噪声使信噪比为-10dB，脉冲压缩后的结果如图4所示，可以看出回波完全淹没在噪声中。采用本发明提出的方法进行转动参数搜索，获得的估计值如表2所示。重建后的小行星三维形状如图5所示，390个散射点的高度值估计偏差如图6所示，其中估计偏差小于5米的散射点有316个，占总散射点数的81％，估计偏差小于10米的散射点有341个，占总散射点数的87％。
表2：转动参数估计结果
a
b
ω
真实值
0
0
2.8
估计值
0.0083315
0.005391
2.7952
偏差值
0.0083315
0.005391
0.0048
本发明可以准确地估计天体的转动参数与三维形状，通过仿真验证了本发明的有效性和实用性。
