Loading [MathJax]/jax/element/mml/optable/MathOperators.js
高级搜索

留言板

尊敬的读者、作者、审稿人, 关于本刊的投稿、审稿、编辑和出版的任何问题, 您可以本页添加留言。我们将尽快给您答复。谢谢您的支持!

姓名
邮箱
手机号码
标题
留言内容
验证码

基于MIMO阵列的综合射频系统鲁棒波形设计算法研究

唐波 吴文俊 史英春 王旭阳 李达

唐波, 吴文俊, 史英春, 王旭阳, 李达. 基于MIMO阵列的综合射频系统鲁棒波形设计算法研究[J]. 电子与信息学报, 2023, 45(11): 3918-3926. doi: 10.11999/JEIT220969
引用本文: 唐波, 吴文俊, 史英春, 王旭阳, 李达. 基于MIMO阵列的综合射频系统鲁棒波形设计算法研究[J]. 电子与信息学报, 2023, 45(11): 3918-3926. doi: 10.11999/JEIT220969
TANG Bo, WU Wenjun, SHI Yingchun, WANG Xuyang, LI Da. Robust Waveform Design Based on MIMO Multi-Function Radio Frequency Systems under Angle Uncertainties[J]. Journal of Electronics & Information Technology, 2023, 45(11): 3918-3926. doi: 10.11999/JEIT220969
Citation: TANG Bo, WU Wenjun, SHI Yingchun, WANG Xuyang, LI Da. Robust Waveform Design Based on MIMO Multi-Function Radio Frequency Systems under Angle Uncertainties[J]. Journal of Electronics & Information Technology, 2023, 45(11): 3918-3926. doi: 10.11999/JEIT220969

基于MIMO阵列的综合射频系统鲁棒波形设计算法研究

doi: 10.11999/JEIT220969
基金项目: 国家自然科学基金(62171450, 61671453),安徽省杰出青年科学基金(2108085J30)
详细信息
    作者简介:

    唐波:男,博士,副教授,研究方向为雷达信号处理、通信信号处理、MIMO系统和阵列信号处理等

    吴文俊:男,硕士生,研究方向为雷达信号处理、波形设计等

    史英春:男,博士,副教授,研究方向为阵列信号处理、通信信号处理等

    王旭阳:男,硕士生,研究方向为雷达信号处理、波形设计等

    李达:男,博士生,研究方向为雷达信号处理、波形设计等

    通讯作者:

    吴文俊  wuwenjun.1010@nudt.edu.cn

  • 中图分类号: TN959.1

Robust Waveform Design Based on MIMO Multi-Function Radio Frequency Systems under Angle Uncertainties

Funds: The National Natural Science Foundation of China (62171450, 61671453), The Anhui Provincial Natural Science Foundation (2108085J30)
  • 摘要: 基于多输入多输出(MIMO)阵列的综合射频(MFRF)系统通过优化多波形,可以在不同角度辐射不同的信号,进而同时实现多种功能。然而,当存在角度误差时,综合射频系统所合成的信号会出现畸变,这将导致系统性能下降。该文主要考虑存在角度误差时,对综合射频系统鲁棒波形设计算法进行研究。首先,提出基于最小最大化(Min-Max)框架的MIMO阵列多波形优化问题。同时,为减少非线性功放引起的波形失真,对发射波形施加了峰均比(PAPR)约束。为求解此问题,该文提出了基于交替方向乘子法(ADMM)和Majorization-Minimization(MM)方法的波形优化算法。仿真实验结果表明,当存在角度误差时,所提算法能够获得更好的性能。
  • 为适应现代信息化战争的需要,武器平台往往需要配备大量的射频电子设备,这会导致作战平台上天线数目增加、隐身性能下降、电磁兼容问题突出,从而严重影响了作战平台(例如飞机、舰船等)的作战效能[1,2]。为有效缓解上述问题,美国等西方国家提出综合射频(Multi-Function Radio Frequency, MFRF)系统这一概念。综合射频系统利用同一宽带天线孔径,将多种射频传感器进行一体化集成,通过开放式信号处理软件架构同时实现目标探测、数据通信、电子对抗和导航制导等功能。相比于传统的电子信息系统,该系统集成性更高,所需天线数量更少,资源利用率更高,体积更小,重量更轻,成本更低,而且能够有效提高平台的机动能力,极大减少平台的反射截面积,使得平台能够适应更复杂的电磁环境。从20世纪80年代美国海军为减少舰艇顶部天线而研究的舰载先进多功能射频概念(Advanced Multifunction Radio Frequency Concept, AMRFC)[3]开始,到如今正在进行的电磁机动指挥与控制(Electronic Maneuvering Command and Control, EMC2)[4]技术发展计划,综合射频系统逐渐成为研究的热点[5,6]

    现有的综合射频系统主要有3种实现形式:时分复用体制、空分复用体制、波形复用体制。时分复用体制需要进行复杂的时间调度规划,且无法同时实现多种功能。空分复用体制由美国于20个世纪90年代开始研究,经过世界各国多年发展,已经形成了相对完整的技术理论体系,并已在一些武器平台上进行了应用,取得了较好的效果,例如,美国F-35第5代战斗机中使用的基于“宝石台计划”的先进综合射频和电子战系统、我国055型驱逐舰上采用的一体化桅杆设计和综合射频技术。然而,空分复用体制的综合射频系统仍存在天线孔径利用率低、天线增益低、频谱资源占用较多等问题。

    对于波形复用体制,目前的主要研究方向是同时实现双功能,例如雷达通信一体化[7]和雷达干扰一体化[8]等。文献[9]提出可以通过改变雷达脉冲重复频率实现数据通信,但是数据传输效率不高;文献[10]提出基于多载频相位编码信号(Multi-Carrier Phased Coded, MCPC)和P3/P4相位编码信号的一体化波形,并在此基础上提出了用于雷达探测和用户通信的波形设计算法,但是所设计波形存在非恒包络、通信速率较低等不足;文献[11]利用频分准正交多载波Chirp信号设计了一体化波形,文献[12]提出将最小频移键控(Minimum Shift Keying, MSK)信号与线性调频(Linear Frequency Modulation, LFM)相结合的一体化波形。文献[13]通过实验验证了噪声信号可以作为雷达发射波形。

    多输入多输出(Multiple-Input-Multiple-Output, MIMO)阵列的每一个通道均可发射独立波形[14]。MIMO阵列所带来的波形分集能力能够显著提升系统性能,已经成功应用于大规模通信和雷达探测[15,16],随着MIMO阵列的蓬勃发展,一些学者也先后提出了基于MIMO阵列的多功能射频系统。文献[17]通过改变频分正交MIMO雷达中各阵元发射信号的初始频率从而实现数据通信和雷达探测双功能;文献[18-20]提出了既能够匹配雷达期望发射方向图又能够实现数据通信的多波形设计方法。文献[21]通过设计MIMO阵列的多波形,不仅能够与合作方进行通信,而且可以干扰窃听方接收机。文献[22,23]对杂波环境中的雷达通信一体化系统发射波形和接收滤波器联合优化设计进行了研究。

    然而,现有方法大多仅仅考虑实现两种功能,很少研究实现3种或以上功能。此外,在设计系统发射波形时,往往需要使用一些先验信息,例如目标或通信接收机所在的角度等。需要指出的是,角度估计误差会导致所设计的波形性能下降。为了提高存在角度估计误差时的波形性能,本文提出一种综合射频系统鲁棒波形设计算法。仿真结果表明,当存在角度误差时,所提算法合成的信号与期望信号之间的误差极小,能够较好地同时实现雷达探测、数据通信和电子干扰等多种功能。

    本文考虑基于MIMO阵列的综合射频系统,旨在利用MIMO阵列同时实现雷达探测、数据通信和电子干扰等多种功能,如图1所示。

    图 1  基于MIMO阵列的综合射频系统

    设综合射频系统天线阵列的发射天线数为NT,将第(n = 1,2, ···, NT)个天线发射的离散基带信号记为 snCL×1(L为基带信号的码长),则综合射频系统的发射波形矩阵可以表示为 S=[s1,s2,,sNT]T CNT×L。在远场窄带条件下,系统在 θ角度上合成的信号 x(θ)可以表示为

    x(θ)=aH(θ)S=(NTn=1an(θ)sn)T (1)

    其中, a(θ)CNT×1为角度 θ处的发射导引矢量, an(θ)为矢量 aH(θ)中的第 n个元素。可以看出,在 θ角度上合成的信号是发射波形的线性组合,其组合系数是角度 θ的函数。因此,对发射波形 S进行优化设计,可以使得系统在不同方向上合成实现不同功能的信号 x(θ)

    θ方向上的期望信号记为 d(θ)CL×1( d(θ)可以为雷达探测信号、通信信号、电子干扰信号等),那么在该方向上,合成信号与期望信号的匹配误差可以表示为

    e(θ)=||x(θ)dT(θ)||22=||aH(θ)SdT(θ)||22 (2)

    假设感兴趣的方向为 θ1,θ2,,θK,相应的期望信号为 d(θ1),d(θ2),,d(θK)。为最小化合成信号与期望信号之间的匹配误差,考虑波形优化问题为

    minSmaxθk||aH(θk)SdT(θk)||22,SS (3)

    其中, a(θk)=wka(θk),d(θk)=wkd(θk)wk>0为权重, k=1,2,,K, S为发射波形的约束集。值得注意的是,通过在不同方向上引入不同的权重,可以在一定程度上控制不同方向上的匹配误差。权重的设置和应用背景、任务需要密切相关。例如,要降低通信信号的误码率,则需要相应地增加其权重系数,此时可以减少通信信号拟合误差,降低拟合误差导致的信号失真对误码率的影响。

    在实际系统中,发射总能量有限,故对发射波形施加能量约束

    tr(SSH)=eT (4)

    其中, tr()表示矩阵的迹,eT为发射总能量。同时,为避免使用昂贵的线性放大器,通常期望发射波形具有较低的峰值平均功率比(Peak-to-Average-Power-Ratio, PAPR, 简称为峰均比)[24],因此,对发射波形进一步施加峰均比约束

    sHnsn=eT/NT,PAPR(sn)ρ,n=1,2,,NT (5)

    其中, 1ρL

    PAPR(sn)=maxl|sn(l)|21/LLl=1|sn(l)|2,n=1,2,,NT (6)

    其中, sn(l)为基带信号 sn的第 l( l=1,2,,L)个元素。结合式(3)和式(5),发射波形设计问题可以表示为

    minSmaxθk||aH(θk)SdT(θk)||22,s.t.sHnsn=eT/NT,PAPR(sn)ρ,n=1,2,,NT (7)

    在式(7)中,假设了 {θk}Kk=1是精确已知的。然而,由于存在角度估计误差,使用式(7)来设计波形会导致系统性能下降。为了提高存在角度误差时所设计波形的性能,设 θkΘk=[θk,l,θk,u],k=1,2,,K(若 θk,l=θk,u, 则说明角度信息是准确的)。为便于后续优化波形,将角度区间 Θk离散化为 Nk个点,其中第 m个离散角度记为 θk,m(m=1,2,,Nk)。综上,将综合射频系统的鲁棒波形设计问题建模为

    minSmaxθk,mΩk||aH(θk,m)SdT(θk)||22,s.t.sHnsn=eT/NT,PAPR(sn)ρ,n=1,2,,NT (8)

    其中, Ωk=θk,1,θk,2,,θk,Nk表示离散角度点集合, k=1,2,,K

    由于峰均比约束是非凸约束,故而式(8)是一个非凸优化问题。接下来利用交替方向乘子法(Alternating Direction Method of Multipliers, ADMM)[25]来求解该问题。首先,引入辅助变量 z,将式(8)改写为

    minS,zz,s.t.||aHk,mSdTk||22z,θk,mΩk,k=1,2,,K,sHnsn=eT/NT,PAPR(sn)ρ,n=1,2,,NT (9)

    其中, ak,m。然后,引入变量 {{\boldsymbol{y}}_{k,m}}(k = 1,2, \cdots\;,K,m = 1,2, \cdots \;,{N_k}),利用变量分裂法将式(9)转化为

    \begin{split} & \mathop {\min }\limits_{{\boldsymbol{S}},z,\{ {{\boldsymbol{y}}_{k,m}}\} } z, \\ & \;\;\;{\text{s}}{\text{.t}}{\text{.}}\;{{\boldsymbol{y}}_{k,m}} = {{\boldsymbol{S}}^{\text{H}}}{{\boldsymbol{a}}_{k,m}} - {\boldsymbol{d}}_k^{\text{*}}, \\ & \;\;\;\;\;\;\;||{{\boldsymbol{y}}_{k,m}}||_2^2 \le z,\forall {\theta _{k,m}} \in {\varOmega _k},k = 1,2, \cdots \;,K, \\ & \;\;\;\;\;\;\;{\boldsymbol{s}}_n^{\text{H}}{{\boldsymbol{s}}_n} = {e_{\text{T}}}/{N_{\text{T}}},{\text{PAPR(}}{{\boldsymbol{s}}_n}{\text{)}} \le \rho {\text{,}}n = 1,{\text{2,}} \cdots \;,{N_{\text{T}}} \end{split} (10)

    式(10)的增广拉格朗日函数为

    \begin{split} & {L_\mu }({\boldsymbol{S}},\{ {{\boldsymbol{y}}_{k,m}}\} ,z,\{ {\gamma _{k,m}}\} )\\ & \quad = z + \frac{\mu }{2}\sum\limits_{k = 1}^K \sum\limits_{m = 1}^{{N_k}} (||{{\boldsymbol{y}}_{k,m}} - {{\boldsymbol{S}}^{\text{H}}}{{\boldsymbol{a}}_{k,m}} + {\boldsymbol{d}}_k^{\text{*}} \\ & \qquad + {\gamma _{k,m}}||_2^2 - ||{\gamma _{k,m}}||_2^2) \end{split} (11)

    其中, \mu 为惩罚因子, {\gamma _{k,m}} 为等式约束 {{\boldsymbol{y}}_{k,m}} = {{\boldsymbol{S}}^{\text{H}}}{{\boldsymbol{a}}_{k,m}} - {\boldsymbol{d}}_k^{\text{*}}的拉格朗日乘子。在ADMM算法的第 l + 1 次迭代中,需要依次更新下列变量

    {{\boldsymbol{S}}^{(l + 1)}} = \arg \mathop {\min }\limits_{\boldsymbol{S}} {L_\mu }({\boldsymbol{S}},\{ {\boldsymbol{y}}_{k,m}^{(l)}\} ,{z^{(l)}},\{ \gamma _{k,m}^{(l)}\} ) (12)
    \begin{split} \{ {\boldsymbol{y}}_{k,m}^{(l + 1)},{z^{(l + 1)}}\} =& \arg \mathop {\min }\limits_{\{ {{\boldsymbol{y}}_{k,m}}\} ,z} {L_\mu }({{\boldsymbol{S}}^{(l + 1)}},\\ & \{ {{\boldsymbol{y}}_{k,m}}\} ,z,\{ \gamma _{k,m}^{(l)}\} ) \end{split} (13)
    \gamma _{k,m}^{(l + 1)} = \gamma _{k,m}^{(l)} + {\boldsymbol{y}}_{k,m}^{(l + 1)} - {({{\boldsymbol{S}}^{(l + 1)}})^{\text{H}}}{{\boldsymbol{a}}_{k,m}} + {\boldsymbol{d}}_k^{\text{*}} (14)

    下面求解优化问题式(12)和式(13),为使得求解过程清晰简洁,在不引起歧义的情况下,将上标省略。

    {\boldsymbol{A}} = \sum\limits_{k = 1}^K {\sum\limits_{m = 1}^{{N_k}} {{{\boldsymbol{a}}_{k,m}}} } {\boldsymbol{a}}_{k,m}^{\text{H}},\;\;{\boldsymbol{B}} = \sum\limits_{k = 1}^K {\sum\limits_{m = 1}^{{N_k}} {{{\boldsymbol{a}}_{k,m}}} } {\boldsymbol{b}}_{k,m}^{\text{H}} (15)

    其中, {{\boldsymbol{b}}_{k,m}} = {{\boldsymbol{y}}_{k,m}} + {\boldsymbol{d}}_k^{\text{*}} + {\gamma _{k,m}} ,则优化问题式(12)可以表示为

    \begin{gathered} \mathop {\min }\limits_{\boldsymbol{S}} {\text{tr(}}{{\boldsymbol{S}}^{\text{H}}}{\boldsymbol{AS}}{\text{)}} - {\text{2Re}}\left[ {{\text{tr(}}{{\boldsymbol{S}}^{\text{H}}}{\boldsymbol{B}})} \right], \\ \;\;\;{\text{s}}{\text{.t}}{\text{.}}\;{\boldsymbol{s}}_n^{\text{H}}{{\boldsymbol{s}}_n} = {e_{\text{T}}}/{N_{\text{T}}},{\text{PAPR(}}{{\boldsymbol{s}}_n}{\text{)}} \le \rho {\text{,}}n = 1,{\text{2,}} \cdots\;,{N_{\text{T}}} \\ \end{gathered} (16)

    这是一个非凸的2次优化问题,可以利用Majorization-Minimization(MM)方法进行求解[26-28]。MM方法是一种迭代算法,在每次迭代中,构造原优化问题目标函数的上界函数,通过求该上界函数的最小值,来逼近原优化问题的最小值,从而实现对原问题的求解。下面对该方法进行推导,首先令 {{\boldsymbol{A}}_{{\text{neg}}}} = {\boldsymbol{A}} - {\lambda _{\max }}({\boldsymbol{A}}){\boldsymbol{I}} ( {\lambda _{\max }}({\boldsymbol{A}}) 表示矩阵 {\boldsymbol{A}} 的最大特征值),易得 {{\boldsymbol{A}}_{{\text{neg}}}} 是负定的,故式(17)恒成立

    {\text{tr}}\left[ {{{{\text{(}}{\boldsymbol{S}} - {{\boldsymbol{S}}^{(q)}}{\text{)}}}^{\text{H}}}{{\boldsymbol{A}}_{{\text{neg}}}}{\text{(}}{\boldsymbol{S}} - {{\boldsymbol{S}}^{(q)}}{\text{)}}} \right] \le 0 (17)

    其中, {{\boldsymbol{S}}^{(q)}} 表示在MM算法中第 q 次(内层)循环时的发射波形矩阵。将式(17)展开,可以得到

    \begin{split} {\text{tr}}({{\boldsymbol{S}}^{\text{H}}}{\boldsymbol{AS}}) \le & {\text{2Re}}\left[ {{\text{tr(}}{{\boldsymbol{S}}^{\text{H}}}{{\boldsymbol{A}}_{{\text{neg}}}}{{\boldsymbol{S}}^{(q)}})} \right] + 2{\lambda _{\max }}({\boldsymbol{A}}){e_{\text{T}}} \\ & - {\text{tr}}\left[ {{{{\text{(}}{{\boldsymbol{S}}^{(q)}}{\text{)}}}^{\text{H}}}{\boldsymbol{A}}{\text{(}}{{\boldsymbol{S}}^{(q)}}{\text{)}}} \right]\\[-15pt] \end{split} (18)

    不难看出,式(18)中不等号右边的函数可以被视为 {\text{tr}}({{\boldsymbol{S}}^{\text{H}}}{\boldsymbol{AS}}) 的上界函数。因此,式(16)上界函数的最小化问题(忽略常数项)为

    \begin{gathered} \mathop {\min }\limits_{\boldsymbol{S}} \; - {\text{2Re}}\left[ {{\text{tr(}}{{\boldsymbol{S}}^{\text{H}}}{{\tilde {\boldsymbol{B}}}})} \right], \\ \;\;\;{\text{s}}{\text{.t}}{\text{.}}\;{\boldsymbol{s}}_n^{\text{H}}{{\boldsymbol{s}}_n} = {e_{\text{T}}}/{N_{\text{T}}},{\text{PAPR(}}{{\boldsymbol{s}}_n}{\text{)}} \le \rho {\text{,}}n = 1,{\text{2,}} \cdots ,{N_{\text{T}}} \\ \end{gathered} (19)

    其中, {{\tilde {\boldsymbol{B}}}} = {\boldsymbol{B}} - {{\boldsymbol{A}}_{{\text{neg}}}}{{\boldsymbol{S}}^{(q)}} = {[{{{\tilde {\boldsymbol{b}}}}_1},{{{\tilde {\boldsymbol{b}}}}_2}, \cdots \;,{{\boldsymbol{\tilde b}}_{{N_{\text{T}}}}}]^{\text{T}}}。令 {{\tilde {\boldsymbol{b}}}}_n^{\text{T}}为矩阵 {{\tilde {\boldsymbol{B}}}}的第 n ( n = 1,2, \cdots \;,{N_{\text{T}}})行向量,可以得到

    {\text{tr(}}{{\boldsymbol{S}}^{\text{H}}}{{\tilde {\boldsymbol{B}}}}) = \sum\limits_{n = 1}^{{N_{\text{T}}}} {{\text{tr(}}{\boldsymbol{s}}_n^*{{\tilde {\boldsymbol{b}}}}_n^{\text{T}})} = \sum\limits_{n = 1}^{{N_{\text{T}}}} {{{{\text{(}}{{\tilde {\boldsymbol{b}}}}_n^{\text{H}}{{\boldsymbol{s}}_n})}^*}} (20)

    根据式(20),将优化问题式(19)拆解为NT个独立的子问题。其中,第n(n = 1,2, ···, NT)个子问题可以表示为

    \begin{split} & \mathop {\max }\limits_{{{\boldsymbol{s}}_n}} {\text{Re(}}{{\tilde {\boldsymbol{b}}}}_n^{\text{H}}{{\boldsymbol{s}}_n}), \\ &\;\;\;{\text{s}}{\text{.t}}{\text{.}}\;{\boldsymbol{s}}_n^{\text{H}}{{\boldsymbol{s}}_n} = {e_{\text{T}}}/{N_{\text{T}}},{\text{PAPR(}}{{\boldsymbol{s}}_n}{\text{)}} \le \rho \end{split} (21)

    优化问题式(21)可以采用文献[29]中的算法2进行求解。

    {{\boldsymbol{u}}_{k,m}} = {{\boldsymbol{S}}^{\text{H}}}{{\boldsymbol{a}}_{k,m}} - {\boldsymbol{d}}_k^{\text{*}} - {\gamma _{k,m}} (22)

    优化问题式(13)可以表示为

    \begin{split} & \mathop {\min }\limits_{\{ {{\boldsymbol{y}}_{k,m}}\} ,z} z + \frac{\mu }{2}\sum\limits_{k = 1}^K {\sum\limits_{m = 1}^{{N_k}} {(||{{\boldsymbol{y}}_{k,m}} - {{\boldsymbol{u}}_{k,m}}||_2^2)} } , \\ & \;\;\;{\text{s}}{\text{.t}}{\text{.}}\;||{{\boldsymbol{y}}_{k,m}}||_2^2 \le z,k = 1,2, \cdots \;,K,m = 1,2, \cdots \;,{N_k} \end{split} (23)

    这是一个2阶锥规划(Second-Order Cone Programming, SOCP)问题,可以在多项式时间内求得最优解(例如,使用内点法)[30]。但是,利用内点法求解式(23)时,计算复杂度高,运算效率低。下面给出该问题的快速解法。首先,容易推出式(23)中 {{\boldsymbol{y}}_{k,m}} 的最优解为

    {{\boldsymbol{y}}_{k,m}} = \left\{ \begin{array}{lllllllll} {{\boldsymbol{u}}_{k,m}},& ||{{\boldsymbol{u}}_{k,m}}||_2^2 \le z \\ {z^{1/2}} \cdot {{\boldsymbol{u}}_{k,m}}/||{{\boldsymbol{u}}_{k,m}}|{|_2},& ||{{\boldsymbol{u}}_{k,m}}||_2^2 > z \end{array} \right. (24)

    然后,将式(24)代入式(23)中,则式(23)的目标函数可以化简为

    f(z) = z + \frac{\mu }{2}\sum\limits_{k = 1}^K {\sum\limits_{m = 1}^{{N_k}} {{\omega _{k,m}}{{(\sqrt z - ||{{\boldsymbol{u}}_{k,m}}|{|_2})}^2}} } (25)

    其中

    {\omega _{k,m}} = \left\{ \begin{gathered} 0,\;\;\; ||{{\boldsymbol{u}}_{k,m}}||_2^2 \le z \\ 1,\;\;\; ||{{\boldsymbol{u}}_{k,m}}||_2^2 > z \\ \end{gathered} \right. (26)

    可以看出 f(z) (0,\infty ) 上是分段凸函数。将集合 \{ ||{{\boldsymbol{u}}_{k,m}}||_2^2,k = 1,2, \cdots ,K,m = 1,2, \cdots,{N_k}\}按照升序进行排列,记为 \{ {a_k}\} _{k = 1}^{\tilde N} ( \tilde N = \sum\nolimits_{k = 1}^K {{N_k}} )。

    接着,将 f(z) 的定义域划分为 \tilde N + 1 个小的区间,即 (0,{a_1}) \cup [{a_2},{a_3}) \cup \cdots \cup [{a_{\tilde N}},\infty )。要求 f(z) 的最小值,可以先分别求 f(z) 在上述 \tilde N + 1 个区间上的最小值,然后比较它们的大小,最终得出 f(z) 的最小值。当 z \in [{a_{\tilde N}},\infty ) 时,由于 ||{{\boldsymbol{u}}_{k,m}}||_2^2 \le {a_{\tilde N}} 恒成立,此时,最优解为 z = {a_{\tilde N}} ,且 f(z) 的最小值为 {a_{\tilde N}} 。当 z 位于第 k 个区间上,即 z \in [{a_{k - 1}},{a_k}) 时( {a_0} = 0,k = 1, 2, \cdots,\tilde N), f(z) 可以化简为

    f(z) = z + \frac{\mu }{2}\sum\limits_{k' = k}^K {{{(\sqrt z - \sqrt {{a_{k'}}} )}^2}} (27)

    f(z) 求导,得

    f'(z) = {c_0} - {b_0}{z^{ - 1/2}} (28)

    其中

    \left. \begin{aligned} & {{c_0} = 1 + \mu /2(\tilde N - k + 1)} \\ & {{b_0} = \frac{\mu }{2} \cdot \sum\limits_{k' = k}^K {\sqrt {{a_k}} } } \end{aligned} \right\} (29)

    可以看出 f'(z) 是关于 z 的单调递增函数。因此,在区间 [{a_{k - 1}},{a_k}) 上,如果 f'({a_{k - 1}}) > 0 ,则 f(z) 的最优解为 z = {a_{k - 1}} ;如果 f'({a_k}) < 0 ,则 f(z) 的最优解为 z = {a_k} ;其他情况下,最优解为 z = b_0^2/c_0^2

    对上述算法进行归纳总结,形成算法1,其停止迭代的条件为

    算法1 基于ADMM算法的综合射频系统鲁棒波形设计算法
     输入:阵列天线数 {N_{\text{T}}} ,期望信号以及角度误差范围 \{ {{\boldsymbol{d}}_k},{\varOmega _k}\} _{k = 1}^K ,发射总能量 {e_{\text{T}}} ,峰均比 \rho ,惩罚因子 \mu
     输出:发射阵列波形 {\boldsymbol{S}}
     初始化 {\boldsymbol{A}} = \sum\nolimits_{k = 1}^K {\sum\nolimits_{m = 1}^{{N_k}} {{{\boldsymbol{a}}_{k,m}}{\boldsymbol{a}}_{k,m}^{\text{H}}} } ,{{\boldsymbol{A}}_{{\text{neg}}}} = {\boldsymbol{A}} - {\lambda _{\max }}({\boldsymbol{A}}){\boldsymbol{I}},l = 0,{{\boldsymbol{S}}^{(l)}},\{ {\boldsymbol{y}}_{k,m}^{(l)},\gamma _{k,m}^{(l)}\}
     重复:
        q = 0,{{\boldsymbol{S}}^{(l,q)}} = {{\boldsymbol{S}}^{(l)}},{\boldsymbol{b}}_{k,m}^{(l)} = {\boldsymbol{y}}_{k,m}^{(l)} + {\boldsymbol{d}}_k^{\text{*}} + \gamma _{k,m}^{(l)},{{\boldsymbol{B}}^{(l)}} = \sum\nolimits_{k = 1}^K {\sum\nolimits_{m = 1}^{{N_k}} {{{\boldsymbol{a}}_{k,m}}{{({\boldsymbol{b}}_{k,m}^{(l)})}^{\text{H}}}} }
       重复:
          {{{\tilde {\boldsymbol B}}}^{(l,q)}} = {{\boldsymbol{B}}^{(l)}} - {{\boldsymbol{A}}_{{\text{neg}}}}{{\boldsymbol{S}}^{(l,q)}} ,利用文献[29]中的算法2求解优化问题式(21), q = q + 1
       直到收敛
        {{\boldsymbol{S}}^{(l + 1)}} = {{\boldsymbol{S}}^{(l,q)}}
        {\boldsymbol{u}}_{k,m}^{(l + 1)} = {({{\boldsymbol{S}}^{(l + 1)}})^{\text{H}}}{{\boldsymbol{a}}_{k,m}} - {\boldsymbol{d}}_k^{\text{*}} - \gamma _{k,m}^{(l)} ,计算 f(z) 的最小值以及最小值所对应的点 {z^{(l + 1)}} , {\boldsymbol{y}}_{k,m}^{(l + 1)} = {\boldsymbol{u}}_{k,m}^{(l + 1)} \cdot \min \left( {\sqrt {{z^{(l + 1)}}} /||{\boldsymbol{u}}_{k,m}^{(l + 1)}|{|_2},1} \right)
        \gamma _{k,m}^{(l + 1)} = \gamma _{k,m}^{(l)} + {\boldsymbol{y}}_{k,m}^{(l + 1)} - {({{\boldsymbol{S}}^{(l + 1)}})^{\text{H}}}{{\boldsymbol{a}}_{k,m}} - {\boldsymbol{d}}_k^{\text{*}},l = l + 1
     直到收敛
    下载: 导出CSV 
    | 显示表格
    \left| {G({{\boldsymbol{S}}^{(l + 1)}}) - G({{\boldsymbol{S}}^{(l)}})} \right| < \vartheta (30)

    其中, G({\boldsymbol{S}}) = \mathop {\max }\limits_{k = 1,2, \cdot \cdot \cdot \;,\tilde N} ||{\boldsymbol{a}}_k^{\text{H}}{\boldsymbol{S}} - {\boldsymbol{d}}_k^{\text{T}}||_2^2 \vartheta 为控制迭代停止的阈值。

    本节中,综合射频系统所采用的MIMO阵列发射天线数为 {N_{\text{T}}} = 16 ,阵元等间距,且间距为半波长。每个天线发射信号的码长为 L = 256 ,总的发射能量为 {e_{\text{T}}} = 2 。期望生成3组信号,其中在 {\theta _1}{\text{ = 4}}{{\text{0}}^ \circ } 附近生成能量为 {e_{\text{i}}} = 0\;{\text{dB}} 的高斯白噪声信号用于干扰;在 {\theta _2}{\text{ = }}{{\text{0}}^ \circ } 附近生成能量为 {e_{\text{c}}} = 0\;{\text{dB}} 的正交相移键控(Quadrature Phase Shift Keying, QPSK)信号用于通信;在 {\theta _3}{\text{ = }} - {\text{3}}{{\text{0}}^ \circ } 附近生成能量为 {e_{\text{r}}} = 5\;{\text{dB}} 的线性调频信号用于雷达探测。期望信号的角度误差范围分别为 {\varTheta _1} = [{38^ \circ },{42^ \circ }], {\varTheta _2} = [ - {2^ \circ },{2^ \circ }]以及 {\varTheta _3} = [ - {32^ \circ }, - {28^ \circ }],为便于设计波形,将角度误差范围间隔 {0.4^ \circ } 进行离散化,即 {N_1} = {N_2} = {N_3} = 11 。将通信方向权重系数设为 {\omega _{\text{c}}} = 5 ,其余方向权重系数均为1,峰均比为 \rho = 1 ,停止迭代阈值为 \vartheta = {10^{ - 8}}

    图2为所提算法进行波形设计时的目标函数值收敛曲线。从图2结果可以看出,随着迭代次数的增加,目标函数值逐渐减小,并趋于收敛。

    图 2  算法收敛曲线

    图3所示为采用所提算法设计波形时,在各个角度上合成的信号与期望信号之间的误差值。可以看出,如果未考虑角度估计误差(即利用式(7)来设计波形),所合成信号与期望信号之间的误差很大,这将会严重影响综合射频系统的工作性能。本文所设计的鲁棒发射波形,即使存在角度估计误差,在目标探测和电子干扰方向上,所合成信号与期望信号之间的误差数量级约在10–2,明显优于未考虑角度估计误差时所设计波形的性能。而在通信方向上,由于施加了权重系数,所合成的信号与期望通信信号之间的误差更小,几乎可以忽略。

    图 3  算法在各个角度上所合成的信号与期望信号的误差

    图4分析了算法所合成的信号性能,其中图4(a)所示为 {\varTheta _3} = [ - {32^ \circ }, - {28^ \circ }] 区间内各个角度所合成雷达信号的自相关函数,可以看出在不同角度上所合成的波形具有相似的自相关函数,副瓣电平高度基本相同(约为–13.5 dB);图4(b) {\varTheta _2} = [ - {2^ \circ },{2^ \circ }] 区间内所合成的QPSK信号星座图,星座图较为理想,表明合成误差较小;图4(c)采用分位数(Quantile-Quantile, Q-Q)图对 {\varTheta _1} = [{38^ \circ },{42^ \circ }] 区间内合成的干扰信号进行正态性检验,可以看出所合成的干扰信号是遵从正态分布的,表明所合成的干扰信号具有类似于高斯噪声的特性。

    图 4  合成信号性能检验

    下面进一步分析算法所合成信号的性能。假设在 - {31^ \circ } 存在两个目标,与综合射频系统的距离分别为32 km和75 km,目标回波的信噪比分别为5 dB和–5.5 dB。将目标回波进行脉冲压缩处理,所得目标距离像如图5所示。可以看出,在32 km和75 km处均发现了目标。

    图 5  所合成的雷达信号脉压性能

    图6为在各角度所合成的通信信号误码率随信噪比的变化曲线,其中重复实验2 000次。相比于理想的QPSK信号,误码率没有明显上升。

    图 6  所合成的通信信号误码率

    下面验证所合成的干扰信号的干扰效果,使用在 {42^ \circ } 方向所合成的干扰信号分别对QPSK信号和雷达回波信号进行干扰,其中假设对方雷达使用LFM信号作为雷达发射波形。图7(a)为QPSK信号误码率随信干比的变化曲线(信噪比为5 dB),其中重复实验2 000次。当信干比为0 dB时,误码率由干扰前的1.6 × 10–3增大至为0.328 3;图7(b)分析了干扰前后雷达目标检测概率随着信噪比的变化曲线(信干比为10 dB,虚警概率为10–6),可以看出,综合射频系统所释放的干扰信号会显著降低对方雷达的目标检测概率。

    图 7  合成干扰信号的效果

    下面分析峰均比约束对于算法性能的影响。图8为不同峰均比条件下算法的收敛曲线,其中除改变峰均比外,仿真条件与图2相同。图9所示为算法收敛时所设计波形在各个角度所合成的信号与期望信号之间的误差曲线。可以得出,当峰均比增大时,目标函数的收敛值变小,这是因为峰均比的增大会使得波形设计的自由度变大,有更多发射波形可供选择,从而会减小所合成信号与期望信号之间的匹配误差。

    图 8  不同峰均比下的算法迭代曲线
    图 9  不同峰均比下算法在各个角度上所合成的信号与期望信号的误差

    图10分析不同峰均比条件下,算法所设计波形在不同角度上合成的雷达信号的自相关函数的峰值旁瓣比(Peak Side Lobe Ratio, PSLR)。很明显,当增加波形的峰均比时,波形匹配误差变小,峰值旁瓣比的取值趋于更加稳定。最后分析在不同峰均比条件下,在 {0.8^ \circ } 方向上所合成的通信信号误码率随信噪比变化的曲线,如图11所示,其中进行2 000次独立实验。可以看出提高峰均比后,通信信号误码率略有下降。

    图 10  不同峰均比下雷达信号峰值旁瓣比
    图 11  不同峰均比下通信信号误码率

    针对存在角度误差时基于MIMO阵列的综合射频系统性能下降问题,本文研究了鲁棒低峰均比波形设计问题,提出一种基于ADMM和MM方法的迭代优化算法。仿真实验表明,算法所设计的波形能够分区域同时进行雷达探测、数据通信以及电子干扰。因此,系统可以在占用少量频谱资源的条件下同时实现多种功能,提高了频谱的利用效率。与此同时,当存在角度误差时,算法在误差区间内所合成的信号与期望信号之间的匹配误差较小,有助于提升综合射频系统在实际环境中的性能。

  • 图  1  基于MIMO阵列的综合射频系统

    图  2  算法收敛曲线

    图  3  算法在各个角度上所合成的信号与期望信号的误差

    图  4  合成信号性能检验

    图  5  所合成的雷达信号脉压性能

    图  6  所合成的通信信号误码率

    图  7  合成干扰信号的效果

    图  8  不同峰均比下的算法迭代曲线

    图  9  不同峰均比下算法在各个角度上所合成的信号与期望信号的误差

    图  10  不同峰均比下雷达信号峰值旁瓣比

    图  11  不同峰均比下通信信号误码率

    算法1 基于ADMM算法的综合射频系统鲁棒波形设计算法
     输入:阵列天线数 {N_{\text{T}}} ,期望信号以及角度误差范围 \{ {{\boldsymbol{d}}_k},{\varOmega _k}\} _{k = 1}^K ,发射总能量 {e_{\text{T}}} ,峰均比 \rho ,惩罚因子 \mu
     输出:发射阵列波形 {\boldsymbol{S}}
     初始化 {\boldsymbol{A}} = \sum\nolimits_{k = 1}^K {\sum\nolimits_{m = 1}^{{N_k}} {{{\boldsymbol{a}}_{k,m}}{\boldsymbol{a}}_{k,m}^{\text{H}}} } ,{{\boldsymbol{A}}_{{\text{neg}}}} = {\boldsymbol{A}} - {\lambda _{\max }}({\boldsymbol{A}}){\boldsymbol{I}},l = 0,{{\boldsymbol{S}}^{(l)}},\{ {\boldsymbol{y}}_{k,m}^{(l)},\gamma _{k,m}^{(l)}\}
     重复:
        q = 0,{{\boldsymbol{S}}^{(l,q)}} = {{\boldsymbol{S}}^{(l)}},{\boldsymbol{b}}_{k,m}^{(l)} = {\boldsymbol{y}}_{k,m}^{(l)} + {\boldsymbol{d}}_k^{\text{*}} + \gamma _{k,m}^{(l)},{{\boldsymbol{B}}^{(l)}} = \sum\nolimits_{k = 1}^K {\sum\nolimits_{m = 1}^{{N_k}} {{{\boldsymbol{a}}_{k,m}}{{({\boldsymbol{b}}_{k,m}^{(l)})}^{\text{H}}}} }
       重复:
          {{{\tilde {\boldsymbol B}}}^{(l,q)}} = {{\boldsymbol{B}}^{(l)}} - {{\boldsymbol{A}}_{{\text{neg}}}}{{\boldsymbol{S}}^{(l,q)}} ,利用文献[29]中的算法2求解优化问题式(21), q = q + 1
       直到收敛
        {{\boldsymbol{S}}^{(l + 1)}} = {{\boldsymbol{S}}^{(l,q)}}
        {\boldsymbol{u}}_{k,m}^{(l + 1)} = {({{\boldsymbol{S}}^{(l + 1)}})^{\text{H}}}{{\boldsymbol{a}}_{k,m}} - {\boldsymbol{d}}_k^{\text{*}} - \gamma _{k,m}^{(l)} ,计算 f(z) 的最小值以及最小值所对应的点 {z^{(l + 1)}} , {\boldsymbol{y}}_{k,m}^{(l + 1)} = {\boldsymbol{u}}_{k,m}^{(l + 1)} \cdot \min \left( {\sqrt {{z^{(l + 1)}}} /||{\boldsymbol{u}}_{k,m}^{(l + 1)}|{|_2},1} \right)
        \gamma _{k,m}^{(l + 1)} = \gamma _{k,m}^{(l)} + {\boldsymbol{y}}_{k,m}^{(l + 1)} - {({{\boldsymbol{S}}^{(l + 1)}})^{\text{H}}}{{\boldsymbol{a}}_{k,m}} - {\boldsymbol{d}}_k^{\text{*}},l = l + 1
     直到收敛
    下载: 导出CSV
  • [1] 肖博, 霍凯, 刘永祥. 雷达通信一体化研究现状与发展趋势[J]. 电子与信息学报, 2019, 41(3): 739–750. doi: 10.11999/JEIT180515

    XIAO Bo, HUO Kai, and LIU Yongxiang. Development and prospect of radar and communication integration[J]. Journal of Electronics&Information Technology, 2019, 41(3): 739–750. doi: 10.11999/JEIT180515
    [2] ZHANG J A, LIU Fan, MASOUROS C, et al. An overview of signal processing techniques for joint communication and radar sensing[J]. IEEE Journal of Selected Topics in Signal Processing, 2021, 15(6): 1295–1315. doi: 10.1109/JSTSP.2021.3113120
    [3] HUGHES P K and CHOE J Y. Overview of advanced multifunction RF system (AMRFS)[C]. 2000 IEEE International Conference on Phased Array Systems and Technology, Dana Point, USA, 2000.
    [4] SCOTT R. ONR plans EMC2 technology maturation effort[J]. Jane's International Defense Review:IHS Jane's International Defense Review, 2016, 49: 10.
    [5] 唐波, 汤俊, 胡元奎. 基于MIMO阵列的综合射频系统技术研究[J]. 信息对抗技术, 2022, 1(1): 62–72. doi: 10.12399/j.issn.2097-163x.2022.01.006

    TANG Bo, TANG Jun, and HU Yuankui. Multifunction radio frequency systems based on MIMO array[J]. Information Countermeasure Technology, 2022, 1(1): 62–72. doi: 10.12399/j.issn.2097-163x.2022.01.006
    [6] TANG Bo and STOICA P. MIMO multifunction RF systems: Detection performance and waveform design[J]. IEEE Transactions on Signal Processing, 2022, 70: 4381–4394. doi: 10.1109/TSP.2022.3202315
    [7] MA Dingyou, SHLEZINGER N, HUANG Tianyao, et al. Joint radar-communication strategies for autonomous vehicles: Combining two key automotive technologies[J]. IEEE Signal Processing Magazine, 2020, 37(4): 85–97. doi: 10.1109/MSP.2020.2983832
    [8] AXELSSON S R J. Random noise radar/sodar with ultrawideband waveforms[J]. IEEE Transactions on Geoscience and Remote Sensing, 2007, 45(5): 1099–1114. doi: 10.1109/TGRS.2007.893742
    [9] FIDEN W H and CZUBIAK D W. Radar-compatible data link system (U)[P]. USA Patent, US7298313B1, 2007.
    [10] ELLINGER J, ZHANG Zhiping, WU Zhiqiang, et al. Dual-use multicarrier waveform for radar detection and communication[J]. IEEE Transactions on Aerospace and Electronic Systems, 2018, 54(3): 1265–1278. doi: 10.1109/TAES.2017.2780578
    [11] 李晓柏, 杨瑞娟, 程伟. 基于频率调制的多载波Chirp信号雷达通信一体化研究[J]. 电子与信息学报, 2013, 35(2): 406–412. doi: 10.3724/SP.J.1146.2012.00567

    LI Xiaobai, YANG Ruijuan, and CHENG Wei. Integrated radar and communication based on multicarrier frequency modulation chirp signal[J]. Journal of Electronics&Information Technology, 2013, 35(2): 406–412. doi: 10.3724/SP.J.1146.2012.00567
    [12] CHEN Xingbo, WANG Xiaomo, ZHANG Zhao, et al. A communication with LFM carrier modulated by MSK scheme[C]. IEEE ICCP2012, Nanjing, China, 2012: 461–463.
    [13] 刘国岁, 顾红, 苏卫民. 随机信号雷达[M]. 北京: 国防工业出版社, 2005: 45–52.

    LIU Guosui, GU Hong, and SU Weimin. Random Signal Radar[M]. Beijing: National Defense Industry Press, 2005: 45–52.
    [14] LI Jian and STOICA P. MIMO Radar Signal Processing[M]. Hoboken: J. Wiley & Sons, 2009: 1–20.
    [15] 何团, 唐波, 张玉. 基于加权SPICE的MIMO-STAP稀疏恢复算法[J]. 信号处理, 2019, 35(8): 1417–1424. doi: 10.16798/j.issn.1003-0530.2019.08.017

    HE Tuan, TANG Bo, and ZHANG Yu. Sparse recovery algorithm of MIMO-STAP based on weighted SPICE[J]. Journal of Signal Processing, 2019, 35(8): 1417–1424. doi: 10.16798/j.issn.1003-0530.2019.08.017
    [16] WANG Xuyang, TANG Bo, and ZHANG Ming. Optimisation of practically constrained waveforms for Rician target detection with multiple-input-multiple-output radar[J]. IET Radar,Sonar&Navigation, 2022, 16(7): 1116–1130. doi: 10.1049/rsn2.12247
    [17] 刘冰凡, 陈伯孝. 基于OFDM-LFM信号的MIMO雷达通信一体化信号共享设计研究[J]. 电子与信息学报, 2019, 41(4): 801–808. doi: 10.11999/JEIT180547

    LIU Bingfan and CHEN Baixiao. Integration of MIMO radar and communication with OFDM-LFM signals[J]. Journal of Electronics&Information Technology, 2019, 41(4): 801–808. doi: 10.11999/JEIT180547
    [18] TANG Bo, WANG Hai, QIN Lilong, et al. Waveform design for dual-function MIMO radar-communication systems[C]. The IEEE 11th Sensor Array and Multichannel Signal Processing Workshop (SAM), Hangzhou, China, 2020.
    [19] LIU Fan, ZHOU Longfei, MASOUROS C, et al. Toward dual-functional radar-communication systems: Optimal waveform design[J]. IEEE Transactions on Signal Processing, 2018, 66(16): 4264–4279. doi: 10.1109/TSP.2018.2847648
    [20] SHI Shengnan, WANG Zhaoyi, HE Zishu, et al. Constrained waveform design for dual-functional MIMO radar-communication system[J]. Signal Processing, 2020, 171: 107530. doi: 10.1016/j.sigpro.2020.107530
    [21] FAN Linhui, TANG Bo, HUANG Zhongrui, et al. Designing constant-envelope transmissions for secret communications in MISO wiretap channels[J]. IEEE Access, 2019, 7: 17791–17797. doi: 10.1109/ACCESS.2019.2896406
    [22] 吴文俊, 唐波, 汤俊, 等. 杂波环境中雷达通信一体化系统波形设计算法研究[J]. 雷达学报, 2022, 11(4): 570–580. doi: 10.12000/JR22105

    WU Wenjun, TANG Bo, TANG Jun, et al. Waveform design for dual-function radar-communication systems in clutter[J]. Journal of Radars, 2022, 11(4): 570–580. doi: 10.12000/JR22105
    [23] TSINOS C G, ARORA A, CHATZINOTAS S, et al. Joint transmit waveform and receive filter design for dual-function radar-communication systems[J]. IEEE Journal of Selected Topics in Signal Processing, 2021, 15(6): 1378–1392. doi: 10.1109/JSTSP.2021.3112295
    [24] TANG Bo, LIU Jun, WANG Hai, et al. Constrained radar waveform design for range profiling[J]. IEEE Transactions on Signal Processing, 2021, 69: 1924–1937. doi: 10.1109/TSP.2021.3065830
    [25] BOYD S, PARIKH N, CHU E, et al. Distributed optimization and statistical learning via the alternating direction method of multipliers[J]. Foundations and Trends ® in Machine Learning , 2010, 3(1): 1–122. doi: 10.1561/2200000016
    [26] HUNTER D R and LANGE K. A tutorial on MM algorithms[J]. The American Statistician, 2004, 58(1): 30–37. doi: 10.1198/0003130042836
    [27] TANG Bo, TUCK J, and STOICA P. Polyphase waveform design for MIMO radar space time adaptive processing[J]. IEEE Transactions on Signal Processing, 2020, 68: 2143–2154. doi: 10.1109/TSP.2020.2983833
    [28] TANG Bo and STOICA P. Information-theoretic waveform design for MIMO radar detection in range-spread clutter[J]. Signal Processing, 2021, 182: 107961. doi: 10.1016/j.sigpro.2020.107961
    [29] TROPP J A, DHILLON I S, HEATH R W, et al. Designing structured tight frames via an alternating projection method[J]. IEEE Transactions on Information Theory, 2005, 51(1): 188–209. doi: 10.1109/TIT.2004.839492
    [30] BOYD S and VANDENBERGHE L. Convex Optimization[M]. Cambridge: Cambridge University Press, 2004: 21–58.
  • 期刊类型引用(1)

    1. 唐波,吴文俊. 基于MIMO阵列的一体化射频技术研究进展及展望. 信息对抗技术. 2024(01): 3-17 . 百度学术

    其他类型引用(0)

  • 加载中
图(11) / 表(1)
计量
  • 文章访问数:  913
  • HTML全文浏览量:  451
  • PDF下载量:  205
  • 被引次数: 1
出版历程
  • 收稿日期:  2022-07-21
  • 修回日期:  2022-09-17
  • 网络出版日期:  2022-09-22
  • 刊出日期:  2023-11-28

目录

/

返回文章
返回