摘要
面向在轨目标的天基逆合成孔径雷达(ISAR)成像是空间态势感知(SSA)的一项关键支撑技术。传统的二维距离-多普勒(RD)图像虽可部分揭示目标的散射特性,但其本质缺乏预测轨道机动和实现非合作目标识别所必需的三维几何信息。而现有基于图像序列的多视角重构方法在天基应用场景下存在固有局限:天基平台与目标卫星间的相对轨道运动会造成有效观测时间受限与成像投影平面不稳定的问题。为应对上述挑战,提出一种可变观测模式的目标重构算法。首先,该算法直接利用各特显点距离徙动(RM)轨迹中蕴含的目标结构信息,规避了易引入误差的图像配准过程。其次,推导了二维与三维旋转模式下距离徙动轨迹的表征模型,并提取了用于旋转模式分类的判别特征。最后,提出一种基于高阶多普勒系数估计的高精度距离徙动估计技术;对于不同的旋转模式,结合截断核范数正则化的因式分解方法可实现受观测误差影响时的二维或三维目标重构。仿真结果表明,所提空间目标重构算法可有效实现散射点提取以及距离徙动矩阵重建,进一步实现距离徙动矩阵的正则化收敛,从而获得不同旋转状态下空间目标重构结果,验证了所提算法的有效性与灵活性。
Abstract
The advancement of spaceborne inverse synthetic aperture radar (ISAR) imaging for on-orbit target represents a critical technology for space situational awareness (SSA). While conventional two-dimensional (2D) range-Doppler (RD) imaging provides valuable scattering intensity distributions, it inherently is the projection of the target’s three-dimensional (3D) structure, thereby losing critical geometric information essential for predicting orbital maneuvers and enabling non-cooperative target recognition. Current multi-views reconstruction methods based on image sequences face inherent limitations in spaceborne scenarios: the relative orbital motion between the spaceborne platform and target satellite induces limited observation time and unstable imaging projection plane. To address these challenges, a target reconstruction algorithm with variable observation mode is proposed. First, the structural information contained in range migration (RM) trajectories is directly exploited to avoid the error-prone image alignment process. Second, we derive a unified geometric model characterizing the range migration (RM) evolution under both 2D and 3D rotation patterns, and extract discriminative features for rotation pattern classification. Finally, we develop a high-precision RM estimation algorithm based on higher-order Doppler coefficient estimation. For different rotation modes, a truncated nuclear norm (TNN) regularization combined with factorization framework enables the reconstruction of 2D or 3D targets under observation error conditions. Simulation results demonstrate that the proposed target reconstruction algorithm effectively achieves scattering point extraction and RM matrix reconstruction. It further ensures the regularized convergence of the RM matrix, thereby obtaining spatial target reconstruction results under various rotation states. This validates the effectiveness and flexibility of the proposed algorithm.
卫星目标在轨状态分析是空间态势感知(space situation awareness,SSA)领域的一项关键任务[1]。卫星的三维几何特征可用来描述卫星的在轨状态,并为后续监测卫星状态、分析卫星意图和识别卫星类型提供重要信息。高分辨率逆合成孔径雷达(inverse synthetic aperture radar,ISAR)成像技术的进步[2],为实现这些特征的实时精准提取提供了技术支撑。当前,此类特征往往通过解译地基ISAR获取的高分辨率图像获得。然而,传统的地基雷达在对空间目标进行连续跟踪时,往往受限于大气干扰及轨道遮蔽[3-4],存在其固有局限。在当前背景下,天基ISAR凭借其全天候、长航时的观测能力,逐步成为在轨目标特征提取的核心技术手段[5-6]。
尽管传统二维距离-多普勒(range-Doppler,RD)成像结果可以提供有价值的散射强度分布信息,但其本质仍是目标三维结构在成像平面上的投影[7],这将导致三维几何信息严重缺失。为获取此类信息,目标三维重构技术致力于恢复目标三维结构参数。目前可实现三维重构的雷达体制主要分为多天线系统[8]与单天线系统[9]两类。多天线系统理论上可实现瞬时三维成像,然而其在天基平台应用中面临载荷资源的约束与高精度星间信号同步的挑战,在实际部署中面临显著障碍。
单天线系统通常需要依赖目标的自旋以及目标与平台间的相对轨道运动,从而获取多视角下的成像序列来实现三维重构[10]。然而,天基ISAR具有其独特的复杂性:1)目标卫星与天基平台间的相对轨道运动会导致成像投影平面(imaging projection plane,IPP)平稳的假设失效[11];2)受限于有限的成像时间,累积的观测角度不足,难以实现对目标的高精度三维重构[5]; 3)复杂的目标旋转模式会显著改变各散射中心的距离徙动(range migration,RM)轨迹,增加目标重构难度[12]。目前基于因式分解的方法[13]依赖特定的目标运动模型,面对上述复杂场景时往往难以有效适用。
为突破现有方法的局限,近期研究正积极探索新的目标重构技术。Zhou等[14]构建了光学-雷达联合处理系统,实现了对自旋稳定卫星姿态动态的高精度估计。Cao等[9]提出的新型单天线三维成像算法能够提取目标的径向高度差信息,可显著简化系统结构。为避免复杂的成像系统设计以及严格的图像配准要求,本文研究将聚焦于利用距离徙动轨迹实现目标三维重构的方法[15]。然而,当前的利用距离徙动轨迹的方法仍存在局限:一方面对复杂多变的目标旋转模式缺乏普适性;另一方面,距离徙动轨迹估计方法的抗噪声的鲁棒性不足。
针对上述问题,本文首先建立了成像几何模型,该模型能够精确描述目标在二维与三维旋转模式下距离徙动轨迹的变化规律,并揭示了距离徙动矩阵在不同旋转模式下呈现出的显著秩退化现象,并以此作为区分目标不同旋转模式的判别性特征。其次,给出了一种利用高阶系数估计的高精度距离徙动估计算法。最后,引入了一种结合新型截断核范数(truncated nuclear norm,TNN)正则化的矩阵分解方法,用以处理受轨道摄动与系统噪声影响的距离徙动矩阵,提出的TNN分解方法能够自适应地保留信号子空间的奇异向量,从而在不同的距离徙动估计误差下实现目标重构。
1 ISAR成像模型
在天基ISAR的成像模型中,天基平台与观测卫星之间的相对运动可分解为平动运动与转动运动,平动运动主要由天基平台与观测卫星各自的轨道运动引起,转动运动主要由目标的三维转动以及雷达视线变化引起。如图1所示展示了存在三维旋转卫星的ISAR几何成像模型。图1中成像坐标系O-NUW原点位于卫星中心。围绕N、U、W轴的三维旋转角(滚转角、俯仰角和偏航角)分别标记为θN(t)=wNt、θU(t)=wUt、θW(t)=wWt,其中wN、wU、wW为三维旋转角的角速度。标记为红色的雷达视线(radar line of sight,RLOS)矢量可由方位角ηl(t)=θa+wat和俯仰角φl(t)=θe+wet表示,其中θe、θa、we、wa分别为方位角和俯仰角的初始角度和角速度。
图1成像几何模型
Fig.1Imaging geometry model
天基平台和观测卫星的轨道信息可根据两行轨道根数以及窄带跟踪系统,提前求解,这为后续平动运动补偿提供了便利。平动补偿后,目标的三维转动以及雷达视线变化引起的距离徙动可表示为
(1)
式中:sP=[nP,uP,wP]是散射点P在O-NUW坐标系中的位置矢量,RRot为三维旋转矩阵,ilos为RLOS矢量[16]。
将RRot,ilos的具体表达式代入式(1)后,Rp(t)可线性表示为
(2)
为便于后续推导,式(2)可用向量RP(m),m=1,2,···,M表示,其中M为离散时间数,全部散射点的距离徙动可表示为下列矩阵形式[15]:
(3)
其中:
式中P为散射点个数。
上述广义成像模型考虑了目标的三维旋转以及随时间变化的雷达视线。其中等效旋转向量weff(t)会随成像过程发生变化,随着目标三维旋转运动趋于稳定且成像时间缩短时,weff(t)保持恒定,使得向量u(t)=weff(t)×ilos(t)的方向保持稳定。此时成像投影平面趋于平稳,可采用二维旋转模式来描述当前运动。
在二维旋转模式下,散射点P的距离徙动可表示为
(4)
式中:θRot(t)=wefft为等效旋转角,可通过等效旋转角速度weff=weff2计算获得;分别为散射点P在图像投影平面上的投影坐标。
与式(3)类似,二维旋转模式下的距离徙动矩阵可表示为
(5)
其中:
式(3)和式(5)展示了不同旋转模式下距离徙动矩阵的形式。对于三维旋转模式,该矩阵R由CM×3和S3×P获得,此时矩阵R的秩等为3。对于二维旋转模式,矩阵R的秩为2。由此可见,旋转模式与矩阵R的秩存在关联。故可通过对矩阵R进行奇异值分解(singular value decomposition,SVD)判断旋转模式。矩阵的SVD可以表示为
(6)
式中:U,V为两个正交矩阵,Δ=diag(λ1,···,λr)为对角矩阵,r为R的秩,λi为第i个奇异值,有λ1≥λ2≥···≥λr。
当各散射点距离徙动的估计值不存在误差时,其二维或三维旋转模式下距离徙动矩阵的秩分别等于2或3。但距离徙动的估计误差不可避免地会干扰秩的估计。但仍然可以通过λi,i≥2之间的关系来确定旋转模式,判决旋转模式的方法为[12]
(7)
下面分析不同旋转模式下比值η的变化规律。对于二维旋转模式,λ3和λ4均为噪声奇异值,λ2为显著奇异值之一,η应接近零值。对于三维旋转模式,λ4为噪声奇异值,λ2和λ3均为显著奇异值,η应远大于零值。需要注意的是,无法通过设置η的阈值来严格判断二维或三维旋转,只能说η越大,三维旋转的特征越明显[12,17]。图2给出了不同成像时间和旋转角速度下,根据散射点距离徙动矩阵获得比值η。可以看出随着成像时间和角速度的增加,该比值逐渐增大,旋转模式趋向于三维旋转。
图2比值变化情况
Fig.2Ratio versus imaging time and angular velocity
2 目标重构算法
上述分析表明,距离徙动矩阵与散射点坐标之间存在特定关系(式(3)和式(5))。距离徙动矩阵的秩可用来进一步判别相应的旋转模式(式(7))。本文将介绍目标重构算法,分别包含距离徙动矩阵估计、距离徙动矩阵正则化,以及目标重构3个步骤。
2.1 距离徙动矩阵重建
首先对接收到的雷达回波数据进行脉冲压缩,并进一步利用天基平台和观测卫星的轨道信息实现平动运动补偿。此时,散射点剩余的转动距离徙动会导致目标回波能量分散至多个距离单元,影响后续参数估计,各散射点的距离徙动各不相同,因此利用楔形变换(keystone transform,KT)校正距离徙动。KT后,进一步实施方位维压缩可获得粗聚焦RD图像。对于单个散射点,其在RD图像上的回波强度分布通常遵循特定模式。因此,采用自适应模板匹配方法提取单个散射点,例如多尺度高斯拉普拉斯(Lplacian of Gaussian,LOG)检测器[18]。将LOG算子与RD图像进行二维卷积运算,并进一步搜索局部极值的坐标可获取孤立散射点的时延、多普勒频率及尺寸参数。
基于LOG检测器的检测结果,可进一步通过图像域滤波器提取每个散射点对应的方位维信号。散射点P的方位维信号可表示为
(8)
式中:σP为信号幅度,tm为KT后的慢时间维时间,λ=c/fc为雷达波长,fc为载频,c为光速,n(tm)为噪声。
RP(tm)可利用泰勒级数展开,此时式(8)可表示为[19]
(9)
式中K0,K1,K2,K3为RP(tm)的各阶泰勒级数展开系数。
接下来,采用迭代自适应法(iterative adaptive approach,IAA)[20]、积分三次相位函数法(integrated cubic phase function,ICPF)[21]和积分广义三次相位函数法(integrated generalized cubic phase function,IGCPF)[22]来估计中心频率、线性调频率及二次线性调频率,这些方法具有高估计精度和交叉项抑制能力。首先,通过IGCPF估算三阶多普勒系数,构建补偿项将式(9)解调为线性调频信号。随后采用ICPF估算二阶系数,将式(9)解调为单频信号。最后通过IAA估算一阶系数。K0的估计结果可由LOG检测器输出获得。
利用估计的系数,即可实现当前散射点的距离徙动向量构建。
(10)
重复上述操作,即可获得全部提取散射点的距离徙动向量,从而获得距离徙动矩阵。
2.2 距离徙动矩阵正则化
为增强后续空间目标重构时,距离徙动矩阵对噪声的鲁棒性,本文提出基于TNN约束[23]的正则化框架。该方法在自适应抑制噪声主导的奇异值的同时,保留了距离徙动矩阵的内在低秩结构,从而提升目标重构精度。从中恢复具有低秩特性的矩阵,这一过程的数学表达式为
(11)
式中:E为干扰项,为正则化恢复的距离徙动矩阵,为l1范数,为截断核范数,即中min(M,P)-r个最小奇异值之和,r为与旋转模式相关的秩,其对应的判别条件见式(7)。
由于截断核范数是非凸的,直接求解式(11)并不容易。因此借助SVD,将式(11)转化为以下的凸优化问题:
(12)
式中:为核范数,后前r个特征值对应的向量矩阵,λ为正则化因子,决定了重建矩阵的稀疏性。
为提高计算效率,采用交替方向乘子法(alternating direction method of multipliers,ADMM)求解上述凸优化问题[24],通过交替优化,将复杂的凸优化问题转化为一系列简单的子优化问题。式(12)的拉格朗日函数为
(13)
式中:Y为拉格朗日乘子,ρ>0为惩罚参数,为Frobenius范数的平方。
上述凸优化问题被分解为3个子优化问题,并通过固定两个变量,交替迭代的方式求解剩余变量:
(14)
式中k为k次迭代。式(14)的解为
(15)
(16)
(17)
式中:矩阵的SVD分解为,sgn(Q)为矩阵Q的元素符号矩阵。
2.3 可变模式目标重构算法
对估计得到的距离徙动矩阵进行矩阵正则化后得到,接下来将介绍针对不同旋转模式的目标重构算法。以三维旋转模式为例,经SVD后,矩阵可进一步表示为[27]
(18)
式中:
式(18)在形式上与式中的C,S矩阵完全一致。然而,该矩阵的分解并非唯一,仍然存在满足下列条件的三维可逆方阵H,即
(19)
此时,问题已转化为寻找一个合适的矩阵H,使得估计矩阵满足特定约束:的行向量是单位向量。根据约束条件:
(20)
式中,为的第m行的行向量。
是一个对称矩阵,W可表示为
(21)
将式(21)代入式(20),并进行展开可得
(22)
其中:
根据展开式,可得
(23)
利用线性回归,可以得到式(22)的最小二乘法解为
(24)
然后,H可通过对W的SVD结果得到。最后,三维位置坐标矩阵的重构结果为
(25)
对于二维旋转模式,目标重构流程与三维旋转模式类似。需特别注意的是,在将矩阵R奇异值分解后,本文仅取前两个奇异值,Δ=diag(λ1,λ2),且后续约束仍要求行向量为单位向量。基于此约束条件,结合式(18)~(25)所述的矩阵分解过程,可轻松获得二维旋转模式下的目标重构结果。与三维旋转模式下获得的不同,为目标在IPP上的二维投影。完整的目标重构流程见图3。
图3目标重构算法流程图
Fig.3Target reconstruction flowchart
为了便于计算复杂度分析,下面给出各关键步骤的计算成本。M,P,I,K分别为成像积累脉冲数,提取散射点个数,IAA迭代次数,ADMM迭代次数。则IGCPF,ICPF,IAA的计算复杂度分别为O(M3),O(M2),O(IM3),SVD的计算复杂度为O(min(M2P,MP2)),矩阵正则化的计算复杂度为O(K·min(M2P,MP2))。考虑到空间平台计算能力的限制,难以在空间平台上实现图像的精细化处理。因此,空间平台仅进行原始数据提取以及粗聚焦成像。并基于初步处理结果,卫星可通过轨道变更操作重新访问目标区域,并将获取的原始数据传输至地面站进行后续处理,如图像精细处理、三维重构等。
3 仿真结果与分析
天基平台相较于地基平台,可对目标卫星实现更近距离的观测。同时Ka频段具有更短的波长和更大的可用带宽[2],这意味着更高的成像分辨率。而较小的天线尺寸可实现高频段信号传输,这对于对载荷尺寸和质量有严格要求的空间平台而言非常适合。此外,空间环境中的大气衰减可以忽略不计,从而充分发挥Ka频段成像模式的优势。
3.1 目标随轨道变化的重构结果
在本文中,具体的仿真试验将被用来验证提出算法的性能。图4(a)展示了包含8个散射点的三维散射点模型,同时,目标存在围绕偏航方向的三维旋转,旋转速度为0.05 rad/s[16,28]。图4(b)展示了天基平台与目标卫星的相对运行轨道,图4(c)~(e)给出了相对运动参数变化,所提算法将分别在时间段1,时间段2,时间段3,时间段4内验证,各时间段内目标具有相同的三维转动,但随着相对轨道运动,成像视角会发生相应改变。表1给出了雷达系统的发射信号参数。
图4仿真背景
Fig.4Simulation background
表1信号参数
Tab.1Signal parameters
图5展示了各时间段内目标重构算法中预处理步骤获得的结果。包括基于LOG检测器的散射点提取结果,以及脉冲压缩结果,其中红色虚线为估计获得的距离徙动重建结果,这验证了提出算法中距离徙动矩阵重建算法的准确性。
图5距离徙动矩阵重建结果
Fig.5Range matrix reconstruction result
考虑到旋转模式会随成像时间和目标三维旋转运动而变化,后续验证中目标保持恒定的三维旋转运动,通过改变成像时间满足不同的旋转模式。在图3所示仿真背景下,分别选取0.5 s以及2.0 s的成像时间以满足不同的旋转模式。图6展示了时间段1内,不同旋转模式下,距离徙动矩阵正则化方法获得的收敛曲线,验证了距离徙动矩阵正则化方法的可收敛性。图7展示了不同时间段内获得的目标重构结果,其中理想结果由散射点实际的距离徙动矩阵获得。重构结果验证了运动模型的正确性,以及所提算法在不同仿真背景下的有效性。
图6收敛曲线
Fig.6Convergence curve
图7目标重构结果
Fig.7Results of target reconstruction
3.2 目标不同旋转状态下的重构结果
本文采用图4(a)所示的空间目标散射点模型,并将仿真背景中的时间段2,用于后续不同旋转状态目标重构实验验证。对不同旋转状态目标进行目标重构所需的仿真参数不同,见表2,其余信号参数见表1。图8展示了不同旋转状态目标的三维重构结果。在各种情况的仿真中,雷达视线不发生变化,随着目标的三维旋转,所需成像时间也发生相应改变。至此,本实验展示了所提算法对不同旋转状态目标三维重构的有效性。
表2旋转参数
Tab.2Rotational parameters
图8目标重构结果
Fig.8Results of target reconstruction
3.3 目标重构算法鲁棒性分析
为验证所提算法的鲁棒性,采用蒙特卡洛实验进行验证。式(26)~(28)给出了目标二维、三维重构与距离徙动矩阵估计的均方根误差(root-mean-square error,RMSE)计算公式。
(26)
(27)
(28)
式中:SRME为均方根误差,Mc为蒙特卡洛实验次数,为第i次蒙特卡洛实验中的目标重建结果,(nP,uP,wP),(XP,YP)为利用散射点实际距离徙动矩阵R获得的目标参考结果,为带误差的距离徙动矩阵。
本文采用图4(a)所示的空间目标散射点模型,并将仿真背景中的时间段2用于后续蒙特卡洛实验验证。图9展示了成像时间为2.0 s时,本文提出的距离徙动矩阵估计算法与历史相位估计法[12]在不同信噪比下的估计精度,可以看出在较低信噪比下,历史相位估计法的估计精度受噪声影响很大,而本文提出算法仍可保持较高的估计精度,从而实现距离徙动矩阵的高精度估计。
图9距离徙动矩阵重建算法鲁棒性验证情况
Fig.9Robustness verification of RM matrix reconstruction
图10(a)、(b)分别展示了成像时间为0.5 s和2.0 s时,目标重构算法在不同旋转模式下、面对不同SRMERM时的目标重构结果的RMSE。图10(c)展示了在SRMERM为0.03 m的条件下,目标重构算法在不同成像时间下的RMSE。从图10可见,随着SRMERM增大,目标重构精度逐渐下降。相较于传统分解方法,提出的TNN分解法在更高SRMERM条件下仍保持更优的估计精度。与二维重构结果相比,三维重构结果对距离徙动矩阵估计误差更为敏感。随着成像时间的延长,距离徙动矩阵的信息量增加,目标重构的估计精度也相应提高。
图10目标重构算法鲁棒性验证情况
Fig.10Robustness verification of target reconstruction algorithm
4 结论
1)针对不同目标旋转模式,建立天基ISAR几何成像模型,同时揭示了距离徙动矩阵在不同旋转模式下呈现出的显著秩退化现象,并以此作为区分目标不同旋转模式的判别性特征。
2)提出一种基于高阶多普勒系数估计的高精度距离徙动估计技术。
3)有效融合截断核范数正则化与因式分解框架,使算法可以实现存在观测误差下的高精度目标重构。最后,通过仿真试验验证了所提算法在随轨道变化以及不同旋转状态情况下的有效性,相比于传统的距离徙动矩阵估计算法以及目标重构算法,所提算法在低信噪比下具有更强的鲁棒性。

