1 引言
近年来,随着民用通信业务的广泛而深入的发展,雷达对于高分辨率的迫切追求,通信和雷达等业务都在积极扩张自身占据的频带宽度,这就导致射频频谱环境中不同业务频谱交叠现象日益突显,从而使得利用频谱共享技术解决频谱拥堵问题逐渐成为雷达界和通信界的研究热门。针对频谱环境尤其是在低频频段射频频谱日益拥堵的现象,充分分析雷达在频谱资源利用方面面临的挑战,本文从雷达波形设计的角度,实现雷达在某些频段与通信或其他一些设备产生频谱交叠的情况下,共享频谱资源,提高频谱利用率。
频谱环境日趋拥堵,特别是4G以及即将到来的5G商用通信技术的发展,这些进步、发展虽然实现了通信领域的高数据率以及高传输速率,另一方面带来的却是因为占用了更高的射频频段,加剧了整个射频频谱环境的拥挤。因而解决频谱拥堵问题实现频谱共享也引起了学界的广泛关注。一般来说从雷达的角度实现与通信系统的频谱共享方案多样。术语“频谱共享”[1-3]是在十多年前引入的,最近文献[4-5]从宏观角度提出几种缓和频谱拥堵的指导性方案,并指出波形设计实现频谱共享目前是较为热门的理论研究方向。为提高频谱资源利用率,文献[6]较早地为实现频谱共享提供了一种波形设计技术,以使雷达信号具有某些期望的谱槽(凹槽)。该方法在[7]中推广到连续相位波形,在特定频率上放置多个凹槽。文献[8]通过实验分析了上述两种方法的有效性,同时提出了一种超宽带雷达波形,可以显著降低超宽带(UWB)雷达对其他外部射频源的干扰。近年来,文献[9]将解决频谱共存问题与设计雷达发射方向图结合进行波形设计,提出了一种有效的循环算法。而传统意义上针对雷达常用的线性调频信号设置凹槽的方法在文献[10-13]完成了从在频谱中特定频段设置单个到多个频谱凹槽的算法研究。其中文献[12]提出的算法是一种具有低计算复杂度的迭代投影算法。
与上述设置零陷进而产生阻带的方法不同,文献[14]提出了一种“整形”算法,这是一种设计所需频谱形状序列的有效方法。此外该团队还提出了不考虑检测的一种循环方法来合成具有频谱阻带和相关约束的波形[15]。利用凸优化手段是研究波形优化的常用技术手段,文献[16]研究波形设计问题是通过施加频域与时域的现实条件约束来进行波形优化。通过添加不同的约束,建立模型,最终通过凸优化手段得到理想的波形。文献[16]的方法计算操作次数比较高,并且恢复信号的所利用的“秩一分解理论”[16]效果特别理想但却只能针对一些固定的约束问题。文献[17]进一步总结文献[16]并提出了一种基于优化理论的波形设计框架,提高目标检测概率,同时控制特定频段产生的干扰能量和传输波形的某些理想特征。还提出了在信号无关和信号相关干扰的情况下,雷达波形优化的算法。此后,文献[18]又将峰均功率比约束纳入,将信干噪比与频谱兼容性约束加权重新构造目标函数,求解优化问题,并且需要合理调整Pareto参数。 文献[19]针对MIMO雷达设计恒模信号,并且将发射端和接收端联合考虑,进行优化,添加相似性约束并建立了凸优化模型,利用高斯随机化方法恢复恒模信号,但未考虑设计信号实现频谱共享,高斯随机化方法适用性更广一些,但计算复杂度随着随机化次数升高明显。
本文针对时域约束与频域约束,通过优化雷达发射波形实现频谱共享,采用一种更为简单明了的方法,解决模型凸松弛之后寻找优化信号的问题,能够降低计算次数,避免了复杂的数学运算,易于理解与仿真,同时在性能上能够得到保证。最后,针对算法误差做分析,给出了误差表达式并结合仿真实例进行说明。
2 建立模型
假设s∈CM′×1为离散时间信号,表征着雷达发射的探测信号波形,其中M′代表离散信号的长度。令y为M′维的快时间维列向量,代表接收信号,那么,从接收一侧考虑,可以建立模型为:
y=βs+n′
(1)
其中β代表着从感兴趣的目标距离方向角范围内的信道传播和后向散射效应的一个复参数,n′代表着独立信号干扰,包括接收机内部热噪声、干扰信号、交叠频段干扰。总之,把n′建模为一个零均值,协方差矩阵为M的复高斯随机向量。
先对信号从时域上施加约束,对信号施加恒定能量约束,假定约束能量是单位的,即s†s=1。另外,雷达发射探测信号为获得较高的距离分辨率,需要探测信号具有良好的自相关特性,相似性约束通过控制信号与参考信号的相似程度能够很好的约束信号的变化,通常,选择探测信号为具有良好自相关特性的常用雷达信号——线性调频信号。用数学表示为:‖s-s0‖2≤κ[16],其中κ≥0,s0为参考信号。
假定该雷达系统工作在C波段以下频率,将通信作为该工作频段上的主用户,设计雷达探测信号与通信实现频谱共享,需要雷达探测信号尽可能减少在重叠频段上对通信系统的影响。可以假设
其中K为重叠频带个数,这里的k与‖s-s0‖2≤κ需要区分,下同)而
分别表示第k个频带的归一化的下边界和上边界频率。从频域的角度来看,可以将k个频带上发射的能量表示为:
s†RIs≤EI[16]
(2)
其中,
(符号“†”表示共轭转置,下同),rk的不同代表着不同的频段允许相互干扰的能力。根据
为能量谱密度(energy spectral density,简写为ESD),其中
为
(p,q)∈{1,…,M′}2
(3)
以雷达系统的检测性能为优化指标,选择输出信干噪比作为优化对象。假设信号通过线性系统,ω作为系统系数,由(1)知,yout=βω†s+ω†n′,则输出信干噪比:
sinr=|β|2|ω†s|2/E[(ω†n′)2]=
|β|2|ω†s|2/E(ω†n′n′†ω)=|β|2|ω†s|2/ω†Mω
(4)
因为M1/2M-1/2=I
则
利用施瓦兹不等式有
(5)
令R=M-1,有
max
SINR=|β|2s†M-1s=|β|2s†Rs
(6)
为了简便,先假设信号能量为1且以信号的频谱能量做约束,得到:
max s†Rs
s.t. s†s=1,s†RIs≤EI
(7)
进一步添加相似性约束,以便优化得到的信号更加贴近工程实践,得到:
s.t. s†s=1,s†RIs≤EI,‖s-s0‖2≤κ
(8)
很明显问题(7)和(8)的约束条件中能量约束是非凸集合,问题模型(7)、(8)都属于二次约束二次规划问题(QCQP)。但是,问题(7)可以做简单的凸松弛处理,可以很快找到一个解。与此对应的是问题(8)复杂许多,因为除了能量约束以外的两个约束条件会使得可行集合压缩,甚至相互矛盾导致空集,Aubry等人在文献[16]做了详细阐述。本文对于模型(7)和(8),仍利用凸优化方法解决,并在恢复优化信号时以两种不同于文献[16]的方法得到,并作对比。
3 算法构建
本节对模型(7)、(8)进行凸松弛,给出解决问题的不同于参考文献[16]的具体算法。针对模型(7),可以很容易的发现该问题是可解的非凸问题。主要的困难集中在能量约束上,抛掉能量约束s†s=1,作凸松弛处理,并将模型(7)写成以下“迹”的形式
(9)
使用凸优化工具箱,可以得到S。下面提出两种算法计算得到ss†=S。
3.1 随机化方法生成优化信号
利用高斯随机化方法[19],随机生成X组零均值高斯随机向量,即ξi∈CM′×1,i=1,...,X,并且同时满足ξi～N(0,S),每生成一组随机向量ξi 代入目标函数表达式中得到相应的SINR,选取使得SINR最大的随机向量作为最终的信号。详细算法步骤如下表所示。
表1 高斯随机化方法
Tab.1 Gauss randomization method
算法1 高斯随机化方法1.利用凸优化工具求解QEI(7)模型并得到复矩阵S′,同时计算S′的秩。2.若rank(S′)=1,直接对S′进行特征值分解即得到最优解。3.若rank(S′)>1,以均值为零,协方差矩阵为S′随机生成向量X组随机向量,同时代入信干噪比表达式中。4.在X组中选取使得信干噪比最大的随机向量作为优化信号s′。
使用该方法,选取生成的随机化次数越高,便能够很好地得到高度近似解,效果越好。
3.2 利用已知矩阵元素近似得到优化信号
下面提出一种新方法,该方法可以以更高的精度,更加接近最优解,同时也可以用于求解模型(8)。
在求解QEI模型并得到复矩阵S′之后,考虑秩一矩阵S=sis†的某一行或者列,这里假定是第i行,i∈[0,M′],“* ”表示共轭。那么假定S已知,即已经知道S的第i行元素Si,:,若si再已知,自然可以得到s,显然s†=Si,:/si。
选择
因此
可以令φ=2π/Y,其中Y∈Z+为一正整数。因此
在这里要强调的是,这个问题已变为,已知一个赫米特矩阵S′满足S′≈ss†,如果不存在某一行或者列使得到的ss†很接近S′,那么也不存在误差较小的s。至此,问题转换为基于先验信息寻找ss†极其接近S″=s′s′†,其中S″作为与已知的S′最为接近的秩一矩阵。具体算法流程如下:
表2 利用已知矩阵元素近似
Tab.2 Approximation using known matrix elements
算法2 利用矩阵行(列同理)元素得到优化解1.利用凸优化工具求解QEI模型并得到复矩阵S′,同时计算S′的秩。2.若rank(S′)=1,直接对S′进行特征值分解即得到最优解。3.若rank(S′)>1,选择 S′中某一行(列同理),给定Y,满足ϕ=2p/Y,Y∈Z+,s″=(Si,:S′i,iejϕ×d)†,d=1,2,…,Y。计算Dd=‖s″s″†-S′‖m∞。4.选取Dd最小的s″,并计算目标函数值。5.若要求精度比较高时,利用全行列元素求解则步骤2中令i=1,完成2,3,4步骤,储存Dd和对应的s″,i=i+1,直到i=M′。选出M′组数据中使得Dd最小的s″并输出。
选择使用S中的某一行或列利用上述方法是合理的,但需要给出结果的可信程度也就是要对误差作分析说明。
针对误差证明如下(特殊说明:在证明过程中S表示已知的模型松弛求解得到的矩阵,s是算法2得到的优化信号,此处的n代表矩阵维数,其中i, j,l,m,p,q∈[1,n]):
假设A=aa†,其中A为最为接近S的秩一矩阵,利用已知元素得到的信号的误差分析如下:
(10a)
(10b)
对(10a)利用三角不等式
(11a)
对(10b)利用三角不等式
(11b)
这里需要注意的是
(12)
所以有
(13)
利用三角不等式
(14)
进一步放缩
≤
(15)
则
(16)
这里定义一个量 Ψ来说明算法2信号得到信号的可信程度
(17)
其中
表示最准确的秩一矩阵与S的误差,同时 ς 表示与利用已知矩阵所求信号成一定关系。分析 Ψ的难点在于
但可以明确的是,当S接近秩一矩阵时,
(18)
利用某行,或者利用任意行元素 Ψ1只会发生微小变化。且此时,‖S‖m
会远大于
令
为式(19),只需要确定 Ψ的上下界,便可以判断算法2信号的可信度。下界用Ψ
表示,上界用 ΨΔ 表示,并且有
Ψ
(20)
(21)
该算法相对于前一种算法具有其独特的优势,高斯随机化方法,完全依赖于生成随机数来保证算法的有效性,这就导致随机生成数据的次数必须足够大,那么,要保证效果,通常生成随机数的次数较随着向量维数的增加而增加更大。而方法2,选择其中一行(列,同理),带来的误差并没有很大,同时利用旋转次数Y使得信号在一定的计算机精度范围上下略微浮动便可以“更加充分”利用该行信息。利用全部行或列元素(以下简称“利用全行元素”)进行算法2需要遍历所有的行(或列),计算次数也随着信号维度增加而线性升高。
为叙述完整,下面针对模型(8)依据文献[16]进行简单概括。将(8)写成如下形式:
(22)
(23)
其中σ=(1-κ/2)2。
这是因为
(24)
文献[16]指出式(24)与|s†s0
存在某种关系, 即s′满足(24),然后有
这说明s和s′存在一个映射关系,显然当s=s′ej arg(s′†s0)时成立。
总之,和分析模型(7)方法一样,将(8)能量约束抛掉,松弛成可解的SDP问题[16],
(25)
在求解该问题之前,需要分析约束条件所构成的可行区域问题。由q中约束条件分析应取
寻找问题(26)的解,并取
构成可行集。
表示问题
的优化值。
(26)
在这里,仍利用前面所提出的采用高斯随机化方法找寻优化解,并令s=s′ej arg(s′†s0)。
同理,将新提出的算法2应用于解决该问题,进行仿真验证并分析高斯随机化方法与算法2,比较两种方法之间的优劣势,同时比较算法2与文献[16]解决方法的效果。
4 仿真结果
4.1 仿真验证分析
假设相似对象为一线性调频信号满足:
(27)
其中
同时令fs=810 kHz。
假设在归一化频带上,存在以下任意取的几个重叠频带,
对于干扰,假定各干扰源相互独立。干扰信号的协方差矩阵建模为
(28)
在这里σ0=0 dB代表内部热噪声水平。K=4代表交叠频带个数,KJ=2,代表未许可的窄带干扰信号源的个数,σI,k=10 dB,k=1,2,…,K,σJ,1=50 dB,σJ,2=40 dB,
其中, fJ,k 代表着干扰频率,同时fJ,1/fs=0.7, fJ,2/fs=0.75,(仿真图中κ用k表示与频段个数k意义不同)。
对于模型(7)设EI=0.03,分别根据算法1和2利用凸优化工具箱做仿真实验效果如下。
图1 两种算法信号能量谱密度比较
Fig.1 Comparison of ESD of signal numbers obtained by two algorithms
从仿真结果来看,高斯随机化方法效果是很明显的,在本仿真数据条件下,在重叠的频段,信号的能量能够降低30 dB,而方法2效果相较而言,更为明显,这说明在约束条件下,能够近似满足达到了S′=ss†,也就是说S′近似地为秩一矩阵。即使不能满足,至少有一行(列同理)是最为近似地满足,可以对能量损失进行矫正,效果是很明显的。当然,在计算操作次数方面,高斯随机化方法比方法2较高,利用单行元素完全可以做到比高斯随机化低几个数量级。
下面对模型(8)进行实验仿真,根据式
与(26)画出可行区域并在图2中标注出的其中几个点的数据然后将这些数据点进行仿真(EI1=0.03,k1=0.0837)(EI2=0.03, k2=0.143)(EI3=0.03, k3=0.236)。
图2 干扰相似性区域
Fig.2 Interference and similarity achievable region
首先,将对高斯随机化方法得到的信号自相关函数(autocorrelation function,简写为ACF)以及能量谱密度进行仿真。采用随机化次数20000次,Y=10。
图3 高斯随机化方法自相关函数
Fig.3 The ACF by Gauss randomization method
图4 利用单行元素算法2信号自相关函数
Fig.4 The ACF of signal obtained by using single-row element algorithm 2
图5 全行元素算法2自相关函数
Fig.5 The ACF of signal obtained by using all-row element algorithm 2
可以看出,算法2近似方法与高斯随机化得到的信号自相关不尽相同,总体上效果差别不大。另外利用单行元素与全部行列元素进行算法2时得到的信号自相关函数因为计算机精度因素在画图时是一致的,当画在同一张图上是会被相互掩盖,这是因单行元素寻找优化解与利用全行元素相差无几,也恰好说明单行元素得到的优化解精度很高。为更好地展示效果,能量谱密度(ESD)也做类似比较。
另外对以下各点使用高斯随机化方法和算法2并作能量谱密度比较如下(EI1=0.03, k1=0.0837)、(EI2=0.074, k2=0.0837)、(EI3=0.126,k3=0.0837)。
图6 能量谱密度比较效果图
Fig.6 Comparison of ESD
图7 能量谱密度比较效果图
Fig.7 Comparison of ESD
图8 EI=0.03,算法2得到的信干噪比随着k变化趋势
Fig.8 The variation trend of SINR of algorithm 2 with k,when EI=0.03
图9 能量谱密度比较效果图
Fig.9 Comparison of ESD
表3 仿真数据点算法2相关数据
Tab.3 Simulation data Point algorithm 2 related data
Eκ0.030.08380.030.1430.030.2360.0740.08380.1260.0838Y▽0.01213.244e-083.743e-053.632e-091.565e-07YΔ 0.01223.244e-083.743e-053.632e-091.565e-07式(19)3.2622e-041.103e-091.46e-061.046e-104.511e-09
图10 使用全行元素算法2 Ψ上下界及式(19)变化曲线
Fig.10 The upper and lower bounds of Ψ and the variation curve of formula (19) when the all-row element algorithm 2 is used
图11 使用全行元素算法2后式(19)变化曲线
Fig.11 The change curve of formula (19) after using the all-row element algorithm 2
为进一步证明方法2的有效性,将仿真结果与文献[16]秩一分解方法进行结果对比,作仿真图比较如下:
图12 自相关函数比较效果图
Fig.12 Comparison of ACF
可以看到秩一分解理论[16]得到的信号与算法2利用全部行列元素得到的信号几乎一样,如果将两种方法得到的自相关函数和能量谱密度曲线放在同一张图中比较会相互覆盖。虽然理论上秩一分解理论更高,但是,算法2得到的信号效果也非常理想。另一方面利用算法2同高斯随机化方法一样不受约束条件的限制,秩一分解理论精度高但是仅限于求解固定的约束问题,可以说各有利弊。如果将(7)中的tr(S)=1 改换为考虑信号恒定幅度松弛问题的diag(S)=1。此时,不能利用”秩一分解”理论恢复信号,但算法2适用性更广,可以利用已知矩阵元素近似的算法2和随机化方法得到信号,而且前者效果比随机化方法也更有优势。
图13 能量谱密度比较效果图
Fig.13 Comparison of ESD
4.2 误差数据与计算操作次数分析
本节对仿真数据进行分析,并给出两种算法的计算操作次数。主要对算法2进行分析,表3给出了在一些点的公式(19)～(21)对应的数值,从这些点可以看出定义式(17)Ψ对应的上下界基本一致,这说明真实存在的秩一矩阵与SDR方法求解之后的矩阵之间的差距非常小,进一步从该表中给出式(19)对应的数值(忽略维度),可以看出算法2找到的信号有效性可以保证。
与表3分析思路相同,进一步分析当EI保持不变,k发生变化式(19)～(21)变化曲线,图10、11说明在模型(8)框架下,利用凸优化求解并利用算法2得到的信号有效性能够得到保证。
从仿真结果来看,单列元素和全行元素效果基本一致,这从侧面证实了利用凸优化方法建立求解SDP问题的有效性。此外,高斯随机化方法的计算操作次数约为O(LN2)[19],另一方面,可以对算法2进行分析,利用N行元素计算操作次数约为O(N3),Y理论上不会影响算法性能,但是考虑到运算精度,可以利用Y扰动信号精度,从而更好地寻找信号。因为该信号维度为162,一般在进行随机化次数时,万次量级一般能够保证稳定的效果,基于此,算法2利用全行元素较低于高斯随机化方法计算操作次数(取Y为10即可),而且利用单行元素时,操作次数约为O(N2)。除此之外,高斯随机化方法若要优化结果信干噪比稳定上升,而提高随机化次数一般需要将L提高一个数量级,但利用单行元素寻找与全行元素寻找优化解,信干噪比都是稳定上升而且效果很理想。
5 结论
本文为缓解射频频谱环境拥堵问题,从波形设计角度出发,利用高斯随机化方法解决了凸优化模型下频谱共享波形设计问题,并进一步提出一种新的算法,相互对比发现能够降低计算操作次数,效果更好。该算法和高斯随机化方法一样可以对接到SDR问题求解,此外,经过理论证明和仿真分析验证了算法的有效性。算法2利用全行(列)元素可以比高斯随机化更为有效、快速地得到优化解,同时和随机化方法一样适用性范围广。本文工作尚未进一步考虑工程实践中信号幅度不能剧烈变化的问题,设计恒模波形实现频谱能量控制进而实现雷达与通信频谱共享问题是下一步的研究方向。
[1] Peha J M. Approaches to spectrum sharing[J]. Communications Magazine IEEE, 2005, 43(2): 10-12.
[2] Ma L, Han X, Shen C C. Dynamic open spectrum sharing MAC protocol for wireless ad hoc network[C]∥New Frontiers in Dynamic Spectrum Access Networks, 2005. DySPAN 2005. 2005 First IEEE International Symposium on. IEEE, 2005.
[3] Visotsky, Kuffner, Peterson. On collaborative detection of TV transmissions in support of dynamic spectrum sharing[C]∥First IEEE International Symposium on New Frontiers in Dynamic Spectrum Access Networks. IEEE, 2005.
[4] Griffiths H, Blunt S, Cohen L, et al. Challenge problems in spectrum engineering and waveform diversity[C]∥Radar Conference (RADAR), 2013 IEEE. IEEE, 2013.
[5] Griffiths H, Cohen L, Watts S, et al. Radar Spectrum Engineering and Management: Technical and Regulatory Issues[J]. Proceedings of the IEEE, 2015, 103(1): 85-102.
[6] Haykin S. Cognitive Dynamic Systems: Radar, Control, and Radio [Point of View][J]. Proceedings of the IEEE, 2012, 100(7): 2095-2103.
[7] Nunn C, Moyer L R. Spectrally-compliant waveforms for wideband radar[J]. Aerospace & Electronic Systems Magazine IEEE, 2012, 27(8): 11-15.
[8] Gerlach K. Thinned spectrum ultrawideband waveforms using stepped-frequency polyphase codes[J]. IEEE Transactions on Aerospace and Electronic Systems, 2002, 34(4): 1356-1361.
[9] 麦超云, 孙进平, 周锐, 等. MIMO雷达方向图与稀疏频率波形联合设计[J]. 信号处理, 2016, 32(11): 1261-1268.
Mai Chaoyun, Sun Jinping, Zhou Rui, et al. MIMO Radar Beampattern and Sparse Frequency Waveform Design[J]. Journal of Signal Processing, 2016, 32(11): 1261-1268.(in Chinese)
[10] Gerlach K, Frey M R, Steiner M J, et al. Spectral Nulling on Transmit via Nonlinear FM Radar Waveforms[J]. IEEE Transactions on Aerospace and Electronic Systems, 2011, 47(2): 1507-1515.
[11] Selesnick I W, Pillai S U, Zheng R. An iterative algorithm for the construction of notched chirp signals[C]∥IEEE Radar Conference. IEEE, 2010.
[12] Cook M R, Higgins T, Shackelford A K. Thinned spectrum radar waveforms[C]∥Waveform Diversity & Design Conference, 2010.
[13] Selesnick, Ivan W, Pillai, et al. Chirp-like transmit waveforms with multiple frequency-notches[C]∥2011 IEEE Radar Conference. 2001: 1106-1110.
[14] Rowe W, Stoica P, Jian L. Spectrally Constrained Waveform Design [sp Tips & Tricks][J]. IEEE Signal Processing Magazine, 2014, 31(3): 157-162.
[15] He H, Stoica P, Li J. Waveform design with stopband and correlation constraints for cognitive radar[C]∥International Workshop on Cognitive Information Processing. IEEE, 2010.
[16] Aubry A, De Maio A, Piezzo M, et al. Radar waveform design in a spectrally crowded environment via nonconvex quadratic optimization[J]. IEEE Transactions on Aerospace and Electronic Systems, 2014, 50(2): 1138-1152.
[17] Aubry A, Carotenuto V, Maio A D, et al. Optimization theory-based radar waveform design for spectrally dense environments[J]. IEEE Aerospace and Electronic Systems Magazine, 2016, 31(12): 14-25.
[18] 邹鲲, 骆艳卜, 李伟, 等. 频谱拥挤环境中峰均功率比约束的认知雷达发射波形设计[J]. 电子与信息学报, 2018, 40(7): 255-259.
Zou Kun, Luo Yanbo, Li Wei, et al. Cognitive Radar Waveform Design with A Peak to Average Power Ration Constraint for Spectrally Dense Environments[J]. Journal of Electronics & Information Technology, 2018, 40(7): 255-259.(in Chinese)
[19] Cui G, Li H, Rangaswamy M. MIMO Radar Waveform Design With Constant Modulus and Similarity Constraints[J]. IEEE Transactions on Signal Processing, 2014, 62(2): 343-353.
