一种利用两阶段学习模型的水下阵列定位方法

张宇, 江鹏, 郭文飞, 张丹, 韩震

张宇, 江鹏, 郭文飞, 张丹, 韩震. 一种利用两阶段学习模型的水下阵列定位方法[J]. 武汉大学学报 ( 信息科学版), 2021, 46(12): 1889-1899. DOI: 10.13203/j.whugis20210466
引用本文: 张宇, 江鹏, 郭文飞, 张丹, 韩震. 一种利用两阶段学习模型的水下阵列定位方法[J]. 武汉大学学报 ( 信息科学版), 2021, 46(12): 1889-1899. DOI: 10.13203/j.whugis20210466
ZHANG Yu, JIANG Peng, GUO Wenfei, ZHANG Dan, HAN Zhen. An Underwater Array Localization Method Using Two-Stage Learning Model[J]. Geomatics and Information Science of Wuhan University, 2021, 46(12): 1889-1899. DOI: 10.13203/j.whugis20210466
Citation: ZHANG Yu, JIANG Peng, GUO Wenfei, ZHANG Dan, HAN Zhen. An Underwater Array Localization Method Using Two-Stage Learning Model[J]. Geomatics and Information Science of Wuhan University, 2021, 46(12): 1889-1899. DOI: 10.13203/j.whugis20210466

一种利用两阶段学习模型的水下阵列定位方法

基金项目: 

国家重点研发计划 2016YFB0501800

详细信息
    作者简介:

    张宇,博士生,主要从事水下数据处理和水下通信研究。zhangziju125@whu.edu.cn

    通讯作者:

    江鹏,博士,研究员,博士生导师。jiangp@whu.edu.cn

  • 中图分类号: P229.5

An Underwater Array Localization Method Using Two-Stage Learning Model

Funds: 

The National Key Research and Development Program of China 2016YFB0501800

More Information
    Author Bio:

    ZHANG Yu, PhD candidate, specializes in underwater data processing and communication. E-mail: zhangziju125@whu.edu.cn

    Corresponding author:

    JIANG Peng, PhD, professor. E-mail: jiangp@whu.edu.cn

  • 摘要: 水下声学(underwater acoustic, UWA)阵列信号处理是常见的水下定位方式之一。针对噪声影响定位精度的问题,提出一种利用两阶段学习模型的定位方法。首先,分别对接收信号的实部和虚部特征进行训练,建立一个基于多层卷积神经网络的学习模型进行降噪处理;然后,构建一个改进的加权延时求和的波束形成器组模型,利用梯度下降准则对各个通道的权重进行调整,得到最优相对时延和最佳角度估计,再通过几何解算得到较为精确的定位信息。仿真实验结果表明,在-25~10 dB的信噪比环境中,所提方法与传统水下声阵列处理相比,应对信噪比变化的鲁棒性更强,且定位精度更高,湖上实验进一步验证了该方法的有效性。
    Abstract:
      Objectives  In the localization method based on underwater acoustic arrays, noise affects the localization accuracy mainly. A localization method using a two-stage learning model is proposed to minimize noise's influence on the localization results.
      Methods  Firstly, a learning model based on a multilayer convolutional neural network is built for noise reduction by training the real and imaginary features of the received signal separately. Secondly, an improved weighted delay summation beamformer group model is constructed, and the gradient descent criterion adjusts the weights of each channel to obtain the optimal relative time delay and the best angle estimation. Finally, more accurate localization information is received by geometric solving.
      Results  Simulation experimental results show that the two-stage noise reduction model has awe-inspiring noise reduction performance compared to conventional underwater acoustic array processing for angular comparisons of 30°, 100°, and 130° in a SNR (signal-to-noise ratio) environment in the range of -25 dB to 10 dB. Lake tests show the same advantage of the proposed method in terms of RMSE (root mean square error) of positioning results.
      Conclusions  The two-stage learning model has excellent robustness in coping with SNR variations, and the received signal after processing by noise reduction has higher localization accuracy after processing based on an improved time-delay and angle estimation model, which can also be applied in practical applications such as indoor acoustic localization and radio radar arrays.
  • 近年来,随着海洋资源的开发和海洋科学技术的发展,水声定位在海洋科学领域的应用越来越广[1]。基于水下声阵列信号处理的定位方法在水下空间的不同位置布放多个传感器组成阵列,利用阵列接收信号再对其进行处理,目的是估计判断主径上的发射信号方向,抑制无用的干扰和噪声。在解算出发射信号到各个阵元间的时间差后,利用方向角度和时延推算出发射单元的位置,从而达到定位的目的。

    目前广泛应用于水下声阵列定位的方法有水下目标波达角(direction of arrival,DOA)估计[2]、水下波束形成技术[3],均借鉴了无线电和声学的阵列处理方式解算水下位置。常用的DOA估计方法包括多信号分类(multiple signal classification,MUSIC)算法[4-6]、旋转不变技术估计信号参数(estimating signal parameters via rotational invariance techniques,ESPRIT)算法[7]、最大似然估计法[8-9]等。上述算法的核心是将噪声假设为高斯过程,通过分离噪声子空间来估计发射信号角度。水下波束形成对多路接收阵信号进行合并处理,抑制非目标方向的干扰并增强目标方向的信号,常见的波束形成算法包括最小均方误差(minimum mean square error,MMSE)波束形成器[10-11]、线性约束最小方差(linearly constrained minimum variance,LCMV)波束形成器[12]。在时延估计中,常采用广义互相关-相位变换(generalized cross correlation phase transformation,GCC-PHAT)算法[13],其核心是通过对水下阵列间信号进行互相关处理,得到使互相关函数最大的值即信号的时间差。但在实际阵列信号接收中,信道变化是一个复杂的过程,会受到噪声的干扰[14]。上述算法在噪声环境下均存在难以克服的问题,例如水下信道的噪声变化并不都是高斯过程,DOA估计算法并不能完整地分离出噪声子空间来估计波达方向角[15-17];波束形成算法大多聚焦于对主径角度的估算,算法的性能过度依赖目标信号的DOA信息[18],而强噪声环境对其处理精度影响十分明显。综上可知,噪声是制约定位精度的主要原因,为了最大程度降低噪声对接收信号的干扰,可应用自适应滤波进行降噪处理,其依靠线性滤波器滤掉噪声频率成分[19],但收敛速率较慢且计算量大,为了克服以上问题,本文将深度学习引入降噪的处理中。

    近年来,深度神经网络(deep neural network,DNN)被广泛应用于图像降噪和语音识别领域,并取得了突破性进展。依据研究侧重的不同可分为聚焦于数据特征方面[20-22]和网络模型方面[23-26]两种。前者将重点放在对接收信号数据特征的训练处理中,根据训练特征不同可大致分为基于掩码估计[27]和基于谱映射估计[28]两类。后者则通过不同结构的神经网络提取噪声特征进行训练,例如文献[29]将信号的理想二值掩蔽(ideal binary mask,IBM)作为降噪神经网络的输入特征,对比信号时频单元的信噪比,学习并估计噪声;文献[30]提出一种基于卷积神经网络的水声通信降噪方法,利用时-频关系对接收信号进行二维化表述,达到降噪的目的。深度学习模型的本质是通过学习大量纯净信号和含噪信号样本数据构造复杂非线性函数,生成新的降噪处理估计。但上述两种模型均未针对水下声学信号进行适应性设计,大多数降噪方式都是直接将语音信号处理中常见的特征作为输入,忽略了水声信号的短时平稳性特点[31],或者直接套用图像降噪的方式,随着噪声的筛除,水声信号中的关键传输信息也会伴随性丢失[32]

    针对以上问题,本文在DOA定位的基础上,提出一种基于降噪学习模型的定位方法,通过构建一个两阶段学习模型来处理水下声阵列接收数据,进而完成定位。在第一阶段模型中,为了克服水下噪声对定位结果的影响,使用神经网络模型对接收信号完成降噪处理,为了最大限度地挖掘水声信号的特征,在特征预处理阶段对其进行分帧处理,将信号的实部和虚部作为训练特征,充分利用了水声信号的短时平稳特性,为了尽可能降低对信号信息正交性的破坏,采用压缩合并和还原策略处理数据。在第二阶段模型中,使用第一阶段输出的降噪还原信号作为第二阶段的输入,通过改进的基于加权延时求和的波束形成方法来完成多空域滤波权值更新,达到对期望波达方向的信号增强,结合角度估计结果和时延估计结果解算发射端位置。该方法显著提升了定位系统对信噪比变化的鲁棒性,改善了系统定位精度。

    图 1为水下声阵列的发射-接收结构示意图,发射端信号经过水下噪声和反射产生多径效应的影响,到达接收阵列,L为阵元间距,θ为入射信号与阵元夹角,c为信号传播速度。

    图  1  水下声阵列的发射-接收结构示意图
    Figure  1.  Transmit-Receive Structure of the Underwater Acoustic Arrays

    图 2展示了两个参考阵元间的结构关系。假设水下声阵列由M个接收阵元组成,有K个信号s1(t),s2(t)sK(t)到达该阵列,则第m个阵元的接收信号可以表示为:

    图  2  参考阵元间的结构关系
    Figure  2.  Structural Relationship of Reference Array Elements
    xm(t)=k=1Ksk(t-τm,k)+nm(t)

    式中,m=1, 2Mt为传播时间;τm,k为第k个信号到达第m个阵元相对于参考阵元的时延;sk(t-τm,k)表示第m个阵元上接收到的第k个发射信号到达波;nm(t)表示第m个阵元上的加性噪声。阵列接收数据按照向量形式可表示为:

    $$\begin{array}{*{20}{l}} {\mathit{\boldsymbol{x}}(t) = \left[ {\begin{array}{*{20}{c}} {{{\rm{e}}^{ - {\rm{j}}w{\tau _{{\rm{1}}, {\rm{1}}}}}}}&{\begin{array}{*{20}{c}} {{{\rm{e}}^{ - {\rm{j}}w{\tau _{{\rm{1}}, {\rm{2}}}}}}}& \cdots \end{array}}&{{{\rm{e}}^{ - {\rm{j}}w{\tau _{{\rm{1}}, K}}}}}\\ {{{\rm{e}}^{ - {\rm{j}}w{\tau _{{\rm{2}}, {\rm{1}}}}}}}&{\begin{array}{*{20}{c}} {{{\rm{e}}^{ - {\rm{j}}w{\tau _{{\rm{2}}, {\rm{2}}}}}}}& \cdots \end{array}}&{{{\rm{e}}^{ - {\rm{j}}w{\tau _{{\rm{2}}, K}}}}}\\ \vdots &{\begin{array}{*{20}{c}} \vdots & \ddots \end{array}}& \vdots \\ {{{\rm{e}}^{ - {\rm{j}}w{\tau _{M, {\rm{1}}}}}}}&{\begin{array}{*{20}{c}} {{{\rm{e}}^{ - {\rm{j}}w{\tau _{M, {\rm{2}}}}}}}& \cdots \end{array}}&{{{\rm{e}}^{ - {\rm{j}}w{\tau _{M, K}}}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{s}}_{\rm{1}}}(t)}\\ {{\mathit{\boldsymbol{s}}_{\rm{2}}}(t)}\\ \vdots \\ {{\mathit{\boldsymbol{s}}_K}(t)} \end{array}} \right] + }\\ {\; \; \; \; \; \; \; \; \; \; \; {{\left[ {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{n}}_{\rm{1}}}(t)}&{{\mathit{\boldsymbol{n}}_{\rm{2}}}(t)} \end{array}}& \cdots &{{\mathit{\boldsymbol{n}}_M}(t)} \end{array}} \right]}^{\rm{T}}}} \end{array}$$

    式中,w为数字域频率。式(2)可简化为:

    x(t)=As(t)+n(t)

    其中,

    $$x(t) = {[\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{x}}_{\rm{1}}}(\mathit{t})}&{{\mathit{\boldsymbol{x}}_{\rm{2}}}(\mathit{t})} \end{array}}& \cdots &{{\mathit{\boldsymbol{x}}_\mathit{M}}(\mathit{t})} \end{array}]^{\rm{T}}}$$
    $$A = [\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{a}}({{\rm{ \mathit{ θ} }}_{\rm{1}}})}&{\mathit{\boldsymbol{a}}({{\rm{ \mathit{ θ} }}_{\rm{2}}})} \end{array}}& \cdots &{\mathit{\boldsymbol{a}}({{\rm{ \mathit{ θ} }}_{\rm{K}}})} \end{array}]$$
    a(θk)=[ejwτ1,kejwτ2,kejwτM,k]T

    式中,θk为第k个发射信号的方位角。

    根据式(3)得到参考接收信号x1(t)x2(t)为:

    x1(t)=As(t)+n1(t)
    x2(t)=As(t)+n2(t)

    时延τ可通过广义互相关(generalized cross correlation,GCC)算法[33]对两个参考阵元间的接收信号峰值检测得到,计算如下:

    Rx1x2(τ)=x1(w)x2*(w)e-jwτdw
    τ̂=argmaxτRx1x2(τ)

    式中,“*”表示矩阵共轭;τ̂表示估计时延,一般为相关运算峰值最大的时刻。

    图 2可知,估算时延τ̂和发射信号的波达方向估计值θ̂的关系可表示为:

    sinθ̂=cτ̂L

    式(1)中,以第一个阵元为参考阵元,若第m个阵元的位置坐标为(xm,ym),根据式(10)式(11)可得到[34]

    τm,k=1c[xmsinθk+ymcosθk]

    通过几何解算可完成发射端的定位,但是在实际环境中,噪声的存在会导致式(9)中相关运算Rx1x2(τ)的峰值不明显,降低了时延估计的精度,影响定位精度。因此如何准确地估计出θ̂τ̂成为水下声阵列定位系统的关键。

    为了在噪声影响严重的条件下准确地估计角度信息θ̂和时延信息τ̂,本文提出一种基于两阶段学习模型的水下声学定位方法。在第一阶段模型中,构建一个基于深度卷积神经网络(deep convolution neural network,DCNN)的降噪模型[35],将接收端信号进行基于时-频关系的二维矩阵化[36],对其进行提取特征输入,通过该模型在纯净不含噪的发射信号和带噪的接收信号之间建立神经网络函数映射,进而得到具有降噪效果的接收信号;在第二阶段模型中,为了提升角度估计的准确性,设计一种改进的基于加权延时求和的波束形成模型,该模型利用梯度下降准则确定各个通道的最优相对延时,根据实际环境对各个通道的权重进行迭代更新[37],得到最优权值,最终输出估计出的角度信息θ̂和时延τ̂,得到位置估计。图 3为本文两阶段学习模型的结构示意图。

    图  3  两阶段学习模型结构
    Figure  3.  Structure of Two-Stage Learning Model

    第一阶段模型的核心是一个降噪学习网络模型。由于降噪的过程不是简单的线性过程,对实部和虚部进行分离再合并的过程会破坏二者的正交性,为了使这种影响最小化,本文对数据预处理进行了策略优化:(1)在预处理阶段对接收信号数据进行重新组合分帧,在短时平稳区间内,截取数据长度间存在帧重叠,尽量保持正交的关联性;(2)在得到分帧数据后进行实部和虚部分离,将分离后的实部和虚部数据进行函数压缩,保证数据特性的统一;(3)利用归一化函数将降噪后的实部和虚部数据合并。对式(1)接收信号数据进行分帧和加窗处理,窗函数window(n)采用海明窗,α0取值0.538 36,则取窗函数为:

    window(n)=α0-(1-α0)α1cos(2πnN-1)

    式中,0nN-1N表示海明窗取值长度。

    对加窗处理后的帧数据进行重新组合分帧,由于接收的信号具有短时平稳的特点,在20~30 ms表现为较为平稳的幅度特征,分帧时每一帧长度为20 ms,帧间设置重叠长度为10 ms,分帧组合策略如图 4所示。

    图  4  数据预处理的分帧策略
    Figure  4.  Framing Strategy of Data Preprocessing

    首先将每帧信号表示为时间-频率的二维化矩阵,进行短时傅里叶变换(short time Fourier transform,STFT),计算如下:

    xrevSTFT(t,f)=n=0N-1x(t,f)window(t,f)e-j2πknN

    式中,f为信号频率。

    然后进行实部和虚部分离:xrevSTFT(t,f)(Rx,Ix),得到接收信号输入特征。同理,对发射信号s(t)也进行同样的操作,得到sSTFT(t,f)(Rs,Is),RxRs分别是xrevSTFT(t,f)sSTFT(t,f)对应的实部数据;IxIs分别是xrevSTFT(t,f)sSTFT(t,f)对应的虚部数据。

    最后采用双曲正切函数分别对实部和虚部进行压缩,得到对应的实部双曲正切压缩(hyperbolic tangent comprssed real component,TR)和虚部双曲正切压缩(hyperbolic tangent comprssed imaginary component,TI)作为网络的输入输出特征,计算如下:

    TR(t,f)=β(1-e-αZR(t,f))1+e-αZR(t,f)
    TI(t,f)=β(1-e-αZI(t,f))1+e-αZI(t,f)

    式中,αβ为限制因子,取值分别为0.5和10;ZR(t,f)表示RxRs经压缩后的结果;ZI(t,f)表示IxIs经压缩后的结果。图 5为降噪网络模型的神经网络结构,R̂Î分别代表经降噪网络输出后的实部、虚部。

    图  5  降噪网络模型的神经网络结构
    Figure  5.  Neural Network Structure of the Noise-Reduction Network Model

    降噪网络是基于DCNN的结构进行适应性调整得到的,常见的DCNN已包含完整的训练结构,为使其在实部虚部特征训练过程中具有去噪作用,本文引入通道关注机制[38]进行特征运算,同时将金字塔去噪结构分配到每个分支,只聚焦在同一个尺度内的特征提取[39],提取全局信息并保留局部细节,为全面去噪进行初始化准备,使用U-Net来分割池化层降采样的特征[40],激活函数选用Sigmoid采用的交叉熵损失函数,计算如下:

    Loss=-1Bb=1B[nlnn̂+(1-n)ln(1-n̂)]

    式中,B为接收阵元个数;n为预期输出;n̂为神经网络输出,即降噪后的R̂zÎz

    将降噪输出的实部和虚部进行归一化还原,得到:

    R̂(n)=-1αlogβ-R̂z(n)β+R̂z(n)
    Î(n)=-1αlogβ-Îz(n)β+Îz(n)

    降噪还原信号为:

    y(t)=R̂(t)+Î(t)×j

    表 1为该阶段模型的具体参数设置。

    表  1  降噪神经网络的参数设置
    Table  1.  Parameter Setting of the Noise-Reduction DNN
    设置项 参数值
    卷积层 4
    全连接层 2
    卷积层滤波器 5×5、3×3、2×2
    步长 2×2
    激活函数 Sigmoid函数
    池化层 3
    输出层 6×161
    下载: 导出CSV 
    | 显示表格

    第二阶段模型的核心是一个改进的基于加权延时求和的波束形成组,将第一阶段降噪处理的接收信号作为第二阶段的输入,进行角度θ̂估计和时延τ̂估计。

    为了更准确地完成期望波达方向上的角度θ̂估计,需要抑制降低其他方向的干扰,本文提出一种使用多个波束形成器对期望信号估计策略,如图 6所示,水下声阵列接收信号分别与P个波束形成器相连,对感兴趣用户的第p条路径信号进行波束形成,即将阵列波束主瓣对准其到达方向,并在其他路径信号和多址干扰信号到达方向上形成期望方向零陷或较低的增益[41]

    图  6  多个波束形成器组成的角度估计结构
    Figure  6.  Angle Estimation Structure Consisting of Multiple Beamformers

    波束形成器组输出向量η1,η2ηM,各个接收支路降噪接收后波束形成器的输出信号可表示为:

    y˜(k)=p=1Lηpyp(k-τp)=p=1LηpVpHx(k-τp)

    式中,ηp表示每个分支波束形成器中期望方向的加权系数;τp为相对时延;Vp表示阵列加权向量,

    可采用MMSE准则[42]迭代计算得到:

    Vp=[α1Me2πifΔτ1α2Me2πifΔτ2αMMe2πifΔτM]T

    式中,αm为第m个接收阵元上的接收信号经频域变化后的相位调节参数,由迭代计算产生。

    Vp代入波束形成器组,为了使期望角度上的输出功率为最大,从而抑制干扰和噪音完成输出θ̂估计,θ̂的判断依据为:

    $$\hat \theta = {\rm{arg}}\; {\rm{max}}[E|{\rm{d}}(t){|^{\rm{2}}}||V_p^Ha(\theta)|{|^{\rm{2}}} + {\sigma ^{\rm{2}}}||{V_p}|{|^{\rm{2}}}]$$

    计算时延τ̂需对式(2)式(4)进行转换,则第m个接收阵元上接收信号可表示为:

    xm(t)=As(t-τm)+nm(t)

    式中,m=1, 2Mτm为两个参考阵元间的相对时延,可采用互相关求解得到:

    Rxm-1xm(τ)=E(xm-1(t)xm-1(t-τ))

    式(25)可进一步表示为:

    Rxm-1xm(τ)=Am-1AmE(s(t-τm-1)s(t-τm-τ))+Am-1E(s(t-τm-1)nm(t-τ)+AmE(s(t-τm-τ)nm(t))+E(nm-1(t)nm(t-τ))

    当限制条件τ=τm-1-τm成立时,Rxm-1xm(τ)取最大值,则两个参考阵元间的时延可表示为:

    $$\Delta {\tau _{m - {\rm{1}}, m}} = {\rm{arg}}\; {\rm{max}}{\mathit{\Psi }_{{x_{m - {\rm{1}}}}{x_m}}}(\tau)$$

    式中,Ψxm-1xm(τ)为互相关系数[26],由快速傅里叶逆变换(inverse fast Fourier transforms,IFFT)展开得到,对于估计时延τ̂,计算如下:

    Ψxm-1xm(τ̂)=IFFTxm-1(f)xm*(f)|xm-1(f)xm*(f)|

    将时延τ̂求解转换为对Ψxm-1xm(τ̂)的最佳加权估计[43]是一个基于梯度下降准则的白化滤波过程,可逐次迭代权值直至得到最优估计,定义梯度为Wt,前后两次权值的加权平均和为Htδ为学习率,表示矩阵逐元素相乘,则有:

    Wt+1=Wt-δWt
    Ht=δHt-1+(1-δ)WtWt
    Wt+1=Wt-δHtWt
    Ht=δt-0H̑+(1-δ)i=1tδt-iH̑

    式中,H̑=W0W0

    Ψxm-1xm(τ)设置为初始权值,按照以下条件更新:

    δHt||δWtWt||

    最后通过式(28)完成对τ̂的有效估计。

    仿真实验采用Bellhop仿真工具[44]模拟水下阵列的发射-接收场景,如图 7所示。发射端为1个声源,接收端为四元均匀线性阵列(uniform linear array,ULA),声场限定在300 m×300 m×160 m的水域内,阵元间距为半波长d,与阵元中心O的距离R=100 m,其中发射端T1T2T3的入射角发生了角度变化,对应角度分别为θ1=60°θ2=100°θ3=130°,发射信号由线性调频波(linear frequency modulation,LFM)和正交相移键控(quadrature phase shift keying,QPSK)组成,噪声采用高斯白噪声,为全频带噪声,信噪比设置为-25~10 dB,多径数设置为5。

    图  7  发射-接收的仿真示意图
    Figure  7.  Transmit-Receive Simulation Schematic

    第一阶段模型的输出为降噪后的实部和虚部,通过短时傅里叶逆变换(inverse STFT,ISTFT)得到去噪接收信号。当信噪比为5 dB时,分别取3个角度的某一帧抽样显示去噪接收信号的时域-频域二维重构结果,如图 8所示。由图 8可以看出,对于接收信号而言,第一阶段降噪学习网络模型具有良好的降噪效果。

    图  8  信号降噪后的时域-频域二维重构
    Figure  8.  Time -Frequency Domain 2-Dimensional Reconstruction of the Received Signal After Noise Reduction

    图 9为去噪前后第二阶段模型角度估计结果。图 9(a)为未经过降噪网络处理直接采用波束形成进行角度估计的结果,图 9(b)为经过第一阶段降噪网络后再使用波束形成估计的角度。由图 9可以看出.角度估计的准确率有了较大幅度的提升。

    图  9  去噪前后第二阶段模型角度估计对比
    Figure  9.  Angle Estimation of the Second-Stage Model Before and After Noise Reduction

    广义互相关-相位变换方法根据两个参考阵元的互相关函数峰值来估计时延值,在本文仿真模型中,四阵元接收阵列两两可组成一对参考阵元。

    根据式(7)、(8)得到参考阵元接收信号x1x2,将其转换至频域X1X2,则互相关(generalized cross correlation,GCC)过程的互相关函数和互功率谱关系为:

    RX1X2(τ)=0πX1(w)X2*(w)e-jwτdw

    由于噪声的影响,RX1X2(τ)的相关峰值会受影响,这降低了时延的估计精度。为了锐化峰值,在频域内对互功率谱进行加权,从而抑制干扰。则式(34)可变为:

    RX1X2(τ)=0πφ12(w)X1(w)X2*(w)e-jwτdw

    式中,φ12(w)表示频域加权函数。

    加权函数φ12(w)利用PHAT(phase transformation)加权,其对噪声的约束效果较好,但当信号能量较小时,处理误差增大。φ12(w)可展开为:

    φPHAT(w)=1/||Gx1x2||=1/||X1(w)X2*(w)||

    经过PHAT加权的GCC互功率谱近似于单位冲击响应,可突出时延的峰值,进而解算出时延。在对多径的处理中,采用判断策略确定主径,将接收阵列首次到达的信号方向(时间最短)设为主径的入射方向,在GCC-PHAT中,两个参考阵元件间,互相关运算后时间最短的即为相对时延。

    克拉美-罗下界(Cramer-Rao lower bound,CRLB)[45]常用于描述计算理论能达到的最佳估计精度。对于真实值θτ,其CRLB根据定义可表示为:

    Pθ=(JθTRθ-1Jθ)-1
    Pτ=(JτTRτ-1Jτ)-1

    式中,JθJτ分别表示θτ的观测函数关于发射端位置的雅克比(Jacobi)矩阵;R为协方差矩阵。式(12)中第m个阵元位置为(xm,ym),则Jθ中各行元素为:

    Jθm=-ymrm2xmrm2

    Jτ中各行元素为:

    Jτm=-ymrm2+y1r12xmrm2-x1r12

    则联合定位误差的CRLB为:

    P=(JTR-1J)-1

    式(37)式(38)代入式(41)可得:

    P=(JθTRθ-1Jθ+JτTRτ-1Jτ)-1=(Pθ-1+Pτ-1)-1
    Pθ-1=1σθ2m=1Mym2rm4-m=1Mxmymrk4-m=1Mxmymrk4m=1Mxm2rm4
    Pτ-1=(σθ2/στ2)Pθ-1-(1/στ2)aaT
    a=-m=1Mymrm2m=1Mxmrm2

    式中,σθ2στ2分别为θτ的噪声方差。

    应用本文方法和其他两种方法进行仿真实验,设置信噪比为-25~10 dB,多径数设置为5时,分别验证了θ1=60°θ2=100°θ3=130° 3个角度下定位结果的均方根误差(root mean square error,RMSE)和CRLB对比,结果如图 10所示。其中方法1表示直接利用传统算法GCC-PHAT进行定位,方法2表示在第一阶段降噪模型处理后联合GCC-PHAT算法进行处理,方法3为本文方法。由图 10可知,本文方法模型在应对信噪比变化的鲁棒性更高。

    图  10  定位结果的RMSE对比
    Figure  10.  RMSE Comparison of Positioning Results

    同时,本文在仿真实验中还进行了收发拉距测试,以验证方法的适用边界,发现在收发距离超过1 km后,本文方法对于噪声的鲁棒性逐渐降低。

    湖上实验采用2021年7月在安徽省安庆市花亭湖采集的数据,发射信号由LFM和QPSK组成,多帧连续的波形组合生成WAV格式文件,通过功放由发射换能器发出,接收端为四元ULA,实验水域丰水期平均水深为90 m,发射端变换两次位置,角度分别为30°和60°。实验具体参数设置如下:LFM带宽为7~13 kHz,QPSK载波频率为10 kHz,采样率为200 kHz,发射端水深为50 m、86 m,收发距离为100 m。

    实验中使用LFM宽带信号对水下信道的信噪比进行测量估算,平均信噪比为10 dB。图 11为4个接收阵元接收到的水声信号。

    图  11  4个阵元的接收信号
    Figure  11.  Received Signals of the 4 Array Elements

    对角度估计值和位置估计值的计算结果进行100次统计,图1213分别为角度估计、定位结果的对比,其中方法1、方法2、方法3与§3.1相同。

    图  12  角度估计结果对比
    Figure  12.  Comparison of Angle Estimate Results
    图  13  位置估计统计
    Figure  13.  Statistics of Positioning Estimates

    图 12可知,在角度估计中,方法1未经过第一阶段降噪模型处理,角度的估计结果与真实值差距较大;方法2经过第一阶段降噪模型处理,明显改善了噪声对估计结果的影响;方法3即本文方法改善了噪声的影响,准确地估算出期望方向的角度,验证了方法的稳定性。由图 13可知,在定位结果中,方法1定位结果的解算收敛性不佳,误差较大;方法2定位结果有所改善;方法3的定位解算结果具有十分良好的收敛性,定位精度更高,更稳定。

    表 2为湖上实验的结果对比,包括30°和60°角度估计的平均值和均方根误差,以及真实值分别为(20,50)m、(56,86)m中定位估计的平均值和均方根误差。由表 2可知,本文提出的两阶段方法在实测数据中也保持了较好的优势。

    表  2  湖上实验结果对比统计
    Table  2.  Result Comparison on the Lake Test
    方法 角度/(°) 定位结果/m
    真实值 估计平均值 均方根误差 真实值 估计平均值 均方根误差
    方法1 30 22.3 0.52 (20, 50) (15.32, 51.67) 0.36
    方法2 27.4 0.31 (21.02, 53.43) 0.25
    方法3 30.8 0.12 (20.87, 50.24) 0.10
    方法1 60 53.6 0.49 (56, 86) (52.30, 84.28) 0.34
    方法2 56.3 0.33 (54.24, 87.51) 0.29
    方法3 60.5 0.14 (56.38, 86.55) 0.15
    下载: 导出CSV 
    | 显示表格

    水下声阵列信号处理是水下定位研究的重要手段之一。引入降噪处理、加权时延估计和基于梯度下降准则的权值迭代等方法,可提升定位系统对于噪音变化影响的鲁棒性和准确率,本文利用以上策略建立了两阶段模型方法,有效地应对了仿真环境中信噪比为-25~10 dB、多径数为5的干扰场景,并在湖上实验也得到了有效的验证。对接收信号分段处理的两阶段模型结构适用于多数阵列信号处理的定位场景,包括室内声学定位和无线电雷达阵列等领域[46]。在未来的工作中,可持续性发掘深度学习模型中数据特征提取的优势[47],在多个场景下的水声传播数据场景中,构建覆盖范围更广的训练模型库,使得预训练结构具有更佳的范用价值[48],从而达到更佳的定位效果。

  • 图  1   水下声阵列的发射-接收结构示意图

    Figure  1.   Transmit-Receive Structure of the Underwater Acoustic Arrays

    图  2   参考阵元间的结构关系

    Figure  2.   Structural Relationship of Reference Array Elements

    图  3   两阶段学习模型结构

    Figure  3.   Structure of Two-Stage Learning Model

    图  4   数据预处理的分帧策略

    Figure  4.   Framing Strategy of Data Preprocessing

    图  5   降噪网络模型的神经网络结构

    Figure  5.   Neural Network Structure of the Noise-Reduction Network Model

    图  6   多个波束形成器组成的角度估计结构

    Figure  6.   Angle Estimation Structure Consisting of Multiple Beamformers

    图  7   发射-接收的仿真示意图

    Figure  7.   Transmit-Receive Simulation Schematic

    图  8   信号降噪后的时域-频域二维重构

    Figure  8.   Time -Frequency Domain 2-Dimensional Reconstruction of the Received Signal After Noise Reduction

    图  9   去噪前后第二阶段模型角度估计对比

    Figure  9.   Angle Estimation of the Second-Stage Model Before and After Noise Reduction

    图  10   定位结果的RMSE对比

    Figure  10.   RMSE Comparison of Positioning Results

    图  11   4个阵元的接收信号

    Figure  11.   Received Signals of the 4 Array Elements

    图  12   角度估计结果对比

    Figure  12.   Comparison of Angle Estimate Results

    图  13   位置估计统计

    Figure  13.   Statistics of Positioning Estimates

    表  1   降噪神经网络的参数设置

    Table  1   Parameter Setting of the Noise-Reduction DNN

    设置项 参数值
    卷积层 4
    全连接层 2
    卷积层滤波器 5×5、3×3、2×2
    步长 2×2
    激活函数 Sigmoid函数
    池化层 3
    输出层 6×161
    下载: 导出CSV

    表  2   湖上实验结果对比统计

    Table  2   Result Comparison on the Lake Test

    方法 角度/(°) 定位结果/m
    真实值 估计平均值 均方根误差 真实值 估计平均值 均方根误差
    方法1 30 22.3 0.52 (20, 50) (15.32, 51.67) 0.36
    方法2 27.4 0.31 (21.02, 53.43) 0.25
    方法3 30.8 0.12 (20.87, 50.24) 0.10
    方法1 60 53.6 0.49 (56, 86) (52.30, 84.28) 0.34
    方法2 56.3 0.33 (54.24, 87.51) 0.29
    方法3 60.5 0.14 (56.38, 86.55) 0.15
    下载: 导出CSV
  • [1] 刘经南, 陈冠旭, 赵建虎, 等. 海洋时空基准网的进展与趋势[J]. 武汉大学学报·信息科学版, 2019, 44(1): 17-37 doi: 10.13203/j.whugis20180340

    Liu Jingnan, Chen Guanxu, Zhao Jianhu, et al. De velopment and Trends of Marine Space-Time Frame Network[J]. Geomatics and Information Science of Wuhan University, 2019, 44(1): 17-37 doi: 10.13203/j.whugis20180340

    [2]

    Tao J W, Chang W X, Cui W. Vector Field Smoothing for DOA Estimation of Coherent Under water Acoustic Signals in Presence of a Reflecting Boundary[J]. IEEE Sensors Journal, 2007, 7(8): 1 152-1 158 doi: 10.1109/JSEN.2007.900963

    [3]

    Yang T C. Deconvolved Conventional Beamforming for a Horizontal Line Array[J]. IEEE Journal of Oceanic Engineering, 2018, 43(1): 160-172 doi: 10.1109/JOE.2017.2680818

    [4]

    Friedlander B. A Sensitivity Analysis of the MUSIC Algorithm[J]. IEEE Transactions on Acoustics, Speech, and Signal Processing, 1990, 38(10): 1 740-1 751 doi: 10.1109/29.60105

    [5]

    Rahamim D, Tabrikian J, Shavit R. Source Localiza tion Using Vector Sensor Array in a Multipath Envi ronment[J]. IEEE Transactions on Signal Processing, 2004, 52(11): 3 096-3 103 doi: 10.1109/TSP.2004.836456

    [6]

    Paulraj A, Reddy V U, Shan T J, et al. Perfor mance Analysis of the Music Algorithm with Spatial Smoothing in the Presence of Coherent Sources [C]// IEEE Military Communications Conference, Monterey, California, USA, 1986

    [7]

    Roy R, Kailath T. ESPRIT-Estimation of Signal Parameters via Rotational Invariance Techniques [J]. IEEE Transactions on Acoustics, Speech, and Signal Processing, 1989, 37(7): 984-995 doi: 10.1109/29.32276

    [8]

    Li M H, Lu Y L. Maximum Likelihood DOA Esti mation in Unknown Colored Noise Fields[J]. IEEE Transactions on Aerospace and Electronic Systems, 2008, 44(3): 1 079-1 090 doi: 10.1109/TAES.2008.4655365

    [9]

    Stoica P, Gershman A B. Maximum-Likelihood DOA Estimation by Data-Supported Grid Search [J]. IEEE Signal Processing Letters, 1999, 6(10): 273-275 doi: 10.1109/97.789608

    [10]

    Iltis R A, Kim S J, Hoang D A. Noncooperative Iterative MMSE Beamforming Algorithms for Ad Hoc Networks[J]. IEEE Transactions on Com munications, 2006, 54(4): 748-759 doi: 10.1109/TCOMM.2006.873095

    [11]

    Choi S, Choi J, Im H J, et al. A Novel Adaptive Beamforming Algorithm for Antenna Array CDMA Systems with Strong Interferers[J]. IEEE Transac tions on Vehicular Technology, 2002, 51(5): 808-816 doi: 10.1109/TVT.2002.801546

    [12]

    Koutrouvelis A I, Hendriks R C, Heusdens R, et al. Relaxed Binaural LCMV Beamforming[J]. ACM Transactions on Audio, Speech, and Lan guage Processing, 2017, 25(1): 137-152 http://ieeexplore.ieee.org/document/7742910

    [13]

    Blandin C, Ozerov A, Vincent E. Multi-source TDOA Estimation in Reverberant Audio Using An gular Spectra and Clustering[J]. Signal Processing, 2012, 92(8): 1 950-1 960 doi: 10.1016/j.sigpro.2011.09.032

    [14]

    Kilfoyle D B, Baggeroer A B. The State of the Art in Underwater Acoustic Telemetry[J]. IEEE Jour nal of Oceanic Engineering, 2000, 25(1): 4-27 doi: 10.1109/48.820733

    [15]

    Odendaal J W, Barnard E, Pistorius C W I. TwoDimensional Superresolution Radar Imaging Using the MUSIC Algorithm[J]. IEEE Transactions on Antennas and Propagation, 1994, 42(10): 1 386- 1 391 doi: 10.1109/8.320744

    [16] 赵建虎, 严俊, 张红梅, 等. 基于海底底质回波特征的多波束声呐图像角度响应影响消除[J]. 武汉大学学报·信息科学版, 2018, 43(8): 1 228-1 233 doi: 10.13203/j.whugis20160148

    Zhao Jianhu, Yan Jun, Zhang Hongmei, et al. Using Backscatter Characteristic of Seabed Sediment to Remove Angular Response Effect in Multibeam Images[J]. Geomatics and Information Science of Wuhan University, 2018, 43(8): 1 228-1 233 doi: 10.13203/j.whugis20160148

    [17]

    Li J F, Ma P H, Zhang X F, et al. Improved DFT Algorithm for 2D DOA Estimation Based on 1D Nested Array Motion[J]. IEEE Communications Letters, 2020, 24(9): 1 953-1 956 doi: 10.1109/LCOMM.2020.2997030

    [18]

    Van Veen B D, Buckley K M. Beamforming: A Versatile Approach to Spatial Filtering[J]. IEEE ASSP Magazine, 1988, 5(2): 4-24 doi: 10.1109/53.665

    [19]

    Pillai S U, Bar-Ness Y, Haber F. A New Ap proach to Array Geometry for Improved Spatial Spectrum Estimation[J]. Proceedings of the IEEE, 1985, 73(10): 1 522-1 524 doi: 10.1109/PROC.1985.13324

    [20]

    Ren A, Li Z, Ding C W, et al. SC-DCNN: High ly-Scalable Deep Convolutional Neural Network Using Stochastic Computing[C]// The 22nd Inter national Conference on Architectural Support for Pro gramming Languages and Operating Systems, Xian, China, 2017

    [21]

    Huang H J, Yang J, Huang H, et al. Deep Learning for Super-Resolution Channel Estimation and DOA Estimation Based Massive MIMO System[J]. IEEE Transactions on Vehicular Technology, 2018, 67(9): 8 549-8 560 doi: 10.1109/TVT.2018.2851783

    [22]

    Zheng W Q, Zou Y X, Ritz C. Spectral Mask Esti mation Using Deep Neural Networks for Inter-Sen sor Data Ratio Model Based Robust DOA Estimation [C]// IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), South Brisbane, QLD, Australia, 2015

    [23]

    Wang D S, Zou Y X, Wang W W. Learning Soft Mask with DNN and DNN-SVM for Multi-Speaker DOA Estimation Using an Acoustic Vector Sensor [J]. Journal of the Franklin Institute, 2018, 355 (4): 1 692-1 709 doi: 10.1016/j.jfranklin.2017.05.002

    [24]

    Lee H Y, Cho J W, Kim M, et al. DNN-Based Feature Enhancement Using DOA-Constrained ICA for Robust Speech Recognition[J]. IEEE Signal Processing Letters, 2016, 23(8): 1 091-1 095 doi: 10.1109/LSP.2016.2583658

    [25]

    Xiang H H, Chen B X, Yang M L, et al. A Novel Phase Enhancement Method for Low-Angle Estima tion Based on Supervised DNN Learning[J]. IEEE Access, 2019, 7: 82 329-82 336 doi: 10.1109/ACCESS.2019.2924156

    [26]

    Cong J Y, Wang X P, Huang M X, et al. Robust DOA Estimation Method for MIMO Radar via Deep Neural Networks[J]. IEEE Sensors Journal, 2021, 21(6): 7 498-7 507 doi: 10.1109/JSEN.2020.3046291

    [27]

    Rizwan T, Cai Y Z, Ahsan M, et al. Neural Net work Approach for 2-Dimension Person Pose Esti mation with Encoded Mask and Keypoint Detection [J]. IEEE Access, 2020, 8: 107 760-107 771 doi: 10.1109/ACCESS.2020.3001473

    [28]

    Selva J. Computation of Spectral and Root MUSIC Through Real Polynomial Rooting[J]. IEEE Transac tions on Signal Processing, 2005, 53(5): 1 923- 1 927 doi: 10.1109/TSP.2005.845489

    [29]

    Abdullah S, Zamani M, Demosthenous A. To wards more Efficient DNN-Based Speech Enhance ment Using Quantized Correlation Mask[J]. IEEE Access, 2021, 9: 24 350-24 362 doi: 10.1109/ACCESS.2021.3056711

    [30]

    Dong Y D, Wang H X, Yao Y D. Channel Estima tion for One-Bit Multiuser Massive MIMO Using Conditional GAN[J]. IEEE Communications Let ters, 2021, 25(3): 854-858 doi: 10.1109/LCOMM.2020.3035326

    [31]

    Larsson E G, Edfors O, Tufvesson F, et al. Mas sive MIMO for Next Generation Wireless Systems [J]. IEEE Communications Magazine, 2014, 52 (2): 186-195 doi: 10.1109/MCOM.2014.6736761

    [32]

    Xia H J, Yang K D, Ma Y L, et al. Noise Reduc tion Method for Acoustic Sensor Arrays in Underwa ter Noise[J]. IEEE Sensors Journal, 2016, 16 (24): 8 972-8 981 doi: 10.1109/JSEN.2016.2618770

    [33]

    Lee R, Kang M S, Kim B H, et al. Sound Source Localization Based on GCC-PHAT with Diffuse ness Mask in Noisy and Reverberant Environments [J]. IEEE Access, 2020, 8: 7 373-7 382 doi: 10.1109/ACCESS.2019.2963768

    [34]

    Khabbazibasmenj A, Vorobyov S A, Hassanien A. Robust Adaptive Beamforming Based on Steering Vector Estimation with as Little as Possible Prior Information[J]. IEEE Transactions on Signal Pro cessing, 2012, 60(6): 2 974-2 987 doi: 10.1109/TSP.2012.2189389

    [35]

    Wei Z, Chen X D. Deep-Learning Schemes for Full-Wave Nonlinear Inverse Scattering Problems [J]. IEEE Transactions on Geoscience and Remote Sensing, 2019, 57(4): 1 849-1 860 doi: 10.1109/TGRS.2018.2869221

    [36]

    Wang Z Q, Zhang X L, Wang D L. Robust Speaker Localization Guided by Deep Learning-Based TimeFrequency Masking[J]. IEEE/ACM Transactions on Audio, Speech, and Language Processing, 2019, 27(1): 178-188 doi: 10.1109/TASLP.2018.2876169

    [37] 刘经南, 曾文宪, 徐培亮. 整体最小二乘估计的研究进展[J]. 武汉大学学报·信息科学版, 2013, 38(5): 505-512 http://ch.whu.edu.cn/article/id/2641

    Liu Jingnan, Zeng Wenxian, Xu Peiliang. Overview of Total Least Squares Methods[J]. Geomatics and Information Science of Wuhan University, 2013, 38 (5): 505-512 http://ch.whu.edu.cn/article/id/2641

    [38]

    Chadwell D, Bürgmann R, Zhang Bo, et al. The Nonsub Sampled Contourlet Transform[J]. Transactions on Image Processing, 2006, 15(10): 3089-3101 doi: 10.1109/TIP.2006.877507

    [39]

    Morozs N, Gorma W, Henson B T, et al. Channel Modeling for Underwater Acoustic Network Simula tion[J]. IEEE Access, 2020, 8: 136 151-136 175 doi: 10.1109/ACCESS.2020.3011620

    [40]

    Hinton G, Deng L, Yu D, et al. Deep Neural Net works for Acoustic Modeling in Speech Recognition: The Shared Views of Four Research Groups[J]. IEEE Signal Processing Magazine, 2012, 29(6): 82-97 doi: 10.1109/MSP.2012.2205597

    [41]

    Zhang Y, Yan D Q. Path Generation Method for Aero-Engine Free-Form Surface Blade in Laser Solid Forming[J]. IEEE Access, 2019, 7: 140 588- 140 595 doi: 10.1109/ACCESS.2019.2942378

    [42]

    Protter M, Yavneh I, Elad M. Closed-Form MMSE Estimation for Signal Denoising Under Sparse Repre sentation Modeling over a Unitary Dictionary[J]. IEEE Transactions on Signal Processing, 2010, 58 (7): 3 471-3 484 doi: 10.1109/TSP.2010.2046596

    [43]

    Sundar H, Sreenivas T V, Seelamantula C S. TDOA-Based Multiple Acoustic Source Localiza tion Without Association Ambiguity[J]. IEEE/ ACM Transactions on Audio, Speech, and Lan guage Processing, 2018, 26(11): 1 976-1 990 doi: 10.1109/TASLP.2018.2851147

    [44]

    Stamatiou K, Casari P, Zorzi M. The Throughput of Underwater Networks: Analysis and Validation Using a Ray Tracing Simulator[J]. IEEE Transactions on Wireless Communications, 2013, 12(3): 1 108-1 117 doi: 10.1109/TWC.2013.012513.120234

    [45]

    Balkan G O, Gezici S. CRLB Based Optimal Noise Enhanced Parameter Estimation Using Quantized Observations[J]. IEEE Signal Processing Letters, 2010, 17(5): 477-480 doi: 10.1109/LSP.2010.2043787

    [46] 刘经南, 高柯夫. 智能时代测绘与位置服务领域的挑战与机遇[J]. 武汉大学学报·信息科学版, 2017, 42(11): 1 506- 1 517 doi: 10.13203/j.whugis20170324

    Liu Jingnan, Gao Kefu. Challenges and Opportuni ties for M-apping and Surveying and Location Based Service in the Age of Intelligence[J]. Geomatics and Information Science of Wuhan University, 2017, 42(11): 1 506-1 517 doi: 10.13203/j.whugis20170324

    [47] 朱建军, 宋迎春, 胡俊, 等. 测绘大数据时代数据处理理论面临的挑战与发展[J]. 武汉大学学报·信息科学版, 2021, 46(7): 1 025-1 031 doi: 10.13203/j.whugis20210232

    Zhu Jianjun, Song Yingchun, Hu Jun, et al. Chal lenges and Development of Data Processing Theory in the Era of Surveying and Mapping Big Data[J]. Geomatics and Information Science of Wuhan Uni versity, 2021, 46(7): 1 025-1 031 doi: 10.13203/j.whugis20210232

    [48] 刘经南, 方媛, 郭迟, 等. 位置大数据的分析处理研究进展[J]. 武汉大学学报·信息科学版, 2014, 39(4): 379-385 doi: 10.13203/j.whugis20140210

    Liu Jingnan, Fang Yuan, Guo Chi, et al. Research Progress in Location Big Data Analysis and Processing [J]. Geomatics and Information Science of Wuhan University, 2014, 39(4): 379-385 doi: 10.13203/j.whugis20140210

  • 期刊类型引用(3)

    1. 孟怡悦,郭迟,刘经南. 采用注意力机制和奖励塑造的深度强化学习视觉目标导航方法. 武汉大学学报(信息科学版). 2024(07): 1100-1108+1119 . 百度学术
    2. 张银胜,崔志强,王兴涛,孙佳琪,胡宇翔,单慧琳. 基于单目深度估计的智能驾驶路径规划方法. 国外电子测量技术. 2023(08): 71-79 . 百度学术
    3. 李新凯,虎晓诚,马萍,张宏立. 基于改进DDPG的无人驾驶避障跟踪控制. 华南理工大学学报(自然科学版). 2023(11): 44-55 . 百度学术

    其他类型引用(22)

图(13)  /  表(2)
计量
  • 文章访问数:  1243
  • HTML全文浏览量:  333
  • PDF下载量:  86
  • 被引次数: 25
出版历程
  • 收稿日期:  2021-08-24
  • 发布日期:  2021-12-04

目录

/

返回文章
返回