Long-range infrasound source localization algorithms combining spherical pseudolinear estimator and propagation noise estimation
-
摘要:
球面伪线性估计器的提出为解决广域次声源定位问题提供了一个更加稳健和高效的方案, 然而次声远距离传播易受大气横向风的影响, 导致后方位角产生偏离, 使得球面伪线性估计器的定位结果和性能分析不准确。为此, 首先将次声传播引入噪声模型中, 重新分析了两种球面伪线性估计器的偏差和误差协方差, 给出了理论表达式; 然后, 给出了一种基于时间、空间相关的混合大气模型和射线追踪方法的传播噪声估计方法; 最后, 提出了球面伪线性估计器和传播噪声估计相结合的广域次声源定位算法。仿真和实验验证了理论分析的准确性和所提算法定位性能的优势。
Abstract:The spherical pseudolinear estimator offers a more stable and efficient solution to the long-range infrasound source localization. However, the propagation of infrasound over extended distances is vulnerable to atmospheric crosswinds, which cause deviations in the back-azimuth measurements. This phenomenon introduces inaccuracies in the localization results and performance evaluation of the spherical pseudolinear estimator. This research incorporates the effect of infrasound propagation into the noise model, reanalyzes the bias and error covariance performance of two spherical pseudolinear estimators, and provides theoretical expressions. Subsequently, a propagation noise estimation method is proposed that utilizes hybrid time-space dependent atmospheric modeling and ray tracing. Furthermore, long-range infrasound source localization algorithms are presented that integrate the spherical pseudolinear estimator with propagation noise estimation. The accuracy of the theoretical analysis and the performance advantages of the proposed algorithms are validated by simulations and an experiment.
-
Keywords:
- Infrasound /
- Sound source localization /
- Angle of arrival /
- Atmospheric model /
- Propagation noise
-
引言
地震、火山喷发等自然灾害, 爆炸、飞机音爆等人为事件均会产生人耳无法察觉的低频大气振动, 即次声。使用布设在国内和世界各地的监测台站采集次声信号, 可以实现对自然和人为事件的分析和定位。1996年《全面禁止核试验条约》的签署和全面禁止核试验条约组织 (CTBTO) 国际监测系统 (IMS) 次声部分的部署, 进一步促进了研究人员对次声监测技术的研究[1]。
次声源定位技术的研究具有悠久的历史。早在1912年Geiger便提出了一种非线性最小二乘迭代算法用于地震定位, 称为盖革法[2]。但这种方法直到上世纪六十年代才随着计算机的普及成为一种实用算法, 并经由Jordan和Sverdrup[3], Bratt和Bache[4]等的发展逐步应用到次声源定位之中。另一类次声源定位算法是网格搜索法。该方法的原理是将目标区域划分为若干网格, 逐点计算并从中选出与观测一致性最高的格点作为次声源位置。相关的方法包括贝叶斯次声源定位法 (BISL)[5,6]、蛮力网格搜索交叉定位法 (IMS_vASC)[7,8]、相似度法[9,10]、逆向时间迁移法 (RTM)[11]等。国内关于次声源定位的研究起步较晚。对于近场次声源, 吕君等[12]使用波束形成方法定位了北京一次地震前的异常次声波源, 对于远场次声源, 郭泉等[13,14]提出了一种最小方差法, 通过计算台站间时延的方差搜索次声源。此外, 赵久彬等[15]使用匹配场技术实现了对滑坡次声源的定位。最近, 文献[16]将广域次声源定位问题建模为球面波达角度 (AOA) 定位问题, 提出了一类球面伪线性估计器: 球面加权伪线性估计器 (SWPLE)和球面加权工具变量估计器 (SWIVE)。这类算法既解决了传统迭代定位算法的收敛性问题, 又具有极低的计算复杂度, 是解决广域次声源定位问题的有效方法。
次声的传播通道复杂时变, 尤其是在远距离传播过程中大气横向风会使传播路径偏离理论方向, 导致后方位角测量值偏离理论值, 使定位结果出现偏差[17-19]。为提升定位精度, 有许多研究已经将次声传播的影响纳入到定位算法之中。Ceranna等[20]基于美军海军实验室地对空 (NRL-G2S) 大气模型和有风大气声波传播 (WASP-3D) 射线追踪方法验证了修正后方位角测量对于盖革法定位结果的提升。Pilger等[21]基于Tau-P射线追踪方法计算了在基于美军海军实验室质谱仪非相干散射雷达−散逸层−00/水平风模型07 (NRLMSISE-00/HWM07)和欧洲中期天气预报中心 (ECMWF) 数据大气模型下的后方位角偏差, 用于盖革法的定位修正。Blom等[5]使用NRL-G2S大气模型的历史数据和几何声学 (GeoAc) 射线追踪方法计算了目标区域一年四季八个风向的后方位角偏差随声源到台站距离的统计模型, 用于BISL算法的后方位角修正。Arrowsmith等[6]则使用基于ECMWF数据的大气模型和特征射线 (Eigenray) 方法[19]搜索声源到各台站的特征射线, 从而确定次声传播的后方位角偏移。De Negri等[8]使用最新的NRLMSIS 2.0/HWM14和ECMWF数据构建大气剖面模型并使用次声几何声学 (InfraGA) 射线追踪方法迭代计算射线, 得到靠近台站位置的地面到达 (射线与地面的交点), 从而计算后方位角偏差来修正IMS_vASC算法。
目前一方面球面伪线性估计器只考虑了由传感器测量和参数估计产生的测量噪声, 并未考虑次声传播产生的后方位角偏移, 因此该类算法的定位结果和性能分析对于次声信号来说并不准确。另一方面现有方位角修正方法中主要存在以下限制和不足: (1) 部分方法仅使用经验气候学模型 (ECM) 和大气历史模型。(2) 绝大部分方法仅基于分层的时间空间无关的大气模型进行分析。(3) 绝大部分方法缺少对后方位角偏差二阶矩的表征, 因此难以对算法定位性能进行理论分析。针对上述问题, 本研究首先对球面波达角度定位问题的后方位角噪声模型进行了完善, 将由次声传播过程导致的后方位角偏移建模为传播噪声, 分析了在含有传播噪声的情况下SWPLE算法和SWIVE算法的偏差和误差协方差性能, 给出了理论表达式。然后设计了测量噪声和传播噪声的估计方法, 并将球面伪线性估计器和传播噪声估计相结合提出了SWPLE-prop算法和SWIVE-prop算法, 其中“prop”是“propagation noise estimation”的简写, 表示算法结合了传播噪声估计。数值仿真和实测数据证明SWIVE-prop算法可以有效提升次声事件的定位精度。
1. 问题构建
1.1 球面AOA定位模型
将广域次声源定位问题建模为球面AOA定位问题, 如图1所示。不失一般性地将地球建模为半径
R=1 的球体, 并在球面上建立忽略高程的大地坐标系, 图中NP表示北极。声源与台站均位于该球体表面, 声源与台站的坐标由纬度φ 和经度λ 表示。\boldsymbol{s}=\left[\varphi_s,\lambda_s\right]^{\text{T}} 表示未知的声源位置坐标,{{\boldsymbol{r}}_n} = {\left[ {{\varphi _n},\;{\lambda _n}} \right]^{\text{T}}} 表示已知的台站n \in \left\{ {1,2, \cdots ,N} \right\} 的位置坐标, 上标{\left( \cdot \right)^{\text{T}}} 表示转置。{\theta _n} \in \left( { - \pi ,\;\pi } \right] 为声源s 到达台站{{\boldsymbol{r}}_n} 的后方位角。后方位角以正北方向为零, 顺时针方向为正。球模型下的后方位角公式[22]为{\theta _n}\left( {\boldsymbol{s}} \right) = \arctan \left( {\frac{{{x_n}}}{{{y_n}}}} \right), (1) 其中
{x_n} = \cos {\varphi _s}\sin \left( {{\lambda _s} - {\lambda _n}} \right), (2a) {y_n} = \cos {\varphi _n}\sin {\varphi _s} - \sin {\varphi _n}\cos {\varphi _s}\cos \left( {{\lambda _s} - {\lambda _n}} \right), (2b) \arctan \left( \cdot \right) 为四象限反正切函数。1.2 噪声模型
将次声传播过程对后方位角测量的影响添加到噪声模型中:
{\widetilde{\boldsymbol \theta }} = {\boldsymbol{\theta }}\left( {\boldsymbol{s}} \right) + {\boldsymbol{\omega }} = {\boldsymbol{\theta }}\left( {\boldsymbol{s}} \right) + {\boldsymbol{e}} + {\boldsymbol{\varepsilon }}, (3) 其中
{\widetilde{\boldsymbol \theta }} = {\left[ {{{\widetilde \theta }_1},{{\widetilde \theta }_2}, \cdots ,{{\widetilde \theta }_N}} \right]^{\text{T}}}, (4a) {\boldsymbol{\theta }}\left( {\boldsymbol{s}} \right) = {\left[ {{\theta _1}\left( {\boldsymbol{s}} \right),{\theta _2}\left( {\boldsymbol{s}} \right), \cdots ,{\theta _N}\left( {\boldsymbol{s}} \right)} \right]^{\text{T}}}, (4b) {\boldsymbol{e}} = {\left[ {{e_1},{e_2}, \cdots ,{e_N}} \right]^{\text{T}}}, (4c) {\boldsymbol{\varepsilon }} = {\left[ {{\varepsilon _1},{\varepsilon _2}, \cdots ,{\varepsilon _N}} \right]^{\text{T}}}, (4d) {\boldsymbol{\omega }} = {\left[ {{\omega _1},{\omega _2}, \cdots ,{\omega _N}} \right]^{\text{T}}} = {\left[ {{e_1} + {\varepsilon _1},{e_2} + {\varepsilon _2}, \cdots ,{e_N} + {\varepsilon _N}} \right]^{\text{T}}}, (4e) {\widetilde{\boldsymbol \theta }} 表示台站处的后方位角测量值,{\boldsymbol{\omega }} 表示{\widetilde{\boldsymbol \theta }} 中包含的总噪声。{\boldsymbol{\omega }} 由两个分量组成, 其一是由传感器测量以及次声检测算法导致的测量噪声{\boldsymbol{e}} , 其二是次声传播过程中横向风等大气环境的影响, 本文中将其定义为传播噪声{\boldsymbol{\varepsilon }} 。不失一般性地, 本文假设传播噪声和测量噪声均服从相互独立的高斯分布:{\boldsymbol{e}}\sim N\left( {{\boldsymbol{0}},\;{{\boldsymbol{C}}_e}} \right), (5a) {\boldsymbol{\varepsilon }}\sim N\left( {{\boldsymbol{\mu }},\;{{\boldsymbol{C}}_\varepsilon }} \right), (5b) 其中
{\boldsymbol{\mu }} = {\left[ {{\mu _1},\,{\mu _2},\, \cdots ,\,{\mu _N}} \right]^{\text{T}}}, (6a) {{\boldsymbol{C}}_e} = {\text{diag}}\left( {\sigma _{e,1}^2,\sigma _{e,2}^2, \cdots ,\sigma _{e,N}^2} \right), (6b) {{\boldsymbol{C}}_\varepsilon } = {\text{diag}}\left( {\sigma _{\varepsilon ,1}^2,\sigma _{\varepsilon ,2}^2, \cdots ,\sigma _{\varepsilon ,N}^2} \right), (6c) {\boldsymbol{\mu }} 表示后方位角偏差,{{\boldsymbol{C}}_e} 和{{\boldsymbol{C}}_\varepsilon } 分别表示{\boldsymbol{e}} 和{\boldsymbol{\varepsilon }} 的协方差矩阵,{\text{diag}}\left( \cdot \right) 表示对角矩阵。于是, 总噪声服从高斯分布{\boldsymbol{\omega }}\sim N\left( {{\boldsymbol{\mu }},\;{{\boldsymbol{C}}_\omega }} \right), (7) 其中
{{\boldsymbol{C}}_\omega } = {\text{diag}}\left( {\sigma _{\omega ,1}^2,\sigma _{\omega ,2}^2, \cdots ,\sigma _{\omega ,N}^2} \right) = {{\boldsymbol{C}}_e} + {{\boldsymbol{C}}_\varepsilon }. (8) 出于便捷性的考虑, 在进行理论推导和分析时假设噪声参数
{\boldsymbol{\mu }} ,{{\boldsymbol{C}}_e} ,{{\boldsymbol{C}}_\varepsilon } 已知, 噪声参数的估计方法将在3.1节中给出。2. 存在传播噪声时的球面伪线性算法性能分析
由于噪声模型中添加了新的分量, 因此有必要对原有球面伪线性估计器的性能进行重新分析, 从而确定新的噪声分量对算法性能的影响。
2.1 SWPLE算法
当总噪声
{\boldsymbol{\omega }} 满足小噪声假设时, 球面AOA定位问题有如下伪线性关系[16]:{\widetilde{\boldsymbol Aq}} = {\boldsymbol{\eta }}, (9) 其中
{\widetilde{\boldsymbol A}} = {\left[ {{{{\widetilde{\boldsymbol a}}}_1},{{{\widetilde{\boldsymbol a}}}_2}, \cdots ,{{{\widetilde{\boldsymbol a}}}_N}} \right]^{\text{T}}}, (10a) {{\widetilde{\boldsymbol a}}_n} = \left[ {\begin{array}{*{20}{c}} { - \sin {{\widetilde \theta }_n}\sin {\varphi _n}\cos {\lambda _n} + \cos {{\widetilde \theta }_n}\sin {\lambda _n}} \\ { - \sin {{\widetilde \theta }_n}\sin {\varphi _n}\sin {\lambda _n} - \cos {{\widetilde \theta }_n}\cos {\lambda _n}} \\ {\sin {{\widetilde \theta }_n}\cos {\varphi _n}} \end{array}} \right], (10b) {\boldsymbol{q}} = {\left[ {\cos {\varphi _s}\cos {\lambda _s},\cos {\varphi _s}\sin {\lambda _s},\sin {\varphi _s}} \right]^{\text{T}}}, (10c) {\boldsymbol{\eta }} = {\left[ {{\eta _1},{\eta _2}, \cdots ,{\eta _N}} \right]^{\text{T}}}, (10d) {\eta _n} = \sin {\omega _n}\sin {\delta _n} \approx {\omega _n}\sin {\delta _n}, (10e) {\widetilde{\boldsymbol A}} 称为测量矩阵,{\boldsymbol{\eta }} 称为伪线性噪声向量,{\delta _n} 表示声源到台站的大圆弧长。SWPLE算法有目标矩阵:{{\boldsymbol{O}}_{{\text{SWPLE}}}} = {{\widetilde{\boldsymbol A}}^{\text{T}}}{{\boldsymbol{W}}^{ - 1}}{\widetilde{\boldsymbol A}}, (11) 其中, 上标
{\left( \cdot \right)^{ - 1}} 表示矩阵逆运算。由于SWPLE算法没有考虑传播噪声{\boldsymbol{\varepsilon }} , 所以权重矩阵{\boldsymbol{W}} 为{\boldsymbol{W}}{\text{ = }}\;{\text{diag}}\left( {\sigma _{e,1}^2{{\sin }^2}{\delta _1},\sigma _{e,2}^2{{\sin }^2}{\delta _2}, \cdots ,\sigma _{e,N}^2{{\sin }^2}{\delta _N}} \right), (12) 为了后续分析方便, 假设
{\boldsymbol{W}} 已知。2.1.1 偏差分析
当观测数 (探测站数) N足够大时SWPLE算法的偏差
{{\boldsymbol{\zeta }}_{{\text{SWPLE}}}} 为[16]\begin{split}& {{\boldsymbol{\zeta }}_{{\text{SWPLE}}}} = {\text{E}}\left\{ {{{{\boldsymbol{\widehat s}}}_{{\text{SWPLE}}}} - {\boldsymbol{s}}} \right\} = {\text{E}}\left\{ {\varDelta {{\boldsymbol{s}}_{{\text{SWPLE}}}}} \right\} \approx\\& \quad - {\boldsymbol{DU}}{\left( {{{\boldsymbol{U}}^{\text{T}}}{\text{E}}\left\{ {\frac{{{{{\widetilde{\boldsymbol A}}}^{\text{T}}}{{\boldsymbol{W}}^{ - 1}}{\widetilde{\boldsymbol A}}}}{N}} \right\}{\boldsymbol{U}}} \right)^{ - 1}}{{\boldsymbol{U}}^{\text{T}}}{\text{E}}\left\{ {\frac{{{{{\widetilde{\boldsymbol A}}}^{\text{T}}}{{\boldsymbol{W}}^{ - 1}}{\boldsymbol{\eta }}}}{N}} \right\}, \end{split} (13) 其中
{\boldsymbol{U}} = \left[ {\begin{array}{*{20}{c}} { - \sin {\lambda _s}}&{ - \sin {\varphi _s}\cos {\lambda _s}} \\ {\cos {\lambda _s}}&{ - \sin {\varphi _s}\sin {\lambda _s}} \\ 0&{\cos {\varphi _s}} \end{array}} \right], (14) {\boldsymbol{D}} = \left[ {\begin{array}{*{20}{c}} 0&0&{\dfrac{1}{{\cos {\varphi _s}}}} \\ { - \dfrac{{\sin {\lambda _s}}}{{\cos {\varphi _s}}}}&{\dfrac{{\cos {\lambda _s}}}{{\cos {\varphi _s}}}}&0 \end{array}} \right]. (15) 因为传播噪声均值不为零, 并且由于
{\widetilde{\boldsymbol A}} 和{\boldsymbol{\eta }} 都是总噪声{\boldsymbol{\omega }} 的函数,{\widetilde{\boldsymbol A}} 和{\boldsymbol{\eta }} 相关, 所以有{\text{E}}\left\{ {\frac{{{{{\widetilde{\boldsymbol A}}}^{\text{T}}}{{\boldsymbol{W}}^{ - 1}}{\boldsymbol{\eta }}}}{N}} \right\} \ne {\boldsymbol{0}}. (16) \boldsymbol{\zeta}_{\text{SWPLE}}\ne\boldsymbol{0}. (17) 因此SWPLE算法在存在传播噪声的环境下依然是有偏估计。
式(13)中第一项期望可以分解为
{\text{E}}\left\{ {\frac{{{{{\widetilde{\boldsymbol A}}}^{\text{T}}}{{\boldsymbol{W}}^{ - 1}}{\widetilde{\boldsymbol A}}}}{N}} \right\} = \frac{1}{N}\sum\limits_{n = 1}^N {\frac{{{\text{E}}\left\{ {{{{\widetilde{\boldsymbol a}}}_n}{\widetilde{\boldsymbol a}}_n^{\text{T}}} \right\}}}{{\sigma _{e,n}^2{{\sin }^2}{\delta _n}}}} . (18) 在满足小噪声假设时向量
{{\widetilde{\boldsymbol a}}_n} 可以分解为{{\widetilde{\boldsymbol a}}_n} \approx {{\boldsymbol{a}}_n} + {\omega _n}{{\boldsymbol{\xi }}_n}, (19) 其中
{{\boldsymbol{a}}_n} = \left[ {\begin{array}{*{20}{c}} { - \sin {\theta _n}\sin {\varphi _n}\cos {\lambda _n} + \cos {\theta _n}\sin {\lambda _n}} \\ { - \sin {\theta _n}\sin {\varphi _n}\sin {\lambda _n} - \cos {\theta _n}\cos {\lambda _n}} \\ {\sin {\theta _n}\cos {\varphi _n}} \end{array}} \right], (20a) {{\boldsymbol{\xi }}_n} = \left[ {\begin{array}{*{20}{c}} { - \cos {\theta _n}\sin {\varphi _n}\cos {\lambda _n} - \sin {\theta _n}\sin {\lambda _n}} \\ { - \cos {\theta _n}\sin {\varphi _n}\sin {\lambda _n} + \sin {\theta _n}\cos {\lambda _n}} \\ {\cos {\theta _n}\cos {\varphi _n}} \end{array}} \right]. (20b) 所以有
{\text{E}}\left\{ {{{{\widetilde{\boldsymbol a}}}_n}{\widetilde{\boldsymbol a}}_n^{\text{T}}} \right\} \approx {{\boldsymbol{a}}_n}{\boldsymbol{a}}_n^{\text{T}} + {\mu _n}{{\boldsymbol{a}}_n}{\boldsymbol{\xi }}_n^{\text{T}} + {\mu _n}{{\boldsymbol{\xi }}_n}{\boldsymbol{a}}_n^{\text{T}} + {m_{2,n}}{{\boldsymbol{\xi }}_n}{\boldsymbol{\xi }}_n^{\text{T}}, (21) 其中
{m_{2,n}} = {\text{E}}\left\{ {\omega _n^2} \right\} = \sigma _{\omega ,n}^2 + \mu _n^2. (22) \begin{split}& {\text{E}}\left\{ {\frac{{{{{\widetilde{\boldsymbol A}}}^{\text{T}}}{{\boldsymbol{W}}^{ - 1}}{\widetilde{\boldsymbol A}}}}{N}} \right\} \approx \frac{1}{N}\left( {{{\boldsymbol{A}}^{\text{T}}}{{\boldsymbol{W}}^{ - 1}}{\boldsymbol{A}} + {{\boldsymbol{A}}^{\text{T}}}{{\boldsymbol{M}}_1}{{\boldsymbol{W}}^{ - 1}}{\boldsymbol{\varXi }}} \right)+ \\ & \quad \frac{1}{N}\left( {{{\boldsymbol{\varXi }}^{\text{T}}}{{\boldsymbol{M}}_1}{{\boldsymbol{W}}^{ - 1}}{\boldsymbol{A}} + {{\boldsymbol{\varXi }}^{\text{T}}}{{\boldsymbol{M}}_2}{{\boldsymbol{W}}^{ - 1}}{\boldsymbol{\varXi }}} \right), \end{split} (23) 其中
{\boldsymbol{A}} = {\left[ {{{\boldsymbol{a}}_1},{{\boldsymbol{a}}_2}, \cdots ,{{\boldsymbol{a}}_N}} \right]^{\text{T}}}, (24a) {\boldsymbol{\varXi }} = {\left[ {{{\boldsymbol{\xi }}_1},{{\boldsymbol{\xi }}_2}, \cdots ,{{\boldsymbol{\xi }}_N}} \right]^{\text{T}}}, (24b) {{\boldsymbol{M}}_1} = {\text{diag}}\left( {{\mu _1},\,{\mu _2},\, \cdots ,\,{\mu _N}} \right), (24c) {{\boldsymbol{M}}_2} = {\text{diag}}\left( {{m_{2,1}},{m_{2,2}}, \cdots ,{m_{2,N}}} \right), (24d) {\boldsymbol{A}} 为没有噪声 ({\boldsymbol{\omega }} = {\boldsymbol{0}} ) 时的测量矩阵,{\boldsymbol{\varXi }} 称为偏差矩阵。式(13)中第二项期望可以分解为
{\text{E}}\left\{ {\frac{{{{{\widetilde{\boldsymbol A}}}^{\text{T}}}{{\boldsymbol{W}}^{ - 1}}{\boldsymbol{\eta }}}}{N}} \right\} = \frac{1}{N}\sum\limits_{n = 1}^N {\frac{{{\text{E}}\left\{ {{\eta _n}{{{\widetilde{\boldsymbol a}}}_n}} \right\}}}{{\sigma _{e,n}^2{{\sin }^2}{\delta _n}}}} . (25) {\text{E}}\left\{ {{\eta _n}{{{\widetilde{\boldsymbol a}}}_n}} \right\} \approx {\mu _n}\sin {\delta _n}{{\boldsymbol{a}}_n} + {m_{2,n}}\sin {\delta _n}{{\boldsymbol{\xi }}_n}, (26) {\text{E}}\left\{ {\frac{{{{{\widetilde{\boldsymbol A}}}^{\text{T}}}{{\boldsymbol{W}}^{ - 1}}{\boldsymbol{\eta }}}}{N}} \right\} \approx \frac{1}{N}\left( {{{\boldsymbol{A}}^{\text{T}}}\left( {{\boldsymbol{\mu }} \circ {{\boldsymbol{c}}_e} \circ {\boldsymbol{\delta }}} \right) + {{\boldsymbol{\varXi }}^{\text{T}}}\left( {{{\boldsymbol{m}}_2} \circ {{\boldsymbol{c}}_e} \circ {\boldsymbol{\delta }}} \right)} \right), (27) 其中
{\boldsymbol{\delta }} = {\left[ {{{\sin }^{ - 1}}{\delta _1},{{\sin }^{ - 1}}{\delta _2}, \cdots ,{{\sin }^{ - 1}}{\delta _N}} \right]^{\text{T}}}, (28a) {{\boldsymbol{c}}_e} = {\left[ {\sigma _{e,1}^{ - 2},\sigma _{e,2}^{ - 2}, \cdots ,\sigma _{e,N}^{ - 2}} \right]^{\text{T}}}, (28b) {{\boldsymbol{m}}_2} = {\left[ {{m_{2,1}},{m_{2,2}}, \cdots ,{m_{2,N}}} \right]^{\text{T}}}, (28c) 运算符
\circ 表示哈德玛 (Hadamard) 积。所以在满足小高斯噪声和一阶近似且当N足够大时, SWPLE算法的偏差{{\boldsymbol{\zeta }}_{{\text{SWPLE}}}} 可以近似表示为\begin{split} {{\boldsymbol{\zeta }}_{{\text{SWPLE}}}} \approx & - {\boldsymbol{DU}}{\left( {{{\boldsymbol{U}}^{\text{T}}}{{\boldsymbol{B}}_1}{\boldsymbol{U}}} \right)^{ - 1}}{{\boldsymbol{U}}^{\text{T}}}\left( {{\boldsymbol{A}}^{\text{T}}}\left( {{\boldsymbol{\mu }} \circ {{\boldsymbol{c}}_e} \circ {\boldsymbol{\delta }}} \right) + \right.\\ & \left. {{\boldsymbol{\varXi }}^{\text{T}}}\left( {{{\boldsymbol{m}}_2} \circ {{\boldsymbol{c}}_e} \circ {\boldsymbol{\delta }}} \right) \right), \end{split} (29) 其中
\begin{split} {{\boldsymbol{B}}_1} =& {{\boldsymbol{A}}^{\text{T}}}{{\boldsymbol{W}}^{ - 1}}{\boldsymbol{A}} + {{\boldsymbol{A}}^{\text{T}}}{{\boldsymbol{M}}_1}{{\boldsymbol{W}}^{ - 1}}{\boldsymbol{\varXi }} +\\& {{\boldsymbol{\varXi }}^{\text{T}}}{{\boldsymbol{M}}_1}{{\boldsymbol{W}}^{ - 1}}{\boldsymbol{A}} + {{\boldsymbol{\varXi }}^{\text{T}}}{{\boldsymbol{M}}_2}{{\boldsymbol{W}}^{ - 1}}{\boldsymbol{\varXi }}. \end{split} (30) 2.1.2 误差协方差分析
当N足够大时SWPLE算法的误差协方差矩阵为[16]
\begin{split} {{\boldsymbol{C}}_{{\text{SWPLE}}}} \approx &\frac{1}{N}{\boldsymbol{DU}}{\left( {{{\boldsymbol{U}}^{\text{T}}}{\text{E}}\left\{ {\frac{{{{{\widetilde{\boldsymbol A}}}^{\text{T}}}{{\boldsymbol{W}}^{ - 1}}{\widetilde{\boldsymbol A}}}}{N}} \right\}{\boldsymbol{U}}} \right)^{ - 1}}{{\boldsymbol{U}}^{\text{T}}}\cdot \\& \;\; {\text{E}}\left\{ {\frac{{{{{\widetilde{\boldsymbol A}}}^{\text{T}}}{{\boldsymbol{W}}^{ - 1}}{\boldsymbol{\eta }}{{\boldsymbol{\eta }}^{\text{T}}}{{\boldsymbol{W}}^{ - 1}}{\widetilde{\boldsymbol A}}}}{N}} \right\}\cdot \\& \;\; {\boldsymbol{U}}\,{\left( {{{\boldsymbol{U}}^{\text{T}}}{\text{E}}\left\{ {\frac{{{{{\widetilde{\boldsymbol A}}}^{\text{T}}}{{\boldsymbol{W}}^{ - 1}}{\widetilde{\boldsymbol A}}}}{N}} \right\}{\boldsymbol{U}}} \right)^{ - 1}}{{\boldsymbol{U}}^{\text{T}}}{{\boldsymbol{D}}^{\text{T}}}, \end{split} (31) 其中
{\text{E}}\left\{ {\frac{{{{{\widetilde{\boldsymbol A}}}^{\text{T}}}{{\boldsymbol{W}}^{ - 1}}{\boldsymbol{\eta }}{{\boldsymbol{\eta }}^{\text{T}}}{{\boldsymbol{W}}^{ - 1}}{\widetilde{\boldsymbol A}}}}{N}} \right\} = \frac{1}{N}\sum\limits_{n = 1}^N {\sum\limits_{i = 1}^N {\frac{{{\text{E}}\left\{ {{\omega _n}{\omega _i}{{{\widetilde{\boldsymbol a}}}_n}{\widetilde{\boldsymbol a}}_i^{\text{T}}} \right\}}}{{\sigma _{e,n}^2\sigma _{e,i}^2\sin {\delta _n}\sin {\delta _i}}}} } . (32) 当
n = i 时有\begin{split} {\text{E}}\left\{ {{\omega _n}{\omega _i}{{{\widetilde{\boldsymbol a}}}_n}{\widetilde{\boldsymbol a}}_i^{\text{T}}} \right\} \approx & {m_{2,n}}{{\boldsymbol{a}}_n}{\boldsymbol{a}}_n^{\text{T}} + {m_{3,n}}{{\boldsymbol{a}}_n}{\boldsymbol{\xi }}_n^{\text{T}} +\\& {m_{3,n}}{{\boldsymbol{\xi }}_n}{\boldsymbol{a}}_n^{\text{T}} + {m_{4,n}}{{\boldsymbol{\xi }}_n}{\boldsymbol{\xi }}_n^{\text{T}}, \end{split} (33) 其中
{m_{3,n}} = {\text{E}}\left\{ {\omega _n^3} \right\} = \sigma _{\omega ,n}^3 + 3\sigma _{\omega ,n}^2{\mu _n} + \mu _n^3, (34a) {m_{4,n}} = {\text{E}}\left\{ {\omega _n^4} \right\} = 3\sigma _{\omega ,n}^4 + 6\sigma _{\omega ,n}^2\mu _n^2 + \mu _n^4. (34b) 当
n \ne i 时有\begin{split} {\text{E}}\left\{ {{\omega _n}{\omega _i}{{{\widetilde{\boldsymbol a}}}_n}{\widetilde{\boldsymbol a}}_i^{\text{T}}} \right\} \approx & {\mu _n}{\mu _i}{{\boldsymbol{a}}_n}{\boldsymbol{a}}_i^{\text{T}} + {\mu _n}{m_{2,i}}{{\boldsymbol{a}}_n}{\boldsymbol{\xi }}_i^{\text{T}} +\\& {m_{2,n}}{\mu _i}{{\boldsymbol{\xi }}_n}{\boldsymbol{a}}_i^{\text{T}} + {m_{2,n}}{m_{2,i}}{{\boldsymbol{\xi }}_n}{\boldsymbol{\xi }}_i^{\text{T}}. \end{split} (35) \begin{split}& {\text{E}}\left\{ {\frac{{{{{\widetilde{\boldsymbol A}}}^{\text{T}}}{{\boldsymbol{W}}^{ - 1}}{\boldsymbol{\eta }}{{\boldsymbol{\eta }}^{\text{T}}}{{\boldsymbol{W}}^{ - 1}}{\widetilde{\boldsymbol A}}}}{N}} \right\} \approx \frac{1}{N}\left( {{{\boldsymbol{A}}^{\text{T}}}{{\boldsymbol{\varDelta }}_1}{\boldsymbol{A}} + {{\boldsymbol{A}}^{\text{T}}}{{\boldsymbol{\varDelta }}_2}{\boldsymbol{\varXi }}} \right) + \\& \qquad \;\; \frac{1}{N}\left( {{{\boldsymbol{\varXi }}^{\text{T}}}{\boldsymbol{\varDelta }}_2^{\text{T}}{\boldsymbol{A}} + {{\boldsymbol{\varXi }}^{\text{T}}}{{\boldsymbol{\varDelta }}_3}{\boldsymbol{\varXi }}} \right), \end{split} (36) 其中
{{\boldsymbol{\varDelta }}_1} = \left[ {\begin{array}{*{20}{c}} {\dfrac{{{m_{2,1}}}}{{\sigma _{e,1}^4{{\sin }^2}{\delta _1}}}}&{\dfrac{{{\mu _1}{\mu _2}}}{{\sigma _{e,1}^2\sigma _{e,2}^2\sin {\delta _1}\sin {\delta _2}}}}& \cdots &{\dfrac{{{\mu _1}{\mu _N}}}{{\sigma _{e,1}^2\sigma _{e,N}^2\sin {\delta _1}\sin {\delta _N}}}} \\ {\dfrac{{{\mu _2}{\mu _1}}}{{\sigma _{e,2}^2\sigma _{e,1}^2\sin {\delta _2}\sin {\delta _1}}}}&{\dfrac{{{m_{2,2}}}}{{\sigma _{e,2}^4{{\sin }^2}{\delta _2}}}}& \cdots &{\dfrac{{{\mu _2}{\mu _N}}}{{\sigma _{e,2}^2\sigma _{e,N}^2\sin {\delta _2}\sin {\delta _N}}}} \\ \vdots & \vdots & \ddots & \vdots \\ {\dfrac{{{\mu _N}{\mu _1}}}{{\sigma _{e,N}^2\sigma _{e,1}^2\sin {\delta _N}\sin {\delta _1}}}}&{\dfrac{{{\mu _N}{\mu _2}}}{{\sigma _{e,N}^2\sigma _{e,2}^2\sin {\delta _N}\sin {\delta _2}}}}& \cdots &{\dfrac{{{m_{2,N}}}}{{\sigma _{e,N}^4{{\sin }^2}{\delta _N}}}} \end{array}} \right], (37a) {{\boldsymbol{\varDelta }}_2} = \left[ {\begin{array}{*{20}{c}} {\dfrac{{{m_{3,1}}}}{{\sigma _{e,1}^4{{\sin }^2}{\delta _1}}}}&{\dfrac{{{\mu _1}{m_{2,2}}}}{{\sigma _{e,1}^2\sigma _{e,2}^2\sin {\delta _1}\sin {\delta _2}}}}& \cdots &{\dfrac{{{\mu _1}{m_{2,N}}}}{{\sigma _{e,1}^2\sigma _{e,N}^2\sin {\delta _1}\sin {\delta _N}}}} \\ {\dfrac{{{\mu _2}{m_{2,1}}}}{{\sigma _{e,2}^2\sigma _{e,1}^2\sin {\delta _2}\sin {\delta _1}}}}&{\dfrac{{{m_{3,2}}}}{{\sigma _{e,2}^4{{\sin }^2}{\delta _2}}}}& \cdots &{\dfrac{{{\mu _2}{m_{2,N}}}}{{\sigma _{e,2}^2\sigma _{e,N}^2\sin {\delta _2}\sin {\delta _N}}}} \\ \vdots & \vdots & \ddots & \vdots \\ {\dfrac{{{\mu _N}{m_{2,1}}}}{{\sigma _{e,N}^2\sigma _{e,1}^2\sin {\delta _N}\sin {\delta _1}}}}&{\dfrac{{{\mu _N}{m_{2,2}}}}{{\sigma _{e,N}^2\sigma _{e,2}^2\sin {\delta _N}\sin {\delta _2}}}}& \cdots &{\dfrac{{{m_{3,N}}}}{{\sigma _{e,N}^4{{\sin }^2}{\delta _N}}}} \end{array}} \right], (37b) {{\boldsymbol{\varDelta }}_3} = \left[ {\begin{array}{*{20}{c}} {\dfrac{{{m_{4,1}}}}{{\sigma _{e,1}^4{{\sin }^2}{\delta _1}}}}&{\dfrac{{{m_{2,1}}{m_{2,2}}}}{{\sigma _{e,1}^2\sigma _{e,2}^2\sin {\delta _1}\sin {\delta _2}}}}& \cdots &{\dfrac{{{m_{2,1}}{m_{2,N}}}}{{\sigma _{e,1}^2\sigma _{e,N}^2\sin {\delta _1}\sin {\delta _N}}}} \\ {\dfrac{{{m_{2,2}}{m_{2,1}}}}{{\sigma _{e,2}^2\sigma _{e,1}^2\sin {\delta _2}\sin {\delta _1}}}}&{\dfrac{{{m_{4,2}}}}{{\sigma _{e,2}^4{{\sin }^2}{\delta _2}}}}& \cdots &{\dfrac{{{m_{2,2}}{m_{2,N}}}}{{\sigma _{e,2}^2\sigma _{e,N}^2\sin {\delta _2}\sin {\delta _N}}}} \\ \vdots & \vdots & \ddots & \vdots \\ {\dfrac{{{m_{2,N}}{m_{2,1}}}}{{\sigma _{e,N}^2\sigma _{e,1}^2\sin {\delta _N}\sin {\delta _1}}}}&{\dfrac{{{m_{2,N}}{m_{2,2}}}}{{\sigma _{e,N}^2\sigma _{e,2}^2\sin {\delta _N}\sin {\delta _2}}}}& \cdots &{\dfrac{{{m_{4,N}}}}{{\sigma _{e,N}^4{{\sin }^2}{\delta _N}}}} \end{array}} \right]. (37c) 将式(23)和式(36)代入式(31), SWPLE算法的误差协方差矩阵在满足小高斯噪声假设和一阶近似且N足够大时可以近似为
{{\boldsymbol{C}}_{{\text{SWPLE}}}} \approx {\boldsymbol{DU}}{\left( {{{\boldsymbol{U}}^{\text{T}}}{{\boldsymbol{B}}_1}{\boldsymbol{U}}} \right)^{ - 1}}{{\boldsymbol{U}}^{\text{T}}}{{\boldsymbol{K}}_1}{\boldsymbol{U}}{\left( {{{\boldsymbol{U}}^{\text{T}}}{{\boldsymbol{B}}_1}{\boldsymbol{U}}} \right)^{ - 1}}{{\boldsymbol{U}}^{\text{T}}}{{\boldsymbol{D}}^{\text{T}}}, (38) 其中
\begin{split} &\\[-10pt]& {{\boldsymbol{K}}_1} = {{\boldsymbol{A}}^{\text{T}}}{{\boldsymbol{\varDelta }}_1}{\boldsymbol{A}} + {{\boldsymbol{A}}^{\text{T}}}{{\boldsymbol{\varDelta }}_2}{\boldsymbol{\varXi }} + {{\boldsymbol{\varXi }}^{\text{T}}}{\boldsymbol{\varDelta }}_2^{\text{T}}{\boldsymbol{A}} + {{\boldsymbol{\varXi }}^{\text{T}}}{{\boldsymbol{\varDelta }}_3}{\boldsymbol{\varXi }}. \end{split} (39) 2.2 SWIVE算法
SWIVE算法的目标矩阵为[16]
{{\boldsymbol{O}}_{{\text{SWIVE}}}} = {{\boldsymbol{G}}^{\text{T}}}{{\boldsymbol{W}}^{ - 1}}{\widetilde{\boldsymbol A}}. (40) IV矩阵
{\boldsymbol{G}} 由式(41)估计, 其中{\widehat \theta _n} 使用BCSWPLE算法[16]估计, 并采用选择角度测量 (SAM) 策略[23,16]:{\boldsymbol{G}} = {\boldsymbol{A}}\left( {{{\widehat \theta }_n}} \right), (41a) {\widehat \theta _n} = \left\{ {\begin{array}{ll} {{\theta _n}\left( {{{{\boldsymbol{\widehat s}}}_{{\text{BCSWPLE}}}}} \right),}&{\left| {{\theta _n}\left( {{{{\boldsymbol{\widehat s}}}_{{\text{BCSWPLE}}}}} \right) - {{\widetilde \theta }_n}} \right| < {\chi _n},} \\ {{{\widetilde \theta }_n},}&{\left| {{\theta _n}\left( {{{{\boldsymbol{\widehat s}}}_{{\text{BCSWPLE}}}}} \right) - {{\widetilde \theta }_n}} \right| \geqslant {\chi _n},} \end{array}} \right.\; (41b) 其中,
{\chi _n} = \rho {\sigma _n} 为SAM阈值,\rho 为SAM阈值系数。同样为了分析方便, 假设权重矩阵{\boldsymbol{W}} 已知。2.2.1 偏差分析
构造IV矩阵
{\boldsymbol{G}} 的后方位角{\widehat \theta _n} 可以表示为{\widehat \theta _n} = {\theta _n}\left( {\boldsymbol{s}} \right) + {\upsilon _n}. (42) 对于绝大多数情况, 有
{\widehat \theta _n} = {\theta _n}\left( {{{{\boldsymbol{\widehat s}}}_{{\text{BCSWPLE}}}}} \right) , 因此有{\upsilon _n} = {\theta _n}\left( {{{{\boldsymbol{\widehat s}}}_{{\text{BCSWPLE}}}}} \right) - {\theta _n}\left( {\boldsymbol{s}} \right). (43) 对
{\theta _n}\left( {{{{\boldsymbol{\widehat s}}}_{{\text{BCSWPLE}}}}} \right) 在{{\boldsymbol{\widehat s}}_{{\text{BCSWPLE}}}} = {\boldsymbol{s}} 处进行泰勒展开并忽略二次以上项, 可以得到\begin{split} {\upsilon _n} \approx &\frac{{\partial {\theta _n}\left( {\boldsymbol{s}} \right)}}{{\partial {{\boldsymbol{s}}^{\text{T}}}}}\left( {{{{\boldsymbol{\widehat s}}}_{{\text{BCSWPLE}}}} - {\boldsymbol{s}}} \right) \approx\\& \frac{{\partial {\theta _n}\left( {\boldsymbol{s}} \right)}}{{\partial {{\boldsymbol{s}}^{\text{T}}}}}\left[ {\left( {{{{\boldsymbol{\widehat s}}}_{{\text{SWPLE}}}} - {{{\boldsymbol{\zeta '}}}_{{\text{SWPLE}}}}} \right) - {\boldsymbol{s}}} \right] = \\& \frac{{\partial {\theta _n}\left( {\boldsymbol{s}} \right)}}{{\partial {{\boldsymbol{s}}^{\text{T}}}}}\left( {\varDelta {{\boldsymbol{s}}_{{\text{SWPLE}}}} - {{{\boldsymbol{\zeta '}}}_{{\text{SWPLE}}}}} \right), \end{split} (44) 其中
\dfrac{{\partial {\theta _n}\left( {\boldsymbol{s}} \right)}}{{\partial {{\boldsymbol{s}}^{\text{T}}}}} = \left[ {\dfrac{{\dfrac{{\partial {x_n}}}{{\partial {\varphi _s}}}{y_n} - {x_n}\dfrac{{\partial {y_n}}}{{\partial {\varphi _s}}}}}{{{x_n}^2 + {y_n}^2}},\dfrac{{\dfrac{{\partial {x_n}}}{{\partial {\lambda _s}}}{y_n} - {x_n}\dfrac{{\partial {y_n}}}{{\partial {\lambda _s}}}}}{{{x_n}^2 + {y_n}^2}}} \right], (45) \frac{{\partial {x_n}}}{{\partial {\varphi _s}}} = - \sin {\varphi _s}\sin \left( {{\lambda _s} - {\lambda _n}} \right), (46a) \frac{{\partial {y_n}}}{{\partial {\varphi _s}}} = \cos {\varphi _n}\cos {\varphi _s} + \sin {\varphi _n}\sin {\varphi _s}\cos \left( {{\lambda _s} - {\lambda _n}} \right), (46b) \frac{{\partial {x_n}}}{{\partial {\lambda _s}}} = \cos {\varphi _s}\cos \left( {{\lambda _s} - {\lambda _n}} \right), (46c) \frac{{\partial {y_n}}}{{\partial {\lambda _s}}} = \sin {\varphi _n}\cos {\varphi _s}\sin \left( {{\lambda _s} - {\lambda _n}} \right), (46d) {{\boldsymbol{\zeta '}}_{{\text{SWPLE}}}} 表示只存在测量噪声时SWPLE算法的偏差。当N足够大时SWIVE算法的偏差{{\boldsymbol{\zeta }}_{{\text{SWIVE}}}} 为[16]{{\boldsymbol{\zeta }}_{{\text{SWIVE}}}} \approx - {\boldsymbol{DU}}{\left( {{{\boldsymbol{U}}^{\text{T}}}{\text{E}}\left\{ {\frac{{{{\boldsymbol{G}}^{\text{T}}}{{\boldsymbol{W}}^{ - 1}}{\widetilde{\boldsymbol A}}}}{N}} \right\}{\boldsymbol{U}}} \right)^{ - 1}}{{\boldsymbol{U}}^{\text{T}}}{\text{E}}\left\{ {\frac{{{{\boldsymbol{G}}^{\text{T}}}{{\boldsymbol{W}}^{ - 1}}{\boldsymbol{\eta }}}}{N}} \right\}. (47) 由式(44)和式(10e)可以看出虽然根据
{\widehat \theta _n} 得到的IV矩阵{\boldsymbol{G}} 和伪线性噪声向量{\boldsymbol{\eta }} 的相关性很小, 但由于总噪声{\boldsymbol{\omega }} 的均值不为零, 因此依然有{\text{E}}\left\{ {\frac{{{{\boldsymbol{G}}^{\text{T}}}{{\boldsymbol{W}}^{ - 1}}{\boldsymbol{\eta }}}}{N}} \right\} \ne {\boldsymbol{0}}, (48) {{\boldsymbol{\zeta }}_{{\text{SWIVE}}}} \ne {\boldsymbol{0}}, (49) 因此存在传播噪声时SWIVE算法不再具有渐近无偏性, 而成为有偏估计。
式(47)中第一项期望可以分解为
{\text{E}}\left\{ {\frac{{{{\boldsymbol{G}}^{\text{T}}}{{\boldsymbol{W}}^{ - 1}}{\widetilde{\boldsymbol A}}}}{N}} \right\} = \frac{1}{N}\sum\limits_{n = 1}^N {\frac{{{\text{E}}\left\{ {{{\boldsymbol{g}}_n}{\widetilde{\boldsymbol a}}_n^{\text{T}}} \right\}}}{{\sigma _{e,n}^2{{\sin }^2}{\delta _n}}}} . (50) 根据式(42), 在满足小噪声假设时, 向量
{{\boldsymbol{g}}_n} 可以分解为{{\boldsymbol{g}}_n} \approx {{\boldsymbol{a}}_n} + {\upsilon _n}{{\boldsymbol{\xi }}_n}. (51) 而N足够大时,
{\boldsymbol{G}} ,{\boldsymbol{\eta }} 的相关性很小, 故有{\text{E}}\left\{ {{\omega _n}{\upsilon _n}} \right\} \approx {\text{E}}\left\{ {{\omega _n}} \right\}{\text{E}}\left\{ {{\upsilon _n}} \right\} , 则{\text{E}}\left\{ {{{\boldsymbol{g}}_n}{\widetilde{\boldsymbol a}}_n^{\text{T}}} \right\} \approx {{\boldsymbol{a}}_n}{\boldsymbol{a}}_n^{\text{T}} + {\mu _n}{{\boldsymbol{a}}_n}{\boldsymbol{\xi }}_n^{\text{T}} + {\nu _n}{{\boldsymbol{\xi }}_n}{\boldsymbol{a}}_n^{\text{T}} + {\mu _n}{\nu _n}{{\boldsymbol{\xi }}_n}{\boldsymbol{\xi }}_n^{\text{T}}, (52) 其中
{\nu _n} = {\text{E}}\left\{ {{\upsilon _n}} \right\} \approx \frac{{\partial {\theta _n}\left( {\boldsymbol{s}} \right)}}{{\partial {{\boldsymbol{s}}^{\text{T}}}}}{{\boldsymbol{\zeta }}_{{\text{BCSWPLE}}}}, (53) \begin{split}& {{\boldsymbol{\zeta }}_{{\text{BCSWPLE}}}} = {\text{E}}\left\{ {{{{\boldsymbol{\widehat s}}}_{{\text{BCSWPLE}}}} - {\boldsymbol{s}}} \right\} \approx {{\boldsymbol{\zeta }}_{{\text{SWPLE}}}} - {{{\boldsymbol{\zeta '}}}_{{\text{SWPLE}}}} \approx \\&\quad - {\boldsymbol{DU}}{\left( {{{\boldsymbol{U}}^{\text{T}}}{{\boldsymbol{B}}_1}{\boldsymbol{U}}} \right)^{ - 1}}{{\boldsymbol{U}}^{\text{T}}}\left( {{{\boldsymbol{A}}^{\text{T}}}\left( {{\boldsymbol{\mu }} \circ {{\boldsymbol{c}}_e} \circ {\boldsymbol{\delta }}} \right) + {{\boldsymbol{\varXi }}^{\text{T}}}\left( {{{\boldsymbol{m}}_2} \circ {{\boldsymbol{c}}_e} \circ {\boldsymbol{\delta }}} \right)} \right) + \\ &\;\; \;\;{\boldsymbol{DU}}{\left[ {{{\boldsymbol{U}}^{\text{T}}}\left( {{{\boldsymbol{A}}^{\text{T}}}{{\boldsymbol{W}}^{ - 1}}{\boldsymbol{A}} + {{\boldsymbol{\varXi }}^{\text{T}}}{{\boldsymbol{C}}_e}{{\boldsymbol{W}}^{ - 1}}{\boldsymbol{\varXi }}} \right){\boldsymbol{U}}} \right]^{ - 1}}{{\boldsymbol{U}}^{\text{T}}}{{\boldsymbol{\varXi }}^{\text{T}}}{\boldsymbol{\delta }}. \end{split} (54) \begin{split} {\text{E}}\left\{ {\frac{{{{\boldsymbol{G}}^{\text{T}}}{{\boldsymbol{W}}^{ - 1}}{\widetilde{\boldsymbol A}}}}{N}} \right\} \approx & \frac{1}{N}\left( {{{\boldsymbol{A}}^{\text{T}}}{{\boldsymbol{W}}^{ - 1}}{\boldsymbol{A}} + {{\boldsymbol{A}}^{\text{T}}}{{\boldsymbol{M}}_1}{{\boldsymbol{W}}^{ - 1}}{\boldsymbol{\varXi }}} \right) + \\& \frac{1}{N}\left( {{{\boldsymbol{\varXi }}^{\text{T}}}{\boldsymbol{N}}{{\boldsymbol{W}}^{ - 1}}{\boldsymbol{A}} + {{\boldsymbol{\varXi }}^{\text{T}}}{{\boldsymbol{M}}_1}{\boldsymbol{N}}{{\boldsymbol{W}}^{ - 1}}{\boldsymbol{\varXi }}} \right), \end{split} (55) 其中
{\boldsymbol{N}} = {\text{diag}}\left( {{\nu _1},{\nu _2}, \cdots ,{\nu _N}} \right). (56) 式(47)中第二项期望可以分解为
{\text{E}}\left\{ {\frac{{{{\boldsymbol{G}}^{\text{T}}}{{\boldsymbol{W}}^{ - 1}}{\boldsymbol{\eta }}}}{N}} \right\} = \frac{1}{N}\sum\limits_{n = 1}^N {\frac{{{\text{E}}\left\{ {{\eta _n}{{\boldsymbol{g}}_n}} \right\}}}{{\sigma _{e,n}^2{{\sin }^2}{\delta _n}}}} . (57) {\text{E}}\left\{ {{\eta _n}{{\boldsymbol{g}}_n}} \right\} \approx {\mu _n}\sin {\delta _n}{{\boldsymbol{a}}_n} + {\mu _n}{\nu _n}\sin {\delta _n}{{\boldsymbol{\xi }}_n}. (58) {\text{E}}\left\{ {\frac{{{{\boldsymbol{G}}^{\text{T}}}{{\boldsymbol{W}}^{ - 1}}{\boldsymbol{\eta }}}}{N}} \right\} \approx \frac{1}{N}\left( {{{\boldsymbol{A}}^{\text{T}}}\left( {{\boldsymbol{\mu }} \circ {{\boldsymbol{c}}_e} \circ {\boldsymbol{\delta }}} \right) + {{\boldsymbol{\varXi }}^{\text{T}}}\left( {{\boldsymbol{\mu }} \circ {\boldsymbol{\nu }} \circ {{\boldsymbol{c}}_e} \circ {\boldsymbol{\delta }}} \right)} \right), (59) 其中
{\boldsymbol{\nu }} = {\left[ {{\nu _1},{\nu _2}, \cdots ,{\nu _N}} \right]^{\text{T}}}, (60) 所以在满足小高斯噪声假设和一阶近似且N足够大时, SWIVE算法在含有传播噪声的环境中的偏差
{{\boldsymbol{\zeta }}_{{\text{SWIVE}}}} 可以近似表示为\begin{split} {{\boldsymbol{\zeta }}_{{\text{SWIVE}}}} \approx & - {\boldsymbol{DU}}{\left( {{{\boldsymbol{U}}^{\text{T}}}{{\boldsymbol{B}}_2}{\boldsymbol{U}}} \right)^{ - 1}}{{\boldsymbol{U}}^{\text{T}}}\left( {{\boldsymbol{A}}^{\text{T}}}\left( {{\boldsymbol{\mu }} \circ {{\boldsymbol{c}}_e} \circ {\boldsymbol{\delta }}} \right) +\right. \\&\left. {{\boldsymbol{\varXi }}^{\text{T}}}\left( {{\boldsymbol{\mu }} \circ {\boldsymbol{\nu }} \circ {{\boldsymbol{c}}_e} \circ {\boldsymbol{\delta }}} \right) \right), \end{split} (61) 其中
\begin{split} {{\boldsymbol{B}}_2} = &{{\boldsymbol{A}}^{\text{T}}}{{\boldsymbol{W}}^{ - 1}}{\boldsymbol{A}} + {{\boldsymbol{A}}^{\text{T}}}{{\boldsymbol{M}}_1}{{\boldsymbol{W}}^{ - 1}}{\boldsymbol{\varXi }} + \\ &{{\boldsymbol{\varXi }}^{\text{T}}}{\boldsymbol{N}}{{\boldsymbol{W}}^{ - 1}}{\boldsymbol{A}} + {{\boldsymbol{\varXi }}^{\text{T}}}{{\boldsymbol{M}}_1}{\boldsymbol{N}}{{\boldsymbol{W}}^{ - 1}}{\boldsymbol{\varXi }}. \end{split} (62) 2.2.2 误差协方差分析
当N足够大时SWIVE算法的误差协方差矩阵为[16]
\begin{split} {{\boldsymbol{C}}_{{\text{SWIVE}}}} \approx &\frac{1}{N}{\boldsymbol{DU}}{\left( {{{\boldsymbol{U}}^{\text{T}}}{\text{E}}\left\{ {\frac{{{{\boldsymbol{G}}^{\text{T}}}{{\boldsymbol{W}}^{ - 1}}{\widetilde{\boldsymbol A}}}}{N}} \right\}{\boldsymbol{U}}} \right)^{ - 1}}{{\boldsymbol{U}}^{\text{T}}} \cdot \\& {\text{E}}\left\{ {\frac{{{{\boldsymbol{G}}^{\text{T}}}{{\boldsymbol{W}}^{ - 1}}{\boldsymbol{\eta }}{{\boldsymbol{\eta }}^{\text{T}}}{{\boldsymbol{W}}^{ - 1}}{\boldsymbol{G}}}}{N}} \right\} \cdot \\& {\boldsymbol{U}}\,{\left( {{{\boldsymbol{U}}^{\text{T}}}{\text{E}}\left\{ {\frac{{{{{\widetilde{\boldsymbol A}}}^{\text{T}}}{{\boldsymbol{W}}^{ - 1}}{\boldsymbol{G}}}}{N}} \right\}{\boldsymbol{U}}} \right)^{ - 1}}{{\boldsymbol{U}}^{\text{T}}}{{\boldsymbol{D}}^{\text{T}}}. \end{split} (63) 式(63)中的第二项期望可以分解为
{\text{E}}\left\{ {\frac{{{{\boldsymbol{G}}^{\text{T}}}{{\boldsymbol{W}}^{ - 1}}{\boldsymbol{\eta }}{{\boldsymbol{\eta }}^{\text{T}}}{{\boldsymbol{W}}^{ - 1}}{\boldsymbol{G}}}}{N}} \right\} = \frac{1}{N}\sum\limits_{n = 1}^N {\sum\limits_{i = 1}^N {\frac{{{\text{E}}\left\{ {{\omega _n}{\omega _i}{{\boldsymbol{g}}_n}{\boldsymbol{g}}_i^{\text{T}}} \right\}}}{{\sigma _{e,n}^2\sigma _{e,i}^2\sin {\delta _n}\sin {\delta _i}}}} } . (64) 当
n = i 时有\begin{split} {\text{E}}\left\{ {{\omega _n}{\omega _i}{{\boldsymbol{g}}_n}{\boldsymbol{g}}_i^{\text{T}}} \right\} \approx & {m_{2,n}}{{\boldsymbol{a}}_n}{\boldsymbol{a}}_n^{\text{T}} + {m_{2,n}}{\nu _n}{{\boldsymbol{a}}_n}{\boldsymbol{\xi }}_n^{\text{T}} + \\& {m_{2,n}}{\nu _n}{{\boldsymbol{\xi }}_n}{\boldsymbol{a}}_n^{\text{T}} + {m_{2,n}}\nu _n^2{{\boldsymbol{\xi }}_n}{\boldsymbol{\xi }}_n^{\text{T}}, \end{split} (65) 当
n \ne i 时有\begin{split} {\text{E}}\left\{ {{\omega _n}{\omega _i}{{\boldsymbol{g}}_n}{\boldsymbol{g}}_i^{\text{T}}} \right\} \approx &{\mu _n}{\mu _i}{{\boldsymbol{a}}_n}{\boldsymbol{a}}_i^{\text{T}} + {\mu _n}{\mu _i}{\nu _i}{{\boldsymbol{a}}_n}{\boldsymbol{\xi }}_i^{\text{T}}+ \\ & {\mu _n}{\mu _i}{\nu _n}{{\boldsymbol{\xi }}_n}{\boldsymbol{a}}_i^{\text{T}} + {\mu _n}{\mu _i}{\nu _n}{\nu _i}{{\boldsymbol{\xi }}_n}{\boldsymbol{\xi }}_i^{\text{T}}. \end{split} (66) \begin{split} {\text{E}}\left\{ {\frac{{{{\boldsymbol{G}}^{\text{T}}}{{\boldsymbol{W}}^{ - 1}}{\boldsymbol{\eta }}{{\boldsymbol{\eta }}^{\text{T}}}{{\boldsymbol{W}}^{ - 1}}{\boldsymbol{G}}}}{N}} \right\} \approx &\frac{1}{N}\left( {{{\boldsymbol{A}}^{\text{T}}}{{\boldsymbol{\varDelta }}_1}{\boldsymbol{A}} + {{\boldsymbol{A}}^{\text{T}}}{{\boldsymbol{\varDelta }}_4}{\boldsymbol{\varXi }}} \right)+ \\ & \frac{1}{N}\left( {{{\boldsymbol{\varXi }}^{\text{T}}}{\boldsymbol{\varDelta }}_4^{\text{T}}{\boldsymbol{A}} + {{\boldsymbol{\varXi }}^{\text{T}}}{{\boldsymbol{\varDelta }}_5}{\boldsymbol{\varXi }}} \right), \end{split} (67) 其中
{{\boldsymbol{\varDelta }}_4} = \left[ {\begin{array}{*{20}{c}} {\dfrac{{{m_{2,1}}{\nu _1}}}{{\sigma _{e,1}^4{{\sin }^2}{\delta _1}}}}&{\dfrac{{{\mu _1}{\mu _2}{\nu _2}}}{{\sigma _{e,1}^2\sigma _{e,2}^2\sin {\delta _1}\sin {\delta _2}}}}& \cdots &{\dfrac{{{\mu _1}{\mu _N}{\nu _N}}}{{\sigma _{e,1}^2\sigma _{e,N}^2\sin {\delta _1}\sin {\delta _N}}}} \\ {\dfrac{{{\mu _2}{\mu _1}{\nu _1}}}{{\sigma _{e,2}^2\sigma _{e,1}^2\sin {\delta _2}\sin {\delta _1}}}}&{\dfrac{{{m_{2,2}}{\nu _2}}}{{\sigma _{e,2}^4{{\sin }^2}{\delta _2}}}}& \cdots &{\dfrac{{{\mu _2}{\mu _N}{\nu _N}}}{{\sigma _{e,2}^2\sigma _{e,N}^2\sin {\delta _2}\sin {\delta _N}}}} \\ \vdots & \vdots & \ddots & \vdots \\ {\dfrac{{{\mu _N}{\mu _1}{\nu _1}}}{{\sigma _{e,N}^2\sigma _{e,1}^2\sin {\delta _N}\sin {\delta _1}}}}&{\dfrac{{{\mu _N}{\mu _2}{\nu _2}}}{{\sigma _{e,N}^2\sigma _{e,2}^2\sin {\delta _N}\sin {\delta _2}}}}& \cdots &{\dfrac{{{m_{2,N}}{\nu _N}}}{{\sigma _{e,N}^4{{\sin }^2}{\delta _N}}}} \end{array}} \right], (68a) {{\boldsymbol{\varDelta }}_5} = \left[ {\begin{array}{*{20}{c}} {\dfrac{{{m_{2,1}}\nu _1^2}}{{\sigma _{e,1}^4{{\sin }^2}{\delta _1}}}}&{\dfrac{{{\mu _1}{\mu _2}{\nu _1}{\nu _2}}}{{\sigma _{e,1}^2\sigma _{e,2}^2\sin {\delta _1}\sin {\delta _2}}}}& \cdots &{\dfrac{{{\mu _1}{\mu _N}{\nu _1}{\nu _N}}}{{\sigma _{e,1}^2\sigma _{e,N}^2\sin {\delta _1}\sin {\delta _N}}}} \\ {\dfrac{{{\mu _2}{\mu _1}{\nu _2}{\nu _1}}}{{\sigma _{e,2}^2\sigma _{e,1}^2\sin {\delta _2}\sin {\delta _1}}}}&{\dfrac{{{m_{2,2}}\nu _2^2}}{{\sigma _{e,2}^4{{\sin }^2}{\delta _2}}}}& \cdots &{\dfrac{{{\mu _2}{\mu _N}{\nu _2}{\nu _N}}}{{\sigma _{e,2}^2\sigma _{e,N}^2\sin {\delta _2}\sin {\delta _N}}}} \\ \vdots & \vdots & \ddots & \vdots \\ {\dfrac{{{\mu _N}{\mu _1}{\nu _N}{\nu _1}}}{{\sigma _{e,N}^2\sigma _{e,1}^2\sin {\delta _N}\sin {\delta _1}}}}&{\dfrac{{{\mu _N}{\mu _2}{\nu _N}{\nu _2}}}{{\sigma _{e,N}^2\sigma _{e,2}^2\sin {\delta _N}\sin {\delta _2}}}}& \cdots &{\dfrac{{{m_{2,N}}\nu _N^2}}{{\sigma _{e,N}^4{{\sin }^2}{\delta _N}}}} \end{array}} \right]. (68b) 将式(55)和式(67)代入式(63), SWIVE算法的误差协方差矩阵在满足小高斯噪声假设和一阶近似且N足够大时可以近似为
{{\boldsymbol{C}}_{{\text{SWIVE}}}} \approx {\boldsymbol{DU}}{\left( {{{\boldsymbol{U}}^{\text{T}}}{{\boldsymbol{B}}_2}{\boldsymbol{U}}} \right)^{ - 1}}{{\boldsymbol{U}}^{\text{T}}}{{\boldsymbol{K}}_2}{\boldsymbol{U}}{\left( {{{\boldsymbol{U}}^{\text{T}}}{\boldsymbol{B}}_2^{\text{T}}{\boldsymbol{U}}} \right)^{ - 1}}{{\boldsymbol{U}}^{\text{T}}}{{\boldsymbol{D}}^{\text{T}}}, (69) 其中
{{\boldsymbol{K}}_2} = {{\boldsymbol{A}}^{\text{T}}}{{\boldsymbol{\varDelta }}_1}{\boldsymbol{A}} + {{\boldsymbol{A}}^{\text{T}}}{{\boldsymbol{\varDelta }}_4}{\boldsymbol{\varXi }} + {{\boldsymbol{\varXi }}^{\text{T}}}{\boldsymbol{\varDelta }}_4^{\text{T}}{\boldsymbol{A}} + {{\boldsymbol{\varXi }}^{\text{T}}}{{\boldsymbol{\varDelta }}_5}{\boldsymbol{\varXi }}. (70) 与克拉美劳下界 (Cramér–Rao Lower Bound, CRLB) [16]进行比较可以看出, SWIVE算法的误差协方差无法达到CRLB, 因此SWIVE算法在传播噪声下失去渐近有效性。需要注意的是SAM阈值系数
\rho 的取值会影响SWIVE算法的性能, 而文献[16]中对于\rho 最优取值的分析并未考虑传播噪声。因此存在传播噪声时, SWIVE算法可能无法达到理论定位性能。所以式(61)和式(69)实际上是SWIVE算法定位性能的下界。3. 结合传播噪声估计的球面伪线性估计器
为了解决传播噪声对广域次声源定位造成的不利影响, 本节首先给出了噪声的参数估计方法, 然后提出了将球面伪线性估计器与传播噪声估计相结合的广域次声源定位算法。
3.1 噪声参数估计
3.1.1 测量噪声参数估计
次声检测算法使用PMCC算法[24,25], PMCC检测结果为一组或多组PMCC族。PMCC族由多个PMCC像素组成, 每个像素的属性包括时间、频率、后方位角和视速度。将次声到达信号PMCC族中的所有像素作为样本, 即可得到该次声信号的后方位角测量值
{\widetilde \theta _n} 和测量噪声的方差估计\widetilde \sigma _{e,n}^2 :{\widetilde \theta _n} = \frac{1}{{{Q_n}}}\sum\limits_{i = 1}^{{Q_n}} {{\theta _{n,i}}} , (71a) \widetilde \sigma _{e,n}^2 = \frac{1}{{{Q_n} - 1}}\sum\limits_{i = 1}^{{Q_n}} {{{\left( {{\theta _{n,i}} - {{\widetilde \theta }_n}} \right)}^2}} , (71b) 其中,
{Q_n} 表示PMCC族n中像素的个数,{\theta _{n,i}} 表示每个像素的后方位角测量值。以相同的方法可以估计视速度的测量值{\widetilde v_{{\text{t}},n}} 和其方差估计\widetilde \sigma _{{v_{\text{t}}},n}^2 :{\widetilde v_{{\text{t}},n}} = \frac{1}{{{Q_n}}}\sum\limits_{i = 1}^{{Q_n}} {{v_{{\text{t}},n,i}}} , (72a) \widetilde \sigma _{{v_{\text{t}}},n}^2 = \frac{1}{{{Q_n} - 1}}\sum\limits_{i = 1}^{{Q_n}} {{{\left( {{v_{{\text{t}},n,i}} - {{\widetilde v}_{{\text{t}},n}}} \right)}^2}} , (72b) 其中,
{v_{{\text{t}},n,i}} 表示每个像素的视速度测量值。{\widetilde \theta _n} ,\widetilde \sigma _{e,n}^2 ,{\widetilde v_{{\text{t}},n}} ,\widetilde \sigma _{{v_{\text{t}}},n}^2 将用于传播噪声参数的估计。3.1.2 传播噪声参数估计
传播噪声的参数估计需要构建大气剖面模型并使用射线追踪方法确定各次声台站处的后方位角偏差。为了最大限度保证射线追踪结果的准确性, 本研究构建了基于ECM和大气再分析数据混合且时间空间相关的大气剖面模型。大气再分析数据是不均匀的气象观测和基于模型的天气预报经过数据同化得到的全球范围时间空间统一的气象数据。本研究ECM使用NRLMSIS2.0/HWM14经验气候学模型, 大气再分析数据使用ECMWF的ERA5数据[26], 该数据可以提供最高0.01 hPa (约80 km), 时间分辨率1 h, 空间分辨率31 km的大气数据。
首先使用SWIVE算法确定声源的初步位置
{{\boldsymbol{\widehat s}}_{{\text{SWIVE}}}} , 结合次声台站地理位置确定射线追踪的经纬度范围并按\varDelta \varphi 和\varDelta \lambda 的间隔划分网格。然后假设次声传播的平均速度为c , 使用最先接收到信号的台站{{\boldsymbol{r}}_1} 的次声到达时间{t_1} (PMCC族中最早的时间) 估计声源的发生时间{\widehat t_s} 和次声信号到达各个格点的时间{\widehat t_{\varphi ,\lambda }} :{\widehat t_s} = {t_1} - \frac{{H\left( {{{{\boldsymbol{\widehat s}}}_{{\text{SWIVE}}}},{{\boldsymbol{r}}_1},{R_{\text{e}}}} \right)}}{c}, (73) {\widehat t_{\varphi ,\lambda }} = {\widehat t_s} + \frac{{H\left( {{{{\boldsymbol{\widehat s}}}_{{\text{SWIVE}}}},{{\left[ {\varphi ,\,\lambda } \right]}^{\text{T}}},{R_{\text{e}}}} \right)}}{c}, (74) 其中
\begin{split}& H\left( {{{\boldsymbol{r}}_i},{{\boldsymbol{r}}_j},R} \right) = \\& 2R\arcsin \left( {\sqrt {{{\sin }^2}\left( {\frac{{{\varphi _j} - {\varphi _i}}}{2}} \right) + \cos {\varphi _j}\cos {\varphi _i}{{\sin }^2}\left( {\frac{{{\lambda _j} - {\lambda _i}}}{2}} \right)} } \right). \end{split} (75) 函数
H\left( { \cdot ,\; \cdot , \cdot } \right) 是半正矢球面距离公式[27], 球体半径取地球的平均半径{R_{\text{e}}} = 6371393 m。此时就得到了数据下载的时间和空间范围。分别按经纬度网格下载对应的ERA5数据并在目标区域内线性插值得到\vartheta _{\varphi ,\lambda ,t,h}^{{\text{ERA5}}} , 以及读取ECM数据\vartheta _{\varphi ,\lambda ,t,h}^{{\text{ECM}}} , 其中大气剖面数据\vartheta 包括高度、温度、纬向风、经向风、总质量密度和气压。然后将\vartheta _{\varphi ,\lambda ,t,h}^{{\text{ECM}}} 和\vartheta _{\varphi ,\lambda ,t,h}^{{\text{ERA5}}} 进行融合得到融合数据\vartheta _{\varphi ,\lambda ,t,h}^{{\text{Hybrid}}} :\begin{split}& \vartheta _{\varphi ,\lambda ,t,h}^{{\text{Hybrid}}} = \\& \left\{ {\begin{array}{*{20}{l}} {\vartheta _{\varphi ,\lambda ,t,h}^{{\text{ECM}}},}&{{h_{\max }} < h,} \\ {\left( {\dfrac{{{h_{\max }} - h}}{{{h_{\max }} - {h_{\min }}}}} \right)\vartheta _{\varphi ,\lambda ,t,h}^{{\text{ERA5}}} + \left( {\dfrac{{h - {h_{\min }}}}{{{h_{\max }} - {h_{\min }}}}} \right)\vartheta _{\varphi ,\lambda ,t,h}^{{\text{ECM}}},}&{{h_{\min }} < h \leqslant {h_{\max }},} \\ {\vartheta _{\varphi ,\lambda ,t,h}^{{\text{ERA5}}},}&{h \leqslant {h_{\min }}.} \end{array}} \right. \end{split} (76) 此时得到了不同时间的空间相关大气剖面。使用与各个格点处次声到达时间
{\widehat t_{\varphi ,\lambda }} 最相近的两组大气剖面按时间线性插值, 即可得到时间空间相关的混合大气剖面模型。使用InfraGA射线追踪方法计算次声传播过程, 从声源估计位置
{{\boldsymbol{\widehat s}}_{{\text{SWIVE}}}} 向各台站方向{\psi _n} \pm {\psi _{{\text{lim}}}} 的方位角范围和\left[ {{\phi _{\min }},\;{\phi _{\max }}} \right] 的俯仰角范围, 以\varDelta \psi 的方位角间隔和\varDelta \phi 的俯仰角间隔计算次声射线, 获得地面到达信号参数, 参数包括经纬度、后方位角、俯仰角等。根据到达信号与台站到声源的距离差\varDelta d 、后方位角\theta 和视速度{v_{\text{t}}} 比较射线追踪得到的地面到达信号和实际观测信号的相似度。台站n处射线地面到达信号j与观测信号的相似度定义为{\alpha _{n,j}} = \sqrt {{{\left( {\frac{{\varDelta {d_{n,j}}}}{{{\sigma _{\varDelta d,n}}}}} \right)}^2} + {{\left( {\frac{{{\theta _{n,j}} - {{\widetilde \theta }_n}}}{{{{\widetilde \sigma }_{e,n}}}}} \right)}^2} + {{\left( {\frac{{{v_{{\text{t}},n,j}} - {{\widetilde v}_{{\text{t}},n}}}}{{{{\widetilde \sigma }_{{v_{\text{t}}},n}}}}} \right)}^2}} , (77) 其中
\varDelta {d_{n,j}} = H\left( {{{{\boldsymbol{\widehat s}}}_{{\text{SWIVE}}}},{{\left[ {{\varphi _j},{\lambda _j}} \right]}^{\text{T}}},{R_{\text{e}}}} \right) - H\left( {{{{\boldsymbol{\widehat s}}}_{{\text{SWIVE}}}},{{\boldsymbol{r}}_n},{R_{\text{e}}}} \right), (78a) {\sigma _{\varDelta d,n}} = 0.05 \cdot H\left( {{{{\boldsymbol{\widehat s}}}_{{\text{SWIVE}}}},{{\boldsymbol{r}}_n},{R_{\text{e}}}} \right), (78b) {v_{{\text{t}},n,j}} = \frac{{340}}{{\cos {\phi _{n,j}}}}. (78c) 最后选取最相似的P个到达信号作为样本计算传播噪声均值
{\widetilde \mu _n} 和方差\widetilde \sigma _{\varepsilon ,n}^2 :{\widetilde \mu _n} = \frac{1}{P}\sum\limits_{j = 1}^P {\left( {{\theta _{n,j}} - {\theta _n}\left( {{{{\boldsymbol{\widehat s}}}_{{\text{SWIVE}}}}} \right)} \right)} , (79a) \widetilde \sigma _{\varepsilon ,n}^2 = \frac{1}{{P - 1}}\sum\limits_{j = 1}^P {{{\left( {{\theta _{n,j}} - {\theta _n}\left( {{{{\boldsymbol{\widehat s}}}_{{\text{SWIVE}}}}} \right) - {{\widetilde \mu }_n}} \right)}^2}} . (79b) 3.2 SWPLE-prop算法
在加入传播噪声的噪声模型下, 后方位角测量值的似然函数[28]为
\begin{split} p\left( {{\widetilde{\boldsymbol \theta }}\,;\,{\boldsymbol{s}}} \right) = &\frac{1}{{{{\left( {2\pi } \right)}^{{N \mathord{\left/ {\vphantom {N 2}} \right. } 2}}}\det {{\left( {{{\boldsymbol{C}}_\omega }} \right)}^{{1 \mathord{\left/ {\vphantom {1 2}} \right. } 2}}}}}\cdot \\[2mm] & \exp \left( { - \frac{1}{2}{{\left( {\widetilde {\boldsymbol{\theta }}- {\boldsymbol{\theta}} \left( {\boldsymbol{s}} \right) - {\boldsymbol{\mu }}} \right)}^{\text{T}}}C_\omega ^{ - 1}\left( {\widetilde {\boldsymbol{\theta}} - {\boldsymbol{\theta}} \left( {\boldsymbol{s}} \right) - {\boldsymbol{\mu }}} \right)} \right), \end{split} (80) 其中,
\det \left( \cdot \right) 表示矩阵行列式。则MLE目标函数为\begin{split} {J_{{\text{MLE}}}}\left( {\boldsymbol{s}} \right) =& \frac{1}{2}{\left( {\widetilde {\boldsymbol{\theta}} - {\boldsymbol{\theta}} \left( {\boldsymbol{s}} \right) - {\boldsymbol{\mu }}} \right)^{\text{T}}}C_\omega ^{ - 1}\left( {\widetilde {\boldsymbol{\theta}} - {\boldsymbol{\theta }}\left( s \right) - {\boldsymbol{\mu }}} \right) = \\[1mm]& \sum\limits_{n = 1}^N {\frac{1}{{2\sigma _{\omega ,n}^2}}{{\left( {\left( {{{\widetilde \theta }_n} - {\mu _n}} \right) - {\theta _n}\left( s \right)} \right)}^2}} . \end{split} (81) 在小噪声假设下, 根据SWPLE算法的推导过程[16]可以得到新的目标函数:
{J_{{\text{SWPLE-prop}}}}\left( {\boldsymbol{q}} \right) = \frac{1}{2}{\left( {{\widetilde{\boldsymbol A'{\boldsymbol{q}}}}} \right)^{\text{T}}}{{\boldsymbol{W'}}^{ - 1}}\left( {{\widetilde{\boldsymbol A'{\boldsymbol{q}}}}} \right)\quad {\text{s}}{\text{.t}}{\text{.}}\;{{\boldsymbol{q}}^{\text{T}}}{\boldsymbol{q}} = 1, (82) 其中
{\widetilde{\boldsymbol A'}} = {\left[ {{{{\widetilde{\boldsymbol a'}}}_1},{{{\widetilde{\boldsymbol a'}}}_2}, \cdots ,{{{\widetilde{\boldsymbol a'}}}_N}} \right]^{\text{T}}}, (83a) {{\widetilde{\boldsymbol a'}}_n} = \left[ {\begin{array}{*{20}{c}} { - \sin \left( {{{\widetilde \theta }_n} - {\mu _n}} \right)\sin {\varphi _n}\cos {\lambda _n} + \cos \left( {{{\widetilde \theta }_n} - {\mu _n}} \right)\sin {\lambda _n}} \\ { - \sin \left( {{{\widetilde \theta }_n} - {\mu _n}} \right)\sin {\varphi _n}\sin {\lambda _n} - \cos \left( {{{\widetilde \theta }_n} - {\mu _n}} \right)\cos {\lambda _n}} \\ {\sin \left( {{{\widetilde \theta }_n} - {\mu _n}} \right)\cos {\varphi _n}} \end{array}} \right], (83b) {\boldsymbol{W'}} = {\text{diag}}\left( {\sigma _{\omega ,1}^2{{\sin }^2}{\delta _1},\sigma _{\omega ,2}^2{{\sin }^2}{\delta _2}, \cdots ,\sigma _{\omega ,N}^2{{\sin }^2}{\delta _N}} \right). (84) 此时即可使用SWPLE算法进行求解。将该算法命名为SWPLE-prop算法, 其目标矩阵为
{\boldsymbol{O}}_{{\mathrm{SWPLE-prop}}}={\widetilde {\boldsymbol{A}}}^{\mathrm{'T}} {\widehat{\boldsymbol{W}}}'^{-1}{\widetilde{\boldsymbol{A}}}', (85) 其中
{\boldsymbol{\widehat W'}} = {\boldsymbol{W'}}\left( {{{\widehat \delta }_n}\left( {{{\boldsymbol{r}}_n},\;{{{\boldsymbol{\widehat s}}}_{{\text{SPLE - prop}}}}} \right)} \right). (86) {\boldsymbol{\widehat W'}} 的构建使用简化版SPLE-prop算法, 其目标矩阵为{{\boldsymbol{O}}_{{\text{SPLE - prop}}}} = {{\widetilde{\boldsymbol A'}}^{\text{T}}}{\boldsymbol{C}}_\omega ^{ - 1}{\widetilde{\boldsymbol A'}}. (87) 求解球面伪线性估计器首先需要对目标矩阵
{\boldsymbol{O}} 进行特征值分解 (EVD) :{\boldsymbol{\widehat q}} = \frac{{{{\boldsymbol{v}}_{\min }}\left( {\boldsymbol{O}} \right)}}{{\left| {{{\boldsymbol{v}}_{\min }}\left( {\boldsymbol{O}} \right)} \right|}}, (88) 其中,
{{\boldsymbol{v}}_{\min }} 为{{\boldsymbol{O}}_{{\text{SPLE - prop}}}} 最小特征值对应的特征向量。然后再将三维笛卡尔坐标估计结果{\boldsymbol{\widehat q}} 转换为大地坐标{\boldsymbol{\widehat s}} :{\boldsymbol{\widehat s}} = {\left[ {\arcsin {{\widehat q}_3},\arctan \left( {\frac{{{{\widehat q}_2}}}{{{{\widehat q}_1}}}} \right)} \right]^{\text{T}}}. (89) 最后进行对跖点判断[16]:
\;{{\boldsymbol{\widehat s}}_{[{\text{alg}}]}} = \left\{ {\begin{array}{*{20}{c}} {{\boldsymbol{\widehat s}},}&{\left| {{\theta _n}\left( {{\boldsymbol{\widehat s}}} \right) - {{\widetilde \theta }_n}} \right| < {\chi _\varDelta },} \\ {{{{\boldsymbol{\widehat s}}}_{{\text{ap}}}},}&{\left| {{\theta _n}\left( {{\boldsymbol{\widehat s}}} \right) - {{\widetilde \theta }_n}} \right| \geqslant {\chi _\varDelta },} \end{array}} \right. (90) 其中,
{\chi _\varDelta } 表示方向一致性阈值, 通常可设置为{\pi \mathord{\left/ {\vphantom {\pi 2}} \right. } 2} ,{{\boldsymbol{\widehat s}}_{{\text{ap}}}} 表示对跖点, 下标{\left( \cdot \right)_{[{\text{alg}}]}} 表示不同的算法名称。根据文献[16]的分析方法可知SWPLE-prop算法是有偏估计, 当N足够大时算法的偏差为\begin{split}& {{\boldsymbol{\zeta }}_{{\text{SWPLE-prop}}}} \approx \\& - {\boldsymbol{DU}}{\left[ {{{\boldsymbol{U}}^{\text{T}}}\left( {{{\boldsymbol{A}}^{\text{T}}}{{{\boldsymbol{W'}}}^{ - 1}}{\boldsymbol{A}} + {{\boldsymbol{\varXi }}^{\text{T}}}{{\boldsymbol{C}}_\omega }{{{\boldsymbol{W'}}}^{ - 1}}{\boldsymbol{\varXi }}} \right){\boldsymbol{U}}} \right]^{ - 1}}{{\boldsymbol{U}}^{\text{T}}}{{\boldsymbol{\varXi }}^{\text{T}}}{\boldsymbol{\delta }}. \end{split} (91) 算法的误差协方差矩阵为
\begin{split} {{\boldsymbol{C}}_{{\text{SWPLE - prop}}}} \approx & {\boldsymbol{DU}}{\left[ {{{\boldsymbol{U}}^{\text{T}}}\left( {{{\boldsymbol{A}}^{\text{T}}}{{{\boldsymbol{W'}}}^{ - 1}}{\boldsymbol{A}} + {{\boldsymbol{\varXi }}^{\text{T}}}{{\boldsymbol{C}}_\omega }{{{\boldsymbol{W'}}}^{ - 1}}{\boldsymbol{\varXi }}} \right){\boldsymbol{U}}} \right]^{ - 1}}\cdot \\& \;\; {{\boldsymbol{U}}^{\text{T}}}\left( {{{\boldsymbol{A}}^{\text{T}}}{{{\boldsymbol{W'}}}^{ - 1}}{\boldsymbol{A}} + {{\boldsymbol{\varXi }}^{\text{T}}}{\boldsymbol{\varDelta \varXi }}} \right){\boldsymbol{U}} \cdot \\ & \;\; {\left[ {{{\boldsymbol{U}}^{\text{T}}}\left( {{{\boldsymbol{A}}^{\text{T}}}{{{\boldsymbol{W'}}}^{ - 1}}{\boldsymbol{A}} + {{\boldsymbol{\varXi }}^{\text{T}}}{{\boldsymbol{C}}_\omega }{{{\boldsymbol{W'}}}^{ - 1}}{\boldsymbol{\varXi }}} \right){\boldsymbol{U}}} \right]^{ - 1}}{{\boldsymbol{U}}^{\text{T}}}{{\boldsymbol{D}}^{\text{T}}}, \end{split} (92) 其中
{\boldsymbol{\varDelta }} = \left[ {\begin{array}{*{20}{c}} {3{{\sin }^{ - 2}}{\delta _1}} & {{{\sin }^{ - 1}}{\delta _1}{{\sin }^{ - 1}}{\delta _2}}& \cdots & {{{\sin }^{ - 1}}{\delta _1}{{\sin }^{ - 1}}{\delta _N}} \\ {{{\sin }^{ - 1}}{\delta _2}{{\sin }^{ - 1}}{\delta _1}} & {3{{\sin }^{ - 2}}{\delta _2}}& \cdots & {{{\sin }^{ - 1}}{\delta _2}{{\sin }^{ - 1}}{\delta _N}} \\ \vdots & \vdots & \ddots & \vdots \\ {{{\sin }^{ - 1}}{\delta _N}{{\sin }^{ - 1}}{\delta _1}} & {{{\sin }^{ - 1}}{\delta _N}{{\sin }^{ - 1}}{\delta _2}} & \cdots & {3{{\sin }^{ - 2}}{\delta _N}} \end{array}} \right]. (93) SWPLE-prop算法总结见表1。
表 1 SWPLE-prop算法流程获取 {{\boldsymbol{r}}_n}, n \in \left\{ {1,2, \cdots ,N} \right\} 根据3.1.1节测量噪声参数估计方法得到 \widetilde \theta , \widetilde \sigma _{e,n}^2
使用SWIVE算法估计初始位置 {{\boldsymbol{\widehat s}}_{{\text{SWIVE}}}}
根据3.1.2节传播噪声参数估计方法得到{\widetilde \mu _n}, \widetilde \sigma _{\varepsilon ,n}^2
构造 {\widetilde{\boldsymbol A'}} (式(83)), {{\boldsymbol{C}}_\omega } (式(8)), {{\boldsymbol{O}}_{{\text{SPLE - prop}}}}(式(87))
计算 {{\boldsymbol{O}}_{{\text{SPLE - prop}}}}的EVD
根据式(88)—式(90)得到 {{\boldsymbol{\widehat s}}_{{\text{SPLE - prop}}}}
构造 {\boldsymbol{\widehat W'}} (式(86)), {{\boldsymbol{O}}_{{\text{SWPLE - prop}}}}(式(85))
计算 {{\boldsymbol{O}}_{{\text{SWPLE - prop}}}}的EVD
根据式(88)—式(90)得到 {{\boldsymbol{\widehat s}}_{{\text{SWPLE - prop}}}}3.3 SWIVE-prop算法
因为SWPLE-prop算法是有偏估计, 所以根据文献[16], 使用SWIVE-prop算法来解决偏差问题。SWIVE-prop算法的目标矩阵为
{\boldsymbol{O}}_{\mathrm{SWIVE-prop}} ={\boldsymbol{G}'}^{\mathrm{T}} {\widehat {\boldsymbol{W}}}'^{-1}\widetilde {\boldsymbol{A}}', (94) 式中
{\boldsymbol{G'}} = {\boldsymbol{A}}\left( {{{\widehat \theta '}_n}} \right), (95a) {\widehat \theta '_n} = \left\{ {\begin{array}{*{20}{l}} {{\theta _n}\left( {{{{\boldsymbol{\widehat s}}}_{{\text{BCSWPLE-prop}}}}} \right),} & {\left| {{\theta _n}\left( {{{{\boldsymbol{\widehat s}}}_{{\text{BCSWPLE - prop}}}}} \right) - \left( {{{\widetilde \theta }_n} - {\mu _n}} \right)} \right| < {\chi _n},} \\ {{{\widetilde \theta }_n} - {\mu _n},} & {\left| {{\theta _n}\left( {{{{\boldsymbol{\widehat s}}}_{{\text{BCSWPLE - prop}}}}} \right) - \left( {{{\widetilde \theta }_n} - {\mu _n}} \right)} \right| \geqslant {\chi _n},} \end{array}} \right.\; (95b) {\boldsymbol{\widehat W'}} = {\boldsymbol{W}}\left( {{{\widehat \delta }_n}\left( {{{\boldsymbol{r}}_n},\;{{{\boldsymbol{\widehat s}}}_{{\text{BCSWPLE - prop}}}}} \right)} \right), (96) 其中,
{\chi _n} = \rho {\sigma _n} 为SAM阈值,\rho 为SAM阈值系数。{\boldsymbol{\widehat W'}} 和{\boldsymbol{G'}} 的构建需要使用BCSWPLE-prop算法, 其目标矩阵为{{\boldsymbol{O}}_{{{{\mathrm{BCSWPLE-prop}}}}}} = \widetilde {\boldsymbol{A}}'^{\mathrm{T}}\widehat {\boldsymbol{W}}^{\mathrm{'-1}}\widetilde {\boldsymbol{A}}'-\widehat{\boldsymbol{\varXi}} '^{\mathrm{T}}{\boldsymbol{C}}_\omega \widehat {\boldsymbol{W}}^{\mathrm{'-1}}\widehat{\boldsymbol{\varXi}}', (97) 其中
{\boldsymbol{\widehat \varXi '}} = {\boldsymbol{\varXi }}\left( {{\theta _n}\left( {{{{\boldsymbol{\widehat s}}}_{{\text{SPLE - prop}}}}} \right)} \right). (98) 根据文献[16]的分析方法可知SWIVE-prop算法是渐近无偏和渐近有效估计, 即算法偏差为
{{\boldsymbol{\zeta }}_{{\text{SWIVE - prop}}}} = {\boldsymbol{0}}\quad {\text{as}}\;N \to \infty , (99) 误差协方差矩阵为
{{\boldsymbol{C}}_{{\text{SWIVE - prop}}}} \approx {\boldsymbol{DU}}{\left( {{{\boldsymbol{U}}^{\text{T}}}{{\boldsymbol{A}}^{\text{T}}}{{{\boldsymbol{W'}}}^{ - 1}}{\boldsymbol{AU}}} \right)^{ - 1}}{{\boldsymbol{U}}^{\text{T}}}{{\boldsymbol{D}}^{\text{T}}}. (100) SWIVE-prop算法总结见表2。
表 2 SWIVE-prop算法流程获取 {{\boldsymbol{r}}_n}, n \in \left\{ {1,2, \cdots ,N} \right\} 根据3.1.1节测量噪声参数估计方法得到 \widetilde \theta , \widetilde \sigma _{e,n}^2
使用SWIVE算法估计初始位置 {{\boldsymbol{\widehat s}}_{{\text{SWIVE}}}}
根据3.1.2节传播噪声参数估计方法得到{\widetilde \mu _n}, \widetilde \sigma _{\varepsilon ,n}^2
构造 {\widetilde{\boldsymbol A'}} (式(83)), {{\boldsymbol{C}}_\omega } (式(8)), {{\boldsymbol{O}}_{{\text{SPLE - prop}}}}(式(87))
计算 {{\boldsymbol{O}}_{{\text{SPLE - prop}}}}的EVD
根据式(88)—式(90)得到 {{\boldsymbol{\widehat s}}_{{\text{SPLE - prop}}}}
构造 {\boldsymbol{\widehat \varXi '}} (式(98)), {\boldsymbol{\widehat W'}} (式(86)), {{\boldsymbol{O}}_{{\text{BCSWPLE - prop}}}} (式(97))
计算 {{\boldsymbol{O}}_{{\text{BCSWPLE - prop}}}} 的EVD
根据式(88)—式(90)得到 {{\boldsymbol{\widehat s}}_{{\text{BCSWPLE - prop}}}}
构造 {\boldsymbol{G'}}(式(95)), {\boldsymbol{\widehat W'}} (式(96)), {{\boldsymbol{O}}_{{\text{SWIVE - prop}}}} (式(94))
计算 {{\boldsymbol{O}}_{{\text{SWIVE - prop}}}} 的EVD
根据式(88)—式(90)得到 {{\boldsymbol{\widehat s}}_{{\text{SWIVE - prop}}}}4. 仿真分析
本节将提出的SWPLE-prop算法和SWIVE-prop算法与SWPLE算法以及SWIVE算法的定位性能进行了比较。
4.1 试验参数与评价指标
为了较为全面地对算法定位性能进行评价, 仿真在
{\left[ {20,\;40} \right]\text{°} }{\text{N}} 、{\left[ {100,\;120} \right]\text{°} }{\text{E}} 的经纬度范围按均匀分布随机设置台站的空间分布, 每种台站分布的台站数量N = 30 , 仿真总台站分布数M = 50 , 声源位置在{\left[ {10,\;50} \right]\text{°} }{\text{N}} 、{\left[ {90,\;130} \right]\text{°} }{\text{E}} 的经纬度范围同样按均匀分布随机设置, 如图2所示, 其中蓝色五星表示随机声源位置的一次实现, 绿色倒三角表示随机台站位置的一次实现, 浅蓝色和浅绿色区域分别为随机声源和随机台站的位置范围。每种台站分布的蒙特卡罗试验次数为L = 2000 。SWIVE算法和SWIVE-prop算法的SAM阈值系数\rho = - 11.92{N^{ - 1.201}} + 3.445 [16]。仿真分析中假设各台站处的传播噪声来自区域内单一方向的大气平流层风, 传播噪声均值{\mu _n} 和方差\sigma _{\varepsilon ,n}^2 分别满足{\mu _n} = \left| {{\mu _0}} \right|b\left( {{d_n}} \right)\sin \left( {{\theta _{\text{w}}} - {\theta _n}} \right), (101a) \sigma _{\varepsilon ,n}^2 = \sigma _{\varepsilon ,0}^2b\left( {{d_n}} \right), (101b) 其中
b\left( {{d_n}} \right) = \left\{ {\begin{array}{ll} {0.005{d_n},} & {0 \leqslant {d_n} \leqslant 200\,{\text{km,}}} \\ {1,} & {{d_n} > 200\,{\text{km,}}} \end{array}} \right. \setcounter{equation}{102} (102) {\mu _0} 和\sigma _{\varepsilon ,0}^2 分别表示传播噪声的参考均值和方差,{d_n} 表示声源到台站的距离,\theta_{\mathrm{w}} 为大气平流层风的去风向, 仿真中设为\theta\mathrm{_w}=-45\text{°} 。算法定位性能的评价指标主要包括偏差距离
{e_{{\text{bias}}}} 和均方根误差 (RMSE){e_{{\text{RMS}}}} 两部分:{e_{{\text{bias}}}} = \frac{1}{M}\sum\limits_{k = 1}^M {H\left( {{\text{geomean}}\left( {{{{\boldsymbol{\widehat s}}}_{k,1}},\;{{{\boldsymbol{\widehat s}}}_{k,2}},\; \cdots ,\;{{{\boldsymbol{\widehat s}}}_{k,L}}} \right),\;{\boldsymbol{s}},{R_{\text{e}}}} \right)} , (103) {e_{{\text{RMS}}}} = \sqrt {\frac{1}{{ML}}\sum\limits_{k = 1}^M {\sum\limits_{l = 1}^L {{{\left( {H\left( {{{{\boldsymbol{\widehat s}}}_{k,l}},\;{\boldsymbol{s}},{R_{\text{e}}}} \right)} \right)}^2}} } } , (104) 其中,
{\text{geomean}}\left( \cdot \right) 表示地理均值,{{\boldsymbol{\widehat s}}_{k,l}} 表示台站分布k第l 次蒙特卡罗试验的声源位置估计。此外计算了给定条件下SWIVE-prop算法的理论RMSE (SWIVE-prop-th) [16]和声源位置估计RMSE的CRLB[16]作为算法性能参照。4.2 测量噪声水平
本仿真假设传播噪声参数不变, 令
\left| {{\mu _0}} \right| = {3\text{°} } 和{\sigma _{\varepsilon ,0}} = {3\text{°} } , 其余设置与4.1节相同。首先分析了测量噪声服从独立同分布的高斯噪声, 即{\boldsymbol{e}} \sim N\left( {{\boldsymbol{0}},\sigma _e^2{\boldsymbol{I}}} \right) 时各算法的定位性能。图3展示了当测量噪声标准差{\sigma _e} \in \left\{ {{{0.5}\text{°} },{1\text{°} }, \cdots ,{5\text{°} }} \right\} 时各算法的偏差和RMSE性能, 其中“-th”后缀表示算法理论结果。从图中可以看出SWPLE算法和SWIVE算法由于没有考虑传播噪声, 在整个测量噪声仿真范围内定位性能较差, 其中SWPLE算法由于本身的偏差特性导致其偏差和RMSE性能是参与比较算法中最差的。存在传播噪声时SWIVE算法也不再具有渐近无偏和渐近有效的特性, 存在较大的定位偏差和RMSE。由于SWIVE算法SAM阈值{\chi _n} 的选取忽略了{\sigma _\varepsilon } 的影响, 因此在{\sigma _\varepsilon } 一定时{\sigma _e} 越大{\chi _n} 越接近理想取值, 导致SWIVE算法的偏差性能随{\sigma _e} 的增大而有所提升。SWPLE-prop算法考虑了传播噪声的影响, 在测量噪声水平较小时其偏差和RMSE表现均要好于传统算法。但是受限于其偏差特性, 算法始终存在较大的偏差且偏差随{\sigma _e} 的增大而增大。因此当测量噪声水平较大时, SWPLE-prop算法的偏差和RMSE性能均劣于SWIVE算法但仍优于SWPLE算法。得益于渐近无偏和渐近有效性, SWIVE-prop算法取得了最优的定位性能, 在测量噪声仿真范围内近乎无偏, 并且RMSE性能接近CRLB。进一步将测量噪声设置为非同分布高斯噪声, 即
{e_n} \sim N\left( {0,\sigma _{e,n}^2} \right) , 其中\sigma _{e,n}^2 = \sigma _{e,0}^2\left( {{{\delta _n^2\left( {\boldsymbol{s}} \right)} \mathord{\left/ {\vphantom {{\delta _n^2\left( {\boldsymbol{s}} \right)} {\delta _0^2}}} \right. } {\delta _0^2}}} \right) 。图4展示了当测量噪声的参考标准差{\sigma _{e,0}} \in \left\{ {{{0.5}\text{°} },{1\text{°} }, \cdots ,{5\text{°} }} \right\} 时各算法的偏差和RMSE性能。图4中展示的仿真结果与图3十分接近, 进一步证明了SWIVE-prop算法在存在传播噪声的环境中可以取得最优的定位性能。从图3和图4中还可以看出SWIVE-prop算法的理论RMSE性能可以达到CRLB。4.3 传播噪声水平
本节假设测量噪声参数不变, 即服从4.2节设定的非同分布高斯噪声, 令
{\sigma _{e,0}} = {3\text{°} } 。首先分析了传播噪声均值变化对各算法定位性能的影响, 图5展示了当{\sigma _{\varepsilon ,0}} = {3\text{°} } ,\left| {{\mu _0}} \right| \in \left\{ {{0\text{°} },{{0.5}\text{°} }, \cdots ,{5\text{°} }} \right\} 时各算法的偏差和RMSE性能。从图中可以看出SWPLE算法和SWIVE算法的定位性能受\left| {{\mu _0}} \right| 的影响显著, 当\left| {{\mu _0}} \right| 增大时这两种算法的偏差和RMSE性能均显著下降。同时可以发现SWIVE算法的偏差在\left| {{\mu _0}} \right| = {0\text{°} } 时较小, 说明存在传播噪声时SWIVE算法的偏差主要来自噪声的均值。SWPLE-prop算法和SWIVE-prop算法由于考虑了传播噪声, 定位性能稳定。SWPLE-prop算法的大偏差和大RMSE源自其偏差特性, 而SWIVE-prop算法在\left| {{\mu _0}} \right| 的仿真范围内均表现出近乎无偏的偏差性能以及接近CRLB的RMSE性能。进一步设置
\left| {{\mu _0}} \right| = {3\text{°} } ,{\sigma _{\varepsilon ,0}} = {3\text{°} } , 通过改变平流层风的去风向\theta\mathrm{_w} 来改变各台站处传播噪声均值{\mu _n} 的大小和正负, 从而进一步考察各算法的稳定性。图6展示了\theta_{\mathrm{w}}\in\left\{-135\text{°},-90\text{°},\cdots,180\text{°}\right\} 时各算法的偏差和RMSE性能。SWPLE算法和SWIVE算法的定位性能随风向的改变并不稳定, 其中SWPLE算法的偏差和RMSE性能起伏较大。而SWPLE-prop算法和SWIVE-prop算法对于不同风向的定位性能比较稳定。最后设置
\left| {{\mu _0}} \right| = {3\text{°} } ,\theta_{\mathrm{w}}=-45\text{°} , 在{\sigma _{\varepsilon ,0}} \in \left\{ {{0.5}\text{°} },{1\text{°} }, \cdots , {5\text{°} } \right\} 时分析各算法的定位性能, 结果如图7所示。可以看出仿真结果与4.2节结果相似, SWIVE-prop算法的定位性能优于其他算法。SWIVE算法在{\sigma _{\varepsilon ,0}} 较大时的偏差增加同样源于忽略传播噪声导致的SAM阈值选取不合理。4.4 观测数
本仿真展示了各算法定位性能随观测数
N \in \left\{ {20,30, \cdots ,100} \right\} 的变化情况, 假设{\sigma _{e,0}} = {3\text{°} } ,\left| {{\mu _0}} \right| = {3\text{°} } ,{\sigma _{\varepsilon ,0}} = {3\text{°} } , 其余设置与4.1节相同, 结果如图8所示。从图中可以看出SWIVE-prop算法得益于其渐近无偏和渐近有效的特性, 在N的仿真范围内始终具有较小的偏差和RMSE, 同时随N的增大偏差逐渐趋于0, RMSE逐渐接近CRLB。其余算法均存在较大的定位偏差和RMSE, 其中定位偏差随N的增加没有明显减小, RMSE随N的增加虽然减小, 但减小幅度不明显, 距离CRLB较远。4.5 算法理论与仿真性能对比
本仿真将各算法的理论偏差距离和RMSE与仿真结果进行对比。由于各个算法定位性能理论值要求观测数N足够大, 所以设置
N = 150 。另外由于SWIVE算法SAM阈值选取会影响其定位性能, 妨碍与理论值的比较, 所以本节中SWIVE算法SAM阈值设置为最优取值。其余仿真设置与4.2节非同分布测量噪声仿真保持一致, 结果如图9所示。可以看出各算法偏差和RMSE性能理论与仿真结果对应良好, 证明了本文中算法理论分析的正确性。为进一步分析观测数N对算法定位性能理论准确性的影响, 对比了
N \in \left\{ {20,40, \cdots ,200} \right\} 时各算法理论定位性能和仿真结果, 如图10所示。可以看出随着N的增加, 各算法的理论定位性能逐渐接近仿真结果。在N较小时SWPLE和SWPLE-prop算法的理论定位性能与仿真结果有一定的差异, 而SWIVE和SWIVE-prop算法的理论定位性能和仿真结果在整个N的仿真范围内均表现出良好的契合度, 因此在实际应用中SWIVE-prop算法的理论误差协方差可以用作定位性能的评价指标。5. 实验验证
英国邦斯菲尔德当地的一个油库于2005年12月11日发生爆炸[20], 爆炸位置 (51.78°N, 0.43°W)。德国IGA、I26次声台站, 瑞典UPP次声台站均接收到了该事件产生的次声信号。该事件的特点是单个台站接收到了多个到达信号, 并且UPP台站和IGA台站处收到的信号受到大气横向风的影响, 传播方向发生了明显偏移。这导致台站处测得的后方位角偏离了真实方位, 因此该事件十分适合用于实验验证。事件的声源台站地理位置关系如图11所示, 其中红色五星表示真实事件位置, 蓝色倒三角表示次声台站位置, 黑色线条表示示向线。
使用PMCC算法[24-25]和3.1.1节描述的测量噪声参数估计方法得到的次声后方位角测量值和测量噪声标准差估计结果如表3所示, 可知IGA和UPP台站的后方位角测量受到了大气横向风影响而明显偏离真实方位。本研究分别在只使用ECM和使用基于ECM和ERA5大气再分析数据的混合大气模型 (Hybrid) 的条件下估计了传播噪声均值和标准差。两种大气剖面模型如图12所示, 混合大气模型相较于ECM在大气中间层及以下展现出更丰富的细节, 更有利于准确估计传播噪声。传播噪声参数估计的参数设置如表4所示, 传播噪声估计结果如表5所示, 使用混合大气模型估计的传播噪声相较于只使用ECM估计更加接近真实情况。
表 3 邦斯菲尔德爆炸事件次声信号检测和测量噪声参数估计结果次声
台站地理位置 后方位角
真值 (°)编号 后方位角
测量值
(°)标准差估
计 (°)IGA 53.26°N, 8.69°E −101.28 I1 −112.91 0.52 I2 −104.09 0.30 I3 −111.60 0.96 I4 −106.61 0.78 I5 −107.38 1.12 I26 48.85°N, 13.71°E −66.66 I1 −65.33 1.68 I2 −67.19 0.38 I3 −67.35 0.30 I4 −67.44 1.04 I5 −66.34 0.23 I6 −68.53 0.18 UPP 59.85°N, 17.61°E −120.77 I1 −132.14 4.83 I2 −124.72 8.84 I3 −130.69 8.24 I4 −124.55 4.87 I5 −131.13 5.90 表 4 传播噪声估计的参数设置参数 取值 平均声速c 300 m/s 纬度间隔\varDelta \varphi 1° 经度间隔\varDelta \lambda 1° 最小融合高度 {h_{\min }} 60 km 最大融合高度 {h_{\max }} 80 km 射线追踪方位角范围{\psi _{{\text{lim}}}} 10° 射线追踪方位角步长\varDelta \psi 1° 射线追踪俯仰角范围\left[ {{\phi _{\min }},{\phi _{\max }}} \right] \left[ {1\text{°},\;45\text{°}} \right] 射线追踪俯仰角步长\varDelta \phi 1° 到达信号样本数P 50 表 5 邦斯菲尔德爆炸事件传播噪声参数估计结果编号 期望估计
(ECM)
(°)标准差估计
(ECM)
(°)期望估计
(Hybrid)
(°)标准差
估计 (Hybrid)
(°)IGAI1 −3.19 0.37 −5.71 1.22 IGAI2 −1.71 0.20 −3.50 0.36 IGAI3 −3.25 0.40 −6.29 1.11 IGAI4 −2.18 0.19 −3.91 0.32 IGAI5 −2.47 0.20 −4.39 0.54 I26I1 0.71 0.44 −0.11 1.56 I26I2 0.52 0.36 1.56 0.18 I26I3 0.59 0.36 1.56 0.18 I26I4 1.20 0.18 −1.56 0.14 I26I5 1.61 0.15 −1.56 0.19 I26I6 1.70 0.19 −1.77 0.16 UPPI1 −4.03 1.70 −6.77 2.09 UPPI2 −3.46 1.42 −5.50 1.67 UPPI3 −2.89 0.54 −5.59 1.63 UPPI4 −4.96 1.21 −8.18 2.68 UPPI5 −2.87 0.55 −5.37 1.56 分别使用SWIVE算法和SWIVE-prop算法对事件进行定位, 结果如图13和表6所示。图13中圆点表示事件位置估计, 椭圆表示99%误差椭圆, 其中蓝色实线、短横线和点线分别表示未考虑传播噪声、使用ECM和混合大气模型计算的SWIVE算法误差椭圆。从实验结果中可以看出由于SWIVE算法未考虑传播噪声, 事件定位结果与真实位置相差较大, SWIVE-prop算法的定位结果更加靠近事件真实位置。由于传统的SWIVE算法性能分析中并未将次声传播的影响考虑其中, 因此图13中蓝色实线表示的未考虑传播噪声计算的误差椭圆低估了SWIVE算法的不确定度, 并不能准确反映算法定位性能。而使用本文提出的考虑传播噪声的SWIVE分析方法并使用混合大气模型计算的误差椭圆更符合算法的真实定位性能, 同时对比只使用ECM (蓝色短横线) 和使用混合大气模型 (蓝色点线) 计算的误差椭圆可以看出大气剖面模型的准确性对于算法不确定度的计算至关重要, 同时可以验证混合大气模型相较于只使用ECM对真实大气的描述更准确。式(105)定义了变量
{x_i} 与变量{x_j} 相比的百分比变化量p , 根据该式可得使用混合大气模型计算的SWIVE-prop算法定位误差比SWIVE算法减少79.6%, SWIVE-prop算法误差椭圆面积比使用混合大气模型计算的SWIVE算法误差椭圆面积减少89.4%。表 6 算法定位性能比较算法 定位误差 (km) 99%误差椭圆面积 (km2) SWIVE 92.35 365.28 6383.99 (ECM) 8836.15 (Hybrid) SWIVE-prop 39.88 772.36 (ECM) 18.88 933.71 (Hybrid) p = \frac{{{x_i} - {x_j}}}{{{x_j}}} \times 100{\text{% }}{\text{.}} (105) 6. 结论
本研究首次在球面伪线性估计器中考虑了次声传播过程的影响, 将其建模为传播噪声加入噪声模型, 并给出了SWPLE算法和SWIVE算法偏差和误差协方差的解析表达式。通过理论分析可以发现, 球面伪线性估计器的定位性能均受传播噪声的影响而明显下降, 其中SWIVE算法失去了渐近无偏性和渐近有效性。为了使算法能更好地解决广域次声源定位问题, 使用基于ECM和大气再分析数据的混合大气剖面模型和射线追踪方法估计传播噪声的均值和方差, 提出了SWPLE-prop算法和SWIVE-prop算法。仿真分析表明, 各算法的理论性能分析准确, 并且SWIVE-prop算法在不同的台站分布、测量噪声标准差、传播噪声均值和标准差以及观测数下均可以达到接近无偏的偏差性能和接近CRLB的RMSE性能。使用邦斯菲尔德爆炸事件验证了SWIVE-prop算法的性能优势, 其中使用混合大气模型的SWIVE-prop算法相较SWIVE算法定位误差减少79.6%, 99%误差椭圆面积减少89.4%。
-
表 1 SWPLE-prop算法流程
获取 {{\boldsymbol{r}}_n}, n \in \left\{ {1,2, \cdots ,N} \right\} 根据3.1.1节测量噪声参数估计方法得到 \widetilde \theta , \widetilde \sigma _{e,n}^2
使用SWIVE算法估计初始位置 {{\boldsymbol{\widehat s}}_{{\text{SWIVE}}}}
根据3.1.2节传播噪声参数估计方法得到{\widetilde \mu _n}, \widetilde \sigma _{\varepsilon ,n}^2
构造 {\widetilde{\boldsymbol A'}} (式(83)), {{\boldsymbol{C}}_\omega } (式(8)), {{\boldsymbol{O}}_{{\text{SPLE - prop}}}}(式(87))
计算 {{\boldsymbol{O}}_{{\text{SPLE - prop}}}}的EVD
根据式(88)—式(90)得到 {{\boldsymbol{\widehat s}}_{{\text{SPLE - prop}}}}
构造 {\boldsymbol{\widehat W'}} (式(86)), {{\boldsymbol{O}}_{{\text{SWPLE - prop}}}}(式(85))
计算 {{\boldsymbol{O}}_{{\text{SWPLE - prop}}}}的EVD
根据式(88)—式(90)得到 {{\boldsymbol{\widehat s}}_{{\text{SWPLE - prop}}}}
下载: 导出CSV
表 2 SWIVE-prop算法流程
获取 {{\boldsymbol{r}}_n}, n \in \left\{ {1,2, \cdots ,N} \right\} 根据3.1.1节测量噪声参数估计方法得到 \widetilde \theta , \widetilde \sigma _{e,n}^2
使用SWIVE算法估计初始位置 {{\boldsymbol{\widehat s}}_{{\text{SWIVE}}}}
根据3.1.2节传播噪声参数估计方法得到{\widetilde \mu _n}, \widetilde \sigma _{\varepsilon ,n}^2
构造 {\widetilde{\boldsymbol A'}} (式(83)), {{\boldsymbol{C}}_\omega } (式(8)), {{\boldsymbol{O}}_{{\text{SPLE - prop}}}}(式(87))
计算 {{\boldsymbol{O}}_{{\text{SPLE - prop}}}}的EVD
根据式(88)—式(90)得到 {{\boldsymbol{\widehat s}}_{{\text{SPLE - prop}}}}
构造 {\boldsymbol{\widehat \varXi '}} (式(98)), {\boldsymbol{\widehat W'}} (式(86)), {{\boldsymbol{O}}_{{\text{BCSWPLE - prop}}}} (式(97))
计算 {{\boldsymbol{O}}_{{\text{BCSWPLE - prop}}}} 的EVD
根据式(88)—式(90)得到 {{\boldsymbol{\widehat s}}_{{\text{BCSWPLE - prop}}}}
构造 {\boldsymbol{G'}}(式(95)), {\boldsymbol{\widehat W'}} (式(96)), {{\boldsymbol{O}}_{{\text{SWIVE - prop}}}} (式(94))
计算 {{\boldsymbol{O}}_{{\text{SWIVE - prop}}}} 的EVD
根据式(88)—式(90)得到 {{\boldsymbol{\widehat s}}_{{\text{SWIVE - prop}}}}
下载: 导出CSV
表 3 邦斯菲尔德爆炸事件次声信号检测和测量噪声参数估计结果
次声
台站地理位置 后方位角
真值 (°)编号 后方位角
测量值
(°)标准差估
计 (°)IGA 53.26°N, 8.69°E −101.28 I1 −112.91 0.52 I2 −104.09 0.30 I3 −111.60 0.96 I4 −106.61 0.78 I5 −107.38 1.12 I26 48.85°N, 13.71°E −66.66 I1 −65.33 1.68 I2 −67.19 0.38 I3 −67.35 0.30 I4 −67.44 1.04 I5 −66.34 0.23 I6 −68.53 0.18 UPP 59.85°N, 17.61°E −120.77 I1 −132.14 4.83 I2 −124.72 8.84 I3 −130.69 8.24 I4 −124.55 4.87 I5 −131.13 5.90
下载: 导出CSV
表 4 传播噪声估计的参数设置
参数 取值 平均声速c 300 m/s 纬度间隔\varDelta \varphi 1° 经度间隔\varDelta \lambda 1° 最小融合高度 {h_{\min }} 60 km 最大融合高度 {h_{\max }} 80 km 射线追踪方位角范围{\psi _{{\text{lim}}}} 10° 射线追踪方位角步长\varDelta \psi 1° 射线追踪俯仰角范围\left[ {{\phi _{\min }},{\phi _{\max }}} \right] \left[ {1\text{°},\;45\text{°}} \right] 射线追踪俯仰角步长\varDelta \phi 1° 到达信号样本数P 50
下载: 导出CSV
表 5 邦斯菲尔德爆炸事件传播噪声参数估计结果
编号 期望估计
(ECM)
(°)标准差估计
(ECM)
(°)期望估计
(Hybrid)
(°)标准差
估计 (Hybrid)
(°)IGAI1 −3.19 0.37 −5.71 1.22 IGAI2 −1.71 0.20 −3.50 0.36 IGAI3 −3.25 0.40 −6.29 1.11 IGAI4 −2.18 0.19 −3.91 0.32 IGAI5 −2.47 0.20 −4.39 0.54 I26I1 0.71 0.44 −0.11 1.56 I26I2 0.52 0.36 1.56 0.18 I26I3 0.59 0.36 1.56 0.18 I26I4 1.20 0.18 −1.56 0.14 I26I5 1.61 0.15 −1.56 0.19 I26I6 1.70 0.19 −1.77 0.16 UPPI1 −4.03 1.70 −6.77 2.09 UPPI2 −3.46 1.42 −5.50 1.67 UPPI3 −2.89 0.54 −5.59 1.63 UPPI4 −4.96 1.21 −8.18 2.68 UPPI5 −2.87 0.55 −5.37 1.56
下载: 导出CSV
表 6 算法定位性能比较
算法 定位误差 (km) 99%误差椭圆面积 (km2) SWIVE 92.35 365.28 6383.99 (ECM) 8836.15 (Hybrid) SWIVE-prop 39.88 772.36 (ECM) 18.88 933.71 (Hybrid)
下载: 导出CSV
-
[1] Le Pichon A, Blanc E, Hauchecorne A. Infrasound monitoring for atmospheric studies: Challenges in middle atmosphere dynamics and societal benefit. Cham: Springer International Publishing, 2019
[2] Lee W H K, Stewart S W. Principles and applications of microearthquake networks. New York: Academic Press, 1981
[3] Jordan T H, Sverdrup K A. Teleseismic location techniques and their application to earthquake clusters in the South-Central Pacific. Bull. Seism. Soc. Am., 1981; 71(4): 1105−1130 DOI: 10.1785/BSSA0710041105
[4] Bratt S R, Bache T C. Locating events with a sparse network of regional arrays. Bull. Seism. Soc. Am., 1988; 78(2): 780−798 DOI: 10.1785/bssa0780020780
[5] Blom P S, Marcillo O, Arrowsmith S J. Improved Bayesian infrasonic source localization for regional infrasound. Geophys. J. Int., 2015; 203(3): 1682−1693 DOI: 10.1093/gji/ggv387
[6] Arrowsmith S, Park J, Che I Y, et al. Event location with sparse data: When probabilistic global search is important. Seismol. Res. Lett., 2021; 92(2A): 976−985 DOI: 10.1785/0220200292
[7] Matoza R S, Green D N, Le Pichon A, et al. Automated detection and cataloging of global explosive volcanism using the International Monitoring System infrasound network. J. Geophys. Res.: Solid Earth, 2017; 122(4): 2946−2971 DOI: 10.1002/2016JB013356
[8] De Negri R, Matoza R S. Rapid location of remote volcanic infrasound using 3D ray tracing and empirical climatologies: Application to the 2011 Cordón Caulle and 2015 Calbuco eruptions, Chile. J. Geophys. Res.: Solid Earth, 2023; 128(3): 1−18 DOI: 10.1029/2022JB025735
[9] Almendros J, Chouet B. Performance of the radial semblance method for the location of very long period volcanic signals. Bull. Seism. Soc. Am., 2003; 93(5): 1890−1903 DOI: 10.1785/0120020143
[10] Mckee K, Fee D, Rowell C, et al. Network-based evaluation of the infrasonic source location at Sakurajima volcano, Japan. Seismol. Res. Lett., 2014; 85(6): 1200−1211 DOI: 10.1785/0220140119
[11] Fee D, Toney L, Kim K, et al. Local explosion detection and infrasound localization by reverse time migration using 3-D finite-difference wave propagation. Front. Earth Sci., 2021; 9: 1−14 DOI: 10.3389/feart.2021.620813
[12] 吕君, 郭泉, 冯浩楠, 等. 北京地震前的异常次声波. 地球物理学报, 2012; 55(10): 3379−3385 DOI: 10.6038/j.issn.0001-5733.2012.10.020 [13] 郭泉, 杨亦春, 吕君, 等. 基于广域次声传感器网络的地震本地次声波监测. 地球科学(中国地质大学学报), 2014; 39(12): 1807−1817 DOI: 10.3799/dqkx.2014.164 [14] 杨亦春, 郭泉, 吕君, 等. 大地震前出现的异常次声波观测研究. 物理学报, 2014; 63(13): 224−237 DOI: 10.7498/aps.63.134302 [15] 赵久彬, 刘元雪, 杨骏堂, 等. 滑坡次声信号简正波模型的匹配场超前定位算法. 岩土力学, 2020; 41(12): 4116−4126 DOI: 10.16285/j.rsm.2020.0345 [16] Zhang T, Teng P, Lyu J, et al. Quasi-closed-form algorithms for spherical angle-of-arrival source localization. IEEE Trans. Signal Process., 2024; 72: 432−448 DOI: 10.1109/TSP.2023.3340879
[17] Diamond M. Cross wind effect on sound propagation. J. Appl. Meteorol. Clim., 1964; 3(2): 208−210 DOI: 10.1175/1520-0450(1964)003<0208:CWEOSP>2.0.CO;2
[18] Smets P S M, Assink J D, Le Pichon A, et al. ECMWF SSW forecast evaluation using infrasound. J. Geophys. Res.: Atmos., 2016; 121(9): 4637−4650 DOI: 10.1002/2015JD024251
[19] Blom P, Waxler R. Modeling and observations of an elevated, moving infrasonic source: Eigenray methods. J. Acoust. Soc. Am., 2017; 141: 2681−2692 DOI: 10.1121/1.4980096
[20] Ceranna L, Le Pichon A, Green D N, et al. The Buncefield explosion: A benchmark for infrasound analysis across Central Europe. Geophys. J. Int., 2009; 177(2): 491−508 DOI: 10.1111/j.1365-246X.2008.03998.x
[21] Pilger C, Ceranna L, Ross J O, et al. The European infrasound bulletin. Pure Appl. Geophys., 2018; 175(10): 3619−3638 DOI: 10.1007/s00024-018-1900-3
[22] Vincenty T. Direct and inverse solutions of geodesics on the ellipsoid with application of nested equations. Surv. Rev., 1975; 23(176): 88−93 DOI: 10.1179/sre.1975.23.176.88
[23] Doǧançay K. 3D pseudolinear target motion analysis from angle measurements. IEEE Trans. Signal Process., 2015; 63(6): 1570−1580 DOI: 10.1109/TSP.2015.2399869
[24] Cansi Y. An automatic seismic event processing for detection and location: The P.M.C.C. method. Geophys. Res. Lett., 1995; 22(9): 1021−1024 DOI: 10.1029/95GL00468
[25] Havelock D, Kuwano S, Vorländer M. Handbook of signal processing in acoustics. New York: Springer, 2008: 1425−1435
[26] Hersbach H, Bell B, Berrisford P, et al. The ERA5 global reanalysis. Q. J. Roy. Meteor. Soc., 2020; 146(730): 1999−2049 DOI: 10.1002/qj.3803
[27] Korn G A, Korn T M. Mathematical handbook for scientists and engineers: Definitions, theorems, and formulas for reference and review. Dover ed. Mineola, N. Y.: Dover Publications, 2000
[28] Torrieri D J. Statistical theory of passive location systems. IEEE Trans. Aerosp. Electron. Syst., 1984; 20(2): 183−198 DOI: 10.1109/TAES.1984.310439
计量
- 文章访问数: 46
- HTML全文浏览量: 5
- PDF下载量: 16
