Three-dimensional Imaging of Tomographic SAR Based on Adaptive Elevation Constraint
-
摘要: 层析合成孔径雷达成像(TomoSAR)是2010年以来SAR成像领域尤其是城市三维成像的热门研究方向。但在TomoSAR三维重建中,相位缠绕会引起高程散射剖面的周期性谱峰,并导致散射体高程向位置的错误估计和三维成像结果中建筑点云的分层,即高程模糊。该文针对这一现象,提出一种自适应调整高程搜索范围的方法,以提升散射体高程估计的准确度,并改善高程模糊。该方法首先进行场景的高度预估计,然后根据高度预估计构建高程采样中心线性函数并计算搜索半径,从而确定并更新各像素的高程搜索范围,保留真实谱峰并隔离模糊峰值。机载和星载的实测数据实验表明所提方法明显改善了高程模糊和伪影问题,提高了三维点云的空间集中度和连续性。
-
关键词:
- 层析合成孔径雷达成像(TomoSAR) /
- 相位缠绕 /
- 高程模糊 /
- 散射体 /
- 高程搜索范围
Abstract: Synthetic Aperture Radar Tomography (TomoSAR) has emerged as a hot research topic in the field of SAR imaging, particularly for three-dimensional (3D) urban imaging in recent years. However, in TomoSAR 3D reconstruction, due to the phase unwrapping difficulty, periodic spectral peaks appear in the reconstruction results of the reflectivity profile along the elevation. This results in errors in estimating the elevation locations of the scatterers and causing layering effects in 3D imaging results, which is the elevation ambiguity. In light of this phenomenon observed in TomoSAR, a method for the adaptive adjustment of the elevation search range is proposed to improve the accuracy of the elevation estimation and reduce elevation ambiguity. In this method, the height of the scene is first estimated, the linear function of the elevation sampling center is subsequently constructed based on the height pre-estimations, and the search radius is finally calculated. Thereafter, the elevation search range of each pixel in the SAR image is determined and updated, preserving the true spectral peaks while isolating the ambiguity peaks. The experimental results for airborne and spaceborne measured data demonstrate that the proposed method significantly improves elevation ambiguity and artifacts-related issues while also improving the spatial concentration and continuity of 3D point clouds.-
Key words:
- TomoSAR /
- Phase unwrapping /
- Elevation ambiguity /
- Scatterers /
- Elevation sample range
-
1. 引言
相较于传统的相控阵雷达,多输入多输出(Multi-Input Multi-Output, MIMO)雷达能够发射多种波形信号,这些信号可能相互关联或不相关。得益于波形的多样性,MIMO雷达的目标检测能力、抗干扰能力、目标识别和分类能力得到了显著增强,发射波束设计的灵活性也得到了很大的提升[1,2]。雷达与通信系统在频谱拥挤环境下的共存问题是近年来学者关注的一个重要问题。 当两种系统占用相同的频段时,就需要这种兼容性。例如,在宽带传输中,雷达系统占用很大的带宽,因此可能会与周边无线通信网络的频谱发生重叠[3]。目前关于频谱共存问题的研究大多数是从雷达角度出发,通过对发射波形设计的优化来实现频谱共存。
在MIMO雷达系统的收发联合设计中,可以考虑不同的约束条件来设计和优化发射波形[4-9]。例如,通过施加恒模(Constant Modulus, CM)约束或者峰均功率比(Peak-to-Average Ratio, PAR)约束使得雷达发射机能在饱和工作状态将性能发挥到最大,通过施加相似性约束来确保优化后的发射波形具有良好特性,通过施加频谱约束控制发射波形的频谱能量来实现频谱共存。存在信号相关干扰时,为了检测扩展目标,在发射波形的能量约束下进行收发联合设计,以最大化信干噪比(Signal to Interference pulse Noise Ratio, SINR)[10]。Cui等人[11]将雷达波形优化问题转化为纯相位约束问题,提出了坐标下降(Block Coordinate Descent, BCD)的迭代算法。Liang等人[12]设计了保持方向图形状的相位阵列优化方案,提出了基于交替方向乘子法(Alternating Direction Method of Multipliers, ADMM)的相位恢复算法。Palomar团队[13]考虑了PAR约束下的恒模波形设计问题,提出了优化最小(Majorization-Minimization, MM)算法用以设计雷达发射编码。文献[14]针对信号相关杂波环境下的运动目标检测问题,研究了一种稳健的慢时间发射波形和接收滤波器联合设计方法,并利用Dinkelbach 程序进行求解。
针对雷达和通信系统频谱共存问题,De Maio和Abury等人[15]以输出SINR最大化为准则,研究了在能量约束和相似性约束下的雷达频谱兼容波形设计方法,并利用半正定规划与秩一分解(Rank-one decomposition)理论对雷达发射序列进行优化。虽然基于半正定规划的秩一分解技术能够获得好的近似解,但求解半正定规划问题的计算复杂度较高,不太利于工程的实现。在此基础上,Wu等人[16]考虑频谱约束条件,提出了在MM算法框架下的改进(Successive Convex Approximation, SCA)算法,降低了每次迭代的计算复杂度。Tang等人[17]提出了一种基于ADMM的快速算法,数值结果证明该算法在计算时间上优于半正定规划算法,所设计的波形也能折中SINR、频谱共存和自相关函数的特性。随后,Fan等人[18]考虑了宽频带MIMO雷达在频谱密集环境下的发射波束合成问题,提出了一种改进的ADMM (Majorization-ADMM, M-ADMM)算法,通过利用一个光滑的凸函数对原始非光滑的非凸目标函数进行近似,然后利用ADMM算法的更新规则并行求解约束优化问题,最终降低了求解的计算时间。
针对频谱约束条件下MIMO雷达系统的联合收发设计问题,已有的方法包括顺序优化算法(Sequential Optimization Algorithm, SOA)[19]、逐次二次约束二次规划(Quadratically Constrained Quadratic Programming, QCQP)求解法[20]、BCD法[11,21]和MM求解法[13,22]。在文献[23]中考虑了PAR和相似度约束,并利用乘子块逐次上界最小化法求解了文献[19-22]中的优化问题。为保证与授权辐射器的频谱兼容性,在频谱、相似性和能量约束下的联合设计问题也通过文献[19]中的MM方法进行求解。在文献[24]中,作者利用半定松弛(Semie-Definite Relaxation, SDR)方法[25]将上述方法扩展到包括频谱、相似性和能量约束的问题中。然而,对于多约束优化问题[20],文献[24]中的可行集并不能总是得到保证。在文献[26]中,通过在相似度和通信速率约束下联合设计MIMO雷达收发器的空时协方差矩阵,确保MIMO雷达和MIMO通信系统的共存。文献[27]设计了发射信号和接收滤波器组的迭代优化程序。它单调地改善了最坏情况的SINR,并收敛到一个平稳点。该方法还能够通过根据规定的时间或复杂度要求适当地选择迭代次数权衡SINR性能和算法复杂度。
ADMM算法通过将原问题分解为更小的子问题进行求解,更好地适用于大规模的约束凸问题。然而,对于一般的非凸问题,ADMM的收敛性分析仍然是一个开放性问题。迭代分块连续上界最小化方法(Block Successive Upper-bound minimization Method, BSUM)[28-32]通过引入原始目标的上界函数,可以将求解的非凸问题转换为依次优化凸函数问题。本文以输出SINR最大化为准则,构建了频谱约束和相似性约束条件下的发射波形与接收滤波器的联合设计问题,结合ADMM和BSUM算法的优点,设计了BSUM算法的变体(Alternating BSUM, ABSUM)算法,提出了一个原始对偶框架,该算法可以收敛到一个静态稳定点,使算法能够有效地解决约束非凸问题,数值仿真实验证明了该算法能够在合理的计算时间内取得更好的SINR,得到更优的目标检测性能。
符号说明:
(⋅)T ,(⋅)H ,(⋅)∗ ,||⋅|| ,|⋅| 和E{⋅} 分别表示转置、共轭转置、共轭、范数、取绝对值和统计期望值;IN×N 表示N×N 的单位矩阵;CN×1 为包含N×1 向量的复数集合;CN×N 为包含N×N 矩阵的复数集合;⊗ 表示Kronerker积;diag{⋅} 表示矩阵的对角元素;vec{⋅} 表示将矩阵按列排列形成列向量;ℜ{⋅} 和ℑ{⋅} 分别表示取实部和取虚部操作。2. 信号模型
假设MIMO雷达系统有
Nt 个发射天线和Nr 个接收天线,接收阵列和发射阵列都采用间距为半波长的均匀线阵。每个发射天线发射不同的雷达编码ci(l) ,其中i=1,2,⋯,Nt ,l=1,2,⋯,Ns ,Ns 是发射波形的采样总数。c(l)=[c1(l),c2(l),⋯,cNt(l)]T 为Nt 个发射天线在第l个采样时刻的发射波形,C=[c(1),c(2),⋯,c(Ns)]∈CNt×Ns 为发射波形的矩阵形式。NrNs 维接收信号的复向量形式可以表示为[33]y=α0H(θ0)c+K∑k=1αkH(θk)c+v (1) 其中,
c=vec(C) ,α0 和αk 分别表示目标和第k个干扰源的幅值。向量n既包括内部热噪声,也包括未知的、未经许可的、可能是敌对干扰器造成的干扰信号,以及与目标雷达共存相同频率的授权电信网络。此外,向量n被建模为一个复的、零均值的、圆对称的高斯随机向量,协方差矩阵E{nn†}=M 。θ0 和θk 分别为目标和第k个干扰源的角度,且θk≠θ0 。H(θ) 为均匀线性阵列天线的转向矩阵,H(θ)=INs⊗[ar(θ)aTt(θ)] 。以输出SINR作为性能指标,联合优化发射波形与接收滤波器来最大化输出SINR。在实际场景中,可能无法获得目标角度和干扰的准确信息,因此,可以合理地假设θ0 和θk 是不相关的均匀分布随机变量,分别具有已知的平均值˜θ0 和˜θk ,即θ0∼H(˜θ0−δt2,˜θ0+δt2)θk∼H(˜θk−δk2,˜θk+δk2),k=1,2,⋯,K (2) 假设接收信号y被接收滤波器
w∈CNrNs 过滤,则SINR(c,w) 可以表示为SINR(c,w)=|α0|2|wHH(θ0)c|2wHΣI(c)w+wHMw (3) 其中,
wHΣI(c)w 表示为杂波能量,ΣI(c)=∑Kk=1σ2kH(θk)ccHHH(θk) 。其中αk,k=1,2,⋯,K 被假设为独立的零均值随机变量,其方差为σ2k={|αk|2} 。为了方便计算,将SINR重新表示为以下两个等价表达式:SINR(c,w)=wHΨ(c)wwHΨin(c)w=cHΠ(w)ccHΠin(w)c (4) 其中,
Π(w)=|α0|2(WT⊗INt)Oδt˜θt(WT⊗INt)H ,Πin(w)=∑Kk=1σ2k(WT⊗INt)Oδk˜θk(WT⊗INt)H+||w||2M ,Ψ(c)=|α0|2(CT⊗INr)HPδt˜θt(CT⊗INr)H 和Ψin(c)=∑Kk=1σ2k(CT⊗INr)Pδk˜θk(CT⊗INr)H+M ,Pδ˜θ 和Oδ˜θ 表示为Pδ˜θ=∫˜θ+δ2˜θ−δ2(at(θ)aHt(θ))⊗(ar(θ)aHr(θ))dθ ,Oδ˜θ=∫˜θ+δ2˜θ−δ2(a∗r(θ)aTr(θ))⊗(a∗t(θ)aTt(θ))dθ ,W是w的矩阵形式,为了使得雷达系统与周边的授权通信系统频谱共存,在优化模型中加入频谱兼容性约束,频谱兼容性约束为[15]cHRc≤EI ,其中cHRc 表示雷达波形在M个通信系统频带上频谱能量总和,EI 是最大允许干扰量,R是频谱兼容性矩阵,可以表示为R=M∑i=1ωi(Γi⊗Mi) (5) 其中,
Γi 和Mi 被定义为Mi(p,q)={sinθi2−sinθi1, p=qejπsinθi2(p−q)−ejπsinθi1(p−q)jπ(p−q), p≠qΓi(m,l)={fiupper−filower, m=lej2πfiupper(m−l)−ej2πfilower(m−l)j2π(m−l), m≠l (6) 其中,
ωi≥0,i=1,2,⋯,M 是对应于第i个共存无线网络的权重系数,θi1 和θi2 分别表示第i个周围电磁辐射器的上下夹角,每个共存辐射器被认为位于一个标准化的空间Θi=[sinθi1,sinθi2] 内。fiupper 和filower 分别表示第i个无线网络较高和较低的归一化频率。频谱兼容约束意味着在这M个频段上所有Nt 个发射波形的总能量不能大于一个阈值,使得雷达能够与通信设备频谱共存。此外,任意波形会导致较大的模数变化、较差的多普勒分辨率和较高的峰值旁瓣电平,因此还需对发射波形施加相似性约束,这个约束为||c−c0||∞≤ε ,c0 为参考波形,ε 是用户定义的参数,用于控制相似度。该式可以表示为非齐次二次不等式约束,即(c−c0)HEn(c−c0)≤ε2,n=1,2,⋯,NtNs,En(i,j)={1, i=n,j=n0,其他 (7) 基于SINR最大化的优化准则,在能量、相似性和频谱兼容约束条件下,雷达发射波形和接收滤波器联合设计问题可以表述为
maxw,c SINR(c,w)=wHΨ(c)wwHΨin(c)ws.t. cHc=1cHRc≤EI(c−c0)HEn(c−c0)≤ε2,n=1,2,⋯,NtNs (8) 从式(8)可以明显看出,此问题为非凸问题,其包含非凸目标函数、非凸二次等式约束和非齐次二次不等式约束,是NP-hard问题[34],很难直接求得上述问题的全局最优解。因此,需要设计一种新的迭代优化算法来求解问题(8)。
3. 波形设计算法
ADMM算法可以将原来复杂的问题分解成为更小的子问题进行求解,更好地适用于大规模的约束问题。也被广泛应用于求解雷达波形设计问题[35-37]。针对凸问题,ADMM算法能够保证收敛到全局最优解,然而针对非凸问题,其收敛性仍是一个开放性问题。结合ADMM算法和BSUM算法的特点,本文提出一种新的优化方法,称其为ABSUM算法,该算法可以收敛到一个静态稳定点。其主要思想是在ADMM算法框架下,使用BSUM算法引入目标函数的上界函数,然后应用ADMM框架更新规则求解约束优化问题。所提出ABSUM算法兼顾了ADMM算法和BSUM算法的优点,具有更好的收敛性。首先将式(8)中的问题改写为实值形式:
maxw,cr cTrΠinr(w)crcTrΠr(w)crs.t. cTrcr=1cTrRrcr≤EI(cr−c0r)HErn(cr−c0r)≤ε2,n=1,2,⋯,NtNs (9) 其中,
c0r ,cr ,Πinr(w) ,Πr(w) ,Rr 和Ern 定义为c0r=(ℜ{c0}ℑ{c0}),cr=(ℜ{c}ℑ{c}),Πinr(w)=(ℜ{Πin(w)}−ℑ{Πin(w)}ℑ{Πin(w)} ℜ{Πin(w)}),Πr(w)=(ℜ{Π(w)}−ℑ{Π(w)}ℑ{Π(w)} ℜ{Π(w)}),Rr=(ℜ{R}−ℑ{R}ℑ{R} ℜ{R}),Ern=(ℜ{En}−ℑ{En}ℑ{En} ℜ{En})。 此外,
cTrcr=1 表示能量约束,用来控制发射信号的能量,cTrRrcr≤EI 为频谱约束,使得雷达发射波形能与周围通信设备达到频谱共存的目的,(cr−c0r)HErn(cr−c0r)≤ε2,n=1,2,⋯,NtNs 为相似性约束,无约束优化会导致信号具有显著的模量变化、较差的距离分辨率、较高的峰旁瓣电平,以及不理想的模糊函数响应,这些缺点可以通过强制解决方案与已知波形c0 相似来部分规避。针对优化问题(9),引入辅助变量tr 和一个线性约束cr−tr=0 ,将问题(9)修正为maxw,cr,tr cTrΠinr(w)trcTrΠr(w)trs.t. cr−tr = 0 cTrtr=1cTrRrtr≤EI(cr−c0r)TErn(tr−c0r)≤ε2,n=1,2,⋯,NtNs (10) 因为在问题(9)中包括非凸目标函数和二次等式约束,其为一个非凸问题,直接求解比较困难。通过引入辅助变量可以把问题(9)中的非凸目标函数转换为关于每个变量的双拟凸函数。此外,问题(9)中的非凸二次等式约束可以修改为双仿射等式约束,即
cr 和tr 的联合仿射。换句话说,对于固定tr ,它在cr 上是拟凸的,对于固定cr ,它在tr 上是拟凸的。因为目标函数是关于ck 的拟凸函数,拟凸函数的上镜图也是凸函数,因此可以通过最小化F(wk+1,cr,tkr) 的上界函数来求解ck ,同理也可求解tk 。因此,利用ADMM算法可以求解问题(10)。在ADMM框架下,将等式约束置于原函数的增广拉格朗日函数中,而不是始终保持等式约束。更确切地说,问题(10)中的增广拉格朗日量为L(w,cr,tr,u,v)=F(w,cr,tr)+uT(cr−tr) + ρ12||cr−tr||2+v(cTrtr−1) + ρ22||cTrtr−1||2 (11) 其中,
F(w,cr,tr)≜ ,{\boldsymbol{u}} \in {\mathbb{R}^{2{N_{\rm{t}}}{N_{\rm{s}}}}} 和v \in \mathbb{R} 是对偶变量,{\rho _1},{\rho _2} > 0 为惩罚参数。根据ADMM更新规则,在第(k + 1) 次迭代中,ADMM算法由以下迭代组成:\left. \begin{split} & {{\boldsymbol{w}}^{k + 1}} = \arg \mathop {\min }\limits_{\boldsymbol{w}} {\text{ }}\mathfrak{L}\left( {{\boldsymbol{w}},{\boldsymbol{c}}_{\rm r}^k,{\boldsymbol{t}}_{\rm r}^k,{{\boldsymbol{u}}^k},{v^k}} \right) \\ & {\boldsymbol{t}}_{\rm r}^{k + 1} = \arg \mathop {\min }\limits_{{{\boldsymbol{t}}_{\rm r}} \in {\mathcal{D}_{\rm{t}}}} {\text{ }}\mathfrak{L}\left( {{{\boldsymbol{w}}^{k + 1}},{\boldsymbol{c}}_{\rm r}^{k + 1},{{\boldsymbol{t}}_{\rm r}},{{\boldsymbol{u}}^k},{v^k}} \right)\\ & {\boldsymbol{c}}_{\rm r}^{k + 1} = \arg \mathop {\min }\limits_{{{\boldsymbol{c}}_{\rm r}} \in {\mathcal{D}_{\rm{c}}}} {\text{ }}\mathfrak{L}\left( {{{\boldsymbol{w}}^{k + 1}},{{\boldsymbol{c}}_{\rm r}},{\boldsymbol{t}}_{\rm r}^k,{{\boldsymbol{u}}^k},{v^k}} \right)\\ & {{\boldsymbol{u}}^{k + 1}} = {{\boldsymbol{u}}^k} + {\rho _1}({\boldsymbol{c}}_{\rm r}^{k + 1} - {\boldsymbol{t}}_{\rm r}^{k + 1})\\ & {v^{k + 1}} = {v^k} + {\rho _2}\left(\left({\boldsymbol{c}}_{\rm r}^{k + 1}\right)^{{\text{T}}}{\boldsymbol{t}}_{\rm r}^{k + 1} - 1\right) \end{split}\right\} (12) {\mathcal{D}_{\rm{t}}} 和{\mathcal{D}_{\rm{c}}} 为两个集合,即{\mathcal{D}_{\bf{t}}} \triangleq \left\{ {{{\boldsymbol{t}}_{\rm{r}}}\left| \begin{gathered} {\left( {{\boldsymbol{c}}_{\rm{r}}^{k + 1}} \right)^{\text{T}}}{{\boldsymbol{R}}_{\rm{r}}}{{\boldsymbol{t}}_{\rm{r}}} \le {E_I},{\left( {{\boldsymbol{c}}_{\rm{r}}^{k + 1} - {{\boldsymbol{c}}_{0{\rm{r}}}}} \right)^{\text{T}}}{\boldsymbol{E}}_n^{\rm{r}}\left( {{\boldsymbol{t}}_{\rm{r}}^{} - {{\boldsymbol{c}}_{0{\rm{r}}}}} \right) \le {\varepsilon ^2} \\ \qquad n = 1,2,\cdots,{N_{\rm{t}}}{N_{\rm{s}}} \\ \end{gathered} \right.} \right\} 和{\mathcal{D}_{\rm{c}}} \triangleq \left\{ {{{\boldsymbol{c}}_{\rm{r}}} \left| \begin{gathered} {\boldsymbol{c}}_{\rm{r}}^{\text{T}}{{\boldsymbol{R}}_{\rm{r}}}{\boldsymbol{t}}_{\rm{r}}^k \le {E_I},{\left( {{{\boldsymbol{c}}_{\rm{r}}} - {{\boldsymbol{c}}_{0{\rm{r}}}}} \right)^{\text{T}}}{\boldsymbol{E}}_n^r\left( {{\boldsymbol{t}}_{\rm{r}}^k - {{\boldsymbol{c}}_{0{\rm{r}}}}} \right) \le {\varepsilon ^2}, \\ \qquad n = 1,2,\cdots,{N_{\rm{t}}}{N_{\rm{s}}} \\ \end{gathered} \right.} \right\} 。3.1 接收滤波器优化
当
{{\boldsymbol{c}}^k} 和{{\boldsymbol{t}}^k} 固定时,接收滤波器向量w的优化问题可以写为\mathop {\max }\limits_{\boldsymbol{w}} {\text{ }}\frac{{{{\boldsymbol{w}}^{\text{H}}}{\boldsymbol{\varPsi }}({{\boldsymbol{c}}^k},{{\boldsymbol{t}}^k}){\boldsymbol{w}}}}{{{{\boldsymbol{w}}^{\text{H}}}{{\boldsymbol{\varPsi }}_{{\text{in}}}}({{\boldsymbol{c}}^k},{{\boldsymbol{t}}^k}){\boldsymbol{w}}}} (13) 其中,
{\boldsymbol{\varPsi }}({\boldsymbol{c}},{\boldsymbol{t}}) = {\text{SNR}}({{\boldsymbol{C}}^{\text{T}}} \otimes {{\boldsymbol{I}}_{{N_{\rm{r}}}}}){\boldsymbol{{P}}}_{{{\tilde \theta }_{\rm{t}}}}^{{\delta _{\rm{t}}}}{({{\boldsymbol{T}}^{\text{T}}} \otimes {{\boldsymbol{I}}_{{N_{\rm{r}}}}})^{\text{H}}} ,{{\boldsymbol{\varPsi }}_{{\text{in}}}}({\boldsymbol{c}},{\boldsymbol{t}}) = \displaystyle\sum\nolimits_{k = 1}^K {{\text{INR}}({{\boldsymbol{C}}^{\text{T}}} \otimes {{\boldsymbol{I}}_{{N_{\rm{r}}}}}){\boldsymbol{{P}}}_{{{\tilde \theta }_k}}^{{\delta _k}}{{( {{\boldsymbol{C}}^{\text{T}}} \otimes {{\boldsymbol{I}}_{{N_{\rm{r}}}}} )}^{\text{H}}} + {{\boldsymbol{I}}_{{N_{\rm{r}}}{N_{\rm{s}}}}}} ,C和T分别表示c和t的矩阵形式。式(13)是一个无约束最大化问题,其最优解可以由矩阵{\boldsymbol{\varPsi }}({{\boldsymbol{c}}^k},{{\boldsymbol{t}}^k}) 和{{\boldsymbol{\varPsi }}_{{\text{in}}}}({{\boldsymbol{c}}^k},{{\boldsymbol{t}}^k}) 的最大特征值所表示。3.2 发射编码优化
固定
{{\boldsymbol{w}}^{k + 1}} 和{\boldsymbol{t}}_{\rm{r}}^k ,求解{{\boldsymbol{c}}_{\rm{r}}} 。式(11)由F({{\boldsymbol{w}}^{w + 1}},{{\boldsymbol{c}}_{\rm{r}}},{{\boldsymbol{t}}_{\rm{r}}}) 和二次函数项两部分组成。由于F({{\boldsymbol{w}}^{w + 1}},{{\boldsymbol{c}}_{\rm{r}}},{\boldsymbol{t}}_{\rm{r}}^k) 是关于{{\boldsymbol{c}}_{\rm{r}}} 的拟凸函数[34],直接求解问题(12)是具有挑战性的。拟凸函数的上镜图是凸的,因此可以通过最小化F({{\boldsymbol{w}}^{k + 1}},{{\boldsymbol{c}}_{\rm{r}}},{\boldsymbol{t}}_{\rm{r}}^k) 的上界函数来求解{{\boldsymbol{c}}_{\rm{r}}} 。换句话说,可以通过使用BSUM算法来近似求解[30],用{u_{{c}}}({{\boldsymbol{c}}_{\rm{r}}};{{\boldsymbol{w}}^{k + 1}},{\boldsymbol{c}}_{\rm{r}}^k,{\boldsymbol{t}}_{\rm{r}}^k) 来表示F({{\boldsymbol{w}}^{w + 1}},{{\boldsymbol{c}}_{\rm{r}}},{\boldsymbol{t}}_{\rm{r}}^k) 的上界函数,即\begin{split} & {u_c}({{\boldsymbol{c}}_{\rm{r}}};{{\boldsymbol{w}}^{k + 1}},{\boldsymbol{c}}_{\rm{r}}^k,{\boldsymbol{t}}_{\rm{r}}^k) \\ & \quad = F({{\boldsymbol{w}}^{k + 1}},{\boldsymbol{c}}_{\rm{r}}^k,{\boldsymbol{t}}_{\rm{r}}^k) + \nabla _{{{\boldsymbol{c}}_{\rm{r}}}}^{\text{T}}F({{\boldsymbol{w}}^{k + 1}},{{\boldsymbol{c}}_{\rm{r}}},{\boldsymbol{t}}_{\rm{r}}^k){|_{{{\boldsymbol{c}}_{\rm{r}}} = {\boldsymbol{c}}_{\rm{r}}^k}}\\ & \qquad \cdot({{\boldsymbol{c}}_{\rm{r}}} - {\boldsymbol{c}}_{\rm{r}}^k) + \frac{1}{2}{({{\boldsymbol{c}}_{\rm{r}}} - {\boldsymbol{c}}_{\rm{r}}^k)^{\text{T}}}{{{{Z}}}_{{c}}}({{\boldsymbol{c}}_{\rm{r}}} - {\boldsymbol{c}}_{\rm{r}}^k)\\[-15pt] \end{split} (14) 其中
\begin{split} \nabla _{{{\boldsymbol{c}}_{\rm{r}}}}^{}F({{\boldsymbol{w}}^{k + 1}},{{\boldsymbol{c}}_{\rm{r}}},{\boldsymbol{t}}_{\rm{r}}^k) =& \frac{{{{\boldsymbol{\varPi }}_{{\text{inr}}}}({{\boldsymbol{w}}^{k + 1}}){\boldsymbol{t}}_{\rm{r}}^k}}{{{\boldsymbol{c}}_{\rm{r}}^{\text{T}}{{\boldsymbol{\varPi }}_{\rm{r}}}({{\boldsymbol{w}}^{k + 1}}){\boldsymbol{t}}_{\rm{r}}^k}} \\ & - \frac{{F({{\boldsymbol{w}}^{k + 1}},{{\boldsymbol{c}}_{\rm{r}}},{\boldsymbol{t}}_{\rm{r}}^k){{\boldsymbol{\varPi }}_{\rm{r}}}({{\boldsymbol{w}}^{k + 1}}){\boldsymbol{t}}_{\rm{r}}^k}}{{{\boldsymbol{c}}_{\rm{r}}^T{{\boldsymbol{\varPi }}_{\rm{r}}}({{\boldsymbol{w}}^{k + 1}}){\boldsymbol{t}}_{\rm{r}}^k}} \end{split} (15) 从式(15)可以看出
{u_c}({\boldsymbol{c}}_{\rm{r}}^k;{{\boldsymbol{w}}^{k + 1}},{\boldsymbol{c}}_{\rm{r}}^k,{\boldsymbol{t}}_{\rm{r}}^k) = F({{\boldsymbol{w}}^{k + 1}}, {\boldsymbol{c}}_{\rm{r}}^k,{\boldsymbol{t}}_r^k) ,令{\boldsymbol{T}}_c 为F({\boldsymbol{w}}^{k+1},{\boldsymbol{c}}_{\rm{r}},{\boldsymbol{t}}^k_{\rm{r}}) 关于{\boldsymbol{c}}_{\rm{r}} 的海森矩阵,以及{\boldsymbol{Z}}_c=\lambda_{\max}({\boldsymbol{T}}_c){\boldsymbol{I}}_{2N_{\rm{t}} N_{\rm{s}}} , 根据文献[30]可得{\boldsymbol{Z}}_c 和{\boldsymbol{Z}}_c-{\boldsymbol{T}}_c 为半正定矩阵,因此{u_c}({{\boldsymbol{c}}_{\rm{r}}};{{\boldsymbol{w}}^{k + 1}},{\boldsymbol{c}}_{\rm{r}}^k,{\boldsymbol{t}}_{\rm{r}}^k) \ge F({{\boldsymbol{w}}^{k + 1}},{{\boldsymbol{c}}_{\rm{r}}},{\boldsymbol{t}}_{\rm{r}}^k),\forall {{\boldsymbol{c}}_{\rm{r}}} \in {{\cal D}_c} 。基于以上分析,可以通过求解以下问题来更新{{\boldsymbol{c}}_{\rm{r}}} ,即:\begin{split} & \mathop {\min }\limits_{{{\boldsymbol{c}}_{\rm{r}}}} {\text{ }}{u_c}({{\boldsymbol{c}}_{\rm{r}}};{{\boldsymbol{w}}^{k + 1}},{\boldsymbol{c}}_{\rm{r}}^k,{\boldsymbol{t}}_{\rm{r}}^k) + {g_c}({{\boldsymbol{c}}_{\rm{r}}}) \\ & {\text{s}}{\text{.t}}{\text{. }}{\boldsymbol{c}}_{\rm{r}}^{\text{T}}{{\boldsymbol{R}}_{\rm{r}}}{\boldsymbol{t}}_{\rm{r}}^k \le {E_I}, \\ & \quad\;\; {\text{(}}{{\boldsymbol{c}}_{\rm{r}}} - {{\boldsymbol{c}}_{0r}}{)^{\text{T}}}{\boldsymbol{E}}_n^r({\boldsymbol{t}}_{\rm{r}}^k - {{\boldsymbol{c}}_{0{\rm{r}}}}) \le {\varepsilon ^2},n = 1,2,\cdots,{N_{\rm{t}}}{N_{\rm{s}}} \end{split} (16) 其中
\begin{split} {g_c}({{\boldsymbol{c}}_{\rm{r}}}) =& {({{\boldsymbol{u}}^k})^{\text{T}}}({{\boldsymbol{c}}_{\rm{r}}} - {\boldsymbol{t}}_{\rm{r}}^k) + \frac{{{\rho _1}}}{2}||{{\boldsymbol{c}}_{\rm{r}}} - {\boldsymbol{t}}_{\rm{r}}^k|{|^2}\\ & {\text{ + }}{v^k}({\boldsymbol{c}}_{\rm{r}}^{\text{T}}{\boldsymbol{t}}_{\rm{r}}^k - 1) + \frac{{{\rho _2}}}{2}||{\boldsymbol{c}}_{\rm{r}}^{\text{T}}{\boldsymbol{t}}_{\rm{r}}^k - 1|{|^2} \end{split} (17) 优化问题(16)是一个经典二次规划(Quadratic Programming, QP)问题,可以通过MATLAB CVX工具箱直接求解[38]。由于
{{\boldsymbol{t}}_{\rm{r}}} 和{{\boldsymbol{c}}_{\rm{r}}} 的对称性,同样可以通过MATLAB CVX工具箱求解以下最小化问题来更新{{\boldsymbol{t}}_{\rm{r}}} ,即\begin{split} & \mathop {\min }\limits_{{{\boldsymbol{t}}_{\rm{r}}}} {\text{ }}{u_t}({{\boldsymbol{t}}_{\rm{r}}};{{\boldsymbol{w}}^{k + 1}},{\boldsymbol{c}}_{\rm{r}}^{k + 1},{\boldsymbol{t}}_{\rm{r}}^k) + {g_t}({{\boldsymbol{t}}_{\rm{r}}}) \\ & {\text{s}}{\text{.t}}{\text{. }}{\boldsymbol{c}}{_{\rm{r}}^{{k + 1}^{{\text{T}}}}}{{\boldsymbol{R}}_{\rm{r}}}{{\boldsymbol{t}}_{\rm{r}}} \le {E_I}{\text{, }} \\ & \quad\;\; {\text{ (}}{\boldsymbol{c}}_{\rm{r}}^{k + 1} - {{\boldsymbol{c}}_{0{\rm{r}}}}{)^{\text{T}}}{\boldsymbol{E}}_n^r({{\boldsymbol{t}}_{\rm{r}}} - {{\boldsymbol{c}}_{0{\rm{r}}}}) \le {\varepsilon ^2},n = 1,2,\cdots,{N_{\rm{t}}}{N_{\rm{s}}} \end{split} (18) 其中,
{\text{ }}{u_t}({{\boldsymbol{t}}_{\rm{r}}};{{\boldsymbol{w}}^{k + 1}},{\boldsymbol{c}}_{\rm{r}}^{k + 1},{\boldsymbol{t}}_{\rm{r}}^k) 是F({{\boldsymbol{w}}^{k + 1}},{\boldsymbol{c}}_{\rm{r}}^{k + 1},{{\boldsymbol{t}}_{\rm{r}}}) 在{\boldsymbol{t}}_{\rm{r}}^k 处的上界函数,且{g_t}({{\boldsymbol{t}}_{\rm{r}}}) 是二次函数项。为了处理具有非凸二次等式约束和二次不等式约束的问题,引入一个辅助变量进行问题修正。因此,对于每个原始变量,非凸目标函数可以转换为双拟凸函数。更重要的是,将非凸二次等式约束转化为双仿射等式约束。随后,利用块逐次上界最小化法和乘法器交替方向法,设计一种原对偶算法来求解改进后的问题。在该框架下,首先通过ADMM算法消除等式约束,然后通过最小化增广拉格朗日函数的上界函数来更新原变量。由于原始变量的每次更新都涉及一个二次规划问题,该问题的最优解可以用内点法求得。内点法的特点是将构造的无约束目标函数定义在可行域内,并在可行域求解极值点,即求解时探索点始终保持在可行域内。这样,在求解优化问题的过程中,所求得的解总是可行解,从而在可行域内部逐步逼近原约束优化问题的最优解。最后,根据上述步骤,总结了本文提出的ABSUM算法。
算法1 基于ABSUM的发射和接收联合设计算法 Alg. 1 ABSUM algorithm for solving transmit-receive design 输入:k = 0,初始化{\boldsymbol{c} }_{\rm{r}}^k, {\boldsymbol{t} }_{\rm{r}}^k, {{\boldsymbol{u}}^k}, {v^k}和收敛参数 {\epsilon}^{{\rm{abs}}} , {\epsilon}^{{\rm{rel}}} 1:重复 2:k = k + 1 3:通过求{\boldsymbol{\varPsi } }({ {\boldsymbol{c} }^k})和{ {\boldsymbol{\varPsi } }_{{\rm{in}}} }({ {\boldsymbol{c} }^k})的最大广义特征值来更新{{\boldsymbol{w}}^k}。 4:使用内点法求解问题(16)和问题(18)更新{\boldsymbol{c} }_{\rm{r}}^k和{\boldsymbol{t} }_{\rm{r}}^k。 5:求解问题(12)更新{{\boldsymbol{u}}^k}和{v^k} 6:如果满足收敛的终止条件,则算法停止迭代。 4. 数值仿真与分析
本节提供了几组数值模拟实验来评估所提算法的性能。在理想情况下,即目标和干扰的角度是已知的情况,设计具有频谱约束和相似性约束的波形。假设MIMO雷达系统的发射天线数
{N_{\rm{t}}} = 4 ,接收天线数{N_{\rm{r}}} = 8 ,每个脉冲的采样数{N_{\rm{s}}} = 64 。线性调频(Linear Frequency Modulation, LFM)波形在脉冲压缩和模糊度方面有很好的特性,它是区分点目标和成像的好候选者。此外,具有正交发射波形的MIMO雷达在目标检测性能、角度测量能力和动态范围特性等方面优于传统雷达。参考波形{{\boldsymbol{c}}_0} 选取正交LFM信号,其形式为[33]\begin{split} & {{\boldsymbol{C}}_0}(n,l) = \frac{{\exp \{ {\rm{j}}2\pi n(l - 1)/{N_{\rm{s}}}\} \exp \{ {\rm{j}}\pi {{(l - 1)}^2}/{N_{\rm{s}}}\} }}{{\sqrt {{N_{\rm{t}}}{N_{\rm{s}}}} }},\\ & \qquad n = 1,2,\cdots,{N_{\rm{t}}},l = 1,2,\cdots,{N_{\rm{s}}} \\[-10pt] \end{split} (19) 波形向量
{{\boldsymbol{c}}_0} = {\text{vec}}({{\boldsymbol{C}}_0}) 。此外,假设目标角度的均值为{\tilde \theta _0} = {15^ \circ } ,功率为\sigma _{\rm{t}}^2 = 20{\text{ dB}} ,干扰角的均值为{\tilde \theta _1} = - {50^ \circ } ,{\tilde \theta _2} = - {10^ \circ } 和{\tilde \theta _3} = {40^ \circ } ,每个干扰功率是\sigma _i^2 = 30{\text{ dB}},i = 1,2,3 ,干扰白噪声的方差为\sigma _n^2 = 0{\text{ dB}} 。对于干扰,假设它是由未授权的窄带连续干扰器、白噪声干扰和与目标雷达频谱重叠的无线通信网络组成的。具体而言,将扰动协方差矩阵建模为{\boldsymbol{M}} = {\sigma _0}{\boldsymbol{I}} + \sum\limits_{k = 1}^K {\frac{{{\sigma _{I,k}}}}{{\Delta {f_k}}}{\boldsymbol{R}}_{}^k} + \sum\limits_{k = 1}^{{K_J}} {{\sigma _{J,k}}{\boldsymbol{R}}_J^k} (20) 其中,
{\sigma _0} = 0{\text{ dB}} 是热噪声;K = 7 是授权辐射器的数目;{\sigma _{I,k}} 表示在归一化频段\left[ {f_2^k,f_1^k} \right] ({\sigma _{I,k}} = 10{\text{ dB}}, k = 1,2,\cdots,K )上第k个共存无线网络的能量;\Delta {f_k} = f_2^k - f_1^k,k = 1,2,\cdots,K 是第k个授权辐射器的相关带宽;{K_J} = 2 表示未授权窄带干扰器的数量;{\sigma _{J,k}}, k = 1,2,\cdots,{K_J} 是第k个主动干扰机的能量({\sigma _{J,1}} = 50{\text{ dB}}, {\sigma _{J,2}} = 40{\text{ dB}} );{\boldsymbol{R}}_J^k,k = 1,2,\cdots,{K_J} 是第k个有源未授权干扰机的归一化干扰协方差矩阵,其定义为{\boldsymbol{R}}_J^k = {{\boldsymbol{r}}_{J,k}}{\boldsymbol{r}}_{J,k}^\dagger ,其中{r_{J,k}}\left( n \right) = {{\rm{e}}^{{\rm{j}}2\pi {f_{J,k}}n/{f_{\rm{s}}}}} ,{f_{J,k}} 表示第k个干扰器的多普勒频移({f_{J,1}}/{f_{\rm{s}}} = 0.7 和{f_{J,2}}/{f_{\rm{s}}} = 0.75 )。此外,假设授权的共存通信系统具有相同的相关性,即{\omega _i} = 1 时,i = 1,2,\cdots,7 。基于假设的二个阻带[f_{{\text{lower}}}^{\text{1}},f_{{\text{upper}}}^{\text{1}}] = [0.2,0.3] 和[f_{{\text{lower}}}^{\text{2}},f_{{\text{upper}}}^{\text{2}}] = [0.75,0.85] 以及权值{\omega _i} ,计算矩阵{\boldsymbol{R}}_{}^k (定义在式(5)),并对发射雷达波形施加通信干扰能量约束。关于所提出算法的参数,设{\boldsymbol{t}}_{\rm{r}}^0 = {\boldsymbol{c}}_{\rm{r}}^0 ,{{\boldsymbol{u}}^0} = 0 和{v^0} = 0 ,惩罚参数{\rho _1} 和{\rho _2} 被设置为{\rho _1} = {\rho _2} = 1 。所有数值示例均使用MATLAB R2016a版本进行运算,并在标准PC上执行(CPU核i5 4.4 GHz和8 GB RAM)。表1列出了主要的仿真参数。表 1 仿真参数Table 1. Simulation parameter参数 数值 参数 数值 发射天线{N_{\rm{t}}} 4 接收天线{N_{\rm{r}}} 8 采样个数{N_{\rm{s}}} 64 干扰源角度{\theta _k} –50°, –10°, 40° 雷达目标角度{\theta _0} {\text{1}}{5^ \circ } 无线网络归一化频率f_{{\rm{lower}}}^i, f_{{\rm{upper}}}^i [0.2,0.3]
[0.75,0.85]最大干扰量{E_I} {10}^{-2}, {10}^{-\text{3} }, {10}^{-\text{5} } 相似性参数\varepsilon 0.5,0.8,1.0,2.0 首先,在确定性情况下设计具有能量、相似性和频谱约束的雷达发射波形,即
{\theta _0} = {\tilde \theta _0} ,{\theta _k} = {\tilde \theta _k}, k = 1,2,3 。在频谱约束{E_I} = {10^{ - 4}} 和相似性约束\varepsilon = 1/\sqrt {{N_{\rm{t}}}{N_{\rm{s}}}} 的条件下,为了简化符号,定义{\varepsilon '} = \varepsilon \sqrt {{N_{\rm{t}}}{N_{\rm{s}}}} 。图1(a)显示了不同相似性参数下ABSUM算法得到的波形相位,图1(b)显示不同约束条件下ABSUM算法得到的波形幅值。4.1 收敛性分析
由于QCQP[20], MM[22], BCD[21]和SQR方法[8]不能解决本文所考虑的问题,因此,在接下来的分析中,我们将所提出的ABSUM算法与GFA[27]和SOA方法[19]进行比较,首先对3种优化波形的SINR迭代曲线进行比较,相比于GFA和SOA算法,ABSUM算法利用拉格朗日乘子将原问题分解成多个可求解的子问题,然后进行交替更新求解,使得优化问题能够获得更好的解。设频谱约束值为
{E_I} = {10^{ - 4}} 时,图2描绘了不同相似性参数下输出SINR值与迭代次数的关系。从图中可以看出ABSUM, GFA和SOA算法的SINR值都随相似性约束的增加而增加,通过增加相似度,即增加可用的自由度,可以实现更好的有用能量分布。此外,在{\varepsilon '} = 0.5, 1.0,1.3 时,ABSUM算法均显著优于GFA和SOA算法,具体而言,相似性约束值为{\varepsilon '} = 1.0 时,所提出的算法与对比算法的增益约为2.3 dB和3.1 dB。假设频谱约束值为
{E_I} = {10^{ - 4}} ,表2分析了ABSUM, GFA和SOA算法在不同相似性约束条件下的SINR性能。结果表明,ABSUM算法在SINR值方面优于GFA和SOA算法。此外,ABSUM算法的全局计算时间明显小于GFA和SOA算法的全局计算时间。表 2 ABSUM, GFA和SOA算法的优化SINR值(dB)和全局计算时间(s)Table 2. SINR value (dB) and global computational times (s) for ABSUM, GFA and SOA{\varepsilon'} ABSUM GFA SOA SINR Time SINR Time SINR Time 0.5 17.2 15.3 15.4 135.8 13.6 637.9 1.0 18.8 19.1 16.5 141.5 14.8 853.1 1.3 19.2 21.5 17.5 154.4 16.5 963.1 2.0 19.3 19.8 19.3 169.3 19.3 812.1 4.2 能量和相似性参数的影响
图3显示了当
{E_I} = {10^{ - 5}},{10^{ - 4}},{10^{ - 2}} 时,ABSUM, GFA和SOA算法的优化SINR值与相似性参数{\varepsilon'} 的关系,结果再次表明ABSUM, GFA和SOA算法的优化SINR值随着相似度值的增加而增加,对于相同的频谱约束值,ABSUM可以获得比GFA和SOA更高的SINR值。此外,可以观察到,当{\varepsilon '} 达到一定值时SINR不再变化,说明当相似性参数超过一定值时,相似性约束将消失。假设相似度参数为
{\varepsilon '} = 0.5,1.0,2.0 时,ABSUM, GFA和SOA算法的优化SINR值与频谱约束值之间的关系如图4所示,可以看出3种方法的优化SINR值随着{E_I} 的增加而增加,如所预料的那样,在相同的相似性约束条件下,ABSUM相比于GFA和SOA能获得更好的SINR性能。此外,可以看出当相似性参数{\varepsilon '} 越大时,频谱约束对SINR的影响越小。特别是,当{\varepsilon '} = 20 时,频谱约束对SINR几乎没有影响。比较图3和图4可以发现,相似性约束对于SINR的影响要明显比频谱约束大得多。因此,设计现代雷达系统时,应该充分考虑相似性约束和频谱约束的权衡。图5显示了ABSUM算法和GFA算法的优化SINR值与不同频谱约束
{E_I} 和相似性约束{\varepsilon'} 的关系。正如预期那样,对于固定的相似性约束和频谱约束,ABSUM具有更好的SINR。另外,两种算法的输出SINR均与输入SNR呈线性关系。这意味着输入SNR对输出SINR损耗没有影响。接下来,讨论目标角
{\theta _0} 和干扰源角度{\theta _k},k = 1, 2,\cdots,K 不确定时SINR值的变化情况,首先假设{\delta _k} = 0 ,图6(a)表明当目标角度{\theta _0} 的不确定性越大时,也就是{\delta _t} 越大时,将导致更大的SINR损耗。同样假设{\delta _t} = 0 ,由图6(b)可得,随着{\delta _k} 的增大SINR的损耗也越大。从图6可以看出,对于不同的{E_I} 和{\varepsilon '} ,所提出的ABSUM算法比GFA算法具有更好的鲁棒性。4.3 能量谱密度分析
关于频谱兼容约束,本文考虑两个共存的频率区间,第1个频率区间为
[f_{{\text{lower}}}^{\text{1}},f_{{\text{upper}}}^{\text{1}}] = [0.2,0.3] ,第2个频率区间为[f_{{\text{lower}}}^{\text{2}},f_{{\text{upper}}}^{\text{2}}] = [0.75,0.85] 。首先假设相似性参数为{\varepsilon '} = 1.0 ,图7提供了优化波形的能量谱密度(Energy Spectral Densities, ESD)和归一化频率的关系图。图7表明对应3个不同的{E_I} 值所优化的波形均在限制频段内产生了一定程度的频谱凹陷。限制频段内的波形功率谱主要受频谱约束{E_I} 的影响,当{E_I} = {10^{ - 4}},{10^{ - 3}} 时,所设计的波形在2个限制频段内均形成较深的频谱凹陷,具有良好的频谱兼容性。对于较小的{E_I} 值,频谱凹陷更深,因为减少{E_I} 意味着在相应的阻带中传输更少的能量。但当{E_I} 增大到{10^{ - 2}} 时,限制频段内的频谱凹陷程度显著减小,这是因为频谱约束实质上是约束限制频段上的总能量,当频谱约束放松时,即{E_I} 值增加,所设计的波形仅需满足所有限制频段内的总能量不超过{E_I} 即可,在这种条件下可能会出现部分频段频谱凹陷消失的现象,但波形能通过更多的自由度来提升SINR值。图8分析了在给定频谱兼容约束
{E_I} = {10^{ - 4}} 的条件下,不同相似性参数{\varepsilon '} 优化波形的功率谱情况。对比图7可以发现,限制频段内的频谱凹陷变化情况受相似性约束的影响较小。相似性参数主要影响非限制频段区的功率谱形状,当相似性参数越小时,波形功率谱在非限制频段更加接近参考信号功率谱。最后,由图7和图8可得,限制频段内的功率谱分布主要受参数{E_I} 的影响,限制频段外的功率谱分布主要受参数{\varepsilon '} 的影响。给定相似性参数
{\varepsilon '} = 1.0 下,表3对比了不同算法在不同频谱约束E_I^{} 的SINR值和全局计算时间。结果再次证实,所提出的ABSUM算法相比于SOA算法和GFA算法的SINR值有明显的提升,特别是当{E_I} = {10^{ - 2}} 时,相比于SOA和GFA的增益分别为2.5 dB和2.4 dB。由于ABSUM算法是将原问题分解为多个子问题并行求解,能够更快地收敛到最优解,从表3中也可以看出全局计算时间明显小于SOA和GFA。相比于MM算法,全局计算时间略慢一些,在SINR增益方面有所提升,当{E_I} = {10^{ - 2}} 时,提升约为0.1 dB。表 3 ABSUM、BCD、MM和GFA的SINR值(dB)和全局计算时间(s)Table 3. SINR value (dB) and global computational times (s) for ABSUM, BCD, MM and GFA{E_I} ABSUM SOA MM GFA SINR time SINR time SINR time SINR time {10^{ - 2}} 19.2 14.3 16.7 45.3 19.1 10.2 16.8 120.8 {10^{ - 3}} 19.0 16.2 16.8 55.2 18.9 13.4 16.7 134.5 {10^{ - 4}} 18.8 19.1 16.8 66.8 18.6 15.7 16.5 141.5 4.4 波束图和模糊度函数分析
假设频谱兼容约束为
{E_I} = {10^{ - 4}} ,图9比较了不同相似值{\varepsilon '} = 0.5,1.0,2.0 下的ABSUM算法和GFA算法的波束图。波束图P\left( \theta \right) 的计算公式表示为P\left( \theta \right) = |{{{\dot {\boldsymbol{w}}}}^{\rm{H}}}{\boldsymbol{H}}(\theta ){\boldsymbol{c}}{|^2} ,其中{{\dot {\boldsymbol{w}}}} 是归一化滤波权重,即{{\dot {\boldsymbol{w}}}} = {{\boldsymbol{w}}}/{{\left\| {\boldsymbol{w}} \right\|}} 。从图9可以看出,对于所有参数{\varepsilon '} = 0.5, 1.0,1.2 ,本文提出的ABSUM算和GFA算法都具有良好的干扰抑制性能。然而,基于ABSUM算法波束图的主瓣增益要优于GFA算法。正如预期的那样,当{\varepsilon '} = 2.0 时,ABSUM算法和GFA算法表现出几乎相同的主瓣增益。最后,本文利用互模糊函数(Cross Ambiguity Function, CAF)[19,39]来分析已优化雷达系统的联合距离和方位特征。图10(a)—图10(d)为不同频谱和相似性约束下CAF的距离-方位切割等高线图。假设频谱兼容约束为
{E_I} = {10^{ - 4}} ,图10(a)和图10(b)分别为算法1在迭代次数为0时(仅为接收机设计)和收敛时的模糊函数特征。距离-方位角图表明,与仅自适应接收不同,利用算法1优化的发射-接收对在收敛时,可以在黑色矩形中发现一个明显的峰值,它代表感兴趣的雷达目标。假设相似性约束为{\varepsilon '} = 0.5 ,图10(c)和图10(d)分别显示了算法1在迭代次数为0时和收敛时的模糊函数特征。可以观察到相同的结果。这一现象表明,在频谱和相似性约束下,ABSUM算法都产生了合适的CAF,从而抑制了信号相关的杂波干扰。5. 结语
针对MIMO雷达系统和周边通信服务网络频谱兼容的问题,本文以发射波形能量、相似性和频谱兼容为约束条件,在信号相关杂波背景下,建立了MIMO雷达发射波形和接收滤波器的联合设计模型。为求解有非凸二次等式约束和二次不等式约束的优化问题,在ADMM框架下设计了一种基于BSUM算法的迭代原对偶算法,所提的方法结合了ADMM算法和BSUM算法的优点,通过数值仿真验证了该方法在输出SINR、频谱特征、波束图、计算复杂度和模糊特性等方面的性能。结果表明:(1)该算法以较低的计算代价实现了优于GFA算法的输出SINR。频谱兼容参数越小,输出SINR越小,相似度参数越小,输出SINR越小;此外,限制频段内的功率谱分布主要受频谱兼容参数的影响,限制频段外的功率谱分布主要受相似度参数的影响。因此,在实践中可以在输出SINR、相似性值和频谱兼容值之间进行适当的权衡。(2)与GFA算法相比,该算法对目标角度的不确定性和干扰具有更好的鲁棒性。此外,不确定性越大,导致的SINR损失越大。(3)在频谱和相似性约束下,该算法产生了合适的CAF,从而抑制了信号相关杂波的干扰。
-
表 1 仿真实验1参数
Table 1. Experimental parameters for the first simulation data
参数 数值 参数 数值 通道数 11 波长 0.02 m 基线跨度 1 m 下视角 45° 最小基线间隔 0.1 m 高程模糊间隔 100 m 斜距 1000 m 高程瑞利分辨率 10 m 表 2 仿真实验2参数
Table 2. Experimental parameters for the second simulation data
参数 数值 参数 数值 航过数 11 波长 0.03 m 基线跨度 450 m 下视角 45° 最小基线间隔 21 m 最大不模糊高程 428.6 m 斜距 600 km 高程瑞利分辨率 20 m 表 3 峨眉数据参数
Table 3. Parameters of Emei data
参数 数值 载波频率 14.5 GHz 最小基线间隔 0.1115 m 最大基线长度 1.1274 m 载机航线海拔 2157 m 场景海拔 420 m 中心斜距 2040.1 m 中心下视角 31.6° 距离向像素尺寸 0.1362 m 方位向像素尺寸 0.1051 m 表 4 峨眉数据的平均邻域高度差(m)
Table 4.
\Delta {{\boldsymbol{h}}_{\bf{E}}} of Emei data (m)未采用自适应高程搜索范围\Delta {h_{\text{E}}} 采用自适应高程搜索范围\Delta {h_{\text{E}}} 12.6512 7.5453 表 5 巴塞罗那数据参数
Table 5. Parameters of Barcelona data
参数 数值 载波频率 9.65 GHz 最小基线间隔 7.98 m 最大基线长度 246.4 m 中心斜距 621.6 km 中心下视角 35.7° 距离向像素尺寸 0.91 m 方位向像素尺寸 1.88 m 表 6 巴塞罗那数据的平均邻域高度差(m)
Table 6.
\Delta {{\boldsymbol{h}}_{\mathbf{E}}} of Barcelona data (m)未采用自适应高程搜索范围\Delta {h_{\text{E}}} 采用自适应高程搜索范围\Delta {h_{\text{E}}} 22.5561 14.3084 -
[1] MUNSON D C, O’BRIEN J D, and JENKINS W K. A tomographic formulation of spotlight-mode synthetic aperture radar[J]. Proceedings of the IEEE, 1983, 71(8): 917–925. doi: 10.1109/PROC.1983.12698 [2] 张斌, 韦立登, 胡庆荣, 等. 基于四阶累积量的机载多基线SAR谱估计解叠掩方法[J]. 雷达学报, 2018, 7(6): 740–749. doi: 10.12000/JR18087ZHANG Bin, WEI Lideng, HU Qingrong, et al. Solution to layover problemin airborne multi-baseline SAR based on spectrum estimation with fourth-order cumulant[J]. Journal of Radars, 2018, 7(6): 740–749. doi: 10.12000/JR18087 [3] 丁赤飚, 仇晓兰, 徐丰, 等. 合成孔径雷达三维成像——从层析、阵列到微波视觉[J]. 雷达学报, 2019, 8(6): 693–709. doi: 10.12000/JR19090DING 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 [4] FORNARO G, SERAFINO F, and SOLDOVIERI F. Three-dimensional focusing with multipass SAR data[J]. IEEE Transactions on Geoscience and Remote Sensing, 2003, 41(3): 507–517. doi: 10.1109/TGRS.2003.809934 [5] 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 [6] FORNARO G, LOMBARDINI F, PAUCIULLO A, et al. Tomographic processing of interferometric SAR data: Developments, applications, and future research perspectives[J]. IEEE Signal Processing Magazine, 2014, 31(4): 41–50. doi: 10.1109/MSP.2014.2312073 [7] 张红, 江凯, 王超, 等. SAR层析技术的研究与应用[J]. 遥感技术与应用, 2010, 25(2): 282–287. doi: 10.11873/j.issn.1004-0323.2010.2.282ZHANG Hong, JIANG Kai, WANG Chao, et al. The current status of SAR tomography[J]. Remote Sensing Technology and Application, 2010, 25(2): 282–287. doi: 10.11873/j.issn.1004-0323.2010.2.282 [8] EL MOUSSAWI I, MINH D H T, BAGHDADI N, et al. L-band UAVSAR tomographic imaging in dense forests: Gabon forests[J]. Remote Sensing, 2019, 11(5): 475. doi: 10.3390/rs11050475 [9] EL MOUSSAWI I, MINH D H T, BAGHDADI N, et al. Monitoring tropical forest structure using SAR tomography at L- and P-band[J]. Remote Sensing, 2019, 11(16): 1934. doi: 10.3390/rs11161934 [10] PIAU P, BRUNIQUEL J, CAEL J C, et al. Analysis of the resolution of a multitemporal SAR System[C]. IEEE International Geoscience and Remote Sensing Symposium, Tokyo, Japan, 1993: 1196–1199. [11] ZHU Xiaoxiang and BAMLER R. Tomographic SAR inversion by L1 norm regularization—the compressive sensing approach[J]. IEEE Transactions on Geoscience and Remote Sensing, 2010, 48(10): 3839–3846. doi: 10.1109/TGRS.2010.2048117 [12] 魏恋欢, 廖明生, BALZ T, 等. 高分辨率SAR层析成像建筑物叠掩散射体提取[J]. 武汉大学学报: 信息科学版, 2014, 39(5): 536–540. doi: 10.13203/j.whugis20120460WEI Lianhuan, LIAO Mingsheng, BALZ T, et al. Layover building scatterers extraction via high-resolution spaceborne SAR tomography[J]. Geomatics and Information Science of Wuhan University, 2014, 39(5): 536–540. doi: 10.13203/j.whugis20120460 [13] 毕辉, 金双, 王潇, 等. 基于高分三号SAR数据的城市建筑高分辨率高维成像[J]. 雷达学报, 2022, 11(1): 40–51. doi: 10.12000/JR21113BI Hui, JIN Shuang, WANG Xiao, et al. High-resolution high-dimensional imaging of urban building based on GaoFen-3 SAR data[J]. Journal of Radars, 2022, 11(1): 40–51. doi: 10.12000/JR21113 [14] LOMBARDINI F and REIGBER A. Adaptive spectral estimation for multibaseline SAR tomography with airborne L-band data[C]. 2003 IEEE International Geoscience and Remote Sensing Symposium, Toulouse, France, 2003: 2014–2016. [15] SCHMIDT R. Multiple emitter location and signal parameter estimation[J]. IEEE Transactions on Antennas and Propagation, 1986, 34(3): 276–280. doi: 10.1109/TAP.1986.1143830 [16] 廖明生, 魏恋欢, 汪紫芸, 等. 压缩感知在城区高分辨率SAR层析成像中的应用[J]. 雷达学报, 2015, 4(2): 123–129. doi: 10.12000/JR15031LIAO Mingsheng, WEI Lianhuan, WANG Ziyun, et al. Compressive sensing in high-resolution 3D SAR tomography of urban scenarios[J]. Journal of Radars, 2015, 4(2): 123–129. doi: 10.12000/JR15031 [17] 仇晓兰, 焦泽坤, 彭凌霄, 等. SARMV3D-1.0: SAR微波视觉三维成像数据集[J]. 雷达学报, 2021, 10(4): 485–498. doi: 10.12000/JR21112QIU Xiaolan, JIAO Zekun, PENG Lingxiao, et al. SARMV3D-1.0: Synthetic aperture radar microwave vision 3D imaging dataset[J]. Journal of Radars, 2021, 10(4): 485–498. doi: 10.12000/JR21112 [18] 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 [19] 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 [20] ZHANG Bangjie, XU Gang, YU Hangwen, et al. Array 3-D SAR tomography using robust gridless compressed sensing[J]. IEEE Transactions on Geoscience and Remote Sensing, 2023, 61: 5205013. doi: 10.1109/TGRS.2023.3259980 [21] 林赟, 张琳, 韦立登, 等. 无先验模型复杂结构设施SAR全方位三维成像方法研究[J]. 雷达学报, 2022, 11(5): 909–919. doi: 10.12000/JR22148LIN Yun, ZHANG Lin, WEI Lideng, et al. Research on full-aspect three-dimensional SAR imaging method for complex structural facilities without prior model[J]. Journal of Radars, 2022, 11(5): 909–919. doi: 10.12000/JR22148 [22] REN Yexian, XIAO Aoran, HU Fengming, et al. Coprime sensing for airborne array interferometric SAR tomography[J]. IEEE Transactions on Geoscience and Remote Sensing, 2022, 60: 5229615. doi: 10.1109/TGRS.2022.3182980 [23] LI Xiaowan, ZHANG Fubo, LI Yanlei, et al. An elevation ambiguity resolution method based on segmentation and reorganization of TomoSAR point cloud in 3D mountain reconstruction[J]. Remote Sensing, 2021, 13(24): 5118. doi: 10.3390/rs13245118 [24] 仇晓兰, 焦泽坤, 杨振礼, 等. 微波视觉三维SAR关键技术及实验系统初步进展[J]. 雷达学报, 2022, 11(1): 1–19. doi: 10.12000/JR22027QIU 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 [25] ZHU Xiaoxiang and BAMLER R. Sparse reconstruction techniques for SAR tomography[C]. 17th International Conference on Digital Signal Processing, Corfu, Greece, 2011: 1–8. [26] 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 [27] WANG Xiao and XU Feng. Tomographic SAR inversion by atomic-norm minimization—the gridless compressive sensing approach[J]. IEEE Transactions on Geoscience and Remote Sensing, 2022, 60: 5239113. doi: 10.1109/TGRS.2022.3223524 [28] DONOHO D L. Compressed sensing[J]. IEEE Transactions on Information Theory, 2006, 52(4): 1289–1306. doi: 10.1109/TIT.2006.871582 [29] CHEN S S, DONOHO D L, and SAUNDERS M A. Atomic decomposition by basis pursuit[J]. SIAM Review, 2001, 43(1): 129–159. doi: 10.1137/S003614450037906X [30] 王金峰, 皮亦鸣, 曹宗杰. 一种机载SAR层析三维成像算法[J]. 电子与信息学报, 2010, 32(5): 1029–1033. doi: 10.3724/SP.J.1146.2009.00737WANG Jinfeng, PI Yiming, and CAO Zongjie. An algorithm for airborne SAR tomography 3D imaging[J]. Journal of Electronics &Information Technology, 2010, 32(5): 1029–1033. doi: 10.3724/SP.J.1146.2009.00737 [31] GUO Rui, WANG Fan, ZANG Bo, et al. High-rise building 3D reconstruction with the wrapped interferometric phase[J]. Sensors, 2019, 19(6): 1439. doi: 10.3390/s19061439 [32] GUO Rui, GAO Yuxin, ZHANG Zhao, et al. Efficient tomographic inversion based on refined scatterer pre-estimation[J]. IEEE Geoscience and Remote Sensing Letters, 2022, 19: 4513305. doi: 10.1109/LGRS.2022.3203037 期刊类型引用(1)
1. 王健,时晨光,周建江,汪飞. 频谱共存下基于OFDM-LFM的机载组网雷达射频隐身波形优化算法. 战术导弹技术. 2024(05): 111-121 . 百度学术
其他类型引用(0)
-