Processing math: 0%
  • 中文核心期刊要目总览
  • 中国科技核心期刊
  • 中国科学引文数据库(CSCD)
  • 中国科技论文与引文数据库(CSTPCD)
  • 中国学术期刊文摘数据库(CSAD)
  • 中国学术期刊(网络版)(CNKI)
  • 中文科技期刊数据库
  • 万方数据知识服务平台
  • 中国超星期刊域出版平台
  • 国家科技学术期刊开放平台
  • 荷兰文摘与引文数据库(SCOPUS)
  • 日本科学技术振兴机构数据库(JST)
微信公众号

微信公众号

基于块稀疏恢复的高频雷达空间角谱和极化谱联合估计方法

许彬, 董英凝, 毛兴鹏, 赵春雷, 刘昱含

许彬, 董英凝, 毛兴鹏, 赵春雷, 刘昱含. 基于块稀疏恢复的高频雷达空间角谱和极化谱联合估计方法[J]. 电波科学学报, 2019, 34(6): 741-750. doi: 10.13443/j.cjors.2019042901
引用格式: 许彬, 董英凝, 毛兴鹏, 赵春雷, 刘昱含. 基于块稀疏恢复的高频雷达空间角谱和极化谱联合估计方法[J]. 电波科学学报, 2019, 34(6): 741-750. doi: 10.13443/j.cjors.2019042901
XU Bin, DONG Yingning, MAO Xingpeng, ZHAO Chunlei, LIU Yuhan. Joint estimation of spatial angular spectrum and polarization spectrum of high frequency radar based on block sparse recovery[J]. CHINESE JOURNAL OF RADIO SCIENCE, 2019, 34(6): 741-750. doi: 10.13443/j.cjors.2019042901
Reference format: XU Bin, DONG Yingning, MAO Xingpeng, ZHAO Chunlei, LIU Yuhan. Joint estimation of spatial angular spectrum and polarization spectrum of high frequency radar based on block sparse recovery[J]. CHINESE JOURNAL OF RADIO SCIENCE, 2019, 34(6): 741-750. doi: 10.13443/j.cjors.2019042901

基于块稀疏恢复的高频雷达空间角谱和极化谱联合估计方法

基金项目: 

国家自然科学基金重点项目 61831009

详细信息
    作者简介:

    许彬  (1995-), 男, 内蒙古人, 哈尔滨工业大学电子与信息工程学院在读硕士研究生, 主要研究方向为阵列信号处理、极化抗干扰技术

    董英凝  (1969-), 女, 黑龙江人, 哈尔滨工业大学电子与信息工程学院副教授、硕士生导师, 主要研究方向为新体制雷达技术和现代信号处理

    毛兴鹏  (1972-), 男, 辽宁人, 哈尔滨工业大学电子与信息工程学院教授、博士生导师, 主要研究方向为电子侦察与电子对抗、弱信号检测、雷达信号处理

    赵春雷  (1992-), 男, 黑龙江人, 哈尔滨工业大学电子与信息工程学院在读博士研究生, 主要研究方向为阵列信号处理、高频地波雷达抗干扰

    刘昱含  (1998-), 女, 黑龙江人, 东南大学信息科学与工程学院在读本科生, 主要研究方向为信号处理

    通信作者:

    毛兴鹏 E-mail:mxp@hit.edu.cn

  • 中图分类号: TN957

Joint estimation of spatial angular spectrum and polarization spectrum of high frequency radar based on block sparse recovery

  • 摘要: 高频地波雷达(high-frequency surface wave radar,HFSWR)中电离层杂波的存在,会极大地影响雷达系统的性能,降低对目标的探测能力.为了精确获得杂波参数,从而更好地抑制电离层杂波,提出了一种基于压缩感知(compressive sensing,CS)的单快拍参数估计方法,用于对电离层杂波的空域和极化域参数估计.该方法基于极化敏感阵列的块稀疏估计模型,应用块正交匹配追踪(block orthogonal matching pursuit,BOMP)算法实现距离-多普勒域的单快拍空间角度和极化参数联合估计,并进一步获得目标和杂波的空间角谱和极化谱.该方法适用于任意极化敏感阵列,在距离-多普勒域单快拍条件下,其估计性能优于传统方法,且计算复杂度极低,可以实现实时处理.仿真结果和某HFSWR系统实测数据处理结果表明了参数估计方法的有效性.
    Abstract: The presence of ionospheric clutter in high-frequency surface wave radar(HFSWR) may extremely affect the performance of radar system and degrade the capability of detecting target. In order to obtain clutter parameters accurately and suppress ionospheric clutter better, a single snapshot parameter estimation method based on compressive sensing(CS) is proposed for estimating the spatial domain and polarization domain parameters of ionospheric clutter. Based on the block sparse estimation model of polarization sensitive array, this method uses block orthogonal matching pursuit(BOMP) algorithm to realize joint estimation of angle and polarization parameters of single snapshot in range-Doppler domain, and further obtains the spatial angular spectrum and polarization spectrum of target and clutter. For arbitrary polarization-sensitive arrays, the performance of estimation is better than that of traditional methods under the condition of single snapshot in range-Doppler domain, and the computational complexity is very low, so real-time processing can be realized. The simulation results and the data processing results of a HFSWR system show that the parameter estimation method is effective.
  • 高频地波雷达(high-frequency surface wave radar, HFSWR)能够以较少的能量损失对海上目标和低空目标实现超视距检测.高频(3~30 MHz)垂直极化波可以在海面完成低损失、远距离传播, HFSWR正是利用这一特性对船只、飞行器甚至巡航导弹进行超视距检测, 从而提供目标的角度、距离和径向速度等信息[1].

    虽然HFSWR拥有超视距探测这一优势, 但不可避免地会受到电离层杂波、海杂波、电台干扰等其他干扰的影响, 这降低了HFSWR的探测性能, 阻碍了它的发展.其中, 以能量强、占时长、覆盖多个连续距离单元和多普勒单元等多种特点的电离层杂波为主, 极大地影响了高频雷达系统的工作性能[2].由于实际高频雷达系统所采用的是非理想的发射天线, 导致部分能量向上扩散而非水平传播, 这部分能量经电离层反射后形成电离层杂波.由于等离子体的漂移和电子浓度的时刻变化, 以及散射信号的存在, 造成电离层杂波的多普勒频移和多普勒展宽[3].文献[1]详细介绍了每一层的分布位置和杂波特性, 同时提出扩散Es层的概念, 并给出了该层的杂波特性.

    极化是继时域、频域、空域信息后, 雷达可使用的另一种重要信息, 根据以往的研究, 目标和电离层杂波除了在俯仰角上存在差异外, 在极化上也表现出明显的区别[4-5], 因此将其用于估计电离层参数甚至抑制电离层杂波对雷达系统的影响方面具有巨大的研究价值[6-7].由于信号的极化信息和方位信息耦合在阵列流型中, 因此当前极化参数估计的过程可以与信号的波达方向(direction of arrival, DOA)估计联合进行.自20世纪90年代, 先后有学者设法将多重信号分类(multiple signal classification, MUSIC)[8]、矩阵束MUSIC[9]、ESPRIT[10-11]等经典的高分辨参数估计方法拓展到极化敏感阵列上来, 从而实现DOA和极化参数的联合估计.近些年, 也有研究者引入四元数[12]、几何代数[13]、张量代数[14]等高维代数理论来提高参数估计的鲁棒性.文献[15]使用稀疏表示和压缩感知(compressive sensing, CS)方法来实现算法适用性的推广和估计性能上的提高, 但该方法仍然要求足够的时间快拍.

    具体到实际高频雷达系统中, 要实现精确的DOA和极化参数估计, 仍面临以下几方面特殊的挑战:

    首先, 传统HFSWR系统多采用垂直极化天线进行单极化接收, 而考虑到尽可能降低各方面成本, 近年来在对原系统进行多极化改造时, 一般选择仅简单地在原接收阵列旁增设水平极化天线[5], 因此往往无法设定某些特殊的阵列排布以适应算法的应用需求; 且由于当前天线工艺的限制, 共点多极化天线仍很难达到令人满意的低互耦和高极化隔离度要求[16].因此, 诸如仅适用于共点多极化阵元及空间旋转不变性阵列的极化ESPRIT算法, 以及要求各阵元均为相同的共点多极化天线的四元数、张量代数等参数估计方法, 仍难以得到实际工程应用.

    其次, 由于接收天线阵元数往往有限, 但雷达回波内包含的目标及干扰个数极多, 直接基于回波时域数据进行目标和干扰的参数估计并不可行.因此, 参数估计通常在距离-多普勒域(range-Doppler, RD)处理后进行, 一方面将目标在距离多普勒域区分开后, 单个单元内目标个数少, 另一方面相干积累后信噪比(signal to noise ratio, SNR)得到明显提升[17].但随之而来的问题是, 积累后可获得的距离多普勒域快拍数量极少(甚至往往仅为1), 而只有在较大快拍数的条件下, 基于阵列协方差矩阵的传统参数估计方法(极化MUSIC等)性能才有保证.文献[18]在距离-多普勒域上使用了基于CS的方法获取空域参数, 能够实现单快拍DOA估计, 而极化估计在测角之后独立进行, 但由于测角时仅使用了阵列中的部分垂直极化天线, 因此性能有所损失.

    最后, 由于参数估计要在每个感兴趣的距离-多普勒单元上分别进行, 且待估计参数包括二维DOA及极化等共四维参数, 因此对算法的实时性也有很高的要求.

    综上, 发展一种适用于任意极化敏感阵列形式、在少快拍或者单快拍条件下性能较高,且复杂度较低的DOA和极化参数联合估计方法十分必要.

    为了精确获取电离层的杂波参数, 本文对电离层杂波的角度信息和极化信息进行了联合估计.首先对雷达回波信号进行距离变换和速度变换, 得到包含各接收通道的相位差信息的距离-多普勒域处理结果(本文中以线性调频断续波为例进行说明), 接着在距离-多普勒域上基于块稀疏重建分别对各单元进行DOA和极化参数联合估计, 从而获得空间角谱和极化谱.该方法可适用于任意空间和极化排布的极化敏感阵列, 仅需距离-多普勒域的单次快拍即可快速实现相应单元的DOA和极化参数联合估计, 算法复杂度极低, 能够实时实现大数据量的电离层杂波参数估计.通过仿真分析, 验证了在任意极化敏感阵列形式、单快拍条件下, 本文方法较极化MUSIC方法在参数估计性能上的优势.最后通过对实测数据的处理, 基于电离层的空域和极化域参数估计结果, 对电离层杂波的方向特性和极化特性进行了分析.

    极化敏感阵列是指极化敏感阵元按一定排列方式在空间放置所构成的阵列系统, 其通过极化敏感阵元获知空间电磁信号的极化信息, 利用阵列几何结构完成空域采样获取信号的空域信息[19].这里以下文实测数据对应的HFSWR实际接收阵列为例对接收信号模型进行说明, 需要注意的是, 下文所述模型及参数估计算法并不局限于此阵列, 而是适用于任意极化敏感阵列.极化敏感接收阵列模型如图 1所示, 阵列为16个阵元组成的L型三极化敏感阵列.其中, 平行于海岸线的8个阵元为主阵, 垂直于海岸线的8个阵元为辅阵, 主阵的阵元间距为d1, 辅阵的阵元间距为d2, 主阵与辅阵间距为d3. 图 1中主阵阵元全部为垂直极化, 用圆圈表示, 辅阵中2、3、5号阵元为水平极化, 且2号阵元垂直于x轴摆放, 用矩形表示; 3号和5号阵元平行于x轴摆放, 用三角形表示; 其余阵元为垂直极化.此外, 图 1中的θ∈[0°, 90°]表示俯仰角, 指向天顶时为0°; φ∈[0°, 360°]表示方位角, 指向海面正前方(主阵法线方向)时为180°

    图  1  极化敏感阵列模型
    Fig.  1  Polarization sensitive array model

    对于完全极化波, 阵元接收的信号的极化矢量可表示为[19]

    {\mathit{\boldsymbol{s}}_{\rm{p}}} = \left[ {\begin{array}{*{20}{c}} {{E_x}}\\ {{E_y}}\\ {{E_z}}\\ {{H_x}}\\ {{H_y}}\\ {{H_z}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{u}}&\mathit{\boldsymbol{v}}\\ \mathit{\boldsymbol{v}}&{ - \mathit{\boldsymbol{u}}} \end{array}} \right] \cdot \left[ {\begin{array}{*{20}{c}} {\cos \gamma }\\ {{{\rm{e}}^{{\rm{j}}\eta }}\sin \gamma } \end{array}} \right]. (1)

    式中:矢量uv

    \left\{ \begin{array}{l} \mathit{\boldsymbol{u}} = {\left[ {\begin{array}{*{20}{c}} { - \sin \varphi }&{\cos \varphi }&0 \end{array}} \right]^{\rm{T}}}\\ \mathit{\boldsymbol{v}} = {\left[ {\begin{array}{*{20}{c}} {\cos \theta \cos \varphi }&{\cos \theta \sin \varphi }&{ - \sin \theta } \end{array}} \right]^{\rm{T}}} \end{array} \right.; (2)

    γ∈[0°, 90°]表示极化角; η∈[0°, 360°]表示极化相差, 参数对(γ, η)与电磁波的极化状态一一对应.

    假设信号是窄带远场信号, 某时刻极化敏感阵列接收数据的极化域-空域联合导向矢量表示为

    {\mathit{\boldsymbol{s}}_{{\rm{s}},{\rm{p}}}} = {\mathit{\boldsymbol{s}}_{\rm{p}}} \otimes {\mathit{\boldsymbol{s}}_{\rm{s}}}. (3)

    式中:sssp分别为空域导向矢量和极化域导向矢量; ⊗表示矩阵的kronecker积.

    由于任意部分极化波(或完全极化波)均可分解为极化状态正交且不相关的两个完全极化分量[19], 显然, 任意极化状态为(γ, η)的完全极化波的极化矢量也可改写为水平极化导向矢量sph(γ=0°)和垂直极化导向矢量spv(γ=90°, 假设η=0°)的线性组合:

    \begin{array}{l} {\mathit{\boldsymbol{s}}_{\rm{p}}}\left( {\theta ,\varphi ,\gamma ,\eta } \right) = {\mathit{\boldsymbol{s}}_{{\rm{ph}}}}\left( {\theta ,\varphi } \right)\cos \gamma + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\mathit{\boldsymbol{s}}_{{\rm{pv}}}}\left( {\theta ,\varphi } \right){{\rm{e}}^{{\rm{j}}\eta }}\sin \gamma , \end{array} (4)
    \left\{ \begin{array}{l} {\mathit{\boldsymbol{s}}_{{\rm{ph}}}}\left( {\theta ,\varphi } \right) = {\left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{u}}&\mathit{\boldsymbol{v}} \end{array}} \right]^{\rm{T}}}\\ {\mathit{\boldsymbol{s}}_{{\rm{pv}}}}\left( {\theta ,\varphi } \right) = {\left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{v}}&{ - \mathit{\boldsymbol{u}}} \end{array}} \right]^{\rm{T}}} \end{array} \right.. (5)

    综上, 任意完全极化波的极化域-空域联合导向矢量ss, p(简记为s)可改写成水平矢量和垂直矢量的线性组合:

    \begin{array}{l} \mathit{\boldsymbol{s}} = {\mathit{\boldsymbol{s}}_{{\rm{ph}}}} \otimes {\mathit{\boldsymbol{s}}_{\rm{s}}} \cdot \cos \gamma + {\mathit{\boldsymbol{s}}_{{\rm{pv}}}} \otimes {\mathit{\boldsymbol{s}}_{\rm{s}}} \cdot {{\rm{e}}^{{\rm{j}}\eta }}\sin \gamma \\ \;\; = {\mathit{\boldsymbol{s}}_{{\rm{s}},{\rm{ph}}}} \cdot \cos \gamma + {\mathit{\boldsymbol{s}}_{{\rm{s}},{\rm{pv}}}} \cdot {{\rm{e}}^{{\rm{j}}\eta }}\sin \gamma \\ \;\; = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{s}}_{\rm{h}}}}&{{\mathit{\boldsymbol{s}}_{\rm{v}}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\cos \gamma }\\ {{{\rm{e}}^{{\rm{j}}\eta }}\sin \gamma } \end{array}} \right]. \end{array} (6)

    式中:将水平矢量ss, ph简记为sh; 垂直矢量ss, pv简记为sv, 下文将沿用这种表示.

    对于高频雷达系统而言, 其时域回波SNR低、目标及杂波分量复杂且数目众多, 因此需预先进行距离-多普勒处理.目前雷达发射信号以线性调频类信号的应用最为广泛, 其中常用的调频断续波有两种:均匀脉冲和伪随机码[20-21], 本文中以前者为例进行说明.在一个调频周期内均匀脉冲截断的线性调频断续波的信号形式为[17]

    \mu \left( t \right) = g\left( t \right){{\rm{e}}^{{\rm{j}}2{\rm{ \mathsf{ π} }}\left( {{f_0}t + \frac{1}{2}k{t^2}} \right)}},0 \le t < {T_{{\rm{sweep}}}}. (7)

    式中:f0为载波频率; k为调频斜率; Tsweep为调频周期; g(t)为均匀截断脉冲, 其表达式为

    g\left( t \right) = \sum\limits_{n = 1}^{{N_{\rm{p}}}} {{\rm{rect}}\left( {\frac{{t - n{T_{{\rm{pp}}}} - \frac{{{T_{{\rm{pw}}}}}}{2}}}{{{T_{{\rm{pw}}}}}}} \right)} . (8)

    式中:Np为单个调频周期内脉冲个数; Tpp为脉冲重复周期; Tpw为脉冲宽度.

    对于阵元数为M的平面阵, 目标数为N, 则到达第m个阵元的第n个目标回波smn(t)(不考虑幅度衰减)可表示为

    {s_{mn}}\left( t \right) = s\left( {t - \Delta {\tau _n} - {\tau _{mn}}} \right). (9)

    式中:m=1, 2, …, M; n=1, 2, …, N; Δτn为该目标回波到达参考阵元时相对于发射信号的时间延迟, τn=2(Rn-υnt)/c, Rn, υn分别为第n个目标的距离和径向速度, c为光速; τmn为回波信号到达第m个阵元的位置(x, y)处时相对于参考阵元的延迟, 可表示为

    {\tau _{mn}} = \left( {x\sin {\theta _n}\cos {\varphi _n} + y\sin {\theta _n}\sin {\varphi _n}} \right)/c. (10)

    式中:θnφn分别为第n个目标的俯仰角和方位角.

    在窄带远场信号源的假设下, 有

    {s_{mn}}\left( t \right) = s\left( {t - \Delta {\tau _n}} \right){{\rm{e}}^{ - {\rm{j}}2{\rm{ \mathsf{ π} }}{f_0}{\tau _{mn}}}}. (11)

    经过混频、低通滤波处理后, 解调信号可表示为

    {{s'}_{mn}}\left( t \right) = {{s'}_n}\left( t \right){{\rm{e}}^{ - {\rm{j}}2{\rm{ \mathsf{ π} }}{f_0}{\tau _{mn}}}}, (12)
    {{s'}_n}\left( t \right) = g\left( {t - \Delta {\tau _n}} \right){{\rm{e}}^{{\rm{j}}2{\rm{ \mathsf{ π} }}\left( {{f_0}\Delta {\tau _n} + kt\Delta {\tau _n} - k\Delta \tau _n^2/2} \right)}}. (13)

    则接收信号在第m个阵元通道可表示为

    {x_m}\left( t \right) = \sum\limits_{n = 1}^N {{{s'}_n}\left( t \right){{\rm{e}}^{ - {\rm{j}}2{\rm{ \mathsf{ π} }}{f_0}{\tau _{mn}}}} + {n_m}\left( t \right)} . (14)

    式中:nm(t)为第m个阵元通道中的噪声.

    对采样后的离散数据xm(i)进行距离变换(记符号为R)和速度变换(记符号为D), 并忽略其中的噪声部分, 有

    \begin{array}{l} {{X'}_m}\left( {r,d} \right) = DFT{\left[ {DFT{{\left[ {{x_m}\left( i \right)} \right]}_r}} \right]_d}\\ \;\;\;\;\;\;\;\;\;\;\;\;\; = \sum\limits_{n = 1}^N {{\rm{D}}{{\left[ {{\rm{R}}{{\left[ {{{s'}_n}\left( t \right)} \right]}_r}} \right]}_d}{{\rm{e}}^{ - {\rm{j}}2{\rm{ \mathsf{ π} }}{f_0}{\tau _{mn}}}}} . \end{array} (15)

    式中:r为距离单元数; d为多普勒单元数.

    由公式(15)可知, 距离多普勒变换过程具有线性性质, 得到的RD域数据中仍包含各接收通道的相位差信息e-j2πf0τmn, 因此对各阵元分别进行RD处理后提取相应RD单元构成数据快拍, 即可用于回波参数估计.

    显然, 对于一段固定时长的回波数据而言, 既可直接对其进行RD处理获得单次快拍, 也可以分段后进行多次RD处理获得多个快拍.而文献[17]证明, 由于相干积累的SNR改善更为显著, 在RD域进行单快拍参数估计可以达到最高速度分辨率和最佳参数估计性能的统一.但由于传统DOA和极化参数估计方法多针对特殊的阵列形式(如空间排布满足旋转不变性或阵元为同种共点多极化天线)且要求较多的快拍数, 在本文针对的高频雷达参数估计中实用性较差, 下文中我们提出了一种适用于任意极化敏感阵列的基于块稀疏重建的单快拍DOA和极化参数估计方法.

    由于RD处理后每个RD单元内目标或干扰个数极少, 在整个感兴趣的空域范围内仅有少数方向上存在信号入射, 因此对于单个RD单元而言, 其中包含的入射信号在空域上存在一定的稀疏性, 鉴于此, 可以建立如下基于RD域数据快拍的块稀疏表示模型.首先, 将整个空间范围的俯仰角、方位角分别划分为NθNφ份, 对应的组合(θ, φ)可分为Nθ, φ=Nθ×Nφ份, 构造冗余字典

    \mathit{\boldsymbol{s'}} = \left[ {\begin{array}{*{20}{l}} {{{\mathit{\boldsymbol{s'}}}_1}}&{{{\mathit{\boldsymbol{s'}}}_2}}& \cdots &{{{\mathit{\boldsymbol{s'}}}_{{N_{\theta ,\varphi }}}}} \end{array}} \right] \in {{\bf{C}}^{M \times 2{N_{\theta ,\varphi }}}}. (16)

    式中, s'n=[sh(θ, φ)n  sv(θ, φ)n]∈CM×2称为字典的块, n=1, 2, …, Nθ, φ. (θ, φ)所构成的垂直极化导向矢量和水平极化导向矢量都分别对应字典中的某一块, 且该块内两个极化矢量的线性组合可以表示该角度对应的任意完全极化信号.结合式(6), 可得(r, d)单元对应的RD域快拍数据块稀疏表示模型如下:

    \mathit{\boldsymbol{x}} = \mathit{\boldsymbol{s'\alpha }} + \mathit{\boldsymbol{n}}. (17)

    式中:α=[α1  α2  …  αNθ, φ]TC2Nθ, φ×1为待恢复的稀疏矢量, 其中各块元素为

    \begin{array}{l} {\mathit{\boldsymbol{\alpha }}_n} = \left[ {\begin{array}{*{20}{c}} {{a_{n,{\rm{h}}}}}&{{a_{n,{\rm{v}}}}} \end{array}} \right]\\ \;\;\;\; = \left[ {\begin{array}{*{20}{c}} {{\alpha _n}\cos {\gamma _n}}&{{\alpha _n}{{\rm{e}}^{{\rm{j}}{\eta _n}}}\sin {\gamma _n}} \end{array}} \right], \end{array} (18)

    αn=D[R[s'n(i)]r]d为距离-多普勒变换后的信号分量; n为噪声.

    显然, 如果第n个角度网格上的信号分量αn为零, 则矢量α中的第nαn必然整体为零, 因此, α是一个块稀疏矢量.如果信号存在, 且位于第n个俯仰角、方位角的组合(θ, φ)n处, 则与s'n对应的αn为非零块.如果信源数为K, 则α中会有K个非零块.因此只要块稀疏矢量α能够精确恢复, 那么就可以通过其中非零块的位置确定信号的入射方向, 通过其中元素的幅相信息根据公式(18)确定信号的极化状态, 进而实现空域和极化域参数的联合估计.

    显然, 以上块稀疏表示模型适用于任意极化敏感阵列形式, 而基于该模型, 采用适当的块稀疏重建算法即可实现目标和杂波的参数估计.由于目标及杂波可能出现在大量的RD单元上, 而参数估计必须对每个单元分别进行, 因此要尽可能满足处理的实时性, 对算法的复杂度要求也较高.故而本文采用复杂度极低的块正交匹配追踪(block orthogonal matching pursuit, BOMP)算法[22]进行块稀疏重建, 以实现空域和极化域参数估计.

    该算法首先将残差初始化为R(0)=x, 第k(k≥1)次迭代选择与残差R(k-1)最匹配的块s'i(k)的位置, 其位置索引为

    {i^{\left( k \right)}} = \mathop {\arg \max }\limits_{1 \le i \le {N_{\theta ,\varphi }}} \sum\limits_{l - 1}^L {{{\left\| {\mathit{\boldsymbol{s'}}_i^{\rm{H}}{\mathit{\boldsymbol{R}}^{\left( {k - 1} \right)}}{\mathit{\boldsymbol{e}}_l}} \right\|}_2}} . (19)

    式中, el是单位矩阵IL×L的第l列.然后, 在得到位置索引后, 即可求出该次迭代下对信号的估计值:

    {{\tilde \alpha }^{\left( k \right)}} = \arg \mathop {\min }\limits_\mathit{\boldsymbol{\alpha }} {\left\| {\mathit{\boldsymbol{x}} - {\mathit{\boldsymbol{s}}^{\left( k \right)}}\mathit{\boldsymbol{\alpha }}} \right\|_2}. (20)

    s(k)满足列满秩矩阵条件时, 以上最小二乘问题的解可由下式给出:

    {{\mathit{\boldsymbol{\tilde \alpha }}}^{\left( k \right)}} = {\left( {{\mathit{\boldsymbol{s}}^{\left( k \right)}}} \right)^\dagger } \cdot \mathit{\boldsymbol{x}}. (21)

    式中, †表示矩阵的广义逆.

    最后, 对残差值进行更新:

    {\mathit{\boldsymbol{R}}^{\left( k \right)}} = \mathit{\boldsymbol{x}} - {\mathit{\boldsymbol{s}}^{\left( k \right)}}{{\mathit{\boldsymbol{\tilde \alpha }}}^{\left( k \right)}}. (22)

    表 1给出BOMP算法的主要步骤.

    表  1  BOMP算法步骤
    Tab.  1  The steps of BOMP algorithm
      输入:观测数据x, 冗余字典s', 块稀疏度K
      输出:非零块的估计值\mathit{\boldsymbol{\widetilde \alpha }}, 残差R(k)
      初始化:迭代次数k=0, s(0)为空矩阵, R(0)=x, 索引集P(0)
      步骤1:第k(k≥1)次迭代, 记录与残差R(k-1)最匹配的块s'i(k)的位置, 其对应的索引可由公式(19)来表示;
      步骤2:更新原子块索引集P(k)=P(k-1)∪{i(k)}以及已选原子块集合s(k)=[s(k-1)  s'i(k)];
      步骤3:根据公式(20)求解最小二乘问题, 得到该次迭代下对信号的估计值{\mathit{\boldsymbol{\widetilde \alpha }}^{\left( k \right)}};
      步骤4:利用公式(22)更新残差值R(k);   
      步骤5:若k=K, 则结束迭代, 估计值\mathit{\boldsymbol{\widetilde \alpha }} = {\mathit{\boldsymbol{\widetilde \alpha }}^{\left( k \right)}}; 否则, k=k+1, 重复步骤1-5
    下载: 导出CSV 
    | 显示表格

    利用Pk中对应的位置索引及已知的字典构造方式就可以明确信号对应的俯仰方位组合, 从而实现DOA估计.

    根据公式(18), 通过重建的块稀疏信号的估计值\mathit{\boldsymbol{\widetilde \alpha }}, 即可由式(23)、(24)进一步实现极化角和极化相差的估计:

    {\gamma _k} = \arctan \left( {\left| {\frac{{{{\tilde \alpha }_{k,{\rm{v}}}}}}{{{{\tilde \alpha }_{k,{\rm{h}}}}}}} \right|} \right), (23)
    {\eta _k} = \angle \left( {\frac{{{{\tilde \alpha }_{k,{\rm{v}}}}}}{{{{\tilde \alpha }_{k,{\rm{h}}}}}}} \right). (24)

    式中:k=1, …, K; ∠(β)表示求β的相位.

    对于RD谱内所有感兴趣的距离-多普勒单元, 分别利用上述方法进行DOA和极化参数联合估计, 将得到的估计结果显示在相应RD单元上, 即可获得高频雷达空间角谱和极化谱.

    由于实际阵列为L型非共点三极化敏感阵列, 而传统方法中仅极化MUSIC方法可不受阵列形式限制, 与本文所提单快拍BOMP算法一样可用于此阵列, 下面给出这两种算法在高斯白噪声下的仿真实验.其中, 所用仿真信号是频率为4.6 MHz、俯仰角为60°、方位角为110°、极化角为25°、极化相差为130°的单快拍单个信号源, 因HFSWR回波信号经相参积累后SNR一般可达到10 dB以上, 所以取SNR的范围为0~25 dB, 阵列形式与图 1类似, 但阵元间距均为半波长, 进行500次蒙特卡洛循环.其中, RMSE的定义为

    {\rm{RMSE}} = \sqrt {\frac{1}{N}\sum\limits_{n = 1}^N {{{\left( {{{\mathit{\boldsymbol{\tilde \alpha }}}_n} - \mathit{\boldsymbol{\xi }}} \right)}^2}} } . (25)

    式中:N=500为蒙特卡洛次数; ξ为待估计参数, 包括俯仰角、方位角、极化角、极化相差; {\mathit{\boldsymbol{\widetilde \alpha }}_n}为第n次的参数估计值.

    不同SNR下, 极化MUSIC和单快拍BOMP算法的参数估计性能与克拉美·罗界(Cramer-Rao bound, CRB)[23]的对比如图 2图 3所示, 二者的运行时间对比如图 4所示.

    图  2  不同SNR下算法的角度估计均方根误差
    Fig.  2  RMSE of angle estimation vs. SNR
    图  3  不同SNR下算法的极化参数估计均方根误差
    Fig.  3  RMSE of polarization parameter vs. SNR
    图  4  不同SNR下算法的运行时间
    Fig.  4  Running time of algorithms vs. SNR

    图 2可以看出:在SNR为0~10 dB的范围内, 单快拍BOMP算法的角度估计性能均优于极化MUSIC算法; 随着SNR的提高, 二者的方位角估计性能基本一致, 但在较高SNR条件(>15 dB)下, BOMP俯仰角的RMSE曲线略有偏离CRB, 其主要原因是随着SNR的提高, 为减小因网格划分引起的角度量化误差, 在极化MUSIC谱峰搜索和BOMP构造字典时采用了更细的角度网格, 导致字典的相关性增强, 并进一步限制了BOMP算法这类稀疏重构类方法的性能提升[24], 而由于俯仰角维度对字典相关性的影响更大, 因此其RMSE曲线在高SNR下不能更好地贴近CRB.

    图 3容易看出, 单快拍BOMP算法的极化角估计性能明显优于极化MUSIC, 基本上贴近于CRB, 而对于极化相差的估计性能则略差于极化MUSIC, 但我们在解决实际问题时基本不会用到极化相差, 因此可以不做考虑.

    图 4可见, 接收信号的SNR的大小并不会影响算法的运行时间, 但是单快拍BOMP算法明显要比极化MUSIC算法耗时小一个数量级, 而由于极化MUSIC算法需要对协方差矩阵进行奇异值分解, 因此在更大阵元数条件下算法的运算效率对比将更为明显, 非常有利于对大批量的电离层杂波数据进行处理, 在较短的时间内获得较好的参数估计结果.

    在3.1.1的仿真中, 阵元间距为半波长31.9 m, 但实际阵列受场地大小等因素的影响, 阵元间距较小, 约为\frac{1}{6}\lambda (λ为波长), 这会在一定程度上影响参数估计的性能, 所以这里给出不同阵元间距的L型非共点三极化敏感阵列对参数估计性能影响的仿真实验, 用以更好地贴近实际雷达接收阵列, 对后续实测数据的处理结果进行解释.仿真中SNR为10 dB, 阵元间距为\frac{1}{7}\lambda \sim \frac{1}{2}\lambda , 其他参数与3.1.1中的仿真实验一致, 这里不再赘述.仿真结果如图 5图 6所示.

    图  5  不同阵元间隔下算法的角度估计均方根误差
    Fig.  5  RMSE of angle estimation vs. element spacing
    图  6  不同阵元间隔下算法的极化参数估计均方根误差
    Fig.  6  RMSE of polarization parameter vs. element spacing

    图 5易知, 随着阵元间距由较小的间隔逐渐接近半波长, 两种算法的估计精度都不断提高, 而且对于方位角而言单快拍BOMP算法要比极化MUSIC算法的估计精度高一些, 对于俯仰角则二者性能接近.

    图 6可以看出, 由于阵列的极化排布不变, 而阵元空间排布对极化参数的估计性能影响较小, 随着阵元间距的增大, 极化角和极化相差的估计精度没有明显提高, 但是单快拍BOMP算法对极化角的估计性能明显优于极化MUSIC算法, 对极化相差的估计性能略差, 这也与前面的仿真结果相一致.

    为验证基于CS的单快拍测角测极化联合估计算法用于电离层杂波的有效性, 现对沿海某HFSWR的实测数据进行处理.首先对实测数据进行距离-多普勒变换, 得到该数据的某个垂直极化通道形成的原始RD谱如图 7所示, 某个水平极化通道形成的原始RD谱如图 8所示.图中横坐标表示速度门(即多普勒单元), 纵坐标表示距离门(即距离单元), 各个距离-多普勒单元的功率根据其数值大小用不同的颜色表征.从图中可以看出, 电离层杂波主要位于第150距离门到第180距离门之间, 并且占据连续的速度门, 具有明显的特征.而目标信号则在RD谱内通常表现为占据不超过5个RD单元的十字形式, 其余为海杂波形成的Bragg峰和其他电台干扰等.

    图  7  某垂直极化通道的RD谱
    Fig.  7  The range-Doppler spectrum of a vertical polarization channel
    图  8  某水平极化通道的RD谱
    Fig.  8  The range-Doppler spectrum of a horizontal polarization channel

    为减少数据计算量, 在图 7图 8距离-多普勒谱的基础上, 设定一个能量门限, 对大于该门限的RD单元做基于CS的空域-极化域参数联合估计, 得到该单元的角度估计值和极化估计值.然后我们将得到的每一个RD单元上信号的俯仰角和方位角根据其数值用不同的颜色来表示, 就得到信号的方位角谱和俯仰角谱, 分别如图 9图 10所示.其中, 图 9中0°对应主阵列法线方向, 指向与海面相反的方向, 即海面正前方对应180°. 图 10中俯仰角0°对应天顶方向, 90°对应水平方向.

    图  9  方位角谱
    Fig.  9  The azimuth spectrum
    图  10  俯仰角谱
    Fig.  10  The pitch angle spectrum

    图 9的方位角谱可以看到, 目标信号和海杂波的方位角主要位于150°~220°, 也就是海面方向, 而电离层杂波的方位角则主要分布于50°~100°和230°~300°, 每个区域中的电离层杂波方位基本一致, 但不同区域的电离层杂波方位角有所不同.其原因是不同区域的电离层杂波可能是由不同的电离层反射产生的, 且各个电离层对应的方位角不同.

    图 10所示的俯仰角谱容易看到, 电离层杂波主要来自天顶方向和近天顶方向, 而海杂波和目标信号则主要来自俯仰角很大的海面方向.显然电离层杂波与目标在俯仰角上的差异更为明显, 因此, 与仅可获得目标一维空间角度的传统雷达线阵相比, 图 10获得的俯仰角信息对于利用空域滤波更好地抑制电离层杂波大有助益.

    类似方位角谱和俯仰角谱, 同样我们将每个RD单元上信号的极化角和极化相差根据数值用不同的颜色表征, 就得到信号在距离-多普勒单元上的极化角谱和极化相差谱, 分别如图 11图 12所示.

    图  11  极化角谱
    Fig.  11  The spectrum of polarization angle
    图  12  极化相差谱
    Fig.  12  The spectrum of polarization phase delay

    图 11可以看到, 电离层杂波的极化角主要分布在50°~60°, 而海杂波和可能存在的海上目标信号的极化角则为近90°, 二者差异较为明显.其原因是, 结合图 10的俯仰角谱可见, 电离层杂波俯仰角较小而且目标回波来自于接近海面方向, 由于水平极化波经海面传播损耗较大, 目标回波中的水平极化分量经过远距离传播后已经被大幅衰减.

    图 12所示的极化相差谱可以看到, 不同层的电离层杂波的极化相差也存在一定的差异.但由于其数值在实测数据中未表现出明显的规律, 因此在杂波抑制处理中仅可作为适当的参考.

    需要指出的是, 在用实测数据进行角度和极化信息的联合估计中, 估计电离层杂波的方位角和俯仰角时, 极化MUSIC对噪声的鲁棒性更好, 往往能得到较好的处理结果, 而在估计极化角时, 由于每个距离多普勒单元上仅采用单个快拍, 因此BOMP算法的性能更佳.

    本文提出了一种基于CS的单快拍参数估计方法, 对电离层杂波的空域和极化域参数进行联合估计.该方法可以适用于任意极化敏感阵列, 仿真结果验证了在距离-多普勒域单快拍条件下, 其估计性能优于传统方法, 且计算复杂度极低, 可以实现实时处理.

    通过对HFSWR回波的实测数据进行处理和分析, 本文对电离层杂波的角度信息和极化信息进行了研究.在距离-多普勒域上将基于CS的单快拍空间角谱和极化谱联合估计方法用于电离层杂波, 验证了距离多普勒上不同区域的电离层杂波方位角不同, 可能来自不同的电离层反射; 而电离层杂波的俯仰角主要来自于近天顶方向, 与目标和海杂波差异显著; 另外, 电离层杂波的极化角与目标和海杂波之间也表现出明显的差异, 且不同电离层的极化相差也存在一定的差异.基于本文方法获得的空间角谱和极化角谱将有助于分析高频雷达回波与电离层杂波在空域和极化域上的差异, 为研究改善电离层杂波抑制方法奠定基础.

  • 图  1   极化敏感阵列模型

    Fig.  1   Polarization sensitive array model

    图  2   不同SNR下算法的角度估计均方根误差

    Fig.  2   RMSE of angle estimation vs. SNR

    图  3   不同SNR下算法的极化参数估计均方根误差

    Fig.  3   RMSE of polarization parameter vs. SNR

    图  4   不同SNR下算法的运行时间

    Fig.  4   Running time of algorithms vs. SNR

    图  5   不同阵元间隔下算法的角度估计均方根误差

    Fig.  5   RMSE of angle estimation vs. element spacing

    图  6   不同阵元间隔下算法的极化参数估计均方根误差

    Fig.  6   RMSE of polarization parameter vs. element spacing

    图  7   某垂直极化通道的RD谱

    Fig.  7   The range-Doppler spectrum of a vertical polarization channel

    图  8   某水平极化通道的RD谱

    Fig.  8   The range-Doppler spectrum of a horizontal polarization channel

    图  9   方位角谱

    Fig.  9   The azimuth spectrum

    图  10   俯仰角谱

    Fig.  10   The pitch angle spectrum

    图  11   极化角谱

    Fig.  11   The spectrum of polarization angle

    图  12   极化相差谱

    Fig.  12   The spectrum of polarization phase delay

    表  1   BOMP算法步骤

    Tab.  1   The steps of BOMP algorithm

      输入:观测数据x, 冗余字典s', 块稀疏度K
      输出:非零块的估计值\mathit{\boldsymbol{\widetilde \alpha }}, 残差R(k)
      初始化:迭代次数k=0, s(0)为空矩阵, R(0)=x, 索引集P(0)
      步骤1:第k(k≥1)次迭代, 记录与残差R(k-1)最匹配的块s'i(k)的位置, 其对应的索引可由公式(19)来表示;
      步骤2:更新原子块索引集P(k)=P(k-1)∪{i(k)}以及已选原子块集合s(k)=[s(k-1)  s'i(k)];
      步骤3:根据公式(20)求解最小二乘问题, 得到该次迭代下对信号的估计值{\mathit{\boldsymbol{\widetilde \alpha }}^{\left( k \right)}};
      步骤4:利用公式(22)更新残差值R(k);   
      步骤5:若k=K, 则结束迭代, 估计值\mathit{\boldsymbol{\widetilde \alpha }} = {\mathit{\boldsymbol{\widetilde \alpha }}^{\left( k \right)}}; 否则, k=k+1, 重复步骤1-5
    下载: 导出CSV
  • [1]

    JIANG W, DENG W B. Characteristic research of ionospheric clutter in over the horizon HFSWR[C]//2009 IET International Radar Conference, 2009.

    [2]

    CHAN H C. Characterization of ionospheric clutter in HF surface-wave radar[R]. Ottawa: Defence R&D, 2003.

    [3] 尚尚, 张宁, 李杨.高频地波雷达电离层杂波统计特性研究[J].电波科学学报, 2011, 26(3):521-527. http://www.cjors.cn/CN/abstract/abstract1389.shtml

    SHANG S, ZHANG N, LI Y. Ionospheric clutter statistical properties in HFSWR[J]. Chinese journal of radio science, 2011, 26(3):521-527.(in Chinese) http://www.cjors.cn/CN/abstract/abstract1389.shtml

    [4] 刘爱军.基于极化信息的高频地波雷达干扰抑制方法研究[D].哈尔滨: 哈尔滨工业大学, 2011.

    LIU A J. Research on interference mitigation methods based on polarization information for high frequency surface wave radar[D]. Harbin: Harbin Institute of Technology, 2011.(in Chinese)

    [5] 洪泓.高频地波雷达多域协同抗电离层杂波干扰方法研究[D].哈尔滨: 哈尔滨工业大学, 2014.

    HONG H. Research on multidomain collaborative ionosphere clutter mitigation methods for high frequency surface wave radar[D]. Harbin: Harbin Institute of Technology, 2014.(in Chinese)

    [6] 代大海, 廖斌, 肖顺平, 等.雷达极化信息获取与处理的研究进展[J].雷达学报, 2016, 5(2):143-155. http://d.old.wanfangdata.com.cn/Periodical/ldxb201602003

    DAI D H, LIAO B, XIAO S P, et al. Advancements on radar polarization information acquisition and processing[J]. Journal of radars, 2016, 5(2):143-155.(in Chinese) http://d.old.wanfangdata.com.cn/Periodical/ldxb201602003

    [7] 赵春雷, 王亚梁, 阳云龙, 等.雷达极化信息获取及极化信号处理技术研究综述[J].雷达学报, 2016, 5(6):620-638. http://d.old.wanfangdata.com.cn/Periodical/ldxb201606004

    ZHAO C L, WANG Y L, YANG Y L, et al. Review of radar polarization information acquisition and polarimetric signal processing techniques[J]. Journal of radars, 2016, 5(6):620-638.(in Chinese) http://d.old.wanfangdata.com.cn/Periodical/ldxb201606004

    [8]

    WONG K T, ZOLTOWSKI M D. Self-initiating MUSIC-based direction finding and polarization estimation in spatio-polarizational beamspace[J]. IEEE transactions on antennas and propagation, 2000, 48(8):1235-1245. doi: 10.1109/8.884492

    [9]

    HUA Y. A pencil-MUSIC algorithm for finding two-dimensional angles and polarizations using crossed dipoles[J]. IEEE transactions on antennas and propagation, 1993, 41(3):370-376. doi: 10.1109/8.233122

    [10]

    LI J, COMPTON Jr R T. Angle and polarization estimation using ESPRIT with a polarization sensitive array[J]. IEEE transactions on antennas and propagation, 1991, 39(9):1376-1383. doi: 10.1109/8.99047

    [11] 徐友根, 刘志文.基于累积量的极化敏感阵列信号DOA和极化参数的同时估计[J].电子学报, 2004, 32(12):1962-1966. doi: 10.3321/j.issn:0372-2112.2004.12.007

    XU Y G, LIU Z W. Cumulant-based two-dimensional DOA and polarization estimation with a polarization sensitive array comprising a spatially stretched tripole[J]. Acta electronica sinica, 2004, 32(12):1962-1966.(in Chinese) doi: 10.3321/j.issn:0372-2112.2004.12.007

    [12]

    MIRON S, LE BIHAN N, MARS J I. Quaternion-MUSIC for vector-sensor array processing[J]. IEEE transactions on signal processing, 2006, 54(4):1218-1229. doi: 10.1109/TSP.2006.870630

    [13]

    XIAO H K, ZOU L, XU B G, et al. Direction and polarization estimation with modified quad-quaternion music for vector sensor arrays[C]//IEEE International Conference on Signal Processing, 2014: 352-357.

    [14]

    GONG X F, LIU Z W, XU Y G. Regularised parallel factor analysis for the estimation of direction-of-arrival and polarization with a single electromagnetic vector-sensor[J]. IET signal processing, 2011, 5(4):390-396. doi: 10.1049/iet-spr.2009.0221

    [15]

    TIAN Y, SUN X, ZHAO S. Sparse-reconstruction-based direction of arrival, polarization and power estimation using a cross-dipole array[J]. IET radar, sonar & navigation, 2015, 9(6):727-731. http://cn.bing.com/academic/profile?id=e51d4517ee4781beb007dd166c1ca33f&encoded=0&v=paper_preview&mkt=zh-cn

    [16]

    PIAO D, YANG L, GUO Q, et al. Measurement-based performance comparison of colocated tripolarized loop and dipole antennas[J]. IEEE transactions on antennas and propagation, 63(8):3371-3379. doi: 10.1109/TAP.2015.2435058

    [17] 赵春雷, 王亚梁, 毛兴鹏, 等.基于压缩感知的高频地波雷达二维DOA估计[J].系统工程与电子技术, 2017, 39(4):733-741. http://d.old.wanfangdata.com.cn/Periodical/xtgcydzjs201704007

    ZHAO C L, WANG Y L, MAO X P, et al. Compressive sensing based two-dimensional DOA estimation for high frequency surface wave radar[J]. Systems engineering and electronics, 2017, 39(4):733-741.(in Chinese) http://d.old.wanfangdata.com.cn/Periodical/xtgcydzjs201704007

    [18]

    YANG Y L, MAO X P, DONG Y N, et al. Space-polarization collaborative suppression method for ionospheric clutter in HFSWR[J]. Journal of radars, 2016, 5(6):673-680. http://d.old.wanfangdata.com.cn/Periodical/pre_1d700613-6309-413f-979b-c9891fe41158

    [19] 徐振海.极化敏感阵列信号处理的研究[D].长沙: 国防科学技术大学, 2004.

    XU Z H. Signal processing based on polarization sensitive array[D]. Changsha: National University of Defense Technology, 2004.(in Chinese)

    [20] 王磊.高频地波雷达信号波形分析与设计[D].哈尔滨: 哈尔滨工业大学, 2007.

    WANG L. Waveform analysis and design for high frequency ground wave radar[D]. Harbin: Harbin Institute of Technology, 2007.(in Chinese)

    [21] 朱朋志.高频雷达干扰信号产生与处理[D].哈尔滨: 哈尔滨工业大学, 2010.

    ZHU P Z. The generation and processing of high frequency radar jamming signals[D]. Harbin: Harbin Institute of Technology, 2010.(in Chinese)

    [22]

    ELDAR Y C, KUPPINGER P, BOLCSKEI H. Block-sparse signals:uncertainty relations and efficient recovery[J]. IEEE transactions on signal processing, 2010, 58(6):3042-3054. doi: 10.1109/TSP.2010.2044837

    [23]

    ZHAO C L, MAO X P, YANG Y L, et al. Deterministic CRB and manifold ambiguity analysis for polarization sensitive arrays[C]//2018 IEEE Radar Conference (RadarConf18), OKC, USA, April 2018: 759-764.

    [24]

    STOICA P, BABU P. Sparse estimation of spectral lines:grid selection problems and their solutions[J]. IEEE transactions on signal processing, 2011, 60(2):962-967. doi: 10.1109/TSP.2011.2175222

  • 期刊类型引用(2)

    1. 王龙岗,岳显昌,吴雄斌,易先洲. 基于奇异值分解的空域海杂波抑制算法. 电波科学学报. 2021(04): 645-652 . 本站查看
    2. 吴双,袁野,吴微微,袁乃昌. 一种宽带相干信源的无网格超分辨DOA估计方法. 电波科学学报. 2020(05): 648-655 . 本站查看

    其他类型引用(3)

图(12)  /  表(1)
计量
  • 文章访问数:  246
  • HTML全文浏览量:  46
  • PDF下载量:  38
  • 被引次数: 5
出版历程
  • 收稿日期:  2019-04-28
  • 网络出版日期:  2020-12-30
  • 发布日期:  2019-12-29
  • 刊出日期:  2019-12-29

目录

/

返回文章
返回