CN102890265A - 一种基于水下声传感器网络的目标被动定位方法 - Google Patents

Abstract

本发明提供了一种基于水下声传感器网络的目标被动定位方法，随机布放多个声传感器节点和融合中心形成一种并行的网络拓扑结构，即水下声传感器网络；声传感器节点采集目标辐射噪声数据，求取该节点接收到目标的能量，各个声传感器节点将求解的能量分别传输至融合中心进行融合处理，依据能量衰减与距目标距离之间关系，对目标进行定位。发明提高了被动定位目标的精度；在传感器节点数目少及低信噪比环境中，定位鲁棒性好；基于RSS无线传感器网络定位装置实施简单。

Description

一种基于水下声传感器网络的目标被动定位方法

技术领域

本发明是一种水下目标被动定位方法，涉及信号处理、应用数学等领域。

背景技术

无线传感器网络技术(Wireless Sensor Networks,WSN)的发展和微机电技术(micro-electro-mechanical systems,MEMS)的进步大大的促进了水下无线传感器网络技术的发展，而源目标的定位问题不仅是声信号处理的重要研究内容，而且是水下无线传感器网络的一个重要应用。无线传感器网络是由大量分布在被监测区域内的传感器节点通过自组网形成的网络，已广泛应用于环境监测、健康卫生、灾难预警等方面。声源定位作为水下无线传感器网络的一个重要应用，已被学者广泛研究，其基本原则是：水下传感器节点将采集到的数据经过一定处理，向并行拓扑结构中的融合中心(Fusion Center,FC)传输处理结果，融合中心依据各个传感器节点传输的数据融合处理，得出目标的位置估计。

依据利用的目标信息不同，可以将定位方法分为：基于到达时间的定位(Time ofarrival,TOA)、基于到达时延的定位(Time delay of arrival,TDOA)、基于波达方向的定位(Direction of arrival,DOA)和基于信号强度的定位(Received signal strength,RSS)。TDOA、TOA需要精确的时间同步，对于水下无线传感器网络而言，由于水声信道带宽和传感器节点能量的限制，精确的时间同步是异常困难的；TDOA和TOA需要检测信号到达时间，这一任务对于被动的目标定位而言仅仅依据目标辐射噪声是困难的，DOA需要接收设备配置测向装置，实现复杂；基于RSS的水下目标定位采集水下声信号能量，与原始数据相比，可以采用较低的采样率获取RSS的数据。

发明专利(专利申请号：CN200810099378.6)公开了一种基于RSSI(Received signalstrength indicator)的无线传感器网络节点定位方法，利用已知坐标节点(锚节点)与未知坐标节点之间的RSSI信息对，调节对未知节点位置的权重，从而获取节点位置信息。该发明在一定程度上说明了可以利用接收到RSSI(或RSS)信息中包含的距离信息对未知节点的位置(或目标位置)进行估计。发明专利(专利申请号：CN200810150500.8)公开了一种利用无线能量强度进行移动目标定位的方法，该方法利用移动节点(目标)在移动过程中，轮换发射不同能量的信号，以在固定的接收节点处获得不同发射功率或能量下的接收信号强度与距离关系，进而获取大量距离信息进行定位。该发明要求目标配合整个定位系统工作，如：煤矿井下人员定位时需要目标(待定位的人员)携带射频信号发射器。因此该发明不适应这种不具有配合能力的目标，如海洋环境中鱼群、海豚等目标定位，普适性不强。发明专利(专利申请号：200810118829.6)公开了一种声源定位的无线测量方法，通过利用在一定区域布放无线声音传感器，在确知各个传感器坐标情况下，依据各个传感器采集到的能量信息(RSS)进行声源定位。汇聚节点选取能量较大的W=3~5个节点参与最终目标定位，将求解问题转换为优化问题，采取最速下降法对目标位置进行搜索，从而实现目标定位。该发明仅仅利用3~5个接收能量较大的传感器节点进行目标定位，没有充分利用无线传感器网络的全部信息；该发明适用于较小规模的无线传感器网络下(节点数目N=6~20)的二维目标定位，实际应用中具有一定局限性。

发明内容

为了克服现有技术不能应用于目标不配合情况下目标定位、适用于小规模无线传感器网络目标定位、不能充分利用无线传感器网络全部信息等方面的不足，本发明提供一种基于水下声传感器网络的水下目标被动定位方法，提高了被动定位目标的精度；在传感器节点数目少及低信噪比环境中，定位鲁棒性好；基于RSS无线传感器网络定位装置实施简单。

本发明解决其技术问题所采用的技术方案包括以下步骤：

步骤一、随机布放多个声传感器节点（节点数目N≥2）和融合中心形成一种并行的网络拓扑结构，即水下声传感器网络；

步骤二、声传感器节点采集目标辐射噪声数据，求取该节点接收到目标的能量(RSS)，包括以下步骤：

(1)信号衰减模型：第i个传感器节点接收到的信号为

s(n)表示距离目标1m处的声强，γ_i是第i个传感器的接收增益，r_s表示目标位置，r_i表示第i个传感器位置，r_s和r_i维数为p×1，p=2、3，

N(x,y)表示均值为x，方差为y的正态分布，传感器节点位置可以通过自定位算法求解。假设每个传感器观测噪声方差相同，即d_i＝||r_s-r_i‖是目标与第i个传感器之间的欧几里得距离。假设s(n)与ω_i(n)相互独立，即E[s(n)ω_i(n)]＝0，E[·]表示期望；

(2)在步骤(1)建立的信号基础上，传感器节点接收到目标的能量为

实际求解过程中，期望是通过在一个窗内加权平均实现的，即 $h_i=g_i\frac{S^{\′}}{{\left|\right|r_s-r_i\left|\right|}^2}+\frac{1}{L}\underset{n=0}{\overset{L-1}{\mathrm{\Σ}}}{\mathrm{\ω}}_i^2\left(n\right)=g_i\frac{S^{\′}}{{\left|\right|r_s-r_i\left|\right|}^2}+{\mathrm{\θ}}_i,$ 其中第一项目是信号能量，第二项是噪声能量，S'＝s²(n)，

L是窗函数长度，即采样点数，θ_i服从

的卡方分布，依据中心极限定理，当L＞＞30时，θ_i~N(σ²,2σ²L)，即近似服从正态分布，假设各个传感器节点接收增益相同，即g₁＝g₂＝…＝g_N，则第i个传感器节点接收到目标的能量可以表示为其中S＝S′×g_i。

(3)步骤(2)得出的第i个传感器节点接收到目标的能量

表达式两边减去随机变量θ_i的均值σ²，得到第i个传感器节点接收到的目标能量演化式

其中

N个传感器节点接收到的目标能量演化式子组成能量矩阵E＝[y₁,y₂,…y_N]；

步骤三、各个声传感器节点将求解的能量分别传输至融合中心进行融合处理，依据能量衰减与距目标距离之间关系，利用多个节点获得的大量距离信息对目标进行定位。

所述的融合处理是指：每个传感器将计算得到的y_i通过水声信道传输至融合中心，融合中心将目标被动定位问题转化为目标位置的最大似然估计问题，进而利用SDR算法将目标位置的最大似然估计问题转化为凸优化问题，利用凸优化工具箱求解目标位置。

融合处理过程中所述的将目标被动定位问题建模为最大似然估计问题的步骤为：

(1)选取矩阵E中最大能量值对应节点为参考节点1，假设其能量值为y₁，联合第i个传感器节点接收到的目标能量y_i(i＝2,3,…N，N为传感器节点数目)和参考节点接收到的目标能量y₁，两式子相除，消去S，得到

其中，

${\mathrm{\ρ}}_i^2=\frac{{\left|\right|r_i-r_1\left|\right|}^2}{{(y_i-y_1)}^2}{y}_i{y}_1,$ ${\mathrm{\υ}}_i=\frac{{\mathrm{\ϵ}}_i{d}_i^2-{\mathrm{\ϵ}}_1{d}_1^2}{y_i-y_1}\≈S\frac{{\mathrm{\ϵ}}_i{y}_1-{\mathrm{\ϵ}}_1{y}_i}{y_i{y}_1(y_i-y_1)}.$ ${\mathrm{\υ}}_i~N(0,{\mathrm{\σ}}_e^2),$ ${\mathrm{\σ}}_e^2=\frac{2{\mathrm{\σ}}_i^2{S}^2}{L}\frac{y_i^2+y_1^2}{y_1^2{y}_i^2{(y_i-y_1)}^2}(i=\mathrm{2,3},...N),$ 这样的等式个数为N-1；

(2)假设g＝[‖r_s-c₂||²,‖r_s-c₃||²,…||r_s-c_N‖²]^T，V＝[-υ₂,-υ₃,…-υ_N]^T，步骤(1)得到的等式写成矩阵形式，有P＝g+V，P的联合概率密度函数可以表示为 $f\left(P\right|r_s)={\left(2\mathrm{\π}\right)}^{-(N-1)/2}{\left|\mathrm{\Σ}\right|}^{-1/2}\exp \{-\frac{1}{2}{(P-g)}^T{\mathrm{\Σ}}^{-1}(P-g)\},$ 其中，Σ＝E[VV^T]是协方差矩阵，元素可以表示为 ${\left[\mathrm{\Σ}\right]}_{\mathrm{ij}}=\left\{\begin{array}{cc}\frac{2{\mathrm{\σ}}^2{S}^2}{L}\frac{y_i^2+y_1^2}{y_i^2{y}_1^2{(y_i-y_1)}^2}& i=j\\ \frac{2{\mathrm{\σ}}^2{S}^2}{L}\frac{y_i{y}_j}{y_1^2{y}_i{y}_j(y_i-y_1)(y_j-y_1)}& i\≠j\end{array}\right.;$

(3)目标位置的最大似然估计为

融合处理过程中所述的将目标位置的最大似然估计问题转化为凸优化问题(SDR算法)并估计出目标位置的步骤是：

(1)目标位置的最大似然估计可以表示为s.t.g_i＝‖r_s-c_i||²，i＝2,3,…N，s.t.表示约束条件；

(2)步骤(1)中等式可以重新写作 $\underset{r_s,g}{\min }\left\{\begin{array}{c}\left[\begin{array}{cc}{g}^T& 1\end{array}\right]\end{array}\left[\begin{array}{cc}{\mathrm{\Σ}}^{-1}& -{\mathrm{\Σ}}^{-1}P\\ -P^T{\mathrm{\Σ}}^{-1}& P^T{\mathrm{\Σ}}^{-1}P\end{array}\right]\left[\begin{array}{c}g\\ 1\end{array}\right]\right\},$

$s.t.g_i=\left[\begin{array}{cc}{r_s}^T& 1\end{array}\right]\left[\begin{array}{cc}I& -c_i\\ {c_i}^T& {\left|\right|c_i\left|\right|}^2\end{array}\right]\left[\begin{array}{c}{r}_s\\ 1\end{array}\right];$

(3)利用矩阵等式X^TAX＝tr(XX^TA)，其中tr(·)表示矩阵的迹，将步骤(2)中等式写作 $\underset{r_s,g,G}{\min }\left\{\begin{array}{c}\mathrm{tr}\left(\begin{array}{c}\left[\begin{array}{cc}G& g\\ g^T& 1\end{array}\right]\left[\begin{array}{cc}{\mathrm{\Σ}}^{-1}& -{\mathrm{\Σ}}^{-1}P\\ -P^T{\mathrm{\Σ}}^{-1}& P^T{\mathrm{\Σ}}^{-1}P\end{array}\right]\end{array}\right)\end{array}\right\},$ s.t. $g_i=\mathrm{tr}\left(\begin{array}{c}\left[\begin{array}{cc}{R}_s& r_s\\ {r_s}^T& 1\end{array}\right]\end{array}\left[\begin{array}{cc}I& -c_i\\ {c_i}^T& {\left|\right|c_i\left|\right|}^2\end{array}\right]\right),$ R_s＝r_sr_s ^T，G＝gg^T；

(4)步骤(3)中得到的目标函数式凸的，但是约束条件不是凸的，采用SDR将非凸问题转化为凸问题： $\underset{r_s,g,G}{\min }\left\{\begin{array}{c}\mathrm{tr}\left(\begin{array}{c}\left[\begin{array}{cc}G& g\\ g^T& 1\end{array}\right]\left[\begin{array}{cc}{\mathrm{\Σ}}^{-1}& -{\mathrm{\Σ}}^{-1}P\\ -P^T{\mathrm{\Σ}}^{-1}& P^T{\mathrm{\Σ}}^{-1}P\end{array}\right]\end{array}\right)\end{array}\right\},$ s.t. $g_i=\mathrm{tr}\left(\begin{array}{c}\left[\begin{array}{cc}{R}_s& r_s\\ {r_s}^T& 1\end{array}\right]\end{array}\left[\begin{array}{cc}I& -c_i\\ {c_i}^T& {\left|\right|c_i\left|\right|}^2\end{array}\right]\right),$ $\left[\begin{array}{cc}{R}_s& r_s\\ {r_s}^T& 1\end{array}\right]\≥0,$ $\left[\begin{array}{cc}G& g\\ g^T& 1\end{array}\right]\≥0,$ 其中[·]≥0表示矩阵半正定。需要说明的是：由于S未知，所以协方差阵Σ包含未知参数S，但是在求解优化问题时，常数项

不会影响目标函数的最小化过程，因此在计算目标函数过程中，可以不必理会未知项S²；

(5)步骤(4)得到的等式及约束条件都是凸的，求解得到目标位置。

本发明的有益效果是：在水下声传感器网络的目标被动定位方面，采用基于RSS的目标被动定位模型，克服了基于TDOA、TOA目标定位需要传感器节点间较为精确的时间同步以及需要确切的目标信号到达时间等缺点，克服了基于DOA目标定位需要配备测向装置的缺点；此外，基于RSS目标被动定位模型中传感器节点在采集数据时可以采用与采集原始数据相比较低采样率达到求解目标辐射能量的目的。将目标被动定位问题建模成目标位置的最大似然估计问题，该最大似然估计模型相比于加权最小二乘定位模型，提高了目标的定位精度。将最大似然估计模型通过SDR算法转化为凸优化问题求解目标位置，可以达到比加权最小二乘模型更高的定位精度，即这种模型转化相比于加权最小二乘定位模型而言，提高了定位精度，且当信噪比低(SNR<10dB)、传感器节点数目小(N＜10)，本发明的方法鲁棒性好，适用性强。

附图说明

图1是基于RSS的被动目标定位模型，其中，空心圆圈代表目标，虚线的圆表示目标辐射噪声或能量，实心圆圈代表传感器节点，FC表示融合中心；

图2不同信噪比条件下，本发明方法(CVX_SDP)与加权最小二乘(OS_WLS)目标被动定位性能对比；

图3不同传感器节点数目条件下，本发明方法(CVX_SDP)与加权最小二乘(OS_WLS)目标被动定位性能对比。

具体实施方式

假设N个传感器节点布放于25×25×25三维区域内，评估目标被动定位性能采取的标准是均方误差(Root mean square erro)

其中N_t是传感器节点分布次数，N_m是仿真过程中蒙特卡罗次数，在验证该发明目标被动定位性能时，为了消除不同传感器节点拓扑结构对目标被动定位性能造成的影响，采用N_t次25×25×25三维区域内的随机分布，为了消除相同信噪比下不同噪声对目标被动定位性能造成影响，采用N_m次蒙特卡洛仿真，即相同的性能分析重复做N_m次，不同的是每次观测的高斯分布值不同。信噪比定义为

其中S是距离目标1m处接收到的能量与各个传感器节点增益乘积，σ²是传感器节点观测噪声方差。

实施实例1：

由于仿真过程中，没有各个传感器节点接收到目标的能量值，因此需要依据信号衰减模型对传感器节点接收到的能量值建模，而在实际的定位中，无需此过程。实施实例参数设置：观测长度L＝5000，传感器节点数目N＝10，蒙特卡洛次数N_m=100，拓扑结构次数N_t=100，信噪比SNR=5dB~30dB，观测噪声方差σ²＝1，信号辐射能量与传感器节点增益乘积

仿真实施步骤如下：

(1)N个传感器节点均匀分布在给定的三维区域内，依据公式计算各个传感器节点接收到目标的能量值，获得能量矩阵E＝[y₁,y₂,…y_N]；

(2)选取矩阵E中最大能量值对应节点为参考节点1，联合第i个传感器节点接收到目标的能量y_i(i＝2,3,…N，N为传感器节点数目)和参考节点能量y₁，消去S，得到 ${\left|\right|r_s-c_i\left|\right|}^2={\mathrm{\&rho;}}_i^2+{\mathrm{\&upsi;}}_i,$ 其中， $c_i=\frac{y_i{r}_i-y_1{r}_1}{y_i-y_1},$ ${\mathrm{\&rho;}}_i^2=\frac{{\left|\right|r_i-r_1\left|\right|}^2}{{(y_i-y_1)}^2}{y}_i{y}_1,$ ${\mathrm{\&upsi;}}_i=\frac{{\mathrm{\&epsiv;}}_i{d}_i^2-{\mathrm{\&epsiv;}}_1{d}_1^2}{y_i-y_1}\&ap;S\frac{{\mathrm{\&epsiv;}}_i{y}_1-{\mathrm{\&epsiv;}}_1{y}_i}{y_i{y}_1(y_i-y_1)}.$ ${\mathrm{\&upsi;}}_i~N(0,{\mathrm{\&sigma;}}_e^2),$ ${\mathrm{\&sigma;}}_e^2=\frac{2{\mathrm{\&sigma;}}_i^2{S}^2}{L}\frac{y_i^2+y_1^2}{y_1^2{y}_i^2{(y_i-y_1)}^2}(i=\mathrm{2,3},...N),$ 这样的等式个数为N-1；

(3)假设g＝[‖r_s-c₂‖²,‖r_s-c₃‖²,…||r_s-c_N‖²]^T，

V＝[-υ₂,-υ₃,…-υ_N]^T，步骤(2)得到的等式写成矩阵形式，有P＝g+V，P的联合概率密度函数可以表示为 $f\left(P\right|r_s)={\left(2\mathrm{\&pi;}\right)}^{-(N-1)/2}{\left|\mathrm{\&Sigma;}\right|}^{-1/2}\exp \{-\frac{1}{2}{(P-g)}^T{\mathrm{\&Sigma;}}^{-1}(P-g)\},$ 其中，Σ＝E[VV^T]是协方差矩阵，元素可以表示为 ${\left[\mathrm{\&Sigma;}\right]}_{\mathrm{ij}}=\left\{\begin{array}{cc}\frac{2{\mathrm{\&sigma;}}^2{S}^2}{L}\frac{y_i^2+y_1^2}{y_i^2{y}_1^2{(y_i-y_1)}^2}& i=j\\ \frac{2{\mathrm{\&sigma;}}^2{S}^2}{L}\frac{y_i{y}_j}{y_1^2{y}_i{y}_j(y_i-y_1)(y_j-y_1)}& i\&NotEqual;j\end{array}\right.;$

(4)目标位置的最大似然估计可以表示为s.t.g_i＝‖r_s-c_i‖²，i＝2,3,…N，s.t.表示约束条件；

(5)步骤(4)中等式可以重新写完 $\underset{r_s,g}{\min }\left\{\begin{array}{c}\left[\begin{array}{cc}{g}^T& 1\end{array}\right]\end{array}\left[\begin{array}{cc}{\mathrm{\&Sigma;}}^{-1}& -{\mathrm{\&Sigma;}}^{-1}P\\ -P^T{\mathrm{\&Sigma;}}^{-1}& P^T{\mathrm{\&Sigma;}}^{-1}P\end{array}\right]\left[\begin{array}{c}g\\ 1\end{array}\right]\right\},$ s.t. $g_i=\left[\begin{array}{cc}{r_s}^T& 1\end{array}\right]\left[\begin{array}{cc}I& -c_i\\ {c_i}^T& {\left|\right|c_i\left|\right|}^2\end{array}\right]\left[\begin{array}{c}{r}_s\\ 1\end{array}\right];$

(6)利用矩阵等式X^TAX＝tr(XX^TA)，其中tr(·)表示矩阵的迹，将步骤(5)中等式写作 $\underset{r_s,g,G}{\min }\left\{\begin{array}{c}\mathrm{tr}\left(\begin{array}{c}\left[\begin{array}{cc}G& g\\ g^T& 1\end{array}\right]\left[\begin{array}{cc}{\mathrm{\&Sigma;}}^{-1}& -{\mathrm{\&Sigma;}}^{-1}P\\ -P^T{\mathrm{\&Sigma;}}^{-1}& P^T{\mathrm{\&Sigma;}}^{-1}P\end{array}\right]\end{array}\right)\end{array}\right\},$ s.t. $g_i=\mathrm{tr}\left(\begin{array}{c}\left[\begin{array}{cc}{R}_s& r_s\\ {r_s}^T& 1\end{array}\right]\end{array}\left[\begin{array}{cc}I& -c_i\\ {c_i}^T& {\left|\right|c_i\left|\right|}^2\end{array}\right]\right),$ R_s＝r_sr_s ^T，G＝gg^T；

(7)步骤(6)中得到的目标函数式凸的，但是约束条件不是凸的，采用SDR将非凸问题转化为凸问题： $\underset{r_s,g,G}{\min }\left\{\begin{array}{c}\mathrm{tr}\left(\begin{array}{c}\left[\begin{array}{cc}G& g\\ g^T& 1\end{array}\right]\left[\begin{array}{cc}{\mathrm{\&Sigma;}}^{-1}& -{\mathrm{\&Sigma;}}^{-1}P\\ -P^T{\mathrm{\&Sigma;}}^{-1}& P^T{\mathrm{\&Sigma;}}^{-1}P\end{array}\right]\end{array}\right)\end{array}\right\},$ s.t. $g_i=\mathrm{tr}\left(\begin{array}{c}\left[\begin{array}{cc}{R}_s& r_s\\ {r_s}^T& 1\end{array}\right]\end{array}\left[\begin{array}{cc}I& -c_i\\ {c_i}^T& {\left|\right|c_i\left|\right|}^2\end{array}\right]\right),$ $\left[\begin{array}{cc}{R}_s& r_s\\ {r_s}^T& 1\end{array}\right]\&GreaterEqual;0,$ $\left[\begin{array}{cc}G& g\\ g^T& 1\end{array}\right]\&GreaterEqual;0,$ 其中[·]≥0表示矩阵半正定。需要说明的是：由于S未知，所以协方差阵Σ包含未知参数S，但是在求解优化问题时，常数项

不会影响目标函数的最小化，因此在计算目标函数过程中，可以不必理会未知项S²；

(8)步骤(7)得到的等式及约束条件都是凸的，在MATLAB软件中可以直接采用凸优化工具箱来求解，从而得到目标位置的估计值。凸优化工具箱主要有：Michael Grant和Stephen Boyd的CVX、J.F.Sturm的SeDumi。

为了说明该发明基于水下声传感器网络的目标被动定位效果，分析对比了该发明与加权最小二乘被动目标定位精度，实施步骤如下：

(1)设置拓扑结构索引值Topo_index=1；设置蒙特卡洛仿真次数索引值Moto_index=1，index=1；

(2)设置信噪比SNR=5×index；

(3)设置节点均匀分布在25×25×25三维区域内，传感器节点形成一次拓扑结构；

(4)在相同信噪比不同的观测噪声值情况下，循环执行仿真实施步骤(1)~(8)，得出N_m次蒙特卡洛目标被动定位结果；

(5)重复本步骤(2)~(3)，总共得到N_mN_t个目标被动定位结果，依据公式

求解出一个信噪比SNR条件下的目标被动定位均方误差；

(6)index＝index+1，判断：if index＞6，绘制不同信噪比下目标定位的均方误差曲线，得到结果如图2所示，步骤结束，否则，重复执行步骤(2)~(5)。

实施实例2：

实施实例参数设置：观测长度L＝5000，传感器节点数目N＝5~20，蒙特卡洛次数N_m=100，拓扑结构次数N_t=100，信噪比SNR=5dB，观测噪声方差σ²＝1，信号辐射能量与传感器节点增益乘积

仿真实施步骤如下：

(1)N个传感器节点均匀分布在给定的三维区域内，依据公式

计算各个传感器节点接收到目标的能量值，获得能量矩阵E＝[y₁,y₂,…y_N]；

(2)选取矩阵E中最大能量值对应节点为参考节点l，联合第i个传感器节点接收到目标的能量y_i(i＝2,3,…N，N为传感器节点数目)和参考节点能量y_l，消去S，得到 ${\left|\right|r_s-c_i\left|\right|}^2={\mathrm{\&rho;}}_i^2+{\mathrm{\&upsi;}}_i,$ 其中， $c_i=\frac{y_i{r}_i-y_l{r}_l}{y_i-y_l},$ ${\mathrm{\&rho;}}_i^2=\frac{{\left|\right|r_i-r_l\left|\right|}^2}{{(y_i-y_l)}^2}{y}_i{y}_l,$ ${\mathrm{\&upsi;}}_i=\frac{{\mathrm{\&epsiv;}}_i{d}_i^2-{\mathrm{\&epsiv;}}_1{d}_l^2}{y_i-y_l}\&ap;S\frac{{\mathrm{\&epsiv;}}_i{y}_l-{\mathrm{\&epsiv;}}_l{y}_i}{y_i{y}_l(y_i-y_l)}.$ ${\mathrm{\&upsi;}}_i~N(0,{\mathrm{\&sigma;}}_e^2),$ ${\mathrm{\&sigma;}}_e^2=\frac{2{\mathrm{\&sigma;}}_i^2{S}^2}{L}\frac{y_i^2+y_l^2}{y_l^2{y}_i^2{(y_i-y_l)}^2}(i=\mathrm{2,3},...N),$ 这样的等式个数为N-1；

(3)假设g＝[‖r_s-c₂||²,||r_s-c₃||²,…||r_s-c_N||²]^T，

V＝[-υ₂,-υ₃,…-υ_N]^T，步骤(2)得到的等式写成矩阵形式，有P＝g+V，P的联合概率密度函数可以表示为 $f\left(P\right|r_s)={\left(2\mathrm{\&pi;}\right)}^{-(N-1)/2}{\left|\mathrm{\&Sigma;}\right|}^{-1/2}\exp \{-\frac{1}{2}{(P-g)}^T{\mathrm{\&Sigma;}}^{-1}(P-g)\},$ 其中，Σ＝E[VV^T]是协方差矩阵，元素可以表示为 ${\left[\mathrm{\&Sigma;}\right]}_{\mathrm{ij}}=\left\{\begin{array}{cc}\frac{2{\mathrm{\&sigma;}}^2{S}^2}{L}\frac{y_i^2+y_1^2}{y_i^2{y}_1^2{(y_i-y_1)}^2}& i=j\\ \frac{2{\mathrm{\&sigma;}}^2{S}^2}{L}\frac{y_i{y}_j}{y_1^2{y}_i{y}_j(y_i-y_1)(y_j-y_1)}& i\&NotEqual;j\end{array}\right.;$

(4)目标位置的最大似然估计可以表示为

s.t.g_i＝‖r_s-c_i||²，i＝2,3,…N，s.t.表示约束条件；

(5)步骤(4)中等式可以重新写完 $\underset{r_s,g}{\min }\left\{\begin{array}{c}\left[\begin{array}{cc}{g}^T& 1\end{array}\right]\end{array}\left[\begin{array}{cc}{\mathrm{\&Sigma;}}^{-1}& -{\mathrm{\&Sigma;}}^{-1}P\\ -P^T{\mathrm{\&Sigma;}}^{-1}& P^T{\mathrm{\&Sigma;}}^{-1}P\end{array}\right]\left[\begin{array}{c}g\\ 1\end{array}\right]\right\},$ s.t. $g_i=\left[\begin{array}{cc}{r_s}^T& 1\end{array}\right]\left[\begin{array}{cc}I& -c_i\\ {c_i}^T& {\left|\right|c_i\left|\right|}^2\end{array}\right]\left[\begin{array}{c}{r}_s\\ 1\end{array}\right];$

(6)利用矩阵等式X^TAX＝tr(XX^TA)，其中tr(·)表示矩阵的迹，将步骤(5)中等式写作 $\underset{r_s,g,G}{\min }\left\{\begin{array}{c}\mathrm{tr}\left(\begin{array}{c}\left[\begin{array}{cc}G& g\\ g^T& 1\end{array}\right]\left[\begin{array}{cc}{\mathrm{\&Sigma;}}^{-1}& -{\mathrm{\&Sigma;}}^{-1}P\\ -P^T{\mathrm{\&Sigma;}}^{-1}& P^T{\mathrm{\&Sigma;}}^{-1}P\end{array}\right]\end{array}\right)\end{array}\right\},$ s.t. $g_i=\mathrm{tr}\left(\begin{array}{c}\left[\begin{array}{cc}{R}_s& r_s\\ {r_s}^T& 1\end{array}\right]\end{array}\left[\begin{array}{cc}I& -c_i\\ {c_i}^T& {\left|\right|c_i\left|\right|}^2\end{array}\right]\right),$ R_s＝r_sr_s ^T，G＝gg^T；

(7)步骤(6)中得到的目标函数式凸的，但是约束条件不是凸的，采用SDR将非凸问题转化为凸问题： $\underset{r_s,g,G}{\min }\left\{\begin{array}{c}\mathrm{tr}\left(\begin{array}{c}\left[\begin{array}{cc}G& g\\ g^T& 1\end{array}\right]\left[\begin{array}{cc}{\mathrm{\&Sigma;}}^{-1}& -{\mathrm{\&Sigma;}}^{-1}P\\ -P^T{\mathrm{\&Sigma;}}^{-1}& P^T{\mathrm{\&Sigma;}}^{-1}P\end{array}\right]\end{array}\right)\end{array}\right\},$ s.t. $g_i=\mathrm{tr}\left(\begin{array}{c}\left[\begin{array}{cc}{R}_s& r_s\\ {r_s}^T& 1\end{array}\right]\end{array}\left[\begin{array}{cc}I& -c_i\\ {c_i}^T& {\left|\right|c_i\left|\right|}^2\end{array}\right]\right),$ $\left[\begin{array}{cc}{R}_s& r_s\\ {r_s}^T& 1\end{array}\right]\&GreaterEqual;0,$ $\left[\begin{array}{cc}G& g\\ g^T& 1\end{array}\right]\&GreaterEqual;0,$ 其中[·]≥0表示矩阵半正定。需要说明的是：由于目标辐射能量与增益乘积S未知，所以协方差阵Σ包含未知参数S，但是在求解优化问题时，常数项不会影响目标函数的最小化，因此在计算目标函数过程中，可以不必理会未知项S²；

(8)步骤(7)得到的等式及约束条件都是凸的，在MATLAB软件中可以直接采用凸优化工具箱来求解，从而得到目标位置的估计值。凸优化工具箱主要有：Michael Grant和Stephen Boyd的CVX、J.F.Sturm的SeDumi。

为了说明该发明基于水下声传感器网络的目标被动定位效果，分析对比了该发明与加权最小二乘被动目标被动定位精度，实施步骤如下：

(1)设置拓扑结构索引值Topo_index=1；设置蒙特卡洛仿真次数索引值Moto_index=1，index=1；

(2)设置传感器节点数目SNR=5×index；

(3)设置节点均匀分布在25×25×25三维区域内，传感器节点形成一次拓扑结构；

(4)在相同信噪比不同的观测噪声值情况下，循环执行仿真实施步骤(1)~(8)，得出N_m次蒙特卡洛目标被动定位结果；

(5)重复本步骤(2)~(3)，总共得到N_mN_t个目标被动定位结果，依据公式

求解出一个传感器数目条件下的目标被动定位均方误差；

(6)index＝index+1，判断：if index＞4，绘制不同传感器节点数目条件下目标被动定位的均方误差曲线，得到结果如图3所示，步骤结束，否则，重复执行步骤(2)~(5)。

实施实例中所述的加权最小二乘目标被动定位算法(OS_WLS)实施步骤为：

(1)N个传感器节点均匀分布在给定的三维区域内，依据公式

计算各个传感器节点接收到目标的能量值，获得能量矩阵E＝[y₁,y₂,…y_N]；

(2)选取矩阵E中最大能量值对应节点为参考节点l,联合第i个传感器节点接收到目标的能量y_i(i＝2,3,…N，N为传感器节点数目)和参考节点能量y₁，消去S，得到 ${\left|\right|r_s-c_i\left|\right|}^2={\mathrm{\&rho;}}_i^2+{\mathrm{\&upsi;}}_i,$ 其中， $c_i=\frac{y_i{r}_i-y_l{r}_l}{y_i-y_l},$ ${\mathrm{\&rho;}}_i^2=\frac{{\left|\right|r_i-r_l\left|\right|}^2}{{(y_i-y_l)}^2}{y}_i{y}_l\&ap;S\frac{{\mathrm{\&epsiv;}}_i{y}_l-{\mathrm{\&epsiv;}}_l{y}_i}{y_i{y}_l(y_i-y_l)}.$ ${\mathrm{\&upsi;}}_i~N(0,{\mathrm{\&sigma;}}_e^2),$ ${\mathrm{\&sigma;}}_e^2=\frac{2{\mathrm{\&sigma;}}_i^2{S}^2}{L}\frac{y_i^2+y_l^2}{y_l^2{y}_i^2{(y_i-y_l)}^2}(i=\mathrm{2,3},...N),$ 这样的等式个数为N-1；

(3)将步骤(2)中等式展开得到

其中，1≤m≤N且m≠l，l为最大能量值所对应节点，写成矩阵形式有Gθ_OS-LS＝c+υ，其中 $G=\left[\begin{array}{cc}-2c_2^T& 1\\ .& .\\ .& .\\ .& .\\ -2c_N^T& 1\end{array}\right],$ ${\mathrm{\&theta;}}_{\mathrm{OS}-\mathrm{LS}}=\left[\begin{array}{c}{r}_s\\ {\left|\right|r_s\left|\right|}^2\end{array}\right],$ υ＝[υ₂,υ₃…υ_N]^T；

(4)假设θ_OS-LS中r_s和||r_s||²相互独立，则步骤(3)中最小二乘目标被动定位的结果为θ_OS-LS＝(G^TG)^-1G^Tc，加权最小二乘目标被动定位的结果为θ_OS-WLS＝(G^TΣ^-1G)^-1G^TΣ^-1c，其中协方差阵Σ与该发明的协方差阵相同。

Claims (1)

1.一种基于水下声传感器网络的目标被动定位方法，其特征在于包括下述步骤：

步骤一、随机布放多个声传感器节点和融合中心形成一种并行的网络拓扑结构，即水下声传感器网络；

步骤二、声传感器节点采集目标辐射噪声数据，求取该节点接收到目标的能量，包括以下步骤：

(1)第i个传感器节点接收到的信号

其中，s(n)表示距离目标1m处的声强，γ_i是第i个传感器的接收增益，r_s表示目标位置，r_i表示第i个传感器位置，每个传感器观测噪声方差相同，s(n)与ω_i(n)相互独立；

(2)传感器节点接收到目标的能量为

实际求解过程中， $h_i=g_i\frac{S^{\&prime;}}{{\left|\right|r_s-r_i\left|\right|}^2}+\frac{1}{L}\underset{n=0}{\overset{L-1}{\mathrm{\&Sigma;}}}{\mathrm{\&omega;}}_i^2\left(n\right)=g_i\frac{S^{\&prime;}}{{\left|\right|r_s-r_i\left|\right|}^2}+{\mathrm{\&theta;}}_i,$ 其中

是信号能量，是噪声能量，S'＝s²(n)，

L是窗函数长度，θ_i服从

的卡方分布，当L＞＞30时，θ_i~N(σ²,2σ²/L)，各个传感器节点接收增益相同，则第i个传感器节点接收到目标的能量可以表示为

其中S＝S′×g_i；

(3)得到第i个传感器节点接收到的目标能量演化式

其中，

N个传感器节点接收到的目标能量演化式子组成能量矩阵E＝[y₁,y₂,…y_N]；

步骤三、各个声传感器节点将求解的能量分别传输至融合中心进行融合处理，依据能量衰减与距目标距离之间关系，对目标进行定位，包括以下步骤：

(1)选取矩阵E中最大能量值y₁对应节点为参考节点1，联合第i个传感器节点接收到的目标能量y_i，i＝2,3,…N，N为传感器节点数目，得到

其中， $c_i=\frac{y_i{r}_i-y_1{r}_1}{y_i-y_1},$ ${\mathrm{\&rho;}}_i^2=\frac{{\left|\right|r_i-r_1\left|\right|}^2}{{(y_i-y_1)}^2}{y}_i{y}_1,$ ${\mathrm{\&upsi;}}_i=\frac{{\mathrm{\&epsiv;}}_i{d}_i^2-{\mathrm{\&epsiv;}}_1{d}_1^2}{y_i-y_1}\&ap;S\frac{{\mathrm{\&epsiv;}}_i{y}_1-{\mathrm{\&epsiv;}}_1{y}_i}{y_i{y}_1(y_i-y_1)},$ ${\mathrm{\&upsi;}}_i~N(0,{\mathrm{\&sigma;}}_e^2),$ ${\mathrm{\&sigma;}}_e^2=\frac{2{\mathrm{\&sigma;}}_i^2{S}^2}{L}\frac{y_i^2+y_1^2}{y_1^2{y}_i^2{(y_i-y_1)}^2},$ 这样的等式个数为N-1；

(2)假设g＝[‖r_s-c₂‖²,‖r_s-c₃‖²,…||r_s-c_N‖²]^T，

V＝[-υ₂,-υ₃,…-υ_N]^T，上一步骤得到的等式写成矩阵形式，有P＝g+V，P的联合概率密度函数可以表示为 $f\left(P\right|r_s)={\left(2\mathrm{\&pi;}\right)}^{-(N-1)/2}{\left|\mathrm{\&Sigma;}\right|}^{-1/2}\exp \{-\frac{1}{2}{(P-g)}^T{\mathrm{\&Sigma;}}^{-1}(P-g)\},$ 其中，Σ＝E[VV^T]是协方差矩阵，元素可以表示为 ${\left[\mathrm{\&Sigma;}\right]}_{\mathrm{ij}}=\left\{\begin{array}{cc}\frac{2{\mathrm{\&sigma;}}^2{S}^2}{L}\frac{y_i^2+y_1^2}{y_i^2{y}_1^2{(y_i-y_1)}^2}& i=j\\ \frac{2{\mathrm{\&sigma;}}^2{S}^2}{L}\frac{y_i{y}_j}{y_1^2{y}_i{y}_j(y_i-y_1)(y_j-y_1)}& i\&NotEqual;j\end{array}\right.;$

(3)目标位置的最大似然估计为

(4)目标位置的最大似然估计可以表示为

s.t.g_i＝‖r_s-c_i‖²；

(5)上一步骤中的等式可以重新写作 $\underset{r_s,g}{\min }\left\{\begin{array}{c}\left[\begin{array}{cc}{g}^T& 1\end{array}\right]\end{array}\left[\begin{array}{cc}{\mathrm{\&Sigma;}}^{-1}& -{\mathrm{\&Sigma;}}^{-1}P\\ -P^T{\mathrm{\&Sigma;}}^{-1}& P^T{\mathrm{\&Sigma;}}^{-1}P\end{array}\right]\left[\begin{array}{c}g\\ 1\end{array}\right]\right\},$

$g_i=\left[\begin{array}{cc}{r_s}^T& 1\end{array}\right]\left[\begin{array}{cc}I& -c_i\\ {c_i}^T& {\left|\right|c_i\left|\right|}^2\end{array}\right]\left[\begin{array}{c}{r}_s\\ 1\end{array}\right];$

(6)利用矩阵等式X^TAX＝tr(XX^TA)，其中tr(·)表示矩阵的迹，将上一步骤中等式写作 $\underset{r_s,g,G}{\min }\left\{\begin{array}{c}\mathrm{tr}\left(\begin{array}{c}\left[\begin{array}{cc}G& g\\ g^T& 1\end{array}\right]\left[\begin{array}{cc}{\mathrm{\&Sigma;}}^{-1}& -{\mathrm{\&Sigma;}}^{-1}P\\ -P^T{\mathrm{\&Sigma;}}^{-1}& P^T{\mathrm{\&Sigma;}}^{-1}P\end{array}\right]\end{array}\right)\end{array}\right\},$ s.t. $g_i=\mathrm{tr}\left(\begin{array}{c}\left[\begin{array}{cc}{R}_s& r_s\\ {r_s}^T& 1\end{array}\right]\end{array}\left[\begin{array}{cc}I& -c_i\\ {c_i}^T& {\left|\right|c_i\left|\right|}^2\end{array}\right]\right),$ R_s＝r_sr_s ^T，G＝gg^T；

(7)采用SDR将非凸问题转化为凸问题： $\underset{r_s,g,G}{\min }\left\{\begin{array}{c}\mathrm{tr}\left(\begin{array}{c}\left[\begin{array}{cc}G& g\\ g^T& 1\end{array}\right]\left[\begin{array}{cc}{\mathrm{\&Sigma;}}^{-1}& -{\mathrm{\&Sigma;}}^{-1}P\\ -P^T{\mathrm{\&Sigma;}}^{-1}& P^T{\mathrm{\&Sigma;}}^{-1}P\end{array}\right]\end{array}\right)\end{array}\right\},$ s.t. $g_i=\mathrm{tr}\left(\begin{array}{c}\left[\begin{array}{cc}{R}_s& r_s\\ {r_s}^T& 1\end{array}\right]\end{array}\left[\begin{array}{cc}I& -c_i\\ {c_i}^T& {\left|\right|c_i\left|\right|}^2\end{array}\right]\right),$ $\left[\begin{array}{cc}{R}_s& r_s\\ {r_s}^T& 1\end{array}\right]\&GreaterEqual;0,$ $\left[\begin{array}{cc}G& g\\ g^T& 1\end{array}\right]\&GreaterEqual;0,$ 其中[·]≥0表示矩阵半正定；

(8)求解得到目标位置。

Images