结合盲信号分离算法的局部放电TDOA/DOA混合定位方法

上一篇

下一篇

张泽宇, 汤晓君, 李晓杉, 刘崇智, 周佳泓. 结合盲信号分离算法的局部放电TDOA/DOA混合定位方法[J]. 物理学报, 2025, 74(13): 135207-1. doi: 10.7498/aps.74.20250317
引用本文: 张泽宇, 汤晓君, 李晓杉, 刘崇智, 周佳泓. 结合盲信号分离算法的局部放电TDOA/DOA混合定位方法[J]. 物理学报, 2025, 74(13): 135207-1. doi: 10.7498/aps.74.20250317
Zeyu ZHANG, Xiaojun TANG, Xiaoshan LI, Chongzhi LIU, Jiahong ZHOU. TDOA/DOA hybrid location method of partial discharge combined with blind signal separation algorithm[J]. Acta Physica Sinica, 2025, 74(13): 135207-1. doi: 10.7498/aps.74.20250317
Citation: Zeyu ZHANG, Xiaojun TANG, Xiaoshan LI, Chongzhi LIU, Jiahong ZHOU. TDOA/DOA hybrid location method of partial discharge combined with blind signal separation algorithm[J]. Acta Physica Sinica, 2025, 74(13): 135207-1. doi: 10.7498/aps.74.20250317

结合盲信号分离算法的局部放电TDOA/DOA混合定位方法

    作者简介: 张泽宇.E-mail: 0804160134@163.com .
    通讯作者: E-mail: xiaojun_tang@mail.xjtu.edu.cn.; 
  • 中图分类号: 52.80.Pi, 43.28.Tc, 51.50.+v, 47.50.Cd

TDOA/DOA hybrid location method of partial discharge combined with blind signal separation algorithm

    Corresponding author: E-mail: xiaojun_tang@mail.xjtu.edu.cn.; 
  • MSC: 52.80.Pi, 43.28.Tc, 51.50.+v, 47.50.Cd

  • 摘要: 针对电力设备局部放电(PD)超声波检测中存在的时间-空域特征解耦、硬件成本高及计算量大的技术瓶颈, 提出基于核主成分分析(KPCA)伪白化的改进型非圆FastICA (mnc-FastICA)提取TDOA/DOA参数的混合定位方法. 该方法通过时间-空域特征联合提取与智能优化机制, 实现了小规模传感器阵列下的高精度定位. 本文首先构建KPCA伪白化预处理框架, 利用多项式核函数映射信号非线性升维再降维, 在保留TDOA与DOA特征关联性的同时抑制环境噪声; 其次通过mnc-FastICA算法盲分离超声信号后, 联合广义互相关法GCC与阵列流型解析技术同步提取TDOA/DOA参数; 最后建立融合双参数的最大似然估计模型, 并引入非洲秃鹫优化算法实现全局最优解快速收敛. 实验表明, 在仅配置2个正交阵列(共8个传感器)的小规模硬件架构下, 本方法TDOA估计误差降至2.34%, DOA估计精度优于2°, 定位误差达1.54 cm. 该方法有效解决了PD检测中时间-空域特征联合、硬件成本与定位精度的矛盾, 为电力设备状态监测提供新方案.
  • 加载中
  • 图 1  信号混合与盲信号分离算法的过程

    Figure 1.  The process of signal mixed and blind signal separation.

    图 2  PD信号源与阵列模型图

    Figure 2.  Model diagram of PD source and sensor arrays.

    图 3  TDOA/DOA混合定位算法步骤

    Figure 3.  TDOA/DOA hybrid localization algorithm steps.

    图 4  盲信号分离算法提取信号波形对比

    Figure 4.  Waveform extracted by KPCA-mnc-FastICA compared with original signal.

    图 5  不同数量传感器阵元的DOA估计精度

    Figure 5.  DOA estimation accuracy versus number of sensor array.

    图 6  NCC与分离信号SNR随输入SNR变化情况

    Figure 6.  Variation of NCC and SNR of separated signal with different input SNR.

    图 7  不同SNR的TDOA提取方法性能对比

    Figure 7.  Performance of different TDOA estimation methods with different input SNR.

    图 8  TDOA-chan定位结果

    Figure 8.  Positioning results of TDOA-chan.

    图 9  不同信噪比的DOA估计精度 (a) 方位角; (b) 俯仰角

    Figure 9.  DOA estimation accuracy versus SNR: (a) Azimuth; (b) pitch.

    图 10  不同信噪比下的各种方法定位估计精度

    Figure 10.  RMSE of positioning methods versus SNR.

    图 11  不同阵元数目下的各种方法定位估计精度

    Figure 11.  RMSE of positioning methods versus number of sensor.

    图 12  不同阵列数目下的各种方法定位估计精度

    Figure 12.  RMSE of positioning methods versus number of sensor array.

    图 13  2个正交阵列的定位估计

    Figure 13.  Localization estimation of two orthogonal arrays.

    图 14  试验平台 (a) 示意图; (b) 装置图

    Figure 14.  Schematic diagram of test platform: (a) Schematic diagram; (b) device diagram.

    图 15  3个传感器阵列接收信号的盲信号分离信号 (a) 阵列1; (b) 阵列2; (c) 阵列3

    Figure 15.  The extracted signal obtained by blind signal separation of the received signal from the three sensor arrays: (a) Array 1; (b) array 2; (c) array 3.

    表 1  PD源和3路传感器阵列坐标

    Table 1.  Coordinates of Needle-plate model and three sensor array.

    设备类型坐标/cm
    PD源(20, 10, 10)
    传感器阵列1(0, 30, 15)
    传感器阵列2(40, 0, 5)
    传感器阵列3(60, 20, 20)
    下载: 导出CSV

    表 2  TDOA估计结果

    Table 2.  Results of TDOA estimation of the test.

    时间差 理论时间差/ms 试验时间差/ms 误差/μs
    t21 0.1568 0.1530 –3.8
    t31 0.2275 0.2221 –5.4
    下载: 导出CSV

    表 3  DOA估计结果

    Table 3.  Results of DOA estimation of the test.


    理论方
    位角/(º)
    试验方
    位角/(º)
    方位角
    误差/(º)
    理论俯
    仰角/(º)
    试验俯
    仰角/(º)
    俯仰角
    误差/(º)
    1 –45 –44 1 –10 –10 0
    2 153 155 2 12 11 1
    3 –165 –166 1 –13 –14 1
    下载: 导出CSV

    表 4  局放定位结果

    Table 4.  Localization results of PD source.

    阵列组合阵列序号估计坐标/cm均方根误差/cm
    11和2(20.95, 11.21, 10.09)1.54
    21和3(20.55, 11.79, 9.89)1.88
    32和3(21.25, 12.11, 9.95)2.45
    41, 2和3(21.13, 10.12, 10.14)1.14
    下载: 导出CSV
  • [1] Md-Rashid H, Refaat S, Abu-Rub H 2021 IEEE Ace 9 64587 doi: 10.1109/ACCESS.2021.3075288
    [2] 牛勃, 马飞越, 周秀 2019 高压电器 55 108 doi: 10.13296/j.1001-1609.hva.2019.08.015 Niu B, Ma F Y, Zhou X 2019 High Volt. Appr. 55 108 doi: 10.13296/j.1001-1609.hva.2019.08.015
    [3] Yan B W, Chang D G, Fan Y H, Zhang G J, Zhan J Y, Shuo X J 2017 IEEE Trans. Dielect. El. In. 24 3647 doi: 10.1109/TDEI.2017.006857
    [4] 唐炬, 陈娇, 张晓星, 许中荣 2009 中国电机工程学报 29 125 doi: 10.3321/j.issn:0258-8013.2009.19.019 Tang J, Chen J, Zhang X X, Xu Z R 2009 Procee. CSEE 29 125 doi: 10.3321/j.issn:0258-8013.2009.19.019
    [5] 李保伟, 张兴敢 2020 南京大学学报 56 917 doi: 10.13232/j.cnki.jnju.2020.06.015 Li B W, Zhang X G 2020 J. Nanjing Univ. 56 917 doi: 10.13232/j.cnki.jnju.2020.06.015
    [6] 李易, 任明, 王彬, 席英杰, 刘阳, 董明 2022 电网技术 47 4351 doi: 10.13335/j.1000-3673.pst.2022.1467 Li Y, Ren M, Wang B, Xi Y J, Liu Y, Dong M 2022 Power Syst. Tech. 47 4351 doi: 10.13335/j.1000-3673.pst.2022.1467
    [7] Ghosh G, Chatterjee B, Dalai S 2017 IEEE Trans. Dielect. El. In. 24 237 doi: 10.1109/TDEI.2016.006080
    [8] 关羽, 董明, 席英杰, 李易, 张崇兴 2023 电网技术 48 1721 doi: 10.13335/j.1000-3673.pst.2022.2475 Guan Y, Dong M, Xi Y J, Li Y, Zhang C X 2023 Power Syst. Tech. 48 1721 doi: 10.13335/j.1000-3673.pst.2022.2475
    [9] 王玉龙, 张晓虹, 李丽丽, 高俊国, 郭宁, 程成 2021 物理学报 70 095209 doi: 10.7498/aps.70.20201727 Wang Y L, Zhang X H, Li L L, Gao J G, Guo N, Cheng C 2021 Acta Phys. Sin. 70 095209 doi: 10.7498/aps.70.20201727
    [10] Chang Y T, Wu C L, Cheng H 2014 2014 International Symposium on Computer, Consumer and Control Taichung, China, June 10–12, 2014 p1 doi: 10.1109/IS3C.2014.281
    [11] Thu L N N, Tuan D V, Yoana S 2019 Sensers 19 2121 doi: 10.3390/s19092121
    [12] Jia T Y, Wang H Y, Shen X H, Gao J J, Liu X 2017 An accurate TDOA-AOA Localization Method Using Structured Total Least Squares Aberdeen, UK, June 19–22, 2017 p1 doi: 10.1109/OCEANSE.2017.8084582
    [13] Muhammad U L, Hafiz S M, Amna R, Zakria Q, Abbas Z K, Mahmud M A P 2021 Engrgies 14 3910 doi: 10.3390/en14133910
    [14] 王国利, 杜非, 潘伟, 王婷婷, 罗勇芬 2019 高电压技术 45 2509 doi: 10.13336/j.1003-6520.hve.20190731020 Wang G L, Du F, Pan W, Wang T T, Luo Y F 2019 High Vol. Eng. 45 2509 doi: 10.13336/j.1003-6520.hve.20190731020
    [15] Xu S 2020 IEEE Commun. Lett. 24 1966 doi: 10.1109/LCOMM.2020.2996259
    [16] 谢庆, 张建涛, 陈艺丹, 刘怡, 张莹 2019 高压电器 55 1 doi: 10.13296/j.1001-1609.hva.2019.12.001 Xie Q, Zhang J T, Chen Y D, Liu Y, Zhang Y 2019 High Volt Appr. 55 1 doi: 10.13296/j.1001-1609.hva.2019.12.001
    [17] Visalakshi A, Deepika R S, Chinthaginjala V R, Bagadi K, Mohammad A, Ayman A A 2022 IEEE Access 10 132875 doi: 10.1109/ACCESS.2022.3230661
    [18] Gulia S, Prasad P, Goyal S K, Raakesh K 2020 Atmop. Poll. 99 1588 doi: 10.1016/j.apr.2020.06.016
    [19] Ning S, He Y G, Ali F, Wu Y T, Tong J 2021 IEEE Sensors J. 21 6741 doi: 10.1109/JSEN.2020.3037699
    [20] Fokin G 2019 2019 21st International Conference on Advanced Communication Technology (ICACT) PyeongChang, Korea (South), February 17–20, 2019 p1
    [21] Yue Y G, Cao L, Hu J, Cai S T, Hang B, Wu H 2019 IEEE Access 7 58541 doi: 10.1109/ACCESS.2019.2914924
    [22] Chen T, Wang M X, Huang X S 2018 2018 14th IEEE International Conference on Signal Processing (ICSP) Beijing, China, August 12–16, 2018 p1
    [23] Li X L, Yi X, Zhang Z K 2021 Int. J. Antenn. Propag. 1 5512395 doi: 10.1155/2021/5512395
    [24] Wu P, Su S J, Zuo Z, Guo X J, Sun B, Wen X D 2019 Sensors 19 2554 doi: 10.3390/s19112554
    [25] Jiang F, Zhang Z K 2021 IET Commun. 15 802 doi: 10.1049/cmu2.12122
    [26] Liu Z X, Wang R, Zhao Y J 2018 Sensors 18 3747 doi: 10.3390/s18113747
    [27] 阮宗利, 魏平, 钱国兵, 袁晓垒 2017 电子科技大学学报 46 505 doi: 10.3969/j.issn.1001-0548.2017.03.005 Ruan Z L, Wei P, Qian G B, Yuan X L 2017 J. Univ. Electr. Sci. Tech. 46 505 doi: 10.3969/j.issn.1001-0548.2017.03.005
    [28] 肖剑, 刘经纬, 胡欣, 齐小刚 2023 吉林大学学报(工学版) 54 3558 doi: 10.13229/j.cnki.jdxbgxb.2023007 Xiao J, Liu J W, Hu X, Qi X G 2023 J. Jilin Univ. Eng. Tech. 54 3558 doi: 10.13229/j.cnki.jdxbgxb.2023007
    [29] Bingham E, Hyvärinen A 2000 Int. J. Neural Syst. 10 1 doi: 10.1142/S0129065700000028
    [30] Novey M, Adali T 2008 IEEE Trans. Sig. Proc. 56 2148 doi: 10.1109/TSP.2007.911278
    [31] Zhou L J, Cai J Y, Hu J J, Guo L, Liao W 2021 Sensors 18 3747 doi: 10.1109/TPWRD.2020.3011455
    [32] Pei S T, Liu D W, Ye Z J, Yang J J, Liu Y P 2023 IET Sci. Meas. Technol. 17 11 doi: 10.1049/smt2.12124
  • 加载中
图( 15) 表( 4)
计量
  • 文章访问数:  40
  • HTML全文浏览数:  40
  • PDF下载数:  0
  • 施引文献:  0
出版历程
  • 收稿日期:  2025-03-11
  • 刊出日期:  2025-07-05

结合盲信号分离算法的局部放电TDOA/DOA混合定位方法

    通讯作者: E-mail: xiaojun_tang@mail.xjtu.edu.cn.; 
    作者简介: 张泽宇.E-mail: 0804160134@163.com
  • 1. 西安交通大学电气工程学院, 西安 710049
  • 2. 西安交通大学仪器科学与技术学院, 西安 710049
  • 3. 西安交通大学精密微纳制造技术全国重点实验室, 西安 710049

摘要: 针对电力设备局部放电(PD)超声波检测中存在的时间-空域特征解耦、硬件成本高及计算量大的技术瓶颈, 提出基于核主成分分析(KPCA)伪白化的改进型非圆FastICA (mnc-FastICA)提取TDOA/DOA参数的混合定位方法. 该方法通过时间-空域特征联合提取与智能优化机制, 实现了小规模传感器阵列下的高精度定位. 本文首先构建KPCA伪白化预处理框架, 利用多项式核函数映射信号非线性升维再降维, 在保留TDOA与DOA特征关联性的同时抑制环境噪声; 其次通过mnc-FastICA算法盲分离超声信号后, 联合广义互相关法GCC与阵列流型解析技术同步提取TDOA/DOA参数; 最后建立融合双参数的最大似然估计模型, 并引入非洲秃鹫优化算法实现全局最优解快速收敛. 实验表明, 在仅配置2个正交阵列(共8个传感器)的小规模硬件架构下, 本方法TDOA估计误差降至2.34%, DOA估计精度优于2°, 定位误差达1.54 cm. 该方法有效解决了PD检测中时间-空域特征联合、硬件成本与定位精度的矛盾, 为电力设备状态监测提供新方案.

English Abstract

    • 变压器作为电力系统关键设备, 其绝缘性能直接决定电力系统的安全运行. 局部放电(partial discharge, PD)是指发生在电力设备绝缘介质局部区域的非贯穿性放电现象. 该现象会导致绝缘材料局部劣化、碳化或电蚀, 进而引发绝缘性能的渐进式衰退, 最终可能诱发绝缘系统失效[1]. 因此, 实现PD信号的精确检测与放电源的精确定位, 对提升电力设备运维水平具有重要的工程意义.

      PD信号作为电力设备绝缘状态评估的有效表征手段, 其时间-空域特性为放电源定位提供了重要依据. 目前, 该领域已形成电气检测与非电气检测相结合的多元方法体系, 主要包括特高频检测法、脉冲电流法、超声波检测法及油中溶解气体分析法等[1]. 其中, 超声波检测法因兼具传感器小型化、抗干扰性强以及电磁兼容性优异等特点, 在工程实践中获得广泛应用. 当前超声波信号处理技术主要涵盖傅里叶变换分析、小波变换分析、希尔伯特-黄变换及盲源分离技术; 而超声波定位方法则形成四大技术体系: 到达时间法(time of arrival, TOA)[2]、时差定位法(time difference of arrival, TDOA)[310]、波达方向法(direction of arrival, DOA)[11,12]以及信号强度定位法(received signal strength, RSS)[1315].

      然而, 现有研究普遍将PD超声波信号的特征提取与定位解算视作独立环节, 存在以下技术局限性: 其一, 时域特征与空间方位信息的解耦处理导致现有定位方法仅能依赖单一物理量; 其二, 超声波信号固有频带宽(50—400 kHz)、高频分量显著等特性, 使得DOA定位需配置高密度阵列传感器[16], 而TDOA定位则需四通道同步采集系统, 对传感器布设空间提出严苛要求[1719]. 研究显示, 融合TDOA与DOA的混合定位方法通过多维特征联合优化, 在无线传感网络领域展现出较单参量方法更优的定位精度[2022], 这为电力设备局放定位提供了新的技术路径.

      混合定位算法的构建面临双重挑战: 一方面, 初值敏感且计算量较大, 现有算法如三维闭式解算法[23]虽通过伪线性化处理降低计算量, 但存在参数近似误差; 基于最小二乘(least squares, LS)的定位模型[24]虽能改善估计精度, 仍未能突破算法对初值的敏感性; 另一方面, 误差补偿机制较少, 两步加权最小二乘法(two-stage weighted least squares, TSWLS) [25]虽引入误差补偿机制, 其计算复杂度仍制约工程应用. 相较而言, 最大似然估计法(maximum likelihood estimation, MLE)通过构建传感器位置误差与测量噪声的联合概率模型, 结合偏差补偿机制, 可获得更稳健的估计结果[26].

      针对上述问题, 本研究提出一种TDOA/DOA混合定位方法, 创新性地将核主成分分析(kernel principal component analysis, KPCA)与改进非圆FastICA (modified noncircular FastICA, mnc-FastICA)[27]算法相结合. 首先利用KPCA对非高斯特征的超声信号进行伪白化预处理, 通过mnc-FastICA算法实现信号分离并获得分离矩阵. 继而采用广义互相关(generalized cross-correlation, GCC)[5]法从分离信号中提取TDOA参数, 同时通过分离矩阵的数学变换解析DOA参数. 在此基础上, 构建MLE模型融合多维特征, 并引入非洲秃鹫优化算法(African vulture optimization algorithm)[28]建立适应度函数, 通过群体智能优化机制降低计算复杂度. 仿真与实验表明, 所提出的KPCA-mnc-FastICA混合定位框架能有效协同TDOA与DOA信息, 降低硬件开销, 实现厘米级定位精度.

    • 盲信号分离算法作为信号处理领域的重要分支, 其核心优势在于无需源信号先验信息及混合过程参数, 仅基于源信号与观测噪声间的统计独立性即可实现混合信号的解耦重构[27]. 针对PD检测中源信号未知且混合过程存在复杂时变特性的实际工况, 该算法通过构建多维信号统计学模型, 可有效解析传感器阵列接收信号中的独立源分量, 其强鲁棒性与环境自适应性使其在PD超声波信号处理领域具有独特优势.

      图1展示了信号混合与盲信号分离算法的数学模型. 设$p$个PD源产生时域信号构成的源矩阵$ {\boldsymbol{S}}(t) = \left[ s_{1}(t),\; s_{2}(t),\; \cdots ,\; s_{p}(t) \right]^{\text{T}} $, 经$n \times p$维混合矩阵A作用后生成观测到的接收信号$ {\boldsymbol{X}}(t) = \left[x_{1}(t), x_{2}(t),\;\cdots ,\;x_{n}(t)\right]^{\text{T}} $. 为保证盲信号算法的分离效果, 需满足接收信号$ {\boldsymbol{X}} $的维度$n$大于源信号数量$p$[27]. 实际情况中, 混合矩阵A仅含有传感器几何配置等有限先验知识. 盲信号分离算法通过构建目标函数优化分离矩阵W, 使得分离信号$ \hat {\boldsymbol{S}}(t) = [{\hat s}_{1}(t), {\hat s}_{2}(t),\; \cdots ,\;{\hat s}_p(t) ]^{\text{T}} $能够最大化保留源信号的统计独立性.

      为实现TDOA与DOA信息的高效协同, 本研究构建了分阶段特征融合的混合定位框架. 第1部分在第2.2节构建了基于KPCA-mnc-FastICA的信号提取框架; 第2部分在第3节构建了TDOA和DOA特征提取框架, 并提出基于最大似然估计MLE的TDOA/DOA混合定位模型. 两阶段处理有效实现了时间-空域特征参数的联合优化.

    • 在窄带信号假设下, 超声波传播可简化为平面波模型. 如图2所示, PD源发射的超声波信号$ {\boldsymbol{s}}(t) $经空间传播到达两个正交分布的阵列O1O2. 以阵列O1为例, 设信号$ {\boldsymbol{s}}(t) $从PD源到阵列O1$ j $个传感器阵元的时间延迟为$ {\tau _j} $, 阵元接收信号的时域复包络形式表示为

      式中, $ {\boldsymbol{\alpha}} (t) $为信号$ {\boldsymbol{s}}(t) $的时变包络幅度, $ \omega $为信号中心频率, $ t $为时间变量, $ {\tau _j} $为信号$ {\boldsymbol{s}}(t) $从PD源到阵列O1j个阵元的传播时延.

      基于窄带信号假设$ {\boldsymbol{s}}(t + {\tau _j}) \approx {\boldsymbol{s}}(t) $, 从(1)式可推导出:

      则阵列O1中第j个阵元接收到的信号可表示为

      式中, $ {{\boldsymbol{x}}_j}(t) $为阵列O1中第j个阵元接收到的信号; $ {{\boldsymbol{n}}_j}(t) $为第j个阵元接收到的噪声.

      若将n个传感器排成一排列矢量, 按(3)式计算接收信号, 可得:

      式中, 当$ {\boldsymbol{s}}(t) $长度为$ N $时, 接收信号为$ n \times N $维数据矩阵, 命名为接收矩阵$ {\boldsymbol{X}}(t) $; 向量$ \left[\exp ( - {\text{i}}\omega {\tau _1}), \exp ( - {\text{i}}\omega {\tau _2}),\; \cdots ,\; \exp ( - {\text{i}}\omega {\tau _n}) \right]^{\text{T}} $的维度为$ n \times 1 $, 命名为导向矢量a; 接收的噪声为$ n \times N $维矩阵, 命名为噪声矩阵$ {\boldsymbol{N}}(t) $.

      由(4)式可明确导向矢量a的物理内涵, 该矢量本质为信号空间传播特性的参数化表征, 其元素对应信号入射角作用下阵列各阵元的相位延迟量. 如图1所示, 当存在$p$个PD源信号$ {{\boldsymbol{s}}_i}(t), (1 \leqslant i \leqslant p) $时, 阵列将对应$p$组时延参数, 则阵列O1的流型矩阵可构建为

      式中, 第$i$个列向量$ {{\boldsymbol{a}}_i} $由相应PD源信号$ {{\boldsymbol{s}}_i}(t) $的时延信息计算而成.

      基于图2所示双阵列拓扑结构, 对于$p$个PD源构成的源信号矩阵$ {\boldsymbol{S}}(t) $, 阵列O1O2的接收信号矩阵依照(4)式的方式进行建模:

      式中, $ m = 1, 2 $, 代表图2所示阵列的序号; $ {{\boldsymbol{X}}_m}(t) $, $ {{\boldsymbol{A}}_m} $$ {{\boldsymbol{N}}_m}(t) $分别为阵列Om接收到的信号、阵列流型矩阵和噪声.

      由(4)式和(5)式可知, 导向矢量$ {{\boldsymbol{a}}_i} $与信号的时延信息相关. 设方形阵列O1$ n \times n $个传感器阵元, 阵元间隔为$ d $($ d = \lambda /2 $), 以阵列左上角阵元为参考点, 信号$ s(t) $以方位角$ \varphi $与俯仰角$ \theta $入射到的天线阵列上, 方位角$ \varphi $表示信号源PD到参考点间的连线在x-y平面上投影与x轴夹角, 俯仰角度$ \theta $是信号源PD到参考点间的连线与x-y平面的夹角, 则信号传播到第$ j $个阵元与参考阵元的时延为

      式中, $ {x_j} $$ {y_j} $表示第$ j $个阵元距离x轴、y轴的距离, $ j = 1, 2, \cdots , n \times n $; c表示信号传播速度.

      $ \omega = 2{{\pi}}c/\lambda $与(7)式代入(4)式中的导向矢量$ {{\boldsymbol{a}}_i} $, 得

      其中$ \lambda $为信号波长, 该式具体推导过程请见补充材料(online). 结合(5)式可知, 导向矢量$ {{\boldsymbol{a}}_i} $的时延信息转换为DOA信息, $ {{\boldsymbol{a}}_i} $的表达式可变为$ {{\boldsymbol{a}}_i} = {{\boldsymbol{a}}_i}({\varphi _i}, {\theta _i}) $, 该矢量描述了阵列的空间响应特性.

      结合(4)式、(6)式和(8)式, 对阵列流型矩阵$ {\boldsymbol{A}_m} $重构为

      结合(9)式与(6)式可知, 源信号矩阵$ {\boldsymbol{S}}(t) $经由阵列流型矩阵$ {{\boldsymbol{A}}_m} $混合计算得到接收矩阵$ {{\boldsymbol{X}}_m}(t) $, 则依照$ {{\boldsymbol{A}}_m} $的功能, 也可将之称为混合矩阵. 当$ {\boldsymbol{S}}(t) $维度为$ p \times N $时, $ {{\boldsymbol{A}}_m} $维度为$ (n \times n) \times p $, 则$ {{\boldsymbol{X}}_m}(t) $维度为$ (n \times n) \times N $.

    • 盲信号分离作用是分离接收信号, 并重构入射信号. 按照统计独立原则将多通道的观测信号建立代价函数, 通过优化算法将其分解为多个独立分量, 以此实现源信号复原. 结合图1所示信号混合与盲信号分离算法过程与(6)式, 分离模型公式可表示为

      式中, $ m = 1, {\text{ }}2 $, 代表图2所示阵列的序号; $ {{\boldsymbol{X}}_m} $为(6)式的接收矩阵$ {{\boldsymbol{X}}_m}(t) $; $ {{\boldsymbol{W}}_m} $是盲信号分离矩阵; $ {\hat {\boldsymbol{S}}_m} $是盲信号分离算法从接收矩阵$ {{\boldsymbol{X}}_m} $分离出的信号, $ {\hat {\boldsymbol{S}}_m} $是(6)式中源信号矩阵$ {\boldsymbol{S}}(t) $的估计.

      为更好地实现盲信号分离效果, 快速不动点算法(fixed-pointed algorithm, FastICA)因其优异性能备受科研人员关注, 其扩展算法有复数 (circular FastICA, c-FastICA)[29], 非圆信号(noncircular FastICA, nc-FastICA)[30]以及去噪的非圆复信号 (modified noncircular FastICA, mnc-FastICA)[27]. mnc-FastICA对于去除噪声有着极佳的效果, 故本文提出KPCA伪白化处理法, 再使用mnc-FastICA来计算分离信号和矩阵, 可以达到更好的数据处理效果.

      KPCA是一种处理非线性数据的有效方法, 利用非线性映射函数将数据从低维空间升维到高维空间, 准确区分不同来源的非线性数据. 考虑到(10)式中$ {{\boldsymbol{X}}_m} $包含的超声波信号为带噪声的非高斯信源信号的特点[31], 噪声和超声波信号难以线性分离, 因此使用KPCA的非线性映射方法将$ {{\boldsymbol{X}}_m} $映射到高维空间, 记为$ \phi :{R^2} \to {R^{\text{D}}}, D \gg 2 $, 这个高维空间称为特征空间$ \mathcal{F} $. KPCA的映射升维基于列向量开展, 故将$ (n \times n) \times N $维接收矩阵$ {{\boldsymbol{X}}_m} $用列向量表示出来:

      将(10)式的所有样本都映射到特征空间$ \mathcal{F} $中, 得到了$ D \times N $新矩阵$ \phi {(}{{\boldsymbol{X}}_m}{)} $:

      接着在特征空间$ \mathcal{F} $中对$ \phi {(}{{\boldsymbol{X}}_m}{)} $做PCA降维. 首先对$ \phi {(}{{\boldsymbol{X}}_m}{)} $做中心化处理, 即$ \displaystyle\sum\nolimits_{i = 1}^N {\phi {(}{{\boldsymbol{x}}_i}{)}} = 0 $. 计算空间$ \mathcal{F} $$ \phi {(}{{\boldsymbol{X}}_m}{)} $协方差矩阵:

      式中, $ {{\boldsymbol{C}}_m}_\mathcal{F} $$D \times D$矩阵.

      求解协方差矩阵的特征值. 求解特征值公式如下:

      式中, $ {{\boldsymbol{p}}_m} $$D$维列向量, 是特征空间$ \mathcal{F} $的权重向量, 即为主成分方向; ${\lambda _m}$为特征值. 结合(13) 式和(14)式求解特征值, 将(13)式中${1 {/ } N}$简化, 推导出特征值求解公式如下:

      式中, $ \phi ({{\boldsymbol{x}}_i})\phi {({{\boldsymbol{x}}_i})^{\text{T}}} $为低维空间到高维的映射关系, 本文使用多项式核函数$ k(x, y) = \phi \left( x \right)\phi {\left( y \right)^{\text{T}}} $作为映射关系, 其公式如下:

      式中, $ x $$ y $均为(6)式中接收矩阵$ {{\boldsymbol{X}}_m} $的元素, $ b \geqslant 0 $, $ d $为自然数.

      在对$ \phi {(}{{\boldsymbol{X}}_m}{)} $进行PCA降维时, 只需计算较大特征值对应的特征向量, 无需计算特征值0对应的特征向量. 因此, ${\lambda _m} \ne 0$时, 由(15)式推导出特征向量公式如下:

      式中, 等号右侧$ \left[ {\phi {{({{\boldsymbol{x}}_i})}^{\text{T}}}{{\boldsymbol{p}}_m}} \right] $为标量, 则特征向量$ {{\boldsymbol{p}}_m} $可表示为$ \phi ({{\boldsymbol{x}}_i}) $的线性组合, 可表示为

      式中, $ {{\boldsymbol{\beta}} _m} $N维列向量${\left[ {\begin{array}{*{20}{c}} {{\beta _1}}&{{\beta _2}}& \cdots &{{\beta _N}} \end{array}} \right]^{\text{T}}}$.

      将(18)式代入(16)式, 得

      对(19)式等号两侧都左乘矩阵$ \phi {({{\boldsymbol{X}}_m})^{\text{T}}} $, 得

      定义矩阵$ {{\boldsymbol{K}}_m} = \phi {({{\boldsymbol{X}}_m})^{\text{T}}}\phi ({{\boldsymbol{X}}_m}) $, 则$ {{\boldsymbol{K}}_m} $$ N \times N $对称半正定阵, 将$ {{\boldsymbol{K}}_m} $代入(19)式得

      式中, 矩阵$ {{\boldsymbol{K}}_m} $的元素$ {{\boldsymbol{K}}_m}_{ij} = \phi {({x_i})^{\text{T}}}\phi ({x_j}) $由(6)式中接收矩阵${{\boldsymbol{X}}_m}$的低维向量经过(15)式计算得到高维向量后, 计算高维向量间的点积. 这种无需定义映射关系, 只需定义特征空间中向量点积的方式称为核技巧(kernel trick).

      通过上述步骤推导出$ \phi {(}{{\boldsymbol{X}}_m}{)} $的协方差矩阵$ {{\boldsymbol{C}}_m}_\mathcal{F} $, 即为${\boldsymbol{K}}_m^{\text{T}}$. 对$ {{\boldsymbol{C}}_m}_\mathcal{F} $做中心化处理. 定义一个元素值为$ 1/N $$ N \times N $维度矩阵为unit, 中心化处理即$ {\boldsymbol{K}}_m^{\text{T}} - {\text{unit}} \cdot {\boldsymbol{K}}_m^{\text{T}} - {\boldsymbol{K}}_m^{\text{T}} \cdot {\text{unit}} + {\text{unit}} \cdot {\boldsymbol{K}}_m^{\text{T}} \cdot {\text{unit}} $得到中心化处理后的${\boldsymbol{K}}_m^{\text{T}}$.

      使用奇异谱分析法识别接收信号中的超声波源数量. 将${\boldsymbol{K}}_m^{\text{T}}$第一行数据转换为轨迹矩阵$ {{\boldsymbol{Y}}_m}({\boldsymbol{K}}_m^{\text{T}}(:, 1)) $, 构建出汉克矩阵. 对$ {\boldsymbol{Y}}_m^{\text{T}}{{\boldsymbol{Y}}_m} $进行特征值分解, 得到特征值, 计算特征值累积贡献率, 即可确定超声波信号数量$ {{\boldsymbol{p}}_m} $.

      对中心化处理后的信号矩阵${\boldsymbol{K}}_m^{\text{T}}$做特征值分解$ {\boldsymbol{K}}_m^{\text{T}} = {{\boldsymbol{U}}_m}{{\boldsymbol{\varLambda}} _m}{\boldsymbol{V}}_m^{\text{T}} $. 特征值依照降序排列得到$ \left[ {{\lambda _m}_1 \gg {\lambda _m}_2 \geqslant , \cdots , \geqslant {\lambda _m}_p, 0, \cdots 0} \right] $, 对特征谱$ N - p $个较小的值取平均, 作为噪声的近似, 估计噪声方差$ \sigma \approx \displaystyle\sum\nolimits _{p}^{N}{\lambda }_{m}{}_{p}/(N-p) $. 取信号子空间特征值为$ {\sigma }_{m}{}_{i}={\lambda }_{m}{}_{i}-{\sigma }_{m} $, 信号的特征值矩阵和向量阵分别为$ {{\boldsymbol{\varLambda}} _m}_p = {\text{diag\{ }}{\sigma _m}_1{, }{\sigma _m}_2, \cdots , {\sigma _m}_p{\text{\} }} $, $ {{\boldsymbol{U}}_m}_p = [{{\boldsymbol{u}}_m}_1,{{\boldsymbol{u}}_m}_2, \cdots , {{\boldsymbol{u}}_m}_p] $, 构造KPCA伪白化矩阵$ {{\boldsymbol{V}}_m}_s = {\boldsymbol{\varLambda}} _{mp}^{ - 1}{\boldsymbol{U}}_{mp}^{\text{T}} $. 即可得到伪白化处理的数据:

      建立初始化单位阵$ {{\boldsymbol{W}}_m} = {{\boldsymbol{I}}_p} $, 表示为$ p \times p $单位阵. 对于$ {{\boldsymbol{W}}_m} $各列向量$ {{\boldsymbol{w}}_m}_i $, 依次代入以下mnc-FastICA的代价函数进行迭代:

      式中, $ e $为迭代序号; $ {{\boldsymbol{y}}_m} = {\boldsymbol{w}}_{mi}^{\text{T}}{{\boldsymbol{q}}_m} $; $ g $$ G(u) $一阶导, $ G(u) = \sqrt {0.1 + u} $是非线性函数.

      迭代循环(23)式, 即可得到(10)式的分离矩阵$ {\boldsymbol{W}}_m^{} $. 依照(10)式计算得到(6)式的源信号矩阵$ {\boldsymbol{S}}(t) $的估计$ {\hat {\boldsymbol{S}}_m} $.

    • 本文提出的信号提取与定位方法, 提取工作由第2.2节的盲信号分离算法完成, 得到了分离信号$ {\hat {\boldsymbol{S}}_m} $和分离矩阵$ {\boldsymbol{W}}_m^{} $. 本节则对上述结果进行深层次应用, 如3.1节所示TDOA计算步骤, 利用分离信号$ {\hat {\boldsymbol{S}}_m} $, 计算出TDOA; 3.2节所示DOA计算步骤利用分离矩阵$ {\boldsymbol{W}}_m^{} $, 经过矩阵变换后得到了混合矩阵$ {\boldsymbol{A}}_m^{} $, 从混合矩阵$ {\boldsymbol{A}}_m^{} $的元素比值推导出角度信息. 此外, 本节还提出了基于2.2节所提算法的TDOA或DOA单任务定位方法.

    • 2.2节所示KPCA伪白化处理后盲信号分离算法过程, 可以从接收信号矩阵分离出有效信号. 假设图2的两个阵列O1O2接收源信号矩阵$ {\boldsymbol{S}} $的一个信号$ {{{s}}_1} $, $ {{{s}}_1} $维度为$1 \times N$, 阵列得到了接收矩阵$ {{\boldsymbol{X}}_{1}} $$ {{\boldsymbol{X}}_{2}} $, 经盲信号分离算法分离出信号$ {\hat {{s}}_{{1}{{\text{O}}_{1}}}} $$ {\hat {{s}}_{{1}{{\text{O}}_{2}}}} $, 这两个信号是从阵列O1O2的接收矩阵提取出的, 均为源信号$ {{{s}}_1} $的估计. 信号$ {\hat {{s}}_{{1}{{\text{O}}_{1}}}} $传播距离为$ {r_1} $, 信号$ {\hat {{s}}_{{1}{{\text{O}}_{2}}}} $传播距离为$ {r_2} $, 传播距离不同, 则信号到达阵列的时间便有区别, 则通过GCC方法处理$ {\hat {{s}}_{{1}{{\text{O}}_{1}}}} $$ {\hat {{s}}_{{1}{{\text{O}}_{2}}}} $即可得到信号$ {{{s}}_1} $到达阵列O1O2的TDOA.

      分离信号$ {\hat {{s}}_{{1}{O_{1}}}} $, $ {\hat {{s}}_{{1}{O_{2}}}} $与源信号$ {{{s}}_1} $的特性差异可由下式表示:

      式中, $ {\tau _1} $$ {\tau _2} $分别为源信号$ {{{s}}_1} $到达阵列O1O2的时延.

      计算(24)式中两个信号的互相关函数. 互相关函数$ {R_{{{\hat {\boldsymbol{s}}}_{{1}{O_{1}}}}{{\hat {\boldsymbol{s}}}_{{1}{O_{2}}}}}}(\tau ) $定义为

      式中, $ {\tau _1} $$ {\tau _2} $是两个信号的时延, $ {R_{{{\hat {\boldsymbol{s}}}_{{1}{O_{1}}}}{{\hat {\boldsymbol{s}}}_{{1}{O_{2}}}}}}(\tau ) $是互相关函数.

      互相关的频域表示为

      式中, $ {{\boldsymbol{F}}_{{{\hat {\boldsymbol{s}}}_{{1}{{{O}}_{1}}}}}}(f) $$ {\boldsymbol{F}}_{{{\hat {\boldsymbol{s}}}_{{1}{{{O}}_{2}}}}}^{}(f) $分别为$ {\hat {{s}}_{{1}{O_{1}}}} $$ {\hat {{s}}_{{1}{O_{2}}}} $的傅里叶变换; $ {\mathcal{F}^{ - 1}} $表示逆傅里叶变换.

      为提高时间差估计精度, 引入加权函数$\psi (f)$, 对互相关函数加权处理, 加权后的互相关函数$ R_{{{\hat {\boldsymbol{s}}}_{{1}{O_{1}}}}{{\hat {\boldsymbol{s}}}_{{1}{O_{2}}}}}^{{\text{GCC}}}(\tau ) $如下:

      式中, 加权函数使用ROTH加权, 能够抑制噪声影响. 其表达式为

      寻找加权互相关函数$ R_{{{\hat {\boldsymbol{s}}}_{1{O_1}}}{{\hat {\boldsymbol{s}}}_{1{O_2}}}}^{{\text{GCC}}}(\tau ) $的峰值位置, 即为估计的到达时间差$\tau $:

      通过上述步骤, 利用盲信号分离方法获得信号到达不同阵列的TDOA. 由于TDOA定位需要4个不同空间位置的传感器, 则按照上述步骤计算4个阵列接收信号$ {{{s}}_1} $并提取时间差, 并按如下方式建立双曲面方程组:

      式中, PD点坐标$ {\text{PD}}({l_i}, {h_i}, {g_i}) $, 局放超声波信号 向四周传播, 4个传感器阵列布放在不同位置上, 第$j$个阵列的坐标$ {O_j}({l_j}, {h_j}, {g_j}) $, 超声波信号抵达第一个阵列的所需时间为$ {t_0} $, 波速为$ v $. $ {\text{Di}}{{\text{s}}_{ij}} $为PD源坐标与第$j$个传感器阵列之间的距离, 由$ {\text{Di}}{{\text{s}}_{ij}} = {\left( {{l_i} - {l_j}} \right)^2} + {\left( {{h_i} - {h_j}} \right)^2} + {\left( {{g_i} - {g_j}} \right)^2} $计算可得; $ {\Delta }{t_{ij}} $为信号$ {{{s}}_1} $抵达不同阵列的时间差, 不同阵列接收矩阵经过第2.2节所示步骤分离出超声波信号, 再由(24)式—(29)式计算信号TDOA, 即$ \Delta {t_{ij}} $. 求解(30)式即可得到PD源坐标.

    • 如(6)式与(10)式所示, 通过求解分离矩阵$ {\boldsymbol{W}}_m^{} $使得$ {\hat {\boldsymbol{S}}_m} = {\boldsymbol{W}}_m^{\text{T}}{{\boldsymbol{X}}_m} = {\boldsymbol{W}}_m^{\text{T}}[{{\boldsymbol{A}}_m}{\boldsymbol{S}}(t) + {{\boldsymbol{N}}_m}(t)] $, 得到分离信号$ {\hat {\boldsymbol{S}}_m} $. (6)式与(10)式联立可得

      式中, 分离矩阵$ {\boldsymbol{W}}_m^{} $的目的是筛选有效信号并且抑制干扰噪声, 其物理意义为在寻找一个向量在空间信号中筛选出有效信号, 即在空域中寻找一个角度使得有效信号值最大, 此时可认为向量含有角度信息. 此时, 将DOA视为一种特殊的盲信号分离过程.

      考虑到盲信号分离结果存在幅相以及顺序不确定性, 本文分别使用幅相矩阵$ {{\boldsymbol{D}}_m} $和置换矩阵$ {{\boldsymbol{P}}_m} $指代幅相模糊和顺序模糊. 分离信号$ {\hat {\boldsymbol{S}}_m} $与源信号$ {\boldsymbol{S}}(t) $在理想情况下是一致的, 由(31) 式可得$ {\boldsymbol{W}}_m^{\text{T}}{{\boldsymbol{A}}_m} = 1 $.

      矩阵$ {{\boldsymbol{D}}_m} $$ {{\boldsymbol{P}}_m} $发生在盲信号分离过程, 导致分离信号产生误差, 该过程可以表示为

      式中, $ {{\boldsymbol{D}}_m} $代表幅相模糊, 可用对角阵表示$ {{\boldsymbol{D}}_m} = {\text{diag}}(\begin{array}{*{20}{c}} {{\gamma _{m0}}, }&{{\gamma _{m1}}, }&{ \cdots , }&{{\gamma _{mp}})} \end{array} $; 代表顺序模糊的置换矩阵$ {{\boldsymbol{P}}_m} $表示顺序变化, 当顺序发生变化时, $ {{\boldsymbol{P}}_m} $的对应元素值为1, 无变化时, $ {{\boldsymbol{P}}_m} $对应元素值为0. 在不考虑噪声的情形下, (32)式可忽略$ {\boldsymbol{W}}_m^{\text{T}}{{\boldsymbol{N}}_m}(t) $, 得$ {{\boldsymbol{D}}}_{m}{{\boldsymbol{P}}}_{m}={{\boldsymbol{W}}}_{m}^{\text{T}}{{\boldsymbol{A}}}_{m} $.

      由于(23)式得到的分离矩阵$ {\boldsymbol{W}}_m^{} $通常不为方阵, 计算$ {\boldsymbol{W}}_m^{} $的伪逆表达式为

      式中, 置换矩阵$ {{\boldsymbol{P}}_m} $的逆$ {\boldsymbol{P}}_m^{ - 1} $仍为置换矩阵, $ {{\boldsymbol{A}}_m}{\boldsymbol{P}}_m^{ - 1} $表示对混合矩阵$ {{\boldsymbol{A}}_m} $的部分列交换位置, 将对角阵$ {\boldsymbol{D}}_m^{ - 1} $代入(33)式可得

      式中, 可看出$ {{\boldsymbol{B}}_m} $第一行元素为矩阵$ {{\boldsymbol{P}}_m} $作用于$ {{\boldsymbol{D}}_m} $的结果, 即计算分离矩阵$ {\boldsymbol{W}}_m^{} $伪逆并提取第一行可得矩阵$ {{\boldsymbol{D}}_m}{{\boldsymbol{P}}_m} $, 其表达式如下:

      由(33)式与(35)式联立可得混合矩阵$ {{\boldsymbol{A}}_m} $:

      由(9)式可知混合矩阵$ {{\boldsymbol{A}}_m} $的每个导向矢量单元包含方位角$ \varphi $与俯仰角$ \theta $, 将混合矩阵$ {{\boldsymbol{A}}_m} $依照$ n \times n $阵列形状, 由列向量重塑为同维度矩阵. 在第2.1节中模型设定窄带信号以平面波入射到阵列, 则可推测出每个阵元接收到的信号方向是一致的. 结合(9)式, 按(37)式计算矩阵$ {{\boldsymbol{A}}_m} $相邻元素的比值即可计算出信号DOA:

      式中, $ d $为阵列阵元的中心间距, 如图2所示. 求解方程组(37)即可得到超声波信号的方位角$ \varphi $与俯仰角$ \theta $. 在三维空间中, 定位坐标至少需要3个传感器阵列, 即3个方位角与俯仰角, 通过求解3个空间直线的交点, 得到目标位置.

    • 图1中源信号矩阵$ {\boldsymbol{S}} $的信号$ {s_1} $传播到图2所示阵列O1O2为例, 由(30)式得到的超声波$ {s_1} $传递到两个阵列的距离差表达式为

      式中, $ {e_{21}} $为TDOA测量误差导致的距离误差. 该值可以按(39)式推导:

      式中, $ {e_{21}} $与信号源到阵列的TDOA参数有关. TDOA测量值误差服从正态分布, 即$ {e_{21}} $均值为0, 标准差为$ \sigma $.

      除了时间差、方位角$ \varphi $与俯仰角$ \theta $可提供信号源的方向性约束. 则DOA与信号源$ {\text{PD}}({l_i}, {h_i}, {g_i}) $、阵列坐标的几何关系表示为

      式中, $ {e_\varphi } $$ {e_\theta } $均符合均值为0, 标准差分别为$ {\alpha _\varphi } $$ {\alpha _\theta } $正态分布; $ {\varphi _{ij}} $$ {\theta _{ij}} $分别为PD源坐标$ ({l_i}, {h_i}, {g_i}) $到第$j$个阵列$ ({l_j}, {h_j}, {g_j}) $的方位角和俯仰角, $i = 1$; $j = 1, 2$.

      依据(38)式和(39)式的正态分布时间误差, 建立其概率密度函数:

      将方位角和俯仰角信号纳入优化过程, 其概率密度函数如下:

      对(41)式建立似然函数并取对数, 得到公式如下:

      对(42)式建立似然函数并取对数得

      求解位置的最大化对数似然函数, 通常需要最小化目标函数, 综合(43)式和(44)式, 即可以得到TDOA/DOA定位的最大似然估计式:

      由于(45)式非线性方程难以解析. 本文引入非洲秃鹫算法[28]进行求解. 将(45)式视为信号源PD节点的估计函数$ f({P_i}) $, 若是$ p $个PD点即得到$ p $个估计函数值, 使用智能算法做适应度估计表达式如下:

      根据上述适应度函数, 即可求解最优的PD源坐标值.

    • 图3为TDOA/DOA混合定位算法步骤, 本文所提方式的核心目标是通过信号解混、TDOA和DOA估计以及建立联合方程, 以求解PD源坐标. 具体算法步骤如下.

      输入: 接收信号矩阵$ {{\boldsymbol{X}}_m}(t) $, 维度为$ (n \times n) \times N $, 阵列形状$ n \times n $.

      输出: 目标估计位置$ {\text{PD}}({l_i}, {h_i}, {g_i}) $.

      第1阶段:

      1)对信号矩阵$ {{\boldsymbol{X}}_m}(t) $的所有样本映射到高维空间, 依照(12)式方式得到新矩阵$ \phi {(}{{\boldsymbol{X}}_m}{)} $;

      2)对$ \phi {(}{{\boldsymbol{X}}_m}{)} $求解协方差矩阵的特征值, 并中心化处理, 得到${\boldsymbol{K}}_m^{\text{T}}$;

      3)将${\boldsymbol{K}}_m^{\text{T}}$的第1行数据转为轨迹矩阵$ {{\boldsymbol{Y}}_m} ({\boldsymbol{K}}_m^{\text{T}}(:, 1)) $, 构建汉克矩阵, 对$ {\boldsymbol{Y}}_m^{\text{T}}{{\boldsymbol{Y}}_m} $进行特征值计分解, 计算特征值累积贡献率, 确定信号数量;

      4)对${\boldsymbol{K}}_m^{\text{T}}$做特征值分解, 计算KPCA伪白化矩阵$ {{\boldsymbol{V}}_m}_s $, 求解伪白化数据;

      5)迭代运行(23)式, 可得分离矩阵$ {\boldsymbol{W}}_m^{} $, 依据(10)式, 计算分离信号$ {\hat {\boldsymbol{S}}_m} $.

      算法第1阶段使用结合KPCA的mnc-FastICA, 是由于传感器接收到的信号包括多种混合成分, 包括PD超声波、环境噪声. 这些因素使得接收信号具有非线性特征, 给信号的特征提取带来了极大挑战. 而KPCA可以在高维空间挖掘信号的非线性结构, 抑制噪声干扰, 提取出PD超声波主成分. 但降维后的主成分仍包含一定的背景噪声和非目标信号成分. 为进一步提纯PD信号, 需要基于PD信号与其他成分的统计特征差异, 利用mnc-FastICA进行独立信号分离, 以及分离矩阵. 二者结合在局放超声信号处理中可以有效提高信号的分离质量, 为后续TDOA和DOA计算提供更可靠的数据基础.

      第2阶段:

      6)对从第1阶段计算得到两个阵列的分离矩阵做伪逆处理, 计算幅相模糊矩阵;

      7)如(36)式和(37)式计算混合矩阵以及DOA.

      第3阶段:

      8)对第1阶段计算得到两个阵列的分离信号, 计算互相关函数;

      9)加权函数处理对互相关函数;

      10)寻找峰值位置, 估计TDOA. 算法第2阶段利用了mnc-FastICA的分离矩阵计算DOA. 分离矩阵的目的是选择最优的矩阵筛选PD信号同时尽量抑制噪声, 其实质是在空域中选择一个PD信号含量最强的角度, 该角度对应的信号成分最为接近真实PD信号, 因此能够利用分离矩阵计算出DOA. 算法第3阶段利用PD信号的传播时间差异, 计算TDOA. 在传统方法中, 由于传感器接收到的信号往往包含环境噪声, 直接计算TDOA 可能会导致较大的时间偏差. 然而, mnc-FastICA分离的信号更加纯净, 其背景噪声和干扰成分被有效抑制, 因此基于GCC的TDOA计算能够获得更高的精度计算.

      第4阶段:

      11)如(41)式与(42)式建立误差概率密度函数, 并建立似然函数以及取对数;

      12)如(45)式建立最大似然估计函数;

      13)利用非洲秃鹫算法求解最优的PD源坐标.

      在算法第4阶段, 结合TDOA和DOA误差概率密度函数, 建立最大似然估计函数, 以最小化估计误差, 从而在噪声干扰下得到最符合真实信号传播特性的参数估计. 通过最大化似然函数可有效较低噪声和非理想因素对定位结果的影响. 但最大化似然函数通常涉及非线性优化问题, 可能陷入局部最优解, 导致定位误差变大. 因此本文引入非洲秃鹫算法, 以全局搜索策略求解最优坐标.

    • 通过对上述所提算法进行详细分析, 可统计处算法所需计算复杂度. 具体而言, 在第1阶段, 首先使用KPCA算法处理维度为$ M \times N $的数据, KPCA需构建$ N \times N $的核矩阵, 其中矩阵元素由预先选定的多项式核函数计算得到, 核函数计算的时间复杂度为$ O(1) $, 则构造核矩阵的复杂度为$ O({N^2}) $; 其次进行中心化处理, 以消除数据均值的影响, 其复杂度为$ O({N^2}) $; 下一步对中心化后的核矩阵进行特征值分解, 其复杂度为$ O({N^3}) $; 综合可得, KPCA运算阶段的时间复杂度为$ O({N^3}) $. KPCA处理后的数据维度为$ p \times N $, p为样本数, 迭代计算mnc-FastICA代价函数的复杂度为$ O(Np) $, 则总体复杂度为$ O(Np + {N^3}) $. 第2阶段, 通过计算矩阵运算求解DOA, 矩阵伪逆使用奇异值分解, 其复杂度为$ O({p^2}M) $. 第3阶段, 利用GCC方法计算TDOA, 傅里叶变换和逆傅里叶变换的复杂度均为$ O(N\log N) $, 频谱逐点乘法的复杂度为$ O(N) $; 第4阶段, 最大化似然函数的复杂度为$ O(N) $, 非洲秃鹫算法的复杂度为$ O({\text{Thd}}) $, 其中, 搜索空间维度为$ d $, 迭代次数为$ T $, 种群大小为$ h $. 综合上述分析, 本文所提算法的整体计算复杂度可简化为$ O({N^3} + Np + {p^2}M + {\text{Thd}}) $. 经典算法如MUSIC算法的总复杂度为$ O({M^2}N + {M^3}) $, CAPON算法复杂度为$ O({M^2}N + {M^3}) $. 对比经典算法, 本文所提方法的复杂度依赖于信号样本大小$ N $, 在大样本情况下, 本文所提算比MUSIC和CAPON复杂度高.

      理论上, 算法复杂度是评估算法性能的重要指标. 但在实际中, 复杂度忽略了常数因子和低阶项, 并不能完全代表算法的实际运行速度、计算时间和内存占用情况. 在实际定位算法应用中, MUSIC, CAPON算法需要3—4个阵列, 每个阵列使用至少8个传感器[16,32], 带来了庞大的硬件、计算以及存储开销, 况且多通道高速度采集卡的采样带宽有限, 难以满足大规模阵列的数据需求, 进一步限制了这些算法的应用. 鉴于此, 本文所提算法综合考虑定位需求以及硬件开销的基础上, 通过设计计算流程, 实现了更高效的信号处理方案, 保证了较高的定位效果, 更加适用于实际场景.

    • 本研究揭示了局放超声波信号的提取与定位之间的强关联性. 基于第2.2节所提KPCA-mnc-FastICA算法, 通过对接收信号的盲信号分离及分离矩阵解析, 实现了TDOA和DOA参数的协同提取. 仿真输入信号采用双指数衰减振荡函数模拟的超声波信号, 长度为5000, 加入机械振动(振动频率低于500 Hz, 幅值远高于超声波)以及高斯白噪声, 方形阵列接收超声波信号, 阵元间隔小于超声波半波长.

      图4所示, 原始超声波信号(蓝色曲线)在12 μs时起始, 其幅值1.5 mV远低于机械振动的幅值30 mV, 导致混合后的接收信号(荧光绿色曲线)中超声波特征被完全掩蔽. 经所提算法处理后, 分离信号(红色线条)精准还原了超声波起始时间及包络形态, 但非信号区间(0—12 μs及16—50 μs)呈现均匀噪声分布. 此现象源于源信号建模特性: (6)式定义的混合模型中, 非超声波时段(幅值为0)的高斯噪声经线性叠加后, 在接收端形成全域背景噪声基底, 使得分离信号在无超声波区间仍保留噪声分量.

      对于DOA估计, 现有的研究均基于大阵列的DOA估计性能, 如文献[16,32]使用了至少8个阵列, 考虑到超声波信号的宽频带特性, 如若组成多阵列方案则硬件开销极大. 基于上述情况, 本节将使用多种方法对比不同数量传感器阵元情况下的DOA估计精度, 如图5所示. 使用了MUSIC, Capon等经典算法来对比本文3.2节所提方法的DOA估计性能, 并用克拉美-罗界(Cramer-Rao bound, CRB)反映阵列DOA的估计性能. 其中, DOA性能评估使用均方根误差(root mean square error, RMSE), 表达式为

      其中, $ {\theta _k} $为声源$ k $入射俯仰角, RMSE结果在R = 100次蒙特卡罗方法计算中获得.

      图5所示, 随着传感器阵列规模从4阵元扩展至36阵元, 角度估计均方根误差(RMSE)呈现非单调递减趋势. 当阵元数为4时, 传统MUSIC与Capon算法因空间分辨率不足导致估计失效(RMSE > 100°), 而本方法在极小阵列条件下仍保持度级精度(RMSE < 7°). 当阵元数增至9时, 经典算法性能逐步接近实用阈值(RMSE < 15°), 此时本方法精度(RMSE ≈ 2°)仍显著优于对比方法. 值得注意的是, 当阵元数超过16时, 所有方法均逼近CRB下界, 表明阵列规模达到临界点后, 硬件扩展对精度提升贡献有限. 这一现象验证了本方法在小规模阵列(4阵元)场景下的独特优势, 其核心创新在于通过KPCA-mnc-FastICA算法增强时间-空间特征耦合解析能力.

      为进一步验证本文所提方法分别在TDOA和DOA估计上的性能, 图6图8从波形的相似程度(normalized correlation coefficient, NCC)、信噪比(signal-to-noise ratio, SNR)、时间差对比以及定位结果来评估TDOA估计性能, 图9为不同SNR下的角度估计误差来评估DOA估计性能.

      图6所示, 通过NCC量化分离信号与源信号的波形相似性. 实验表明, 在输入信噪比(SNR)从–5 dB增至20 dB时, NCC值由0.92单调上升至0.9999(SNR = 10 dB时已达0.999), 证实本方法具有优异的信号保真特性. 值得关注的是, 分离信号SNR在输入SNR变化时保持稳定, 凸显算法在强噪声干扰下的鲁棒性.

      图7对比不同TDOA提取方法的时延误差, 基于能量累积法的误差高达2 ms (SNR = –5 dB), 而GCC与互相关法的误差分别稳定在0.8 μs与0.9 μs内. 其原因是所提方法的分离信号存在一定的噪声分量, 能量累积法虽然结合滤波器过滤噪声, 但仍然对结果产生了负面影响.

      使用GCC提取时延用于TDOA定位算法, 结合图8的定位结果, 验证GCC方法在分离信号时延提取的适用性. 依照文献[1719]提到TDOA定位法需要4个不同位置的传感器, 本节使用4个传感器阵列进行信号接收、分离以及TDOA计算的仿真. 结合(30)式建立TDOA-chan定位模型, 在1 m × 1 m × 1 m空间中得到的定位误差为1 cm, 信号源位置与估计位置接近, 能够满足定位需求. 综合图6图8的分析结果可见, 本文所提方法分离出与源信号NCC相近的波形, 保留了源信号到达不同阵列的时间信息, 使得GCC计算分离信号的TDOA更为符合实际情况.

      图9揭示方位角与俯仰角的估计的误差变化规律. 依照3.2节所示算法, KPCA-mnc-FastICA计算得到分离矩阵, 从分离矩阵的元素值比计算出角度. 图9图5采用了相同的对比方法, 当阵列为$3 \times 3$时, 对比方位角与俯仰角的估计精度. 由图5的对比结果可得, 在使用$2 \times 2$阵列时, 本文所提方法的在角度估计上更为精确, 而经典方法的估计值完全偏离实际值从而导致无法在小阵列上使用DOA估计方法.

      若使用$3 \times 3$阵列这一发挥经典算法性能的硬件配置, 图9给出了更进一步的算法性能对比. 随着源信号SNR的增长, 3.2节所提算法在方位角(图9(a))和俯仰角(图9(b))的估计精度逐步上升, 而这一上升趋势与CRB曲线走向一致, 表明了其性能优势. Capon算法的估计值误差较大且比较稳定, MUSIC算法的角度估计性能随着源信号SNR增长而增大, 但表现较本文所提算法有所不如. 此外, 在源信号SNR低于5 dB时, 本文所提算法估计能力仍高于其他算法. 图6图9的分析结果更进一步看出本文所提算法在TDOA和DOA信息提取上的优势, 通过盲信号分离算法计算分离信号和分离矩阵, 深度挖掘分离信号和TDOA、分离矩阵和DOA的耦合关系, 实现了超声波信号提取与定位任务的强关联.

      基于上述分析结果, 验证了本文所提方法在TDOA和DOA的估计优势. 为更进一步验证TDOA/DOA混合定位在精度以及硬件开销上的优势, 图10给出了4个$2 \times 2$传感器阵列下本文所提TDOA, DOA和TDOA/DOA混合方法定位精度对比, 随着源信号SNR的增大, TDOA/DOA混合定位性能保持优异稳定, 定位误差小于其他对比方法, 且定位精度毕竟CRB界. TDOA-chan定位效果优于DOA定位, 这主要是和阵列传感器阵元数量有关, 如图5所示, 阵元数越大, 则DOA估计性能越高, 反之则相反. 图中表明, TDOA/DOA混合方法的CRB界最低, 可以佐证该方法的估计准确上限最高. 此外, 随着源信号SNR的变化, 图中各方法的定位性能较为稳定, 证明本文所提方法在TDOA和DOA的估计优势, 能够适应高噪声的定位任务.

      图11给出了在源信号SNR为10 dB, 4个阵列时, 不同阵元数量阵列下的定位精度对比. 可以明显看出, 随着SNR增长, DOA定位精度逐渐上升, 此结果与图10中由于使用了$2 \times 2$传感器阵列导致DOA定位精度低于TDOA-chan的结果一致, 这表明了DOA定位方法发挥性能优势的前提是使用的大阵元数量的阵列, 但超声波信号频率高且带宽较大, 大阵列所需硬件开销及计算资源较多, 导致了成本居高不下, 此现象也从反面论证了本文所提方法在DOA估计上的优势, 能以小阵元数量的阵列发挥较高性能. TDOA-chan定位精度维持稳定的性能表现, 表明本文所提算法在TDOA估计时无需使用大阵元阵列. 而TDOA/DOA混合定位方法综合了TDOA和DOA定位法的优势, 其定位精度即使在阵列数量较大的时候也优于DOA方法.

      图12给出了在源信号SNR为10 dB时, 不同$2 \times 2$阵列数量下的定位精度对比. 阵列均为两两正交, 数量从2—4, TDOA-chan, DOA和TDOA/DOA混合定位方法的定位精度由高到低, 而TDOA/DOA混合定位方法又优于其他两种方法. 在阵列数量为2的时候, TDOA-chan与DOA定位误差较大, TDOA/DOA混合定位法的定位精度较高, 能够满足实际定位需求. 在阵列数量增加为3时, TDOA/DOA混合定位方法的定位精度更加精准.

      图13给出了在源信号SNR为10 dB时, 采用两个正交$2 \times 2$阵列进行定位, 其中定位误差仅为0.8 cm. 结果表明, 在仅使用8个传感器的情况下, 获得的定位位置与真实位置高度吻合, 验证了算法的优越性能. 由此可见, 本文所提TDOA/DOA混合方法能够在以2个$2 \times 2$阵列获取较为准确的定位性能, 这进一步凸显了该方法在降低硬件成本和计算资源上的优势.

    • 本节搭建了局放超声波放电实验平台, 该平台有一个PD源产生超声波信号, 通过3个$2 \times 2$传感器阵列接收声源, PD源为针板电极, 电极间距3 mm, 超声波传感器谐振频率为150 kHz, 频率带宽为100—400 kHz, 试验平台如图14所示, 图14(a)为试验平台示意图, 图14(b)为装置图. PD源坐标和传感器阵列的坐标如表1所示. 该阵列方案可灵活组成4种不同的阵列配置, 包括本文提出使用的2个正交阵列(共8个传感器)方案, 以验证在不同阵列数量条件下的定位精度.

      3个处于不同方位的传感器阵列收集到的PD超声波信号, 经过2.2节所提算法KPCA-mnc-FastICA分离出超声波信号后, 得到了图15(a)(c)所示波形图. 由图15可知, 分离出的波形均存在着明显的噪声, 结合图4图7的分析结果, 分离信号中的噪声并不影响TDOA计算的精准性, 这一结论也可由表2得出.

      以传感器阵列1作为参考点, 使用GCC提取声波到达阵列2和3的时间差, 计算结果如表2所示. 从表2可以看出, 阵列2到阵列1的时间误差t21为–3.8 μs, 阵列3到阵列1的时间误差t31为–5.4 μs, 波动控制在±5 μs, 平均误差为2.34%, 可见本文所提KPCA-mnc-FastICA分离出的超声波信号具有实际环境中, 计算TDOA的能力.

      3.2节所提算法计算超声波DOA, 以空间直角右手坐标系原点为参考点, xoy平面为入射参考面, 计算3个阵列的DOA信息, 其结果如表3所示, 阵列DOA最大误差小于2°.

      结合表2表3的计算结果, 即可按照本文所提TDOA/DOA混合定位方法定位信号源坐标. 依照图12图13所示, 本文所提算法在使用两个阵列时, 可获得较为精准的定位结果, 使用3个阵列的时, 定位结果更加准确, 因此, 本文使用将3个阵列按照不同组合进行定位, 其结果如表4所示. 两个阵列组合的TDOA/DOA混合定位的距离误差分别为1.54 cm, 1.88 cm和2.45 cm, 3个阵列组合的定位误差为1.14 cm, 对比可发现, 3个阵列组合的定位精度相比两个阵列, 其定位精度提升效果有限, 而双阵列组合在减小硬件开销的前提下, 保证了较为精准的定位结果, 而且易于实际环境的部署. 可见, 本文3.1节与3.2节所提基于2.2节KPCA-mnc-FastICA算法的TDOA与DOA估计效果较好, 挖掘了超声波信号所含时间与角度深度关联, 通过KPCA-mnc-FastICA计算出的分离信号和分离矩阵, 实现了TDOA与DOA的同步提取, 在此基础上, 提出了3.3节TDOA/DOA混合定位方法, 能够降低硬件开销且保证定位精度. 可见, 本文所提方法解决了超声波信号提取与定位任务相互独立的现状, 挖掘了信号提取与定位的耦合关系, 更有利于超声波的定位工作.

    • 1)联合特征解析算法创新. 提出基于核主成分分析KPCA-改进型非圆FastICA算法的时间-空域特征协同提取框架. 首先, 提出了KPCA多项式核映射非线性升维及线性降维的伪白化处理方法, 能够保留信号的有效信息同时抑制噪声. 其次, 使用改进型非圆FastICA从信号中分离出超声波信号以及分离矩阵, 并从中解析出超声波的TDOA和DOA参数, 实验结果表明, TDOA平均误差为2.34%, 波动控制在±5 μs, DOA估计误差小于2°.

      2)混合定位理论突破. 结合超声波的TDOA和DOA耦合特征, 建立融合双参数的最大似然估计模型, 联合非洲秃鹫算法, 实现全局最优解快速收敛. 实验表明, 配置2个正交阵列, 共8个传感器时, 便能实现定位功能且误差最小为1.54 cm.

      3)工程优化验证. 揭示传感器阵列规模与定位精度的定量关系. 当采用2个正交方形阵列(总传感器数8)时, 其定位误差对比3阵列方案的1.14 cm误差, 仅增加0.4 cm误差, 但硬件成本降低1/3. 这表明所提方法在小规模阵列(4阵元阵列)场景下具备工程实用价值, 为电力设备紧凑型局放监测装置研发提供了理论支撑.

    参考文献 (32)

目录

/

返回文章
返回