Loading [MathJax]/jax/element/mml/optable/GeneralPunctuation.js

基于机载多通道雷达迭代超分辨估计的前视成像

任凌云 吴迪 朱岱寅 孙伟杰

任凌云, 吴迪, 朱岱寅, 等. 基于机载多通道雷达迭代超分辨估计的前视成像[J]. 雷达学报, 2023, 12(6): 1166–1178. doi: 10.12000/JR23085
引用本文: 任凌云, 吴迪, 朱岱寅, 等. 基于机载多通道雷达迭代超分辨估计的前视成像[J]. 雷达学报, 2023, 12(6): 1166–1178. doi: 10.12000/JR23085
REN Lingyun, WU Di, ZHU Daiyin, et al. Forward-looking imaging via iterative super-resolution estimation in airborne multi-channel radar[J]. Journal of Radars, 2023, 12(6): 1166–1178. doi: 10.12000/JR23085
Citation: REN Lingyun, WU Di, ZHU Daiyin, et al. Forward-looking imaging via iterative super-resolution estimation in airborne multi-channel radar[J]. Journal of Radars, 2023, 12(6): 1166–1178. doi: 10.12000/JR23085

基于机载多通道雷达迭代超分辨估计的前视成像

DOI: 10.12000/JR23085
基金项目: 国家自然科学基金(62271252)
详细信息
    作者简介:

    任凌云,博士生,主要研究方向为雷达前视成像、MIMO雷达信号处理

    吴 迪,教授,主要研究方向为雷达前视成像、雷达信号处理、空时自适应处理技术等

    朱岱寅,教授,主要研究方向为合成孔径雷达/逆合成孔径雷达(SAR/ISAR)成像以及自聚焦算法、MIMO雷达信号处理、干涉SAR成像、SAR地面动目标指示以及机载雷达空中动目标指示技术等

    孙伟杰,硕士生,主要研究方向为雷达前视成像、阵列信号处理

    通讯作者:

    吴迪 wudi82@nuaa.edu.cn

    朱岱寅 zhudy@nuaa.edu.cn

  • 责任主编:李亚超 Corresponding Editor: LI Yachao
  • 中图分类号: TN957.5; TN959.4

Forward-looking Imaging via Iterative Super-resolution Estimation in Airborne Multi-channel Radar

Funds: The National Natural Science Foundation of China (62271252)
More Information
  • 摘要: 波达角估计算法用于机载多通道雷达前视成像时可以突破瑞利极限,实现同一波束主瓣宽度内的多目标分辨,改善成像的方位向分辨率,然而天线波束覆盖有限且其快速扫描使得可用于协方差矩阵估计的数据样本缺乏,导致对目标位置和幅度估计出现误差。该文提出了一种基于单快拍迭代超分辨处理的多通道雷达前视成像算法,通过对单个空域快拍的迭代谱估计可获得目标的准确位置和幅度信息,再通过多个脉冲的非相干累积得到前视方位高分辨成像。仿真和实测数据处理结果表明,所提算法具有分辨多目标的能力,相较于传统前视成像算法显著提高了前视图像的方位分辨率,同时保证了点目标的精确重构和面目标的轮廓重构。

     

  • 机载雷达前视成像因其全天时、全天候等突出优势,在敌情探测、飞机盲降、导弹制导和复杂地形探测等应用中发挥着不可替代的作用[1]。一般情况下,对地对海成像分辨率主要体现在距离向和方位向两个维度,距离向通过发射大时宽带宽积信号,并利用脉冲压缩技术,得到距离向高分辨率。而提高方位分辨率并不容易,其受到天线孔径的限制,并且随着作用距离增大急剧恶化[2]。现有的高分辨雷达成像技术,如合成孔径雷达(Synthetic Aperture Radar, SAR)[3]和多普勒波束锐化(Doppler Beam Sharpening, DBS)[46],均在雷达前视区域形成成像盲区。实孔径雷达(Real Aperture Radar, RAR)利用雷达实波束对前方观测区域进行扫描,积累各个方位的回波进而成像[7,8],但是其方位向分辨率低,已经不能满足人们对于雷达前视分辨率日益提升的需求[9]

    近些年来,各种方位高分辨率及超分辨率算法不断涌现,这些算法可以在不同配置下提供飞行器前方区域的较清晰成像结果,完成对前方目标的有效侦察,已经在各个方面得到推广。目前,主要的前视成像技术包括单脉冲成像技术、双基地SAR前视成像技术、解卷积成像技术、阵列雷达超分辨成像技术等。南京航空航天大学吴迪等人[10,11]对单脉冲成像算法和自聚焦算法进行了研究,并通过实测数据处理验证了所提算法对成像质量的改善,其测角精度高,但性能受限于空域自由度,无法对同一波束内的多目标进行分辨[12]。对于双基地SAR前视成像技术,国内外关于其理论和成像算法的研究在不断完善和成熟中[13,14],其系统的高复杂性,包括特殊几何形状的约束、时频同步、跨平台通信等诸多因素,限制了其在实际应用中的性能[15,16]。解卷积成像技术将方位向回波视为某一距离门上目标沿方位向的分布函数与天线方向图两者的卷积,通过解卷积运算反演目标分布函数[17,18],国内以电子科技大学为代表对其进行了广泛而深入的研究,并基于实测数据验证了所提前视成像算法的角度分辨能力[1922]

    与传统的单天线扫描雷达相比,阵列雷达超分辨成像技术采用阵列接收系统,由于其具备更高的空间分辨能力,可以在波束宽度内分辨多个目标,在对复杂场景成像时能获得更好的分辨率。方位向基于阵列信号处理实现超分辨成像的物理本质在于天线孔径长度在方位向的延伸,所以将阵列数据直接通过线性预测的方法进行扩展可以有效提升方位分辨率,自回归(Autoregressive, AR)模型是外推方法中常用的预测模型[23]。2015年,Cho等人[24]证实了将基于线性预测方法的阵列数据扩展用于阵列雷达前视成像中,可以提高前视探地雷达方位向分辨率。2018年,南京航空航天大学张洁[25]利用Burg递推算法处理机载阵列雷达回波数据,验证了对阵元快拍数据合理外推,相当于增加阵列天线长度,可以有效提升前视方位分辨率。另一种机载阵列雷达前视方位超分辨成像的思路是通过对各阵元接收到的回波信号,进行空间谱估计,例如多重信号分类(Multiple Signal Classification, MUSIC)算法等已应用到机载雷达前视成像中,并通过仿真和实测数据处理验证了其可行性[26]。但是为了对回波协方差矩阵进行准确估计,需要积累大量的快拍,由于载机的运动和天线波束的快速扫描,这在机载雷达实际应用中是难以获得的。

    为了解决经典谱估计算法需要大量快拍,对方位角分辨率提升有限以及受噪声影响较大等问题,本文提出了单快拍迭代超分辨前视成像算法。首先,该算法基于最小均方误差准则建立代价函数并求解[27];然后通过迭代的方式对单个快拍的2阶统计量进行精确估计[28];最后,在获得前视成像区域内所有目标的角度和幅度后,通过非相干累积得到二维高分辨率前视成像结果。仿真和实测数据处理结果表明:相较于传统前视成像算法,单快拍迭代超分辨前视成像算法有效提高了前视成像的方位分辨率,成像效果更好。

    脉冲多普勒雷达固定在运动平台上,采用沿切航迹方向的阵列天线对前视区域进行扫描,扫描角速度为 $ \omega $ ,以均匀水平线性阵列(Uniform Linear Array, ULA)为例,建立机载多通道雷达前视成像模型,如图1所示。将阵列接收天线划分为若干个子阵列,每个子阵列即构成了一个独立接收通道,用M代表通道总数,d代表通道间距,示意图如图2所示。为在保证远探测距离的同时,获得较高的距离分辨率,雷达发射线性调频信号(Linear Frequency Modulation, LFM)。现将前视成像区域进行距离向和方位向的划分,分为L个距离单元,K个方位角度单元。第l个距离门,第k个方位角度单元上的目标的斜距和方位角度分别用 $ {R_l} $ $ {\theta _k} $ 表示,散射系数 $ \sigma ({R_l},{\theta _k}) $ 表示,在本文的讨论中,只考虑前视成像方位向的超分辨,俯仰角度近似一定为 $ \varphi $ 。在扫描整个前视区域 $ \varOmega $ 后,第m个通道回波表示为

    图  1  机载多通道雷达前视成像三维几何结构图
    Figure  1.  3D geometry for airborne multi-channel radar forward-looking imaging
    图  2  收发天线示意图
    Figure  2.  Schematic diagram of transmitting and receiving antennas
    sm(τ,t)=Ll=1Kk=1σ(Rl,θk)h(tθkω)×rect(ττl(t,θk)Tp)×exp{j2\pi fc(ττl(t,θk))}×exp{j\pi Kr[ττl(t,θk)]2} (1)

    其中, $ \tau $ 为距离向快时间域变量,t为方位向慢时间域变量,与雷达运动和天线扫描有关, $ h(t) $ 代表方位向上的双程天线方向图,其随着慢时间变量t变化,代表了方位向的调制作用。 $ {\tau _l}(t,{\theta _k}) $ 为发射信号到达第l个距离单元 $ {\theta _k} $ 方位上目标后经目标反射到第m个通道的传播时延, $ \mathrm{rect}(\cdot) $ 代表矩形窗函数, $ {f_{\text{c}}} $ 为雷达发射信号的载频, $ {K_{\text{r}}} $ 为距离向的线性调频率, $ {T_{\text{p}}} $ 为信号的脉冲宽度。

    传播时延 $ {\tau _l}(t,{\theta _k}) $ 表示为

    τl(t,θk)=2Rl(t,θk)c+(m1)dsinθkcosφc (2)

    其中, $ {\text{c}} $ 代表光速。瞬时斜距 $ {R_l}(t,{\theta _k}) $ 随机载雷达飞行而变化,在任意t时刻近似表示为

    Rl(t,θk)=Rlvtcosθkcosφ (3)

    其中,v为平台的运动速度。由式(3)可见,雷达与目标间存在距离徙动,其在带来多普勒频移的同时,还造成了距离方位向的耦合,为了解耦合需要对回波进行距离徙动校正。由于一般平台的距离徙动会远小于场景中心距离运动平台的初始距离,跨距离门数较少,距离空变的影响不大,且前视区域的目标集中在航迹 $ \pm 10^\circ $ 左右区域,可将方位余弦 $ \cos {\theta _k} $ 近似为1, $ {R_l}(t,{\theta _k}) $ 记为 $ {R_l}(t) $ 。那么在距离频域乘以相位因子以完成距离徙动校正

    HRCM=exp(j2\pi fr2vtcosφc) (4)

    经过下变频处理和距离向匹配滤波和徙动校正后的回波表示为

    sm(τ,t)=Ll=1Kk=1σ(Rl,θk)h(tθkω)×sinc(B(τ2Rlc))×exp(j2π(m1)dλsinθkcosφ)×exp{j4πλRl(t)} (5)

    对某一特定距离单元,将距离 $ {R_l} $ 的双程延时相位补偿之后,得到该距离单元的方位向回波

    sm(t)=Kk=1σ(Rl,θk)h(tθkω)×exp(j4πλvtcosφ)×exp(j2π(m1)dλsinθkcosφ) (6)

    因为前视区域多普勒频率梯度小,近似认为各方位目标的多普勒历史相等,并且本文仅限于空域处理,而时域采用非相干积累,多普勒频率项影响不大,便忽略掉其对回波的影响。式(6)变为

    sm(t)=Kk=1σ(Rl,θk)h(tθkω)×exp(j2π(m1)dλsinθkcosφ) (7)

    可见点目标的方位向回波受到方位向双程天线方向图的调制,导致方位分辨率受到天线波束宽度的限制,波束宽度内的多目标难以分辨。

    xk(t)=σ(Rl,θk)h(tθkω) (8)
    a(θk)=exp(j2π(m1)dλsinθkcosφ) (9)

    那么

    sm(t)=Kk=1xk(t)a(θk)+nm(t) (10)

    其中, $ {x_k}(t) $ 为第k个目标天线方向图加权后的分布函数, $ n_m(t) $ 为第m个阵元上的噪声。

    t时刻,波束中心指向 $ \bar \theta $ ,对M个通道的回波进行采样,得到的空域快拍表示为 ${\boldsymbol{s}} = {[{s_1},{s_2}, \cdots ,{s_M}]^{\text{T}}}$ ,其矩阵形式为

    s=Ax+N (11)

    其中, $ {\boldsymbol{N}} \in {\mathbb{C}_{M \times 1}} $ 为观测噪声向量, ${\boldsymbol{N}} = [{n_1}, {n_2}, \cdots ,{n_M}]^{\text{T}}$ ,服从零均值高斯分布; $ {\boldsymbol{x}} \in {\mathbb{C}_{Q \times 1}} $ $ {\boldsymbol{A}} \in {\mathbb{C}_{M \times Q}} $ 为空间频率不混叠范围内的目标分布函数和空域导引矢量矩阵, $ {\boldsymbol{x}} = {[{x_1},{x_2}, \cdots ,{x_Q}]^{\text{T}}} $ ,空间频率不混叠范围为

     ˉθθb<θ<ˉθ+θb , θb = arcsin(λ2dcosφ) (12)

    在此连续角度范围上进行离散化 $({\theta _1},{\theta _2}, \cdots , {\theta _Q})$ ,那么, $ {\boldsymbol{A}} \in {\mathbb{C}_{M \times Q}} $ 表示为

    A=[as(θ1) as(θ2) as(θQ)]=[111ej2πdλsinθ1cosφej2πdλsinθ2cosφej2πdλsinθQcosφej2πdλ(M1)sinθ1cosφej2πdλ(M1)sinθ2cosφej2πdλ(M1)sinθQcosφ] (13)

    在接收到空域快拍之后,需要对天线方向图加权后的目标分布函数进行精确重构。将估计值表示为空域快拍的线性加权和,令 ${\boldsymbol{w}} = [{w_1},{w_2}, \cdots ,{w_M}]$ ,代表待确定的权系数构成的权矢量,那么,待估计量表示为

    ˆx=wHs=Mi=1wisi (14)

    采用最小均方误差准则(Minimum Mean Square Error, MMSE)建立代价函数

    J{ (15)

    将式(15)对w求偏导,解得极小值点为

    {\boldsymbol{w}} = {\left( {{{\rm{E}}} \left\{ {{\boldsymbol{s}}{{\boldsymbol{s}}^{\text{H}}}} \right\}} \right)^{ - 1}}{{\rm{E}}} \left\{ {{\boldsymbol{s}}{{\boldsymbol{x}}^{\text{H}}}} \right\} (16)

    其中, ${\rm{E}}\left\{ {{\boldsymbol{s}}{{\boldsymbol{s}}^{\text{H}}}} \right\}$ 为多通道接收信号的自相关矩阵R,表示为 ${\boldsymbol{R}} = {\boldsymbol{AP}}{{\boldsymbol{A}}^{\text{H}}} + {{\boldsymbol{R}}_{\text{N}}}$ P是待估计信号的空间功率对角矩阵,其对角元素 ${p_k} = {\left| {{x_k}} \right|^2}, k = 1, 2,\cdots ,K$ ,噪声协方差矩阵 ${{\boldsymbol{R}}_{\text{N}}} = \sigma _{\text{N}}^2{\boldsymbol{I}}$ 。由此式(16)化简得

    \begin{split} {\boldsymbol{w}}{\text{ = }}&{\left( {{\boldsymbol{AP}}{{\boldsymbol{A}}^{\text{H}}} + {{\boldsymbol{R}}_{\rm{N}}}} \right)^{ - 1}}{\rm{E}}\left\{ {({\boldsymbol{Ax}} + {\boldsymbol{N}}){{\boldsymbol{x}}^{\text{H}}}} \right\} \\ {\text{ = }}&{\left( {{\boldsymbol{AP}}{{\boldsymbol{A}}^{\text{H}}} + {{\boldsymbol{R}}_{\rm{N}}}} \right)^{ - 1}}{\boldsymbol{AP}} \end{split} (17)

    实际处理时,信号的空间功率分布矩阵P的对角线元素的指数取1到2之间的某个数,以保持矩阵求逆的收敛性和稳定性。由式(17)可以看出,对x的估计依赖于其本身计算得到的空间功率对角矩阵,定义初始空间功率分布为

    {\hat {\boldsymbol{P}}_1} = \left[{\hat {\boldsymbol{x}}_{{\rm{MF}}}}\hat {\boldsymbol{x}}_{{\rm{MF}}}^{\text{H}}\right] \odot {{\boldsymbol{I}}_{K \times K}} (18)

    其中, $ {\hat {\boldsymbol{x}}_{{\rm{MF}}}} $ 为匹配滤波的结果, $ {\hat {\boldsymbol{x}}_{{\rm{MF}}}} = {{\boldsymbol{A}}^{\text{H}}}{\boldsymbol{s}} $ $ \odot $ 表示哈达玛积(Hadamard Product)。

    将式(17)求得的权矢量代入式(14)中,得到关于目标的角度和幅度估计值,但是这样求得的估计是不够精确的。在前视成像的实际应用中,需要对x进行精确估计,此时可以引入迭代处理的思想,以不断缩小估计值和真实值之间的误差。接下来建立迭代过程,用变量下标代表迭代次数,基于前一次迭代的估计值 $ {\hat {\boldsymbol{x}}_{t - 1}} $ 计算空间功率分布矩阵 $ {{\boldsymbol{P}}_t} $

    {\hat {\boldsymbol{P}}_t} = \left[{\hat {\boldsymbol{x}}_{t - 1}}\hat {\boldsymbol{x}}_{t - 1}^{\text{H}}\right] \odot {{\boldsymbol{I}}_{K \times K}} (19)

    然后更新权值矩阵 $ {{\boldsymbol{w}}_{t + 1}} $

    {{\boldsymbol{w}}_{t + 1}} = {\left( {{\boldsymbol{A}}{{\boldsymbol{P}}_t}{{\boldsymbol{A}}^{\text{H}}} + {{\boldsymbol{R}}_{\rm{N}}}} \right)^{ - 1}}{\boldsymbol{A}}{{\boldsymbol{P}}_t} (20)

    最后利用新的权值矩阵对空域快拍进行自适应加权得到新的估计值,即

    {\hat {\boldsymbol{x}}_{t + 1}} = {\boldsymbol{w}}_{t + 1}^{\text{H}}{\boldsymbol{s}} (21)

    通常每个快拍迭代10次左右可以得到较好重构效果,回波的自相关矩阵R不依赖于大量的快拍数据,而是通过迭代的方式得到,由此解决了传统超分辨谱估计算法中依赖大量快拍对自相关矩阵估计的问题。

    通过算法1可以看出,随着迭代次数的增加,每一次迭代过程中的噪声都会对当前估计结果造成扰动,这种影响被逐步放大,可能导致某次迭代出现病态问题。在阵列信号处理领域,当空域扫描范围较大时,这种噪声产生的扰动能够起到有效正则化的作用,在处理过程中可以抑制噪声的放大,实现有效的超分辨估计[29]。但是在前视成像的实际应用中,通过调整发射天线的发射波束模式,天线波束只对前视区域进行扫描,大概在–10°~10°范围。而整个方位角度内,除扫描范围外的其他位置,即使这些位置目标的反射功率很小,这些回波仍然对自相关矩阵的秩有贡献。随着方位角度划分的网格数的减小,自相关矩阵的条件数会增加到不利的水平,此时噪声扰动的正则化无法改善算法的病态问题。并且实际处理过程中,噪声功率不可知,无法在迭代过程中理想的加载噪声协方差矩阵 $ {{\boldsymbol{R}}_{\rm{N}}} $ 而计算权值矩阵,随着迭代超分辨趋近于局部收敛,回波的自相关矩阵趋近奇异矩阵,其病态性不可避免。针对这个问题,本节采用加载对角矩阵进行正则化,即

    1  迭代超分辨算法流程
    1.  Flow chart of iterative super-resolution algorithm
     根据当前波束中心 $ \bar \theta $,由式(12)计算得空域不混叠范围;在此范
     围上构建空域导引矢量矩阵A
     初始化:
      ${\hat {\boldsymbol{x}}_1} = {{\boldsymbol{A}}^{\text{H} } }{\boldsymbol{s}}$,基于此计算信号的空间功率分布矩阵 ${{\boldsymbol{P}}_1}$
     重复以下操作:
      ${{\boldsymbol{w}}_{t + 1} } = {\left( {{\boldsymbol{A}}{{\boldsymbol{P}}_t}{{\boldsymbol{A}}^{\text{H} } } + {{\boldsymbol{R}}_{\rm{N}}} } \right)^{ - 1} }{\boldsymbol{A}}{{\boldsymbol{P}}_t}$
      ${\hat {\boldsymbol{x} }_{t + 1} } = {\boldsymbol{w} }_{t + 1} ^{\text{H} }{\boldsymbol{s} }$
      ${\hat {\boldsymbol{P}}_{t + 1} } = [{\hat {\boldsymbol{x}}_t}\hat {\boldsymbol{x}}_t^{\text{H} }] \odot {{\boldsymbol{I}}_{K \times K} }$
     直到收敛
    下载: 导出CSV 
    | 显示表格
    {\boldsymbol{R}} = {\boldsymbol{AP}}{{\boldsymbol{A}}^{\text{H}}} + \rho {\boldsymbol{I}} (22)

    其中, $ \rho $ 为加载的对角矩阵值,I为单位矩阵。

    通过矩阵条件数的分析可以证明对角加载对矩阵稳定性的影响。矩阵的条件数定义为

    \kappa {\text{(}}{\boldsymbol{R}}{\text{) = }}\frac{{{\lambda _{\max }}}}{{{\lambda _{\min }}}} (23)

    其中, $ {\lambda _{\max }} $ $ {\lambda _{\min }} $ 分别为R矩阵的最大和最小的特征值。

    经过对角加载后的回波自相关矩阵的条件数变为

    \kappa {\text{(}}{\boldsymbol{R}}{\text{) = }}\frac{{{\lambda _{\max }} + \rho }}{{{\lambda _{\min }} + \rho }} = 1 + \frac{{{\lambda _{\max }} - {\lambda _{\min }}}}{{{\lambda _{\min }} + \rho }} (24)

    可见,加载值 $ \rho $ 取值越大,条件数越小,通过合适的对角加载可以消除矩阵的病态性。另外,要考虑加载对角矩阵对估计结果的影响,加载值过高会降低估计的分辨率。所以,综合考虑,正则化参数即对角加载值只能选择略大于使矩阵病态的临界值。在实际应用中,可以选择先不进行加载,直到迭代过程出现病态警告,这时可以用小的对角矩阵替换[30]。在本节算法中,因为算法引入了迭代处理,每一次迭代均需要加载正则化参数,取一个固定的加载值便显得不合理,此时需要对自适应的对角加载进行研究[31,32]。令M维单位矩阵 $ {{\boldsymbol{I}}_{M \times M}} $ 的第m列用 $ {\boldsymbol{{i}}_m} $ 表示,那么加载的正则化参数为

    \rho = \frac{1}{M}{\sum\limits_{m = 1}^M {\left| {\frac{{{\boldsymbol{i}}_m^{\text{H}}{{\boldsymbol{R}}^{ - 1}}{\boldsymbol{s}}}}{{{\boldsymbol{i}}_m^{\text{H}}{{\boldsymbol{R}}^{ - 1}}{{\boldsymbol{i}}_m}}}} \right|} ^\alpha },\alpha = 2 (25)

    通过自适应的正则化,本节所提的单快拍迭代超分辨前视成像算法在一定程度上避免了病态问题[31,32],但是在实际处理中,根据多次试验的经验,需要在该结果的基础上再进行微调作为最终的对角加载值,通过适当减小式(25)中的指数 $ \alpha $ 在1.5~1.8可以在有效解决病态问题的同时保证估计的分辨率。

    图3所示,在获得了单个距离-脉冲单元的方位谱估计结果之后,下一步是根据当前的波束位置,将其累加到最终的图像矩阵中。由于所有脉冲都是独立处理的,非相干积累在实际处理中具有更强的鲁棒性。处理完所有的距离-脉冲单元,循环累积结束便获得了最终的前视区域二维成像图。

    图  3  迭代超分辨多通道雷达前视成像处理流程图
    Figure  3.  Flow chart of iterative super-resolution forward-looking imaging processing of multi-channel radar

    为了验证所提的单快拍迭代超分辨前视成像算法的成像性能,对前视区域仿真点目标进行成像。

    实波束扫描成像结果如图4(a)所示,由于实波束成像方法的分辨率是由波束宽度决定的,可以看出经过实波束成像之后,同一距离门内方位角度差小于波束宽度的两个相近目标无法区分。但是经过迭代超分辨前视成像技术处理后,5个点目标均清晰可见,成像结果如图4(b)所示。图4(c)为同一距离门(1250 m)处的方位向成像剖面图,由此证明了所提的前视成像技术可以有效区分波束主瓣中的临近多点目标,实现了方位向的超分辨率。

    图  4  点目标前视成像结果
    Figure  4.  Forward-looking imaging results for point targets

    点目标的扩散程度是衡量成像系统质量的一个指标,如图5所示,在表1所列仿真参数下在雷达正前方1000 m处设置一个点目标,采用实波束成像之后点扩散函数的3 dB宽度为1.6°,而迭代超分辨成像之后,点扩散函数的3 dB宽度仅为0.15°,可见点目标锐化能力相较于实波束成像算法提升了10倍。

    图  5  点扩散函数仿真
    Figure  5.  Simulation of Point Spread Function (PSF)
    表  1  点目标仿真实验系统参数
    Table  1.  Point target simulation parameters
    参数名称 参数值
    雷达波长(m) 0.03
    系统带宽(MHz) 20
    采样频率(MHz) 30
    脉冲重复频率(Hz) 1000
    通道数 8
    通道间隔(m) 0.09
    波束方位向主瓣宽度(°) 2.11
    机载平台运动速度(m/s) 100
    波束扫描速度(°/s) 100
    前视扫描范围(°) –10~10
    目标点距离(m) 1000, 1250, 1250, 1250, 1500
    目标点方位角(°) 0, –3, –2, 3, 0
    目标点SNR (dB)
    20, 25, 25, 20, 20
    下载: 导出CSV 
    | 显示表格

    本节继续通过仿真两组实测场景对第3节所提算法进行验证,为了更逼真地模拟地面的回波,采用已有的机载SAR系统的高分辨率SAR复图像作为待成像的地面场景,将其SAR图像中的每个像素视为一个地面散射点,其对应的复值为反射率,注意设置的相邻目标散射点之间的距离需要小于距离和方位分辨率。两组场景仿真均基于机载多通道雷达前视成像模型,信噪比设置为25 dB,系统参数如表2所示。

    表  2  场景仿真实验系统参数
    Table  2.  Parameters of simulation for simulated scenario
    参数名称 参数值
    (场景1)
    参数值
    (场景2)
    雷达波长(m) 0.03 0.03
    系统带宽(MHz) 40 40
    采样频率(MHz) 50 50
    脉冲重复频率(Hz) 1000 1000
    通道数 8 8
    通道间隔(m) 0.045 0.120
    波束方位向主瓣宽度(°) 4.23 1.59
    波束扫描范围(°) –20~20 –8~8
    机载平台运动速度(m/s) 100 80
    下载: 导出CSV 
    | 显示表格

    图6为第1组仿真场景成像结果对比,图6(b)为实波束成像结果,方位向分辨率低,导致无法对点目标进行分辨,并且面目标轮廓特征模糊。图6(c)采用DBS成像[4,5],可见DBS在载机正前方区域无法完成成像,因为大斜视区域的多普勒频率梯度差异不大,可见在大斜视区域DBS成像的分辨率也不高,图6(d)在多普勒域进行加权自适应处理即迭代DBS前视成像算法[4],该方法相较于传统DBS成像在前视区域分辨率提升不明显。图6(e)为单快拍迭代超分辨前视成像结果,可见方位向分辨率相较于实波束成像和DBS成像有了很大的提升,点目标和面目标都变得清晰,并且面目标的轮廓线条更加分明。

    图  6  第1组场景仿真成像结果对比
    Figure  6.  Forward-looking imaging results for the first scenario

    第2组实验场景主要包括机场的一些跑道、建筑等,其SAR图像如图7(a)所示。分别采用实波束成像、单脉冲成像[10],基于空间平滑的MUSIC超分辨成像[26]和所提的迭代超分辨成像算法对其前视回波仿真数据进行处理。图7(b)为实波束扫描成像结果,成像质量较差,目标难以分辨。图7(c)为单脉冲成像结果,因为当波束主瓣范围内中存在多个目标时,单脉冲雷达的角度估计能力就失效了,导致成像处理中出现分辨率差甚至假目标的问题。此外,在处理复杂场景时,方位角估计不准确,会导致最终图像中散斑点的增加。这些问题基本上都是由于和差波束单脉冲处理在空间域上缺乏自由度造成的。单脉冲成像在空域具备两个自由度,而文献[26]和第3节所提的算法都是基于机载多通道雷达回波进行处理的,具备多空域自由度,角度分辨能力更强。图7(d)为基于空间平滑的MUSIC超分辨成像结果,因为慢时间域相邻脉冲采样的空域多快拍之间的相关性强,采用这些样本估计的协方差矩阵秩亏损,故利用空间平滑技术改善这一问题,由成像结果可见,对前视成像分辨率的提升有限。图7(e)采用本文所提的迭代超分辨前视成像算法进行成像,证明了利用单个空域快拍迭代估计超分辨谱的方法有效,其对于点目标提升了方位分辨率,特别是对于面目标的成像效果相较于实波束成像、单脉冲成像和基于空间平滑的MUSIC超分辨成像算法有较大的改善。细节图如图8图9所示。

    图  7  第2组场景仿真成像结果对比
    Figure  7.  Forward-looking imaging results for the second scenario
    图  8  区域1成像结果放大图
    Figure  8.  Enlarged image of zone 1
    图  9  区域2成像结果放大图
    Figure  9.  Enlarged image of zone 2

    为了继续验证所提算法在实际应用中的可行性,对两组机载前视阵列雷达实测数据进行处理,主要系统参数如表3所示。

    表  3  实测数据系统参数
    Table  3.  Parameters of measured data
    参数名称 第1组实测数据 第2组实测数据
    雷达波段 X Ku
    脉冲重复频率(Hz) 417 3125
    通道数 4 4
    通道间隔(m) 0.150 0.037
    波束主瓣宽度(°) 2.42 6.31
    波束扫描范围(°) –14~14 –20~20
    扫描速度(°/s) 80 70
    机载平台速度(m/s) 145 68
    下载: 导出CSV 
    | 显示表格

    第1组数据由某型X波段机载阵列雷达录取,成像场景为山地区域。图10给出了通道幅相误差校正后经过各技术处理的成像结果。可以看出,图10(a)为实波束扫描成像结果,方位分辨率差。如图10(b)所示,传统的单脉冲方法可以在一定程度上获得较好的方位分辨率,特别是对于孤立的、强的点状目标,然而,由于单脉冲成像不能对波束内多目标成像,成像区域的面目标的轮廓信息会被遗漏。图10(c)给出了迭代超分辨前视成像算法的成像结果。相比之下比传统的实波束成像方法和单脉冲方法能获得更高的分辨率,更好地重构真实场景,土地轮廓明显更加清晰,部分放大图如图11所示。

    图  10  实测数据成像结果
    Figure  10.  Forward-looking imaging results for the measured data
    图  11  实测数据成像结果细节图
    Figure  11.  Enlarged imaging results for the measured data

    接下来从图像熵和对比度的角度来验证所提前视成像算法的有效性,图像熵反映了图像中的平均信息量,定义为

    H = \sum\limits_{i = 0}^{255} {{P_i}\log {P_i}} (26)

    其中, $ {P_i} $ 表示图像中灰度值为i的像素所占比例。图像对比度定义为图像中的灰度反差值,通常反映了图像画质的清晰程度,即

    C = \sum\limits_\delta {\delta {{(i,j)}^2}} {P_\delta }(i,j) (27)

    其中, $ \delta (i,j) = \left| {i - j} \right| $ ,即相邻像素间的灰度差, $ {P_\delta }(i,j) $ 表示相邻像素间的灰度差为 $ \delta $ 的分布概率。3种前视成像算法对实测数据的成像图像的熵和对比度如表4所示,可见所提的迭代超分辨成像算法的成像结果清晰度最高,质量最好,分辨率得到进一步提升。

    表  4  地面实测数据成像结果图像熵和对比度对比
    Table  4.  Entropy and contrast of the measured data results
    成像方法 图像熵 对比度
    实波束成像 3.8625 4.42
    单脉冲成像结果 1.6530 8.12
    迭代超分辨成像 1.2543 9.59
    下载: 导出CSV 
    | 显示表格

    第2组数据由Ku波段机载前视阵列雷达录取,扫描海面区域,对海平面强杂波干扰下船只等目标进行成像。图12(a)为实波束扫描成像结果,这些海面点目标可以从不同的距离单元中区分出来,但由于分辨率较低,在方位角方向上明显模糊。图12(b)为采用单脉冲成像技术进行处理的结果,相较于实波束成像结果,方位向分辨率有一定提升。图12(c)为本文所提算法处理结果,可见方位分辨率得到了明显提高,海面船只等点目标清晰可见,并且对于杂波和噪声的抑制效果明显。

    图  12  海面区域实测数据成像结果
    Figure  12.  Forward-looking imaging results for the measured data of sea surface

    本文针对机载多通道雷达前视成像方位向分辨率低和难以获取大量独立同分布数据样本等问题,在已建立多通道模型的基础上,提出了单快拍超分辨前视成像算法。该算法在MMSE准则下建立代价函数,通过递归完成样本的自相关矩阵的精确估计。仿真和实测数据验证了所提算法相较于实波束成像和单脉冲成像以及多普勒波束锐化等,方位分辨率更高,成像效果更好。

  • 图  1  机载多通道雷达前视成像三维几何结构图

    Figure  1.  3D geometry for airborne multi-channel radar forward-looking imaging

    图  2  收发天线示意图

    Figure  2.  Schematic diagram of transmitting and receiving antennas

    图  3  迭代超分辨多通道雷达前视成像处理流程图

    Figure  3.  Flow chart of iterative super-resolution forward-looking imaging processing of multi-channel radar

    图  4  点目标前视成像结果

    Figure  4.  Forward-looking imaging results for point targets

    图  5  点扩散函数仿真

    Figure  5.  Simulation of Point Spread Function (PSF)

    图  6  第1组场景仿真成像结果对比

    Figure  6.  Forward-looking imaging results for the first scenario

    图  7  第2组场景仿真成像结果对比

    Figure  7.  Forward-looking imaging results for the second scenario

    图  8  区域1成像结果放大图

    Figure  8.  Enlarged image of zone 1

    图  9  区域2成像结果放大图

    Figure  9.  Enlarged image of zone 2

    图  10  实测数据成像结果

    Figure  10.  Forward-looking imaging results for the measured data

    图  11  实测数据成像结果细节图

    Figure  11.  Enlarged imaging results for the measured data

    图  12  海面区域实测数据成像结果

    Figure  12.  Forward-looking imaging results for the measured data of sea surface

    1  迭代超分辨算法流程

    1.   Flow chart of iterative super-resolution algorithm

     根据当前波束中心 $ \bar \theta $,由式(12)计算得空域不混叠范围;在此范
     围上构建空域导引矢量矩阵A
     初始化:
      ${\hat {\boldsymbol{x}}_1} = {{\boldsymbol{A}}^{\text{H} } }{\boldsymbol{s}}$,基于此计算信号的空间功率分布矩阵 ${{\boldsymbol{P}}_1}$
     重复以下操作:
      ${{\boldsymbol{w}}_{t + 1} } = {\left( {{\boldsymbol{A}}{{\boldsymbol{P}}_t}{{\boldsymbol{A}}^{\text{H} } } + {{\boldsymbol{R}}_{\rm{N}}} } \right)^{ - 1} }{\boldsymbol{A}}{{\boldsymbol{P}}_t}$
      ${\hat {\boldsymbol{x} }_{t + 1} } = {\boldsymbol{w} }_{t + 1} ^{\text{H} }{\boldsymbol{s} }$
      ${\hat {\boldsymbol{P}}_{t + 1} } = [{\hat {\boldsymbol{x}}_t}\hat {\boldsymbol{x}}_t^{\text{H} }] \odot {{\boldsymbol{I}}_{K \times K} }$
     直到收敛
    下载: 导出CSV

    表  1  点目标仿真实验系统参数

    Table  1.   Point target simulation parameters

    参数名称 参数值
    雷达波长(m) 0.03
    系统带宽(MHz) 20
    采样频率(MHz) 30
    脉冲重复频率(Hz) 1000
    通道数 8
    通道间隔(m) 0.09
    波束方位向主瓣宽度(°) 2.11
    机载平台运动速度(m/s) 100
    波束扫描速度(°/s) 100
    前视扫描范围(°) –10~10
    目标点距离(m) 1000, 1250, 1250, 1250, 1500
    目标点方位角(°) 0, –3, –2, 3, 0
    目标点SNR (dB)
    20, 25, 25, 20, 20
    下载: 导出CSV

    表  2  场景仿真实验系统参数

    Table  2.   Parameters of simulation for simulated scenario

    参数名称 参数值
    (场景1)
    参数值
    (场景2)
    雷达波长(m) 0.03 0.03
    系统带宽(MHz) 40 40
    采样频率(MHz) 50 50
    脉冲重复频率(Hz) 1000 1000
    通道数 8 8
    通道间隔(m) 0.045 0.120
    波束方位向主瓣宽度(°) 4.23 1.59
    波束扫描范围(°) –20~20 –8~8
    机载平台运动速度(m/s) 100 80
    下载: 导出CSV

    表  3  实测数据系统参数

    Table  3.   Parameters of measured data

    参数名称 第1组实测数据 第2组实测数据
    雷达波段 X Ku
    脉冲重复频率(Hz) 417 3125
    通道数 4 4
    通道间隔(m) 0.150 0.037
    波束主瓣宽度(°) 2.42 6.31
    波束扫描范围(°) –14~14 –20~20
    扫描速度(°/s) 80 70
    机载平台速度(m/s) 145 68
    下载: 导出CSV

    表  4  地面实测数据成像结果图像熵和对比度对比

    Table  4.   Entropy and contrast of the measured data results

    成像方法 图像熵 对比度
    实波束成像 3.8625 4.42
    单脉冲成像结果 1.6530 8.12
    迭代超分辨成像 1.2543 9.59
    下载: 导出CSV
  • [1] 杨建宇. 雷达对地成像技术多向演化趋势与规律分析[J]. 雷达学报, 2019, 8(6): 669–692. doi: 10.12000/JR19099

    YANG Jianyu. Multi-directional evolution trend and law analysis of radar ground imaging technology[J]. Journal of Radars, 2019, 8(6): 669–692. doi: 10.12000/JR19099
    [2] 樊晨阳, 贺思三, 郭乾. 雷达前视成像技术的研究现状[J]. 电光与控制, 2021, 28(9): 59–64. doi: 10.3969/j.issn.1671-637X.2021.09.013

    FAN Chenyang, HE Sisan, and GUO Qian. Research status of radar forward-looking imaging technology[J]. Electronics Optics&Control, 2021, 28(9): 59–64. doi: 10.3969/j.issn.1671-637X.2021.09.013
    [3] MOREIRA A, PRATS-IRAOLA P, YOUNIS M, et al. A tutorial on synthetic aperture radar[J]. IEEE Geoscience and Remote Sensing Magazine, 2013, 1(1): 6–43. doi: 10.1109/MGRS.2013.2248301
    [4] 陈洪猛. 机载广域监视雷达高分辨成像方法研究[D]. [博士论文], 西安电子科技大学, 2016.

    CHEN Hongmeng. Study of high resolution imaging for airborne wide area surveillance radar[D]. [Ph. D. dissertation], Xidian University, 2016.
    [5] ONAT E and ÖZKAZANÇ Y. An analysis of Doppler beam sharpening technique used in fighter aircraft[C]. The 26th Signal Processing and Communications Applications Conference, Izmir, Turkey, 2018: 1–4.
    [6] MAO Deqing, ZHANG Yongchao, ZHANG Yin, et al. Doppler beam sharpening using estimated Doppler centroid based on edge detection and fitting[J]. IEEE Access, 2019, 7: 123604–123615. doi: 10.1109/ACCESS.2019.2937992
    [7] 杨志伟, 贺顺, 廖桂生. 机载单通道雷达实波束扫描的前视探测[J]. 航空学报, 2012, 33(12): 2240–2245.

    YANG Zhiwei, HE Shun, and LIAO Guisheng. Forward-looking detection for airborne single-channel radar with beam scanning[J]. Acta Aeronautica et Astronautica Sinica, 2012, 33(12): 2240–2245.
    [8] 温晓杨, 匡纲要, 胡杰民, 等. 基于实波束扫描的相控阵雷达前视成像[J]. 航空学报, 2014, 35(7): 1977–1991. doi: 10.7527/S1000-6893.2013.0545

    WEN Xiaoyang, KUANG Gangyao, HU Jiemin, et al. Forward-looking imaging based on real beam scanning phased array radars[J]. Acta Aeronautica et Astronautica Sinica, 2014, 35(7): 1977–1991. doi: 10.7527/S1000-6893.2013.0545
    [9] SUTOR T, WITTE F, and MOREIRA A. New sector imaging radar for enhanced vision: SIREV[C]. SPIE 3691, Enhanced and Synthetic Vision, Orlando, USA, 1999.
    [10] 吴迪, 朱岱寅, 田斌, 等. 单脉冲成像算法性能分析[J]. 航空学报, 2012, 33(10): 1905–1914.

    WU Di, ZHU Daiyin, TIAN Bin, et al. Performance evaluation for monopulse imaging algorithm[J]. Acta Aeronautica et Astronautica Sinica, 2012, 33(10): 1905–1914.
    [11] 吴迪, 杨成杰, 朱岱寅, 等. 一种用于单脉冲成像的自聚焦算法[J]. 电子学报, 2016, 44(8): 1962–1968. doi: 10.3969/j.issn.0372-2112.2016.08.027

    WU Di, YANG Chengjie, ZHU Daiyin, et al. An autofocusing algorithm for monopulse imaging[J]. Acta Electronica Sinica, 2016, 44(8): 1962–1968. doi: 10.3969/j.issn.0372-2112.2016.08.027
    [12] 李悦丽, 马萌恩, 赵崇辉, 等. 基于单脉冲雷达和差通道多普勒估计的前视成像[J]. 雷达学报, 2021, 10(1): 131–142. doi: 10.12000/JR20111

    LI Yueli, MA Meng’en, ZHAO Chonghui, et al. Forward-looking imaging via Doppler estimates of sum-difference measurements in scanning monopulse radar[J]. Journal of Radars, 2021, 10(1): 131–142. doi: 10.12000/JR20111
    [13] WALTERSCHEID I, ESPETER T, GIERULL C, et al. Results and analysis of hybrid bistatic SAR experiments with spaceborne, airborne and stationary sensors[C]. 2009 IEEE International Geoscience and Remote Sensing Symposium, Cape Town, South Africa, 2009: II-238–II-241.
    [14] YANG Jianyu, HUANG Yulin, YANG Haiguang, et al. A first experiment of airborne bistatic forward-looking SAR - Preliminary results[C]. 2013 IEEE International Geoscience and Remote Sensing Symposium, Melbourne, Australia, 2013: 4202–4204.
    [15] LIU Zhutian, YE Hongda, LI Zhongyu, et al. Optimally matched space-time filtering technique for BFSAR nonstationary clutter suppression[J]. IEEE Transactions on Geoscience and Remote Sensing, 2022, 60: 5210617. doi: 10.1109/TGRS.2021.3090462
    [16] 武俊杰, 孙稚超, 杨建宇, 等. 基于GF-3照射的星机双基SAR成像及试验验证[J]. 雷达科学与技术, 2021, 19(3): 241–247. doi: 10.3969/j.issn.1672-2337.2021.03.002

    WU Junjie, SUN Zhichao, YANG Jianyu, et al. Spaceborne-airborne bistatic SAR using GF-3 illumination—technology and experiment[J]. Radar Science and Technology, 2021, 19(3): 241–247. doi: 10.3969/j.issn.1672-2337.2021.03.002
    [17] ZHANG Yongchao, MAO Deqing, ZHANG Qian, et al. Airborne forward-looking radar super-resolution imaging using iterative adaptive approach[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2019, 12(7): 2044–2054. doi: 10.1109/JSTARS.2019.2920859
    [18] ZHANG Qiping, ZHANG Yin, HUANG Yulin, et al. TV-sparse super-resolution method for radar forward-looking imaging[J]. IEEE Transactions on Geoscience and Remote Sensing, 2020, 58(9): 6534–6549. doi: 10.1109/TGRS.2020.2977719
    [19] TUO Xingyu, ZHANG Yin, HUANG Yulin, et al. Fast sparse-TSVD super-resolution method of real aperture radar forward-looking imaging[J]. IEEE Transactions on Geoscience and Remote Sensing, 2021, 59(8): 6609–6620. doi: 10.1109/TGRS.2020.3027053
    [20] HUO Weibo, TUO Xingyu, ZHANG Yin, et al. Balanced tikhonov and total variation deconvolution approach for radar forward-looking super-resolution imaging[J]. IEEE Geoscience and Remote Sensing Letters, 2022, 19: 3505805. doi: 10.1109/lgrs.2021.3072389
    [21] 毛德庆. 机载雷达扫描波束超分辨成像方法研究[D]. [博士论文], 电子科技大学, 2022.

    MAO Deqing. Research on scanning beam super-resolution imaging methods for airborne radar[D]. [Ph. D. dissertation], University of Electronic Science and Technology of China, 2022.
    [22] ZHANG Yin, SHEN Jiahao, TUO Xingyu, et al. Scanning radar forward-looking superresolution imaging based on the Weibull distribution for a sea-surface target[J]. IEEE Transactions on Geoscience and Remote Sensing, 2022, 60: 5116111. doi: 10.1109/TGRS.2022.3194118
    [23] GUPTA I J, BEALS M J, and MOGHADDAR A. Data extrapolation for high resolution radar imaging[J]. IEEE Transactions on Antennas and Propagation, 1994, 42(11): 1540–1545. doi: 10.1109/8.362783
    [24] CHO B L and SUN S G. Cross-range resolution improvement in forward-looking imaging radar using autoregressive model-based data extrapolation[J]. IET Radar,Sonar&Navigation, 2015, 9(8): 933–941. doi: 10.1049/IET-RSN.2014.0495
    [25] 张洁. 阵列雷达前视成像技术研究[D]. [硕士论文], 南京航空航天大学, 2018.

    ZHANG Jie. Study on array radar forward-looking imaging techniques[D]. [Master dissertation], Nanjing University of Aeronautics and Astronautics, 2018.
    [26] ZHANG Jie, WU Di, ZHU Daiyin, et al. An airborne/missile-borne array radar forward-looking imaging algorithm based on super-resolution method[C]. The 10th International Congress on Image and Signal Processing, BioMedical Engineering and Informatics, Shanghai, China, 2017: 1–5.
    [27] BLUNT S D and GERLACH K. Adaptive pulse compression via MMSE estimation[J]. IEEE Transactions on Aerospace and Electronic Systems, 2006, 42(2): 572–584. doi: 10.1109/TAES.2006.1642573
    [28] BLUNT S D, CHAN T, and GERLACH K. Robust DOA estimation: The reiterative superresolution (RISR) algorithm[J]. IEEE Transactions on Aerospace and Electronic Systems, 2011, 47(1): 332–346. doi: 10.1109/TAES.2011.5705679
    [29] YARDIBI T, LI Jian, STOICA P, et al. Source localization and sensing: A nonparametric iterative adaptive approach based on weighted least squares[J]. IEEE Transactions on Aerospace and Electronic Systems, 2010, 46(1): 425–443. doi: 10.1109/TAES.2010.5417172
    [30] ROBERTS W, STOICA P, LI Jian, et al. Iterative adaptive approaches to MIMO radar imaging[J]. IEEE Journal of Selected Topics in Signal Processing, 2010, 4(1): 5–20. doi: 10.1109/JSTSP.2009.2038964
    [31] 张小飞, 张胜男, 徐大专. 自适应对角线加载的波束形成算法[J]. 中国空间科学技术, 2007, 27(2): 66–71. doi: 10.3321/j.issn:1000-758X.2007.02.011

    ZHANG Xiaofei, ZHANG Shengnan, and XU Dazhuan. Adaptive diagonal loading beamforming algorithm[J]. Chinese Space Science and Technology, 2007, 27(2): 66–71. doi: 10.3321/j.issn:1000-758X.2007.02.011
    [32] CARLSON B D. Covariance matrix estimation errors and diagonal loading in adaptive arrays[J]. IEEE Transactions on Aerospace and Electronic Systems, 1988, 24(4): 397–401. doi: 10.1109/7.7181
  • 期刊类型引用(1)

    1. 毛德庆,杨建宇,杨明杰,张永超,张寅,黄钰林. IAA-Net:一种实孔径扫描雷达迭代自适应角超分辨成像方法. 雷达学报. 2024(05): 1073-1091 . 本站查看

    其他类型引用(0)

  • 加载中
图(12) / 表(5)
计量
  • 文章访问数: 764
  • HTML全文浏览量: 250
  • PDF下载量: 265
  • 被引次数: 1
出版历程
  • 收稿日期:  2023-05-10
  • 修回日期:  2023-07-10
  • 网络出版日期:  2023-07-27
  • 刊出日期:  2023-12-28

目录

/

返回文章
返回