为了解决上述问题，文献［4］提出了通过互质采样得到的数据进行DOA估计的方法，与传统的使用均匀线性采样相比，若使用含有N个阵元的互质阵列进行数据采样，由其构成的差联合阵中的虚拟均匀线阵的大小达到了[ON2]，因此其自由度为[ON2]。近年来，不少学者已针对使用互质阵列采样数据进行DOA估计做了大量且深入的研究。文献［5］中给出了一种将互质阵列与普通的MUSIC算法相结合来估计DOA的方法。它首先对两个阵元数互质的子阵列进行DOA估计，然后将这两个估计值进行对比，给出了DOA估计值存在的唯一性的证明。文献［6］中引入了空间平滑思想来对互质阵列得到的数据进行DOA估计。该算法首先利用空间平滑的思想对互质阵列接收到的数据的协方差矩阵进行重构，由此得到的矩阵的维数相比于实际物理阵列得到的矩阵更大，从而提高了阵列的孔径和自由度。文献［7］提出了一种将协方差矩阵重构、补全的DOA估计算法。该算法首先对协方差矩阵重构，一方面抑制了非均匀噪声对信号的影响，另一方面有效地提高了孔径和自由度，然后使用奇异值阈值算法对重构造成的空洞进行补全，最后使用DOA估计算法进行角度估计。
在理想情况下，信号中的噪声通常被假定为不相关的高斯白噪声，一些基于互质线性阵列的DOA估计算法通过对实际阵列的数据协方差矩阵进行重构操作，然后对其进行特征值分解，从而得到信号和噪声向量，然后利用两者正交求倒寻谱峰的方法进行DOA估计，具有良好的性能。但是，在实际应用中，环境噪声的空间分布通常是不均匀的，若使用传统的基于特征分解的算法时，会导致子空间泄露。考虑到非均匀噪声可能与任何两个不同的传感器相关，在当前的工作中使用线性噪声模型。在这个参数化的线性模型［8-9］中，利用截断傅里叶展开来得到其功率密度谱函数，从而建立非均匀噪声场中的信号处理模型。为了处理DOA估计中的非均匀噪声，一些学者提出了基于压缩感知的低秩矩阵恢复理论［10-12］。文献［13］证明使用半正定规划可以完全恢复低秩矩阵。文献［14］提出了一种基于最小化截断核范数（truncated nuclear norm regularization，TNNR）的矩阵补全算法。文献［15］提出了一种利用阵列协方差矩阵稀疏性的不精确增广拉格朗日乘子（inexact alternating direction multiplier， IALM）算法。文献［16］结合Lp范数和截断核范数重构优化问题模型，再使用交替方向乘子法（alternating direction multiplier method， ADMM）算法和增广拉格朗日乘子法（augmented lagrangian multipliers， ALM）算法进行问题模型求解，从而恢复低秩矩阵。文献［17］提出了一种交替投影（alternating projection algorithm，APA）算法来恢复低秩矩阵。
为了实现互质阵列在色噪声条件下的DOA估计，本文提出了一种基于协方差矩阵补全-重构-再补全的估计方法。本文所使用的噪声模型相比于高斯白噪声和普通的非均匀噪声更为复杂，更近似于真实环境。为了更好地估计DOA，本文将多种矩阵补全算法进行对比，最后得出本文所使用组合算法的最优方法，其在运算精度和效率上相比于其它补全算法有所提高。
1 信号模型
如图1所示，该互质线性阵列由阵元数分别为N和2M的两个均匀线性子阵列组成，第一个均匀子阵列的阵元间距为Md，第二个子阵列的阵元间距为Nd，且两个子阵列的第一个阵元位置重合。其中，整数M和N互质，[d=λ/2]，[λ]为入射信号波长。若以原点为参考点，各阵元位置可记为[C=]
[nMd|0≤n≤N-1?mNd|0≤m≤2M-1]。
<G:\武汉工程大学\2023\第1期\宋 鹏-1.tif>
图1 互质阵列示意图：（a）构造互质线阵的两个均匀子阵列， （b）两个均匀子阵列组成的新阵列
Fig. 1 Schematic diagrams of coprime array：
（a） constructing two uniform subarrays of coprime linear
array， （b） new array composed of two uniform subarrays
当有L个远场源从不同的方向发射信号到达互质阵列时，其到达角可以表示为[θ=(θ1,θ2,][...,θL)]，其中，[θl]表示第[l]个声源相对于阵列质心的到达角，则由互质阵列接收到的数据[xt]可以表示为
[xt=l=1Lαθlslt+nt=Ast+nt] （1）
其中，[A=αθ1,αθ2,...,αθL]为流型矩阵，且[αθl=1,ej2πc2λsinθl,...,ej2πc2M+N-1λsinθlT]为导向矢量，[ci∈C,i=1,2,...,2M+N-1]为互质阵列中各阵元的真实位置。[st=s1t,s2t,...,sLtT]为入射信号矢量。[nt=n1t,n2t,...,n2M+N-1tT]为非均匀白高斯噪声。互质阵列接收到的信号的数据协方差矩阵可表示为
[R=ExtxHt=AEstsHtAH+]
[Q=ARssAH+Q] （2）
其中，[E?]用来求数据的期望，[?H]表示求数据的共轭转置。[Q]由空间不相关的非均匀噪声组成
[Q=EntnHt=diagq1,1,q2,2,...,q2M+N-1, 2M+N-1] （3）
将环境噪声信息用线性噪声模型表示，其噪声协方差矩阵可以表示为
[Q=02παθfθαHθdθ] （4）
其中，[fθ]是周期为[2π]的噪声功率密度函数（power density function，PDF）。因此，它可以用傅里叶级数来表示
[fθ=l=0∞alcoslθ+l=1∞blsinlθ] （5）
其中，
[al=12π02πfθdθ, l=01π02πfθcoslθdθ, l>0] （6）
[bl=1π02πfθsinlθdθ] （7）
当[fθ]足够平滑时，高阶项系数随[θ]增大而减小。将傅里叶级数截断为[K-1]项，可以得到一个有限线性模型
[Q=μI2M+N-1+l=0K-2αlQl] （8）
其中，[μ]是白噪声的功率，它是一个表示未知噪声的傅里叶级数的实值。[Ql]表示噪声的相关矩阵，它满足以下性质
[Ql=02παθcoslθαHθdθ, l是偶数02παθsinlθαHθdθ, l是奇数] （9）
因此，[Q]可以定义为
[Q=q1,1q1,20?00q2,1q2,2q2,3?000q3,2q3,3?00??????000?q2M+N-2,2M+N-2q2M+N-2,2M+N-1000?q2M+N-1,2M+N-2q2M+N-1,2M+N-1]
（10）
在实际应用中，协方差矩阵[R]可以用N个快拍数近似计算，即
[R=1Nt=1NxtxHt] （11）
2 色噪声条件下的互质阵列DOA估计
2.1 数据协方差矩阵去噪
矩阵补全是基于压缩感知理论的算法，通常应用于低秩矩阵补全。从式（10）可以看出，噪声主要对数据协方差矩阵的3条主对角线元素产生影响。因此，本文将含噪信号从数据协方差矩阵中去除，然后使用Lp-TNN低秩补全算法对此矩阵进行恢复，从而得到去噪处理后的数据协方差矩阵。
信号协方差矩阵定义为[R]，观测矩阵为[G]，观测项对应的位置集合为[Ω]，设[PΩ]为正交投影算子，则
[PΩGij=Gij,i, j∈Ω0,其它] （12）
其中，[Gij]表示[G]中的第[i, j]个元素。
很明显，[R]是一个低秩矩阵，矩阵补全的过程可以描述为以下秩最小化问题，即
[minRrankR s.t. PΩR=PΩG] （13）
其中，[rank?]为求秩函数。显然，矩阵秩的优化是一个NP-hard的问题。然而，可以适当地进行凸松弛，使用矩阵核范数来近似秩函数，式（13）可以表示为
[minRR? s.t. RΩ=GΩ] （14）
其中，[RΩ=PΩR, GΩ=PΩG]。
为了更好地匹配观察项的对应项，可以适当放松等式的约束，式（14）可以变换为
[minRR?+γRΩ-GΩ22] （15）
其中，[R?=i=1minm,nσiR]表示核范数，[γ]为惩罚参数。
由于Lp范数可以被看作L2范数的扩展，因此，此处可以使用Lp范数[0<p≤1]来代替L2范数，以便增强对异常值处理的稳健性。同时，由于矩阵的最大的[r]个非零奇异值的值不会影响到该矩阵的秩，所以可以将注意力放在最小化除[r]之外剩余的最小奇异值。因此，式（15）可以变换为
[minRRr+γRΩ-GΩpp] （19）
显然，求解[Rr=i=r+1minm,nσiR]是非凸的问题，很难直接进行求解。为了解决这一问题，式（16）可以变换为
[minRR?-maxAAT=BBT=Ir×rTrARBT+γRΩ-GΩpp]
（17）
其中，[A]、[B]分别是实数域里的[r×m]、[r×n]矩阵。
对于式（17），其拉格朗日函数可以表示为
[LR, μ, Δ, ∑, W, EΩ=R?-TrAWBT+γEΩpp+μ2EΩ-RΩ-GΩ+1μΔ2F+μ2R-W+1μ∑2F] （18）
其中，[W=R, EΩ=RΩ-GΩ]，[Δ]和[Σ]表示比例对偶变量。
使用ADMM算法和ALM算法对计算进行简化，则有以下步骤
步骤 1 对[Rk]进行奇异值分解（SVD）：
[Uk,∑kVk=svdRk] （19）
其中，[Uk=uk1,...,ukm∈Rm×m,][Vk=uk1,...,ukn∈][Rn×n]。
步骤 2 更新[Ak+1]和[Bk+1]：
[Ak+1=uk1,...,ukrT, Bk+1=vk1,...,vkrT] （20）
步骤 3 使用ADMM简化运算，更新[Rk+1]：
[Rk+1=argminRRk?+]
[μ2EkΩ-RkΩ-GkΩ+1μΔk2F+]
[μ2Rk-Wk+1μ∑k2F] （21）
步骤 4 更新[Ek+1Ω]：
[Ek+1Ω=argminEΩγEkΩpp+]
[μ2EkΩ-RkΩ-GkΩ+1μΔk2F] （22）
步骤 5 更新[Wk+1]：
[Wk+1=argminW-TrAk+1WkBTk+1+] [μ2Rk+1-Wk+1μ∑k2F] （23）
步骤 6 使用ALM算法来计算，更新[Δk+1]：
[Δk+1=Δk+μEk+1Ω-Rk+1+GΩ] （24）
步骤 7 更新[Σk+1]：
[∑k+1=∑kRk+1-Wk+1] （25）
步骤 8 更新[μ]：
[μ=ρμ] （26）
其中，[ρ>0]。
在循环迭代时，若满足以下条件，迭代停止，输出补全后的矩阵。
[Rk+1-RkFmaxRk+1F,1≤ε or iteration≥n] （27）
其中，[ε]是停止迭代阈值，iteration是迭代次数，[n]为最大迭代次数。
2.2 协方差矩阵重构
基于子空间DOA估计算法要求阵元间距相同，显然在互质阵列的条件下无法满足，因此，需要使用矩阵重构来使其满足这一条件。首先引入差分阵列的概念，集合[Z]可以表述为
[Z=pi-pj, pi, pj∈C] （28）
其中，[C]为互质阵列中各阵元的实际位置信息集合，[pi]和[pj]分别表示第[i]和第[j]个实际阵元位置。[Z]为[C]中任意的两个元素的差所构成的集合，显然[Z]中有多个相同的元素值，将[Z]中不相同的元素[px]构成的集合定义为[Zx]，将[px]在[Z]中出现的次数记为[ωpx]。
由式（28）可得各阵元位置之差
[Z=±Nmd-Mnd, 0≤n≤N-1,]
[0≤m≤2M-1] （29）
去除其中的重复位置差后可以得到集合[Zx]
[Zx=-2M-1Nd, -2MNd,. .., 2M-1Nd]
（30）
显然，差分阵虚拟阵元位置与接收到的数据协方差矩阵[R]中的各元素存在对应关系，且同一虚拟阵元对应的协方差矩阵中的各个元素的值并不相等，因此需要对这些元素取平均
[Rpx=1ωpxi=1ωpxRipx] （31）
其中，[Ripx]表示同一波程差对应的[R]中的第i个协方差矩阵元素。
由式（29）可知，由此构造的差分阵中存在空洞部分，许多方法仅仅对差分阵中的连续虚拟阵列中的均匀部分进行空间平滑处理，而舍弃了其非均匀部分，导致算法精度有所下降。为了充分利用差分阵中的虚拟阵元，对去噪处理后的数据协方差矩阵进行扩展，将原来的[(2M+N-1)×(2M+N-1)]数据矩阵[R]扩展为[(2MN-N+1)×(2MN-N+1)]的Toeplitz矩阵[Rto]，该矩阵里的元素满足以下性质
[Rtoi, j=R(i-j)d,i-jd∈Zx0,i-jd?Zx] （32）
使用该方法扩展得到的矩阵[Rto]中存在零元素，这些零元素与差联合阵产生的空洞存在着一一对应的关系。为了更精确地估计DOA，此时需要对矩阵[Rto]再次进行补全。
2.3 APA矩阵补全算法
为了获得更精确的估计值，此处使用APA算法进行矩阵填补。根据文献［17］，式（13）重写为
[minRR? s.t. RΩ-GΩF≤δ] （33）
其中，[?F]是矩阵的Frobenius范数，[δ]为大于0的控制拟合误差的比较小的参数。
与其它秩或核范数最小化方法不同的是，该方法在于找到一个矩阵，使得它具有低秩且它的样本集[Ω]上的元素与观察项上的值一致。则该问题可以为：
[find R s.t. rankR≤r, RΩ-GΩF≤ε] （34）
然后该矩阵补全问题可以被表述为寻找一个公共点
[find R∈Sr?Sp] （35）
其中
[Sr:=R|rankR≤r] （36）
[Sp:=R|RΩ-GΩF≤ε] （37）
[Sr]表示矩阵秩约束集，[Sp]表示矩阵保真约束集。
给定一个集合[S]，将一个点[X?S]在集合[S]上的投影记为[S(X)]，其可以被定义为：
[S(X):=argminR∈SR-X2F] （38）
根据Eckart-Young定理，该式可以通过[X]的截断奇异值分解求得：
[S(X)=i=1rσiμiυTi] （39）
其中，[σiri=1]是矩阵[X]的[r]个最大的奇异值，[μiri=1]和[υiri=1]分别对应矩阵[X]对应的左右奇异向量。
为了找到公共点，采用交替投影策略，则有：
[Rk=Sr(Xk)Rk+1=Sp(Rk)] （40）
其中，[X0=GΩ]。当满足迭代终止条件时，交替投影停止并输出[Rk+1]，该矩阵即为重构后的协方差矩阵补全空洞后的矩阵。
2.4 DOA估计
N次快拍接收数据通过对协方差矩阵[R]进行去噪，重构，扩展，填充后得到目标矩阵[Rmc]，[Rmc]中包含了N次快拍的数据信息，可以应用MUSIC算法对其进行DOA估计。首先，对[Rmc]进行特征分解，得到信号子空间[Us]和噪声子空间[Un]
[Rmc=UssUHs+UnnUHn] （41）
根据两者正交的特性，可以得到其空间频谱
[Pθ=1αHθUnUHnαθ] （42）
通过搜索空间谱的峰值，找到的对应于最高峰值的角度即为所求。
3 仿真实验
为了验证所提出的算法的有效性和准确性，将该算法与传统的MUSIC算法、文献［8］所提出的算法、基于TNNR的DOA估计算法、基于IALM的DOA估计算法和克拉美劳下界进行对比。互质阵列的子阵列阵元数分别为[M=3, N=5]，因此阵元的数目为[2M+N-1=10]。均方根误差定义为
[RMSE=1qLq=1ql=1Lθlq-θl] （43）
其中，q为蒙特卡罗实验次数，L为入射进阵列的信号数目，[θl]为给定信号的真实角度，[θlq]为第q次估计角度。
在第一个实验里，将所提算法的性能与传统的MUSIC算法、文献［8］所提算法和基于IALM的算法进行比较。设置信号的数目[L=3]，[SNR=0 dB]，快拍数为500，信源的DOA分别为[5°]、[10°]和[15°]。4种算法的实验结果如图2所示。
图2显示，本文所提算法与基于IALM的算法性能最好，传统的MUSIC算法性能最差。这是因为传统的MUSIC算法的假设是均匀白噪声背景下，而在色噪声情况下，估计性能下降。而文献［8］所提算法虽然在非均匀噪声背景下，但该噪声只对数据协方差矩阵的主对角元素产生影响，而且该算法重构去噪的思想适用于非均匀噪声中的噪声功率中存在的最小值与其它值相差较大的情况，而在本文色噪声条件下，其性能下降。本文通过对数据协方差矩阵使用Lp-TNN算法进行去噪和扩展，有效地抑制了色噪声。
在第二个实验里，将所提算法的性能与传统的MUSIC算法、文献［8］所提算法、基于TNNR的算法、基于IALM的算法和克拉美劳下界进行比较。对不同算法在不同SNR条件下的均方根误差进行对比，设置信号的数目[L=2]，快拍数为500，信源的DOA分别为[10°]和[20°]，信噪比的变化范围为[-10～20 dB]，蒙特卡罗实验次数为300。6种算法的实验结果如图3（a）所示。
图3（a）表明，使用2次补全后的后三种算法相比于其他两种算法性能提升很多且与克拉美劳下界接近。上述所提算法由于使用补全算法（具有优越性），使得其估计性能略高于另外两种算法。因为使用矩阵补全算法对欠定无噪协方差矩阵进行补全，使得该协方差矩阵更适用于差联合阵的扩展，而文献［8］所提算法只是将主对角元素替换为最小值，导致使用差联合阵扩展矩阵时损失的元素增多，因此其性能较差。传统的MUSIC算法无法对色噪声进行去噪，性能最差。
在第三个实验里，对6种不同算法在不同快拍数条件下的均方误差进行对比。设置信号的数目[L=2]，快拍数的变化范围为[100~700]，信源的DOA分别为[10°]和[20°]，信噪比为[0 dB]，蒙特卡罗实验次数为300，实验结果如图3（b）所示。图3（b）表明，在快拍数增加的情况下，所有算法的估计性能变得越来越好。其中，使用2次补全后的后三种算法相比于其他两种算法性能提升很多且与克拉美劳下界接近。本文所提算法估计性能略高于另外2种算法。文献［8］算法的性能次之，传统的MUSIC算法的DOA估计性能最差。
最后1个实验里，对6种不同算法的分辨率进行对比。设置信号的数目[L=2]，快拍数为500，信源的DOA的差值为[5°~50°]，信噪比为[0 dB]，蒙特卡罗实验次数为300。实验结果如图3（c）所示。
由图3（c）可知，本文算法和基于TNNR、IALM的算法性能接近且接近于克拉美劳下界。其中，本文算法分辨精度略高于其他两种算法。这是因为使用APA算法进行补全时，只是单纯地补全了空洞部分，而其它部分的值未改变，从而大大保留了信息的准确性，而其他算法没有这个细节，所以本文算法的分辨精度相比于其他两种略高。文献［8］所提算法性能次之，但是该算法在角度间隔较小（小于等于5°）时性能明显低于其他算法，这是由于该算法的去噪效果低于本文算法，且去噪后丢失的信息比较多导致其分辨率性能降低。传统MUSIC算法性能总体最低。因此，本文算法在不同角度间隔下具有高分辨率的优势。
4 结 论
为了解决色噪声条件下的互质阵列DOA估计性能差问题，本文提出了一种将互质阵列和低秩恢复理论相结合来估计DOA的方法；同时，将本文算法与基于TNNR和IALM的算法进行对比，得到了最优的本文算法。首先，使用Lp-TNN算法对含噪矩阵进行去噪处理，一方面有效抑制了色噪声对信号协方差矩阵的影响，另一方面使得该协方差矩阵更适用于后面差联合阵的扩展；然后，使用APA算法对使用差联合方法扩展后的矩阵的空洞进行填补，该算法能在不改变已有信息的条件下补全空洞，因此大大地提高了DOA估计的分辨率。仿真实验表明，色噪声条件下，本文算法相比于其它算法有更好的估计和分辨性能。