Loading [MathJax]/jax/output/SVG/jax.js

微波视觉三维SAR实验系统及其全极化数据处理方法

仇晓兰 罗一通 宋舒洁 彭凌霄 程遥 颜千程 上官松涛 焦泽坤 张柘 丁赤飚

张嘉翔, 张凯翔, 梁振楠, 等. 一种基于深度强化学习的频率捷变雷达智能频点决策方法[J]. 雷达学报(中英文), 2024, 13(1): 227–239. doi: 10.12000/JR23197
引用本文: 仇晓兰, 罗一通, 宋舒洁, 等. 微波视觉三维SAR实验系统及其全极化数据处理方法[J]. 雷达学报(中英文), 2024, 13(5): 941–954. doi: 10.12000/JR24137
ZHANG Jiaxiang, ZHANG Kaixiang, LIANG Zhennan, et al. An intelligent frequency decision method for a frequency agile radar based on deep reinforcement learning[J]. Journal of Radars, 2024, 13(1): 227–239. doi: 10.12000/JR23197
Citation: QIU Xiaolan, LUO Yitong, SONG Shujie, et al. Microwave vision three-dimensional SAR experimental system and full-polarimetric data processing method[J]. Journal of Radars, 2024, 13(5): 941–954. doi: 10.12000/JR24137

微波视觉三维SAR实验系统及其全极化数据处理方法

DOI: 10.12000/JR24137 CSTR: 32380.14.JR24137
基金项目: 国家自然科学基金(61991420, 61991421, 61991424)
详细信息
    作者简介:

    仇晓兰,研究员,博士生导师,主要研究方向为SAR成像处理、SAR图像解译、新体制SAR

    罗一通,硕士,助理研究员,主要研究方向为SAR图像几何处理、雷达摄影测量、SAR三维成像

    宋舒洁,博士生,主要研究方向为全极化无人机载阵列干涉三维成像

    彭凌霄,硕士,助理研究员,主要研究方向为SAR图像信息处理和应用、SAR图像仿真技术

    程 遥,硕士,助理研究员,主要研究方向为SAR成像处理

    颜千程,博士生,主要研究方向为阵列干涉SAR三维成像处理

    上官松涛,博士,助理研究员,主要研究方向为遥感影像定量化标校与应用

    焦泽坤,博士,副研究员,主要研究方向为SAR微波视觉三维成像、多维度微波成像等

    张 柘,研究员,博士生导师,主要研究方向为新体制微波成像、稀疏信号处理及信号处理与人工智能的结合

    丁赤飚,研究员,博士生导师,主要研究方向为合成孔径雷达、遥感信息处理和应用系统等

    通讯作者:

    仇晓兰 xlqiu@mail.ie.ac.cn

  • 责任主编:王岩 Corresponding Editor: WANG Yan
  • 中图分类号: TN957.52

Microwave Vision Three-dimensional SAR Experimental System and Full-polarimetric Data Processing Method

Funds: The National Natural Science Foundation of China (61991420, 61991421, 61991424)
More Information
  • 摘要: 三维合成孔径雷达在测绘制图、防灾减灾等诸多领域有应用潜力,已经成为SAR的重要研究方向。为减少三维SAR的观测次数或天线阵元数量,推动三维SAR的应用和发展,中国科学院空天信息创新研究院牵头研制了微波视觉三维SAR实验系统,旨在为微波视觉SAR三维成像提供实验平台和数据。该文针对微波视觉三维SAR实验系统及其全极化数据处理方法进行介绍,涵盖了极化校正、极化相干增强、极化约束三维成像、三维融合可视化等全流程的关键步骤。基于发布的SAR微波视觉三维成像全极化数据集,给出了三维成像结果示例,验证了微波视觉三维SAR实验系统的全极化性能以及处理方法的有效性。该文发布的数据集将为SAR三维成像研究提供良好的数据条件。

     

  • 在现代战争中,敌方为了获取电磁频谱优势与战场主动权,通常会发射各种有源干扰破坏雷达作战性能,从而掩护目标完成预定的作战任务[1]。雷达为了应对各种干扰,相应的抗干扰技术在对抗中不断升级[2]。一般来说,抗干扰技术按照雷达处理阶段的不同可以分为主动抗干扰和被动抗干扰[3]。在雷达发射信号阶段,主动抗干扰技术可以通过雷达波形设计降低敌方干扰机对雷达信号的截获概率或识别概率,从而降低干扰机的干扰效能[46]。如果雷达已经接收到了干扰信号,被动抗干扰技术可以通过空、时、频等多个处理域完成目标与干扰的分离,达到对干扰抑制的目的[79]

    随着雷达抗干扰研究的不断深入,被动抗干扰手段日益丰富。然而,挂载在掩护目标上的自卫式干扰机通过发射大功率瞄准干扰,使干扰与目标回波在多处理域重叠,难以分离。频率捷变雷达通过使用自主调节发射信号载频的主动抗干扰手段,使得干扰机难以截获和干扰,为对抗自卫式压制干扰提供了可能[10]。其抗干扰性能主要取决于跳频策略,传统随机跳频策略已经被证明不是最佳选择[11]。如何精准预测干扰机下一时刻将要发射的干扰频点,从而指导雷达信号的频点选择,是频率捷变雷达在与干扰机博弈中取胜的主要难点。

    相比针对静态优化问题设计的启发式搜索算法,强化学习可以让智能体与环境不断交互,获得反馈,从而指导智能体在动态环境下进行决策[12]。基于深度学习模型强大的数据表征能力而衍生出的深度强化学习,能够处理高维数据并完成非线性映射,弥补了传统强化学习算法的不足[13],在认知电子战方面已经得到了一定的研究。如果将干扰信息看作环境状态,抗干扰措施看作雷达动作,抗干扰效能看作即时回报,那么认知抗干扰决策问题可以通过强化学习技术解决。文献[14]针对干扰类型和参数固定的复合干扰场景,分别使用Q学习和SARSA (State-Action-Reward-State-Action)探索了抗干扰措施组合选取问题。文献[15]使用改进的DDPG (Deep Deterministic Policy Gradient)算法对12种抗干扰措施进行选择,以实施抗干扰措施前后干扰威胁度变化作为反馈。文献[16]使用DDPG-MADDPG (Deep Deterministic Policy Gradient and the Multi-Agent Deep Deterministic Policy Gradient)对包含复合干扰在内的12种干扰类型,以抗干扰改善因子作为反馈,进行多处理域抗干扰措施自适应选取。

    在频点决策方面,强化学习主要围绕瞄频或扫频干扰的频率捷变波形设计展开研究[17]。文献[18]首次对雷达脉冲级跳频策略展开研究,分别对比了随机频点选择、Q学习、深度Q网络(Deep Q-Network, DQN)等3种策略,证明了DQN在决策方面具备更好的性能。并在文献[19]中继续深化研究内容,将检测概率作为奖励值,而不是之前论文中的信干噪比,同时优化了DQN模型。文献[20]在文献[18]和文献[19]工作的基础上,考虑了一种具备侦收功能的干扰机,以及子脉冲频率捷变雷达,并基于近端策略优化(Proximal Policy Optimization, PPO)算法完成智能决策。文献[21]考虑了网络化无人机雷达工作系统,使用雷达信息表示理论作为奖励函数,基于双贪婪的改进Q学习算法优化系统抗干扰性能。文献[22]假定干扰机也具备马尔科夫性质,在预测得到干扰策略的基础上选择雷达频点与之对抗。文献[23]考虑了跳频速率会影响相干积分性能和多普勒分辨率,使用Q学习自适应调整雷达发射波形的脉宽和频点以对抗扫频干扰。

    总体来说,上述研究均基于雷达不同的性能指标设计奖励函数,以此优化频点等雷达参数。虽然在对抗成功率方面超过随机频点决策方法,然而缺少对抗干扰策略收敛速度的讨论。应当指出,在现代电子战中,干扰机可能具备多种策略,并根据某种规则在不同策略间切换。因此雷达在进行抗干扰策略学习时,应当尽快收敛到最优策略,从而保持对抗先机。如果雷达还未收敛到最优策略时,干扰机改变策略,那么雷达将陷入被动地位。因此,网络收敛时间或是所需样本量是评价一个智能化算法能够应用于实际作战场景的重要衡量指标。

    受上述研究启发,考虑到现代干扰机具备侦收-瞄准-干扰的基本策略,本文针对频率捷变雷达,设计了一种基于强化学习的雷达子脉冲跳频抗干扰策略。将当前时刻感知到的干扰频点以及上一时刻的雷达频点作为状态,将当前时刻的雷达频点选择策略作为动作,以目标检测结果和信干噪比作为即时奖励函数设计强化学习关键要素,基于DQN完成子脉冲频点选取策略的学习。仿真针对两种不同侦收策略的干扰机,证明了所提方法的有效性以及较高的收敛效率。

    与文献[20]不同的是,本文的主要贡献在于如何通过对强化学习关键要素的设计,从而达到快速收敛到最优解的目的,而不是在于网络设计与修改。具体包括4点:(1)虽然干扰机具备侦干周期,但是我们通过状态空间的合理设计,仅使用单个时间步即可学习到干扰周期性策略,同时不需要使用长短期记忆网络(Long Short-Term Memory, LSTM)等时间记忆网络即可完成最优策略学习,显著降低了收敛时间。(2)在动作设计方面,我们设计了一种子脉冲频点可重复选取的特殊波形,增大了动作空间选取范围。(3)在动作选取方面,我们通过ε-贪婪原则,实现了搜索和利用的有效平衡。在训练初期,以随机搜索为主,减小了收敛到局部最优解的概率。随着训练过程的进行,随机搜索概率逐渐降低,选择网络输出动作的概率逐渐增加,便于收敛。(4)在奖励设计方面,围绕目标检测性能,在单次目标检测结果的基础上,引入了更具差异性的信干噪比指标,缓解了因为采样不充分可能收敛到局部最优解的情况。

    由于现代干扰机可以对接收到的雷达信号进行快速测频与频率引导,对传统雷达具备较大威胁。而频率捷变雷达可以实现子脉冲级的频率调制,为与其对抗提供了可能。作为常用的雷达传输信号波形,基于线性调频(Linear Frequency Modulation, LFM)信号的子脉冲频率捷变波形如图1(a)所示,其时域表达式如下:

    图  1  频率捷变波形示意图
    Figure  1.  Schematic diagram of the frequency agility waveform
    st(t)=Nn=1rect[(tτn)/(tτn)TsubTsub]exp[j2πfn(tτn)]exp[jπKn(tτn)2] (1)

    其中,rect()表示矩形窗函数,N表示子脉冲个数,Tsub表示子脉冲脉宽;τn表示第n个子脉冲的延时,fn表示子脉冲频点,Kn表示第n个子脉冲的调频斜率。频率捷变雷达各可选频点应当去相关从而达到频率抗干扰的目的,即保证si(ω)sj(ω)=0,其中,si(ω)表示子脉冲 i 的频谱,sj(ω)表示子脉冲 j的频谱。

    式(1)所定义的传统频率捷变雷达在进行子脉冲频点选取时,通常会选择不同的雷达频点。为扩充频点选取自由度,增大波形复杂度,本文设计了一种子脉冲频点可重复选取的雷达发射波形,如图1(b)所示。当相邻子脉冲选取重复频点时,则将其合成一个宽脉冲,其脉宽为Tcom=NrepTsub,其中Nrep表示选取相同频点的相邻子脉冲数量。同时保证合成后的宽脉冲带宽不变,即Bcom=Bsub。合成后的脉冲数用Ncom表示。

    强化学习可以由马尔科夫决策过程(Markov Decision Process, MDP)描述,满足马尔科夫性质。强化学习的优化目标为最大化累计回报,定义为

    Gt=rt+γrt+1+γ2rt+2+=k=0γkrt+k (2)

    其中,rt表示智能体在状态st下执行动作at并转移到st+1后得到的回报;γ为折扣因子,是st+1及其之后的奖励权重,取值范围为0~1,表示对未来奖励的重视程度。

    由于MDP是一种随机过程,其随机独立性导致累计回报Gt是一个随机变量,无法定量描述,如图2所示。因此可对累计回报取期望,获得状态值函数Vπ(s)和动作状态值函数Qπ(s,a),将优化问题变成找到一种最优策略π,使任意一个状态的Vπ(s)Qπ(s,a)为最大。而Q学习的优化目标是针对Qπ(s,a),其贝尔曼方程及最优动作状态值函数Q(s,a)定义如下:

    图  2  MDP的随机独立性与强化学习的优化目标
    Figure  2.  The random independence of MDP and the optimization objectives of reinforcement learning
    Qπ(s,a)=sSp(s|s,a)[r(s,a,s)+γaAπ(a|s)Qπ(s,a)] (3)
    Q(s,a)=sSp(s|s,a)[r(s,a,s)+γmaxaQ(s,a)] (4)

    其中,rt=r(s,a)=ap(s|s,a)r(s,a,s)p(s|s,a)为某状态s执行动作a后,转移到下一状态s的概率。

    由于在实际场景中,我们可能不知道环境先验信息p(s|s,a),因此无法获得值函数的解析表示。而Q学习可以通过多次取平均的方式,近似估计得到Q。具体来说,从任意状态开始与环境1个时间步长,利用t时刻的即时回报rt和下一时刻最大的状态动作值函数Q(st+1,at+1)对当前时刻动作状态值函数Q(st,at)进行估计,最后重复上述动作多次取平均。值函数的更新公式为

    Q(st,at)=Q(st,at)+α[rt+γmaxaQ(st+1,at+1)Q(st,at)] (5)

    其中,α为学习率,表示更新的步长。

    Q学习通过不断与环境进行交互来获取并更新Q值,并将Q值存入到由状态和动作组成的Q表中。待智能体学习完成后,根据当前状态的Q值来选取能够获取最大收益的动作。

    雷达子脉冲级频点决策往往对应于指数级增长的动作空间,而传统Q学习基于Q表存储和查找Q值,维护难度巨大。而DQN利用神经网络拟合值函数,替换了传统Q表的存储方式,有效解决了高维状态和动作空间的寻优问题。

    DQN与Q学习的主要区别在于网络部分,其采用目标值网络和估计值网络组成的双网络。估计值Q网络输出Q(st,at;θ),用来评估当前状态动作对的未来累计回报期望。目标值ˆQ网络输出ˆQ(st+1,at+1;θ),并根据贝尔曼最优方程,使用y=rt+γmaxˆQ(st+1,at+1;θ)表示Q函数的优化目标。其网络训练过程如图3所示。

    图  3  DQN网络参数的更新过程
    Figure  3.  The network parameter update process of DQN

    输入当前状态st,通过估计值网络预测得到当前状态st对应的不同动作atQ值,然后通过ε-贪婪原则选择at并转至下一状态st+1,同时获得rt。通过目标值网络计算下一状态st+1的最大ˆQ值,将其与估计值作差更新估计值网络参数θ,表示为

    L=[rt+γmaxaˆQ(st+1,at+1;θ)Q(st,at;θ)] (6)

    其中,ε-贪婪原则以概率1ε选择估计值网络输出的具有最大Q值的频点,以概率ε随机选择频点,并随着训练步数的增加减小ε,从而达到搜索和利用的充分结合。

    上述流程经过一定次数后,基于软更新来更新目标值网络参数θ

    θ=τθ+(1τ)θ (7)

    其中,0<τ1表示软间隔更新系数。由于在一段时间内目标值具有一定稳定性,这能在一定程度上降低估计值Q网络和目标值ˆQ网络之间的耦合性,提升了网络的收敛性和稳定性。

    训练完成后,测试时直接输入当前时刻状态至训练好的模型中,即可获取最优动作。

    上述提及的状态、动作和奖励是强化学习的关键要素,其中状态和奖励是算法的输入,动作是算法的输出。设置如下:

    (1) 状态空间:假设雷达能够通过干扰感知等手段获取干扰频点信息,则状态空间由雷达子脉冲频点和干扰频点组成。

    S=[fR,t1,fJ,t]=[fsub1,t1,fsub2,t1,,fsubN,t1,fJ,t] (8)

    其中,fR,t1=[fsub1,t1,fsub2,t1,,fsubN,t1]fJ,t分别表示t1时刻雷达N个子脉冲的频点选择以及t时刻干扰瞄准频点。fJ,t取值范围为1(N+1)1N表示干扰机发射窄带瞄频干扰的瞄准频点,(N+1)表示干扰机发射宽带阻塞干扰。fsubn,t(1nN)的取值范围为1N,表示第n个子脉冲的频点。

    (2) 动作空间:t时刻雷达N个子脉冲频点选择:

    A=fR,t=[fsub1,t,fsub2,t,,fsubN,t] (9)

    (3) 奖励函数:奖励函数应当围绕雷达作战任务设置,本文以预警雷达为例,采用目标检测结果Fd和信干噪比(Signal-to-Jamming-plus-Noise Ratio, SJNR)作为评价指标。前者直接反映了目标检测能力,而后者的存在加快了最优解的收敛速度,降低收敛到局部最优解的可能,从而最大化目标检测性能。定义如下:

    R=Ncomn=1(Nrep,nFd,nSJNRn/SJNRnNcomNcom) (10)
    SJNRn={(PT,nˉPJN,n)/η,Fd,n=10,Fd,n=1 (11)

    其中,对于目标检测结果Fd,我们可以根据提前获取的战场态势信息预估目标距离波门,在子脉冲脉压后基于单元平均恒虚警率(Cell Average-Constant False Alarm Rate, CA-CFAR)检测判断目标能否被检测到[24]。如果第n个子脉冲检测到目标则Fd,n=1,反之则Fd,n=1。同时可以获取目标平均功率PT,n和干扰噪声平均功率ˉPJN,nη为归一化系数,用来将信干噪比限制在0~1之间,从而提高训练稳定性。

    结合状态、动作和奖励的定义,基于深度Q网络的雷达子脉冲频点决策流程如算法1所示。

    1  基于深度Q网络的雷达子脉冲频点决策
    1.  Radar sub-pulse frequency decision based on Deep Q-Network (DQN)
     Step 1:初始化:
      Step 1-1:使用随机参数θ初始化估计值Q网络
      Step 1-2:使用参数θ=θ初始化目标值ˆQ网络
      Step 1-3:初始化经验池D
      Step 1-4:初始化干扰策略,雷达子脉冲数量及频点,折扣因
      子γ,学习率α,贪婪因子ε,软间隔更新系数τ等参数
     Step 2:每幕:
     Step 2-1:设置初始状态s1=[fR,0,fJ,1]
     Step 2-2:每个时间步:
      Step 2-2-1:使用ε-贪婪原则依据估计值网络的输出结果选择
      各子脉冲频点at=fR,t=[fsub1,t,fsub2,t,,fsubN,t],即以
      1ε概率选择估计值网络输出的最佳的频点或者以ε概率随
      机选择频点
      Step 2-2-2:雷达发射子脉冲频率捷变波形,接收到回波后,感
      知得到下一时刻状态st+1并根据目标检测结果和脉压后的信
      干噪比评估当前时刻奖励rt
      Step 2-2-3:将(st,at,rt,st+1)存储到经验池D中,如果经验池
      中的样本数超出预定数量,则删除早期训练样本数据,以便存
      储并使用最新样本数据
      Step 2-2-4:如果经验池D中保存数量超过起始值,则从D中选
      择批大小(batchsize)个样本作为训练集输入到估计值和目标值
      网络中,分别计算得到Q(st,at;θ)y=rt+γmaxˆQ(st+1,
      at+1;θ),并反向梯度求导使误差函数L(θ)=[yQ(st,at;
      θ)]2趋近0,更新估计值网络参数θ
      Step 2-2-5:每隔一定的时间步软更新目标值网络参数θ
     Step 2-3:结束该时间步
     Step 2-4:降低贪婪概率ε
     Step 3:结束该幕
    下载: 导出CSV 
    | 显示表格
    4.1.1   仿真参数设置

    本文以3个子脉冲和3个可选频点为例,讨论DQN应用于子脉冲频点自适应选取的可行性。为避免子脉冲脉压后出现虚假目标,非相邻子脉冲不能选取重复频点,因此动作总数为336=21。频率捷变信号、干扰、DQN的仿真参数分别如表1表3所示。其中,每幕表示1个相参处理间隔(Coherent Processing Interval, CPI),时间步t表示某个CPI中的第t个脉冲重复周期。

    表  1  频率捷变信号参数设置
    Table  1.  The parameter settings of frequency agile signal
    参数 数值
    子脉冲调制类型 LFM
    子脉冲个数 3
    子脉冲频点 [10 MHz, 30 MHz, 50 MHz]
    子脉冲脉宽 5 μs
    子脉冲带宽 5 MHz
    信噪比 0 dB
    下载: 导出CSV 
    | 显示表格
    表  2  干扰参数设置
    Table  2.  The parameter settings of jamming
    干扰类型 参数 数值
    窄带瞄频 瞄准频点 [10 MHz, 30 MHz, 50 MHz]
    带宽 10 MHz
    干噪比 35 dB
    宽带阻塞 带宽 120 MHz
    干噪比 30 dB
    下载: 导出CSV 
    | 显示表格
    表  3  DQN参数设置
    Table  3.  The parameter settings of DQN
    参数 数值
    批大小 64
    学习率 0.001
    折扣因子 0.99
    缓冲区大小 10000
    起始训练样本量 64
    贪婪因子衰减系数 0.2
    32个时间步
    目标值网络更新周期 4个时间步
    目标值网络软间隔更新系数 0.01
    隐藏层数量 2
    隐藏层神经元个数 64
    归一化系数 80
    下载: 导出CSV 
    | 显示表格

    很重要的一个技巧是,本文在基于贪婪原则随机选取动作时,只考虑所有子脉冲选择相同频点的情况,即脉内不跳频。该处理旨在尽可能提高相参处理增益以及使干扰机侦收到单频信号并诱导其发射窄带瞄频干扰,从而加快最优策略学习。同样出于加速收敛的目的,输入到神经网络的奖励按照子脉冲个数进行了归一化。

    估计值网络和目标值网络的结构相同,均使用4层全连接神经网络,分别为输入层、2个隐藏层和输出层。其中,隐藏层的神经元个数均为64,并使用ReLU作为激活函数,如图4所示。

    图  4  全连接神经网络结构示意图
    Figure  4.  The schematic diagram of fully connected neural network structure
    4.1.2   干扰策略设置

    考虑一个具备侦收功能的干扰机,并根据侦-干时间长短分别设置了脉内侦干和脉间侦干等两种固定干扰策略,分别如图5图6所示。由于切片转发干扰的对抗效果受限于切片宽度、转发次数等参数,灵活的参数变化可能会导致对抗失效,因此本文考虑的干扰类型为压制干扰,包括窄带瞄频和宽带阻塞。其中,窄带瞄频干扰的带宽为雷达子脉冲带宽的2倍,更宽的带宽会使得全部状态的奖励值发生整体偏移,但在归一化后会消除该影响。

    图  5  脉内侦干策略
    Figure  5.  The intra-pulse interception-jamming strategy
    图  6  脉间侦干策略
    Figure  6.  The pulse-to-pulse interception-jamming strategy

    对于脉内侦干策略,假设干扰机侦收到雷达脉冲上升沿及下降沿,立即对其测频,转发对应频点的窄带瞄频干扰。值得注意的是,干扰时长设置略小于1个脉冲重复周期(Pulse Repetition Time, PRT),从而使得在当前PRT会同时受到上一时刻以及当前时刻的干扰。因此,雷达在该干扰策略下的一种较为合适的选择为后续子脉冲发射不同于子脉冲1的雷达频点,并且每个PRT均保持相同的发射策略。由于干扰所在频点在滤波后可能会在邻近频点上存在干扰功率残留,因此最优策略为雷达后续子脉冲跳频到距离子脉冲1所选频点的最远频点上。即雷达最优频点选择为[1,N,N][N,1,1]

    对于脉间侦干策略,假设干扰机从侦收到第1个子脉冲开始持续侦收一段时间,直至没有检测到子脉冲时侦收结束。根据侦收结果发射一段时间长度的干扰,干扰时长在3~4个PRT之间。相比脉内侦干策略,后者不会在某个PRT同时受到两部分干扰。在侦收阶段若只侦收到1个频点,则发射对应频点的窄带瞄频干扰,反之则发射宽带阻塞干扰。雷达需要尽量避免干扰机发射宽带阻塞干扰,为此雷达需要在干扰机侦收阶段时只发射单频信号,而在干扰阶段时选择其余频点。类似地,考虑到滤波引起的干扰功率残留,在干扰机侦收时雷达最优策略为[1,1,1][N,N,N],对应的干扰时雷达最优策略为[N,N,N][1,1,1]

    值得注意的是,脉间侦干策略虽然具备周期性,但当前时刻的干扰动作不完全取决于上一时刻的状态,而是按照固定的时序执行侦收和干扰,因此不具备马尔科夫性。脉间侦干策略寻求的是由4个PRT组成的侦干周期的最大奖励,满足式(5)所示的贝尔曼最优方程的价值迭代原理,因此可以使用强化学习解决。

    此时干扰机侦收到1个子脉冲的上升沿与下降沿后,完成测频并立刻发射干扰,雷达频点对抗的训练结果如图7所示。得分曲线在第4个CPI左右即可收敛,在36分附近波动,如图7(a)所示。图7(b)展示了文献[20]提出的基于PPO与LSTM相结合的频点决策算法,其至少需要30幕的时间才能提升到32分附近震荡,因此策略学习耗时且鲁棒性较差。其本质原因在于PPO为on-policy算法,只能利用神经网络进行动作搜索,导致探索性不足,所以存在收敛速度慢、可能会收敛到局部最优解、得分无法保持等诸多问题。

    图  7  脉内侦干策略的子脉冲频点决策训练结果
    Figure  7.  The training results of sub-pulse frequency decision for the intra-pulse interception-jamming strategy

    根据图7(a)的收敛情况,保存前10个CPI的训练模型,每个模型对抗100幕,对抗成功率如图8所示。根据4.1.2节对脉内侦干策略的分析,雷达应将未被侦收到的子脉冲频点设置为距离侦收频点的最远频点。因此,PRT对抗成功定义为{fR=[1,3,3]&fJ=1}{fR=[3,1,1]&fJ=3},即21个动作中只有2个动作为最优,占比9.5%。CPI对抗成功的判决依据是当前CPI内所有PRT均对抗成功。

    图  8  训练用CPI数量对脉内侦干策略下对抗成功率的影响
    Figure  8.  The impact of the number of CPI used for training on the success rate of confrontation for the intra-pulse interception-jamming strategy

    发现训练所用CPI数量对对抗成功率的影响与收敛情况基本对应,从第3个CPI开始,对抗成功率即可达到100%。

    表4展示了随机频点、PPO-LSTM和DQN的单次对抗(PRT)成功率,单幕(CPI)对抗成功率。随机频点决策的成功率与最优动作占比,即理论值大致相同。基于PPO的频点决策虽然在第2个和第3个子脉冲避开了干扰频点,但是由于其搜索力度不够,有一定概率选取到次优策略。而基于DQN的频点决策算法由于使用了ε-贪婪算法,大大扩展了动作搜索空间,更容易收敛到最优策略。

    表  4  脉内侦干策略的对抗成功率(%)
    Table  4.  The success rate of confrontation for the intra-pulse interception-jamming strategy (%)
    策略PRT对抗成功率CPI对抗成功率
    随机频点9.70
    PPO949
    DQN100100
    下载: 导出CSV 
    | 显示表格

    PPO算法由于可以处理连续动作空间问题,并且可以学习到随机策略,因此是强化学习中受众面最广的基线方法。然而在本文研究的频点决策场景中,不涉及连续动作空间,最优策略也可以由随机策略退化到确定性策略,因此PPO算法优势没有得到充分利用。更为重要的是,由于每幕对抗中次优策略不低于最优策略得分的10%,大大提高了仅依靠神经网络参数进行动作搜索的最优策略收敛难度。

    图9(a)展示了雷达和干扰在4个PRT下的频点选取情况。对于第1个PRT,由于初始状态的随机性,雷达选取频点[1,2,3],干扰瞄准频点1。由于单个子脉冲的信噪比增益有限,因此除被干扰的子脉冲外,另有1个子脉冲未能检测到目标,奖励为负值,如图9(b)所示。在第2, 3, 4个PRT,基于训练好的模型,雷达的第2个和第3个子脉冲均选择离干扰频点1最远的频点3,降低了干扰剩余能量的同时,合成了宽脉冲,提高了信噪比增益。

    图  9  雷达与干扰对抗4个PRT的策略及对抗奖励
    Figure  9.  The strategies and rewards for radar anti-jamming during four PRT periods

    最优动作的时频图及一维距离像如图10所示。当前PRT会同时收到瞄准上一时刻第1个子脉冲以及瞄准当前时刻第1个子脉冲的窄带瞄频干扰,后者会在瞄准后立即发射。因此,第1个子脉冲脉压后,目标尖峰出现在当前时刻产生的大功率噪声干扰边缘,导致漏检。第2个子脉冲由于跳频策略与干扰频域正交,因此脉压后能够检测到目标尖峰,具有较高的信干噪比。

    图  10  雷达执行最优策略的时频图及一维距离像
    Figure  10.  The time-frequency map and the one-dimensional High-Resolution Range Profile (HRRP) for radar executing optimal strategy

    本文围绕目标检测性能,基于单个PRT能否检测到目标以及脉压后的信干噪比两方面评价跳频抗干扰效能。表5展示了蒙特卡洛1000次下,雷达的几个典型频点选取策略的目标检测率、脉压后的信干噪比以及平均得分。为便于分析,假设当前时刻和上一时刻均干扰相同的频点,频点[3,1,1]和[1,3,3]为本文所提模型的策略。可以看出:

    表  5  脉内侦干策略下各种雷达策略对抗1000次结果(fJ=fsub1)
    Table  5.  The results of 1000 confrontations with various radar strategies for the intra-pulse interception-jamming strategy (fJ=fsub1)
    雷达频点选择 目标检测率(%) 信干噪比(dB) 平均得分
    [1,1,1] 0 –3.00
    [1,1,2] 0 11.09 –1.12
    [1,1,3] 0 12.25 –0.96
    [1,2,2] 97.6 15.20 1.09
    [1,2,3] 81.7 12.78 0.78
    [1,3,3] 99.7 16.06 1.19
    [2,1,1] 98.3 15.35 1.12
    [2,1,3] 75.6 12.47 0.64
    [2,3,3] 97.7 15.19 1.10
    [3,1,1] 99.6 16.07 1.18
    注:综合考虑噪声随机性引起的得分波动情况,加粗项为最优策略
    下载: 导出CSV 
    | 显示表格

    (1) 由于在当前PRT能同时受到上一时刻和当前时刻的干扰,因此至少有一个雷达频点会被干扰到。根据式(10)所示的奖励函数计算方式,最大得分始终小于2;

    (2) 当子脉冲2和子脉冲3跳频成功时,两个子脉冲均选择离干扰频点的最远频点时,平均得分最高,为最优策略,即[1,3,3]和[3,1,1];

    (3) 诸如[1,2,3]和[2,1,3]等传统频点选取策略,由于脉压增益有限,导致目标检测率较低;而[1,2,2]和[2,1,1]等选择了干扰频点相邻频点的动作,由于滤波后的干扰能量残余,从而降低了信干噪比,非最优策略;

    (4) 次优策略和最优策略的单次对抗得分仅差0.06,网络能够捕获到细微差异,收敛到最优解。

    针对脉间侦干策略,DQN和PPO的训练曲线如图11所示。DQN在第15幕(CPI)左右即可收敛,得分在37分附近。而PPO的训练过程虽然整体呈现上升-平稳,但是其波动始终较为剧烈,且至少需要400幕左右才能趋于平稳。

    图  11  脉间侦干策略的子脉冲频点决策训练结果
    Figure  11.  The training results of sub-pulse frequency decision for the pulse-to-pulse interception-jamming strategy

    图12展示了训练所用CPI数量对对抗成功率的影响,蒙特卡洛次数为100幕。由于雷达初始频点随机选取,不参与决策,因此去除包含初始状态在内的第1个干扰侦干周期。从第2个周期开始统计,即每幕(CPI)对抗28次。根据4.1.2节对脉间侦干策略的分析,雷达应始终发射单频信号,并在干扰机对当前脉冲侦收干扰后的下个脉冲跳到另一频点,从而诱导干扰机在后续干扰周期内发射窄带瞄频干扰,避免发射宽带阻塞干扰导致跳频手段失效。由于干扰机可以在侦收后立即发射对应频点的干扰,所以每个侦干周期内,无论采取何种手段,至少会存在1个PRT抗干扰失败。因此可以仅针对剩余PRT计算抗干扰成功率,将PRT对抗成功定义为干扰机处于发射干扰阶段时雷达选取到最优策略,即{fJ=3&fR=[1,1,1]}{fJ=1&fR=[3,3,3]};CPI对抗成功的判决依据是当前CPI内所有PRT均对抗成功。

    图  12  训练用CPI数量对脉间侦干策略对抗成功率的影响
    Figure  12.  The impact of the number of CPI used for training on the success rate of confrontation for the pulse-to-pulse interception-jamming strategy

    可以发现,在前20个CPI的训练过程中模型学习到的策略不是一直向好,而是波动变化。在第13个PRT策略出现了明显恶化,这与图11(a)的训练结果相一致。此时模型尚未稳定学习到干扰机的侦干策略,因此仍主要处于试错探索阶段。从第15~20个CPI,模型探索到干扰机策略,并学习到有效对抗策略,保持稳定。

    100次蒙特卡洛仿真下的随机频点、PPO和DQN决策的单次对抗(PRT)成功率,单幕(CPI)对抗成功率如表6所示。由于对抗成功率隐含雷达在干扰机侦-干PRT和干扰PRT均发射不同的单频信号,因此随机频点选择的成功概率极低,仅有0.7%。相比PPO,DQN动作搜索更加充分,使对抗成功率得到有效提高,达到100%。

    表  6  脉间侦干策略的对抗成功率(%)
    Table  6.  The success rate of confrontation for the pulse-to-pulse interception-jamming strategy (%)
    策略 PRT对抗成功率 CPI对抗成功率
    随机频点 0.7 0
    PPO 93.6 31
    DQN 100 100
    下载: 导出CSV 
    | 显示表格

    图13(a)展示了干扰机的3个侦干周期下的雷达子脉冲频点选取和干扰瞄准频点。在第1个侦干周期中,由于雷达初始状态的随机性,3个子脉冲分别选取不同频点,导致干扰机在接下来的3个PRT中发射宽带阻塞干扰,此时无论雷达如何跳频,目标均未被检测到,奖励为负值,如图13(b)所示。在第2个侦干周期的第1次对抗中,雷达3个子脉冲均选择频点1,干扰机侦收到并立刻发射对应频点的干扰,因此第1个PRT的奖励为负值。接下来的3个PRT,干扰机继续发射频点1,而雷达选择离频点1最远的频点3。至此第2个侦干周期结束,雷达频点选取成功。在第3个侦干周期中,雷达和干扰的频点选取对调,雷达仍然能够通过频点决策选择受到干扰最小的频点。

    图  13  对抗3个侦干周期的雷达策略及对抗奖励
    Figure  13.  The strategies and rewards for radar anti-jamming during three interception-jamming periods

    以干扰瞄准频点1为例,蒙特卡洛1000次,统计各种策略对抗的目标检测率、脉压后的信干噪比以及平均得分,如表7所示,其中频点[3,3,3]为本文所提模型的策略。可以看出:

    表  7  脉间侦干策略下各种雷达策略对抗1000次的结果(fJ=1)
    Table  7.  The results of 1000 confrontations with various radar strategies for the pulse-to-pulse interception-jamming strategy (fJ=1)
    雷达频点选择 目标检测率(%) 信干噪比(dB) 平均得分
    [1,1,1] 0 –3.00
    [1,2,3] 81.3 12.74 0.76
    [2,2,2] 99.7 17.08 3.17
    [3,3,3] 100 17.58 3.22
    注:加粗项表示最优策略
    下载: 导出CSV 
    | 显示表格

    (1) 对于传统雷达跳频策略[1,2,3],有1个子脉冲会被干扰到,此时奖励虽然为正值,但是较低;

    (2) 对于[2,2,2],虽然从频点数值上看确实跳频成功,但此时瞄准频点1的干扰功率可能未被全部滤掉,有很少一部分的功率会溢出到频点2,使得其信干噪比略低于频点3;

    (3) 当雷达所有子脉冲均选择频点3时,接收到的干扰平均功率达到最小值,平均得分最高,为最优策略。

    针对瞄准式压制干扰,本文面向频率捷变雷达,提出了一种基于深度强化学习的频点自适应快速选取方法。根据当前时刻干扰状态,以及上一时刻雷达动作,依靠神经网络自适应选取当前时刻最优雷达频点,并基于目标检测结果以及脉压后的信干噪比作为奖励反馈,迭代改进策略。仿真部分考虑了具备侦收-瞄准-干扰功能的干扰机,证明了通过关键要素设计可以以单个时间步长作为输入学习到干扰策略的时序性。同时,所用DQN算法配合贪婪准则实现了搜索-利用的平衡,配合信干噪比的反馈加速最优抗干扰策略收敛,相比PPO算法收敛速度提升至少10倍。考虑到实际场景中,干扰频点在滤波后可能在邻近频点存在能量残余的情况,所提频率捷变波形设计方法允许子脉冲多次重复选取距离干扰频点最远的雷达频点,有效降低了回波中的干扰剩余能量,提高了信干噪比。同时扩展了动作空间,提供了最优动作选取的基础。

    通过本文研究发现,当子脉冲数或脉冲数较多时,增大了网络的搜索和决策空间,使得收敛时间进一步增加,并且提高了最优策略的收敛难度。但这不会影响强化学习的关键要素设计,因此所提方法仍能根据交互数据的反馈结果进行策略优化。另外,考虑到子脉冲间、脉冲间的相位不一致,在积累时会带来一定程度上的增益损失。因此在未来的研究中,考虑将子脉冲以及脉冲间的积累情况纳入到奖励函数中,从而指导策略选取。

  • 图  1  MV3DSAR系统实物图

    Figure  1.  A photo of MV3DSAR system

    图  2  全极化MV3DSAR实验系统框图

    Figure  2.  Full-polarimetric MV3DSAR system block diagram

    图  3  苏研院光学影像及目标分布图(A—E区域分别为角反、几何散射体、树林、地下通道入口、车辆、工地塔吊)

    Figure  3.  Optical image and target distribution map (area A—E are corner reflector, geometric scatterer, forest, underpass entrance, vehicle, tower crane)

    图  4  MV3DSAR飞行航迹设计

    Figure  4.  MV3DSAR flight track design

    图  5  MV3DSAR校飞实验基线设计

    Figure  5.  Baseline design of MV3DSAR

    图  6  观测矩阵互相关特性及空间模糊函数曲线

    Figure  6.  Cross-correlation properties of observation matrix and spatial ambiguity function curve

    图  7  最大不模糊高度及瑞利分辨率

    Figure  7.  Maximum unblurred height and Rayleigh resolution

    图  8  激光点云数据

    Figure  8.  Lidar point cloud data

    图  9  光学倾斜摄影结果

    Figure  9.  Optical oblique photography data

    图  10  全极化MV3DSAR数据处理总体流程

    Figure  10.  The overall flow of full-polarimetric MV3DSAR data processing

    图  11  基于能量最大极化投影的极化相干增强流程图

    Figure  11.  Flow chart of polarization coherent enhancement based on energy maximum polarization projection

    图  12  基于极化相似性约束的三维成像流程

    Figure  12.  Flow chart of 3D imaging based on polarization similarity constraint

    图  13  SAR三维可视化产品生产流程图

    Figure  13.  Flow chart of SAR 3D visual product production

    图  14  解叠掩图制作流程图

    Figure  14.  Flow chart of scatterer-unmixed image production

    图  15  手动选择场景中行道树区域作为极化标校对象

    Figure  15.  The street tree area in the scene is manually selected as the polarization calibration object

    图  16  场景在极化校正前后的Pauli伪彩图

    Figure  16.  Pauli pseudo-color image of the scene before and after polarization correction

    图  17  三面角反射器的极化失真评估

    Figure  17.  Evaluation of polarization distortion of a corner reflector

    图  18  单极化HH与本文提出的PEM方法的相干系数统计图(图例中的下角标代表不同通道组合)

    Figure  18.  Statistical diagram of the coherence coefficients of the single polarization HH and the proposed PEM method (the lower corner marks in the legend represent different channel combinations)

    图  19  R1R4通道能量统计图(图例中的下角标代表通道序号)

    Figure  19.  Energy statistical diagram of channel R1R4 (the lower corner mark in the legend represents the channel number)

    图  20  三维成像结果对比(整体)

    Figure  20.  Overall comparison of SAR 3D point clouds

    图  21  三维成像结果对比(配套楼屋顶细节)

    Figure  21.  Detailed SAR 3D point clouds comparison of Area 1

    图  22  某个飞行航迹下的解叠掩图

    Figure  22.  A set of 2D scatterer-unmixed images under a flight track

    图  23  某个飞行航迹下的三维融合可视化产品

    Figure  23.  A 3D fusion visualization product under a flight track

    表  1  全极化Ku-SAR载荷参数

    Table  1.   Parameters of full-polarimetric Ku-SAR payload

    序号 参数名称 参数值或内容
    1 中心频率 15.2 GHz
    2 信号形式 调频连续波(FMCW)
    3 极化方式 HH/HV/VH/VV
    4 信号带宽 1200 MHz
    5 天线尺寸(单通道) 0.05 m(俯仰)×0.32 m(方位)
    6 每个极化的阵列通道数 4
    7 分辨率 优于0.2 m×0.2 m
    8 天线波束宽度 方位≥4° 俯仰≥24°
    9 天线极化隔离度 优于25 dB
    10 通道相位不平衡稳定度 ±5° (10 min内)
    11 通道幅度不平衡稳定度 ±0.2 dB (10 min内)
    12 中心视角 45°
    13 NESZ 不大于–30 dB (最远
    作用距离3.6 km)
    14 Ku-SAR重量 主机、存储、电池、天线、结构等
    一共5.7 kg
    下载: 导出CSV

    表  2  三维重建精度对比

    Table  2.   Comparison of 3D reconstruction accuracy

    三维成像方法 点云三维重建精度(m)
    HH单极化 1.51
    VV单极化 1.46
    无约束的全极化 1.36
    极化相似性约束的全极化 0.93
    下载: 导出CSV
  • [1] KNAELL K K and CARDILLO G P. Radar tomography for the generation of three-dimensional images[J]. IEE Proceedings-Radar, Sonar and Navigation, 1995, 142(2): 54–60. doi: 10.1049/ip-rsn:19951791.
    [2] ZHU Xiaoxiang and BAMLER R. Superresolving SAR tomography for multidimensional imaging of urban areas: Compressive sensing-based TomoSAR inversion[J]. IEEE Signal Processing Magazine, 2014, 31(4): 51–58. doi: 10.1109/MSP.2014.2312098.
    [3] TEBALDINI S, NAGLER T, ROTT H, et al. Imaging the internal structure of an alpine glacier via L-band airborne SAR tomography[J]. IEEE Transactions on Geoscience and Remote Sensing, 2016, 54(12): 7197–7209. doi: 10.1109/TGRS.2016.2597361.
    [4] HUANG Yue, FERRO-FAMIL L, and REIGBER A. Under-foliage object imaging using SAR tomography and polarimetric spectral estimators[J]. IEEE Transactions on Geoscience and Remote Sensing, 2012, 50(6): 2213–2225. doi: 10.1109/TGRS.2011.2171494.
    [5] FORNARO G, LOMBARDINI F, and SERAFINO F. Three-dimensional multipass SAR focusing: Experiments with long-term spaceborne data[J]. IEEE Transactions on Geoscience and Remote Sensing, 2005, 43(4): 702–714. doi: 10.1109/TGRS.2005.843567.
    [6] 张福博. 阵列干涉SAR三维重建信号处理技术研究[D]. [博士论文], 中国科学院大学, 2015.

    ZHANG Fubo. Research on signal processing of 3-D reconstruction in linear array synthetic aperture radar interferometry[D]. [Ph.D. dissertation], University of Chinese Academy of Sciences, 2015.
    [7] REIGBER A and MOREIRA A. First demonstration of airborne SAR tomography using multibaseline L-band data[J]. IEEE Transactions on Geoscience and Remote Sensing, 2000, 38(5): 2142–2152. doi: 10.1109/36.868873.
    [8] ZHU Xiaoxiang and BAMLER R. Super-resolution power and robustness of compressive sensing for spectral estimation with application to spaceborne tomographic SAR[J]. IEEE Transactions on Geoscience and Remote Sensing, 2012, 50(1): 247–258. doi: 10.1109/TGRS.2011.2160183.
    [9] ZHU Xiaoxiang, GE Nan, and SHAHZAD M. Joint sparsity in SAR tomography for urban mapping[J]. IEEE Journal of Selected Topics in Signal Processing, 2015, 9(8): 1498–1509. doi: 10.1109/JSTSP.2015.2469646.
    [10] 丁赤飚, 仇晓兰, 徐丰, 等. 合成孔径雷达三维成像——从层析、阵列到微波视觉[J]. 雷达学报, 2019, 8(6): 693–709. doi: 10.12000/JR19090.

    DING Chibiao, QIU Xiaolan, XU Feng, et al. Synthetic aperture radar three-dimensional imaging—from TomoSAR and array InSAR to microwave vision[J]. Journal of Radars, 2019, 8(6): 693–709. doi: 10.12000/JR19090.
    [11] 仇晓兰, 焦泽坤, 杨振礼, 等. 微波视觉三维SAR关键技术及实验系统初步进展[J]. 雷达学报, 2022, 11(1): 1–19. doi: 10.12000/JR22027.

    QIU Xiaolan, JIAO Zekun, YANG Zhenli, et al. Key technology and preliminary progress of microwave vision 3D SAR experimental system[J]. Journal of Radars, 2022, 11(1): 1–19. doi: 10.12000/JR22027.
    [12] JIAO Zekun, DING Chibiao, QIU Xiaolan, et al. Urban 3D imaging using airborne TomoSAR: Contextual information-based approach in the statistical way[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2020, 170: 127–141. doi: 10.1016/j.isprsjprs.2020.10.013.
    [13] JIAO Zekun, QIU Xiaolan, DONG Shuhang, et al. Preliminary exploration of geometrical regularized SAR tomography[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2023, 201: 174–192. doi: 10.1016/j.isprsjprs.2023.05.019.
    [14] 朱庆涛, 殷君君, 曾亮, 等. 基于邻域一致性的极化SAR图像仿射配准[J]. 雷达学报, 2021, 10(1): 49–60. doi: 10.12000/JR20120.

    ZHU Qingtao, YIN Junjun, ZENG Liang, et al. Polarimetric SAR image affine registration based on neighborhood consensus[J]. Journal of Radars, 2021, 10(1): 49–60. doi: 10.12000/JR20120.
    [15] 仇晓兰, 罗一通, 程遥, 等. SAR微波视觉三维成像数据集-3.0[EB/OL]. 雷达学报. https://radars.ac.cn/web/data/getData?newsColumnId=1cbc9f2d-f2ee-4748-9972-748c007f697f, 2024.

    QIU Xiaolan, LUO Yitong, CHENG Yao, et al. SAR microwave vision 3D imaging Dataset 3.0[EB/OL]. Journal of Radars. https://radars.ac.cn/web/data/getData?newsColumnId=2f2748db-10ef-4ad0-bcc4-f087ce59b6f8&pageType=en, 2024.
    [16] 徐牧, 王雪松, 肖顺平. 基于改善极化相似性的极化SAR目标增强新方法[J]. 电子与信息学报, 2008, 30(5): 1047–1051. doi: 10.3724/SP.J.1146.2007.00754.

    XU Mu, WANG Xuesong, and XIAO Shunping. Target enhancement in POL-SAR imagery based on the improvement of polarization characteristics similarity[J]. Journal of Electronics & Information Technology, 2008, 30(5): 1047–1051. doi: 10.3724/SP.J.1146.2007.00754.
    [17] JIANG Sha, QIU Xiaolan, HAN Bing, et al. A quality assessment method based on common distributed targets for GF-3 polarimetric SAR data[J]. Sensors, 2018, 18(3): 807. doi: 10.3390/s18030807.
    [18] LEE J S and POTTIER E. Polarimetric Radar Imaging: From Basics to Applications[M]. Boca Raton: CRC Press, 2009. doi: 10.1201/9781420054989.
    [19] SONG Shujie, QIU Xiaolan, and SHANGGUAN Songtao. Study on the three dimension imaging methods of fully-polarised array InSAR[C]. IET International Radar Conference, Chongqing, China, 2023: 2999–3003. doi: 10.1049/icp.2024.1571.
    [20] SONG Shujie and QIU Xiaolan. A sparse Bayesian learning 3d imaging methodology based on polarimetric energy maximum in urban area for pol-array-insar[C]. IGARSS Conference, 2024, Athens, Greece, 2024. doi: 10.1109/IGARSS53475.2024.10641994.
  • 期刊类型引用(1)

    1. 王兴家,王彬,刘岳巍,晏学成,丁峰. 基于元知识转移的认知雷达波形设计. 雷达科学与技术. 2024(04): 443-453 . 百度学术

    其他类型引用(4)

  • 加载中
图(23) / 表(2)
计量
  • 文章访问数: 1276
  • HTML全文浏览量: 397
  • PDF下载量: 656
  • 被引次数: 5
出版历程
  • 收稿日期:  2024-07-05
  • 修回日期:  2024-09-02
  • 网络出版日期:  2024-09-18
  • 刊出日期:  2024-10-28

目录

/

返回文章
返回