技术领域
[0001] 本发明涉及阵列信号处理、无源定位技术、DOA估计、无人机故障检测与隔离领域，特别涉及一种基于无人机机载阵列的干扰源定位方法及系统。
相关背景技术
[0002] 无人机具有造价较低、飞行环境要求低、安全性高、便于携带等优点，因此无人机在定位、电力巡检、影视拍摄、农业喷洒等领域有着广泛的应用。近年来，电磁干扰事件经常发生，对卫星通信与导航造成干扰。针对日益增加的监测区域要求，地基设备灵活性较差，而无人机具有较强的灵活性，且观测视距的影响也远小于地面设备。因此基于无人机平台的定位具有一定现实意义。
[0003] 到达角度(DOA)估计在无源定位中有着至关重要的地位，经典的超分辨率DOA估计算法主要有MUSIC算法，ESPRIT算法，这些都是基于特征结构的子空间算法。与传统技术不同，子空间算法利用接收信号的特征结构，而不是利用收集数据的统计特征，突破了“瑞利限”的限制，这使得估计分辨力显著提高。许多天线阵列都可以用来进行DOA估计，最流行的平面阵列是均匀矩形阵列(URA)，其中，元件间的间距需要在半波长以内以避免空间混叠。然而，传统URA的自由度受到传感器数量的严格限制。此外，还受到阵列元件之间严重的相互耦合效应的影响。
[0004] 稀疏阵列可以有效解决上述问题。互质阵列是稀疏阵列的一种，和传统URA相比，互质阵列能够增加自由度，同时提高算法的估计性能。虽然与嵌套阵列相比，互质阵列的空间自由度较低，但是它减少了相互耦合的影响。
具体实施方式
[0124] 下面结合说明书附图和实施例对本发明作进一步限定，但不限于此。
[0125] 实施例1
[0126] 一种基于无人机载天线阵列的卫星导航干扰源定位方法，如图1和图4所示，包括：
[0127] 步骤1，数个搭载有天线阵列的无人机按照固定轨迹运动，采集无人机在不同的飞行状态下的多种传感器数据信号，多种传感器数据信号包括三轴加速度，三轴角速度，滚转角、俯仰角、偏航角，对采集到的传感器数据信号进行去噪和归一化处理，从而尽可能的提取出信号的有效特征；该处有效特征是指原始数据经过去噪处理后消除了数据中包含的噪声等无关信息，然后经过归一化消除数据的量纲是信号具有统一性，这样做可以增强信号的准确性和可靠性，具体的特征是由后面的神经网络进行学习；
[0128] 步骤2，根据采集数据的种类和模式设置神经网络的参数，通过误差函数确定预测误差，当误差满足需要的精度时结束训练，利用训练好的网络进行无人机故障检测；
[0129] 步骤3，无人机上搭载的天线阵列接收多个干扰信号，得到互质阵列接收信号模型；
[0130] 步骤4，计算所接收的干扰信号的协方差矩阵并将其向量化；将机载互质阵列虚拟化得到虚拟差分阵列，分别计算虚拟阵列正负两部分的协方差矩阵并将其向量化；
[0131] 步骤5，通过空间平滑构建原始阵列和虚拟阵列的综合协方差矩阵；
[0132] 步骤6，利用稀疏表示，构造整个监测空间内的字典，进而构造多测向向量稀疏表示模型；
[0133] 步骤7，将多测向向量稀疏表示模型转化为凸优化问题，构建二阶锥优化问题，通过谱峰搜索得到干扰源的DOA角度估计。
[0134] 实施例2
[0135] 根据实施例1所述的一种基于无人机载天线阵列的卫星导航干扰源定位方法，其区别在于：
[0136] 步骤1中，对采集到的传感器数据信号进行去噪和归一化处理，提取出信号的有效特征；包括：
[0137] 无人机的常见故障按时间特性可分为间隔性故障、瞬时性故障、慢性故障，无人机发生飞行故障时通常会使输出的信号产生各种噪声，因此，使用小波阈值降噪对采集到的传感器数据信号抑制噪声，假设含有噪声的信号公式如式(1)所示：
[0138] f(t)＝x(t)+n(t) (1)
[0139] 式(1)中，f(t)表示含有噪声的信号，x(t)表示目标信号，n(t)表示噪声信号；
[0140] 设置适当的阈值分离出噪声；通过软阈值函数得到更加平滑的信号，软阈值函数表达式如下式所示：
[0141]
[0142] 表示阈值化后的小波系数，Wj,k表示原始的小波系数，sign表示符号函数，T表示阈值，j表示尺度因子，k表示位置因子；
[0143] 硬阈值可以保留信号的有效特征，硬阈值函数表达式如下式所示：
[0144]
[0145] 其中， 表示阈值化后的小波系数，Wj,k表示原始的小波系数，T表示阈值。阈值函数如何选择对降噪效果影响很大，硬阈值可以保留信号中的尖锐点，但信号容易产生间断，而软阈值信号更加平滑，但容易使得信号产生偏离，因此，本发明选择软硬阈值折衷的方法分离出噪声，如式(2)所示：
[0146]
[0147] 式(2)中， 表示阈值化后的小波系数，Wj,k表示原始的小波系数，sign表示符号函数，α表示变化因子，T表示阈值，这种方法能够调整α的值来选择合适的阈值进行降噪。通用的计算阈值T的公式如下所示：
[0148]
[0149] T表示阈值，σ表示噪声的均方差，N表示信号的长度。但该公式得到的阈值不能随j的变化而变化，因此，修改阈值计算公式使得阈值应随j的增大而减小，计算阈值T的公式如式(3)所示：
[0150]
[0151] 用信噪比和均方根误差来验证，信噪比的公式如下式所示：
[0152]
[0153] 其中，Ps表示信号的功率，Pn表示噪声的功率，f(n)表示原始的信号， 表示降噪后的信号。均方根误差RMSE的公式如下式所示：
[0154]
[0155] 其中，f(n)表示原始的信号， 表示降噪后的信号，n表示样本的数量大小。
[0156] 对信号进行归一化处理，才能作为神经网络的输入，归一化公式如式(4)所示：
[0157]
[0158] 式(4)中，xi为输入的第i个值，xmin表示输入的最小值，xmax表示输入的最大值，表示归一化之后的结果。
[0159] 步骤2中，根据采集数据的种类和模式设置神经网络的参数，通过误差函数确定预测误差，当误差满足需要的精度时结束训练，利用训练好的网络进行无人机故障检测；包括：
[0160] BP(Back Propagation)神经网络是一种按照误差逆向传播算法训练的多层前馈神经网络，是应用最广泛的神经网络。BP神经网络包括输入层、隐含层和输出层；设x1、x2、…xn为输入，i为输入的节点数，y1、y2、…ym为输出，k为输出节点，隐含层为单隐层，则隐含层第j个节点的输出pj公式如式(5)所示：
[0161]
[0162] 式(5)中，f1表示激活函数，wji表示输入层和隐含层之间的权值，aj表示隐含层的阈值；则输出层第k个节点的输出yk公式如式(6)所示：
[0163]
[0164] 式(6)中，f2表示激活函数，wkj表示输入层和隐含层之间的权值，bk表示隐含层的阈值；设n为输入样本的数量，则网络误差公式即误差函数如式(7)、式(8)所示：
[0165]
[0166]
[0167] 式(7)、式(8)中，p表示样本的数量，m表示输出节点的数量，dk表示期望输出的结果，yk表示实际输出的结果，Ep表示第p个样本的正向传播结果与期望结果的误差， 表示第p个样本的第k个输出节点的期望输出结果， 表示第p个样本的第k个输出节点的实际输出结果，综合起来，误差函数Ep衡量了每个样本的正向传播结果与期望输出之间的差异，而总误差函数E则是所有样本误差的累加；
[0168] 输出层的权值变化量公式如式(9)所示：
[0169]
[0170] 式(9)中，△wjk表示输出层的权值变化量，η表示学习率；输出层的阈值变化量△bk、隐含层的权值变化量△wji、隐含层的阈值变化量△aj公式分别如式(10)、式(11)、式(12)所示：
[0171]
[0172]
[0173]
[0174] 由上面公式可得，权值和阈值的更新公式如式(13)、式(14)、式(15)、式(16)所示：
[0175] wji(n+1)＝wji(n)+△wji (13)
[0176] wkj(n+1)＝wkj(n)+△wkj (14)
[0177] aj(n+1)＝aj(n)+△aj (15)
[0178] bk(n+1)＝bk(n)+△bk (16)
[0179] 式(13)、式(14)、式(15)、式(16)中，wji(n)、wkj(n)、aj(n)、bk(n)分别表示前面一层的权值和阈值，而wji(n+1)、wkj(n+1)、aj(n+1)、bk(n+1)表示更新之后的权值和阈值；多次训练后可使误差逼近极小值，那么就可以完成故障检测的目的。
[0180] 利用训练好的网络进行无人机故障检测：输入层的数据为一个9维的向量，代表x轴加速度、y轴加速度、z轴加速度、x轴角速度、y轴角速度、z轴角速度、滚转角、俯仰角以及偏航角的信号特征，采集到的传感器数据经过小波去噪得到去噪后的信号，再进行归一化处理作为训练好的网络的输入变量，输出层的数据为一个3维向量，代表三种状态：正常状态、慢变故障、突变故障。
[0181] 步骤3中，无人机上搭载的天线阵列接收多个干扰信号，并得到接收信号模型；包括：
[0182] 空间中有K个干扰源sk(t)，k＝1,2,…,K，其来波方向为θ＝[θ1,θ2,…,θK]；空间中有数个无人机，机上各载有一个包含M＝2线阵和N＝3线阵的互质阵列；阵列示意图如图2所示。该互质阵列的传感器分布位置为{0,Md,2Md,…,(N‑1)Md}和{0,Nd,2Nd,…,(M‑1)Nd}，d＝λ2为单位阵元间距，λ是信号波长；两个线阵的第一个元素重叠，互质阵列的阵元数为M+N‑1，互质阵列的位置集合如式(17)所示：
[0183]
[0184] 每个互质阵列的两个线阵接收到的信号如式(18)、式(19)所示：
[0185]
[0186]
[0187] 式(18)、式(19)中，x1(t)和x2(t)是在t时刻无人机上的互质阵列的两个线阵分别接收到的信号，τ1(θk)和τ2(θk)表示来波角度为θk的信号到达两个线阵的时间延迟，n1(t)和n2(t)为观测噪声，且该噪声均符合零均值高斯分布并和信号sk(i)不相关；
[0188] 对x1(t)和x2(t)进行采样得到离散形式并作L点DFT变换，其频域模型如式(20)、式(21)所示：
[0189] X1＝A1S+N1(20)
[0190] X2＝A2S+N2(21)
[0191] 其中，A1＝[α1(θ1),α1(θ2),…,α1(θK)]；A2＝[α2(θ1),α2(θ2),…,α2(θK)]；S＝[s1T T T(t),s2(t),…,sK(t)]；N1＝[n1,0,n1,1,…n1,M‑1]；N2＝[n2,0,n2,1,…n2,N‑1]；A1和A2表示两个线阵的流形矩阵， 和
表示流形矩阵各列指向向量，S是信号源
向量，N1和N2是两个线阵的噪声向量；将式(20)、式(21)中的两个信号合并为一个信号，得到互质阵列接收信号模型X；如式(22)所示：
[0192]
[0193] 互质阵列接收信号模型X的阵列流形矩阵为
[0194] 计算所接收的干扰信号的协方差矩阵并将其向量化；将机载互质阵列虚拟化得到虚拟差分阵列，分别计算虚拟阵列正负两部分的协方差矩阵并将其向量化；包括：
[0195] 计算互质阵列接收信号模型的协方差矩阵；其协方差矩阵可以用外积的统计平均H值表示，即Cxx＝E[XX]。由于信号和噪声相互独立，如式(23)所示：
[0196] Cxx＝E[XXH]＝ACssAH+W (23)H 2
[0197] 式(23)中， W＝E[NN ]＝σ I； 为信号功率，2
σ为噪声功率，I是大小为(M+N‑1)×(M+N‑1)的单位矩阵；
[0198] 将协方差矩阵向量化，得到向量化协方差矩阵，如式(24)所示：
[0199] V1＝D1p+μ (24)*
[0200] 式(24)中， D1＝A ⊙A，为新的阵列流形矩阵，⊙表示矩阵的Khatri‑Rao积，μ为矩阵W的向量化形式；
[0201] 考虑互质阵列虚拟化，式(24)看成一个阵元位置为集合{xi‑xj,0≤i,j≤M+N‑1}的虚拟线性阵列，其中，xi表示第i个阵元的位置，D1即为虚拟阵列流形矩阵；为了将互质阵列虚拟化，利用其稀疏性，对两线阵进行差分操作，获得虚拟阵列；如图3所示，虚拟阵列的阵元位置由集合 给出；可以看到，虚拟阵列中阵元位置除了中间位置，有正有负。因此，对虚拟阵列的正、负两部分分别作与信号模型相同的操作；即：先计算出正部分与负部分的协方差矩阵，再将协方差矩阵向量化；故对于正位置部分，如式(25)所示：
[0202] C2＝E[XXT]＝ASSTAT+W (25)
[0203] 将其向量化得到式(26)：
[0204] V2＝D2p2+μ (26)
[0205] 其中，D2是正位置阵元的流形矩阵，信号功率为实值，故p2＝p；
[0206] 再对虚拟阵列的负位置部分求协方差矩阵，如式(27)所示：
[0207] C3＝E[XXH]＝ASSHAH+W (27)
[0208] 将其向量化得到式(28)：
[0209] V3＝D3p3+μ (28)
[0210] 其中，D3是负位置阵元的流形矩阵，p3＝p。
[0211] 通过空间平滑构建原始阵列和虚拟阵列的综合协方差矩阵；包括：将式(24)、式(26)、式(28)三个向量模型合并，如式(29)所示：
[0212]
[0213] 相应的阵列流形矩阵为
[0214] 合并模型式(29)由原始互质阵列和虚拟阵列合并构成，故其中的元素必然有诸多冗余，并且会出现乱序。因此需要重构该模型，去掉其中的冗余并重新排序。
[0215] 采用空间平滑技术，首先，从流形矩阵D中抽取一个维度是(2MN+2M+2N‑1)×K的矩阵，相当于从流形矩阵D中抽取2MN+2M+2N‑1个不同的行进行排序，这些不同的行对应互质矩阵的差分后所形成的虚拟阵列的自由度。将虚拟阵列中的2MN+2M+2N‑1个阵元分成MN+M+N个子阵，接收信号模型式(29)重构为式(30)：
[0216]
[0217] 且信号模型的综合协方差矩阵为式(31)：
[0218]2 2
[0219] 其中，ω＝σe，Ω＝σE，e是(2MN+2M+2N‑1)×1的列向量， 是(MN+M+N)×K的流形矩阵。从其协方差矩阵的形式上可以看出它与一个由MN+M+N个天线组成的均匀线性阵列的形式相同。
[0220] 步骤(6)中，利用稀疏表示，构造整个监测空间内的字典，进而构造多测向向量稀疏表示模型；包括：
[0221] 基于阵列和式(31)的稀疏性，将观测的空域角度区间按照一定的角度间隔分成Q个角度值；Q＞＞K，定义集合如式(32)所示：
[0222]
[0223] 式(32)，Θ表示整个观测空域的搜索网格； 表示整个观测空间的第q个网格的角度值，q＝1,2,…Q，流形矩阵表示为 其也可被称为过完备基或者字典。
[0224] 那么式(31)的第j列如式(33)所示：
[0225]
[0226] 式(33)中，1≤j≤MN+M+N，gj∈RQ×1是字典 的稀疏表示系数；若gj的第i行是非零元素，则说明 处存在一个信号源。虚拟阵列中的2MN+2M+2N‑1个阵元分成MN+M+N个子阵，则其协方差矩阵中的每一列代表一个子阵，对应一个稀疏表示系数gj，该系数的非零行代表信号源位置，故其非零行数为K。因为稀疏表示系数gj是单个列稀疏向量，故将MN+M+N个第j列合并为一个矩阵形式，得到多测向向量稀疏表示模型，如式(34)所示：
[0227]
[0228] 式(34)中，G＝[g1,g2,…,gMN+M+N]是一个稀疏系数矩阵；并且矩阵中的每列都有相同的稀疏结构，因此G具有联合稀疏性，稀疏表示模型称为多测向向量稀疏表示模型。若求得唯一的稀疏表示系数G，那么通过其对应的非零行位置即找到干扰信号的DOA。
[0229] 步骤(7)中，将多测向向量稀疏表示模型转化为凸优化问题，构建二阶锥优化问题，通过谱峰搜索得到干扰源的DOA角度估计；包括：
[0230] 定义γ(G)表示G中的非零行索引集合，|γ(G)|＝K表示非零行个数，即稀疏度，那么对应的DOA估计问题表示为约束优化问题，如式(35)所示：
[0231]
[0232] 定义 为矩阵G的估计，其行之间具有稀疏性，列之间不具有稀疏性；定义b＝[b1,Tb2,…,bQ]，其每一个元素代表 中对应行向量的2范数，即 而b的1范数||b||1
代表 在空域Θ上的稀疏约束；因此，优化问题更改为式(36)：
[0233]
[0234] 噪声功率用SR的最小特征值作为其估计，定义为υ，故将式(36)正则化得到正则化模型，如式(37)所示：
[0235]
[0236] 其中， 是SR的估计，δ是正则化参数， 代表Frobenius范数；正则化模型为凸优化问题，采用二阶锥优化(SOCP)对其进行求解，如式(38)所示：
[0237]
[0238] 式(38)中，辅助变量ε和η来线性化目标函数，I为Q×1的单位向量，ξ为1×Q的向量，ξq是ξ的第q个元素；通过式(38)求出稀疏向量b后，通过谱峰搜索即得到干扰源的DOA角度估计；根据多个阵列的DOA估计值，通过交叉定位得到目标源的位置。
[0239] 设置7个待测远场源，其到达角为θ＝[‑35°,‑25°,‑15°,0°,10°,20°,35°]。图5为本发明方法与同阵列MUSIC算法的仿真结果示意图，由图5对比可知：MUSIC算法在‑25°和‑15°两个角度的谱峰过于接近，本发明方法具有更好的谱峰分辨率。图6为本发明方法与线阵MUSIC算法的仿真结果示意图，由图5对比可知：本发明方法无论是在定位准确度还是谱分辨率方面都明显优于线阵MUSIC算法。
[0240] 实施例3
[0241] 一种计算机设备，包括存储器和处理器，存储器存储有计算机程序，处理器执行计算机程序时实现实施例1或2所述的基于无人机载天线阵列的卫星导航干扰源定位方法的步骤。
[0242] 实施例4
[0243] 一种计算机可读存储介质，其上存储有计算机程序，计算机程序被处理器执行时实现实施例1或2所述的基于无人机载天线阵列的卫星导航干扰源定位方法的步骤。
[0244] 实施例5
[0245] 一种基于无人机载天线阵列的卫星导航干扰源定位系统，包括：
[0246] 信号有效特征提取模块，被配置为：数个载有天线阵列的无人机按照固定轨迹运动，采集无人机在不同的飞行状态下的多种传感器数据信号，对采集到的传感器数据信号进行去噪和归一化处理，从而尽可能的提取出信号的有效特征；
[0247] 无人机故障检测模块，被配置为：根据采集数据的种类和模式设置神经网络的参数，通过误差函数确定预测误差，当误差满足需要的精度时结束训练，利用训练好的网络进行无人机故障检测；
[0248] 接收信号模型获取模块，被配置为：无人机上搭载的天线阵列接收多个干扰信号，并得到接收信号模型；
[0249] 多测向向量稀疏表示模型构造模块，被配置为：计算所接收的干扰信号的协方差矩阵并将其向量化；将机载互质阵列虚拟化得到虚拟差分阵列，分别计算虚拟阵列正负两部分的协方差矩阵并将其向量化；通过空间平滑构建原始阵列和虚拟阵列的综合协方差矩阵；利用稀疏表示，构造整个监测空间内的字典，进而构造多测向向量稀疏表示模型；
[0250] DOA角度估计模块，被配置为：将多测向向量稀疏表示模型转化为凸优化问题，构建二阶锥优化问题，通过谱峰搜索得到干扰源的DOA角度估计。
