风是相对于大气表面的空气运动，在航空领域、航海领域、工业领域、农业领域、环境监控领域，风参数的测量都发挥着重要作用^[1]。常用的测风方法有机械式、热式、超声波式等。机械式测风仪以风杯式、螺旋桨式为主^[2]，主要基于机械旋转部件来感知风参数，但其机械磨损损耗大，影响测量精度和使用寿命。热式风速仪^[3-4]不存在机械磨损，适用于高分辨率以及低风速测量场合，但是随着风速的增加，灵敏度下降，而且其中的热量元件较脆弱容易损坏^[5]。超声波测量技术与机械式、热式测量相比，无机械磨损，测量速度快，维护成本低，但测量难点在于对超声波飞行时间(time of flight，TOF)的检测精度会直接影响风参数信息的测量精度^[6-7]。

为了发掘多维空域信息，提高风参数的测量精度，基于阵列信号处理思想，文献[8]设计了一发多收的弧形超声波传感器阵列结构。在此基础上，文献[9-10]分别采用MVDR(minimum variance directionless response)算法和MUSIC(multiple signal classification)算法实现了二维风矢量的测量。但目前所提出的阵列结构大多只存在一个发射阵元，在低信噪比或环境因素较差时，风参数估计精度会有所下降。而且空间谱的构建及谱峰搜索占用了MUSIC算法的绝大多数时间，是影响波达方向估计实时性的主要因素^[11]。文献[12]通过计算阵元接收数据之间的延时相关函数得到新的阵列输出矩阵，利用变尺度混沌优化算法进行谱峰搜索，但是以上算法的实现也非常复杂。文献[13]提出基于传播算子的空间平滑技术避免特征值分解，但仍需要全局谱峰搜索。本文通过二级谱峰搜索的方式提高MUSIC算法的实时性。

由于超声波具有指向性强、遇到障碍物会被反射^[14]的特点，文献[15]提出一种基于反射的超声波风向风速测量装置，通过时差法原理进行风参数的测量；文献[16]证明了纵波在0~45°的入射角度范围内的第1个反射面可以发生全反射。而且超声波在波束角内遇到障碍物，波束角内的超声波都会发生全反射^[17]。因此，本文提出基于反射的多发多收超声波传感器阵列测风方法，形成虚拟阵元，增大阵列孔径，通过CRB理论分析和仿真实验验证了本文所提方法的有效性以及测量精度的优越性。

1 信号模型 1.1 构建阵列结构

阵列结构如图 1所示，其中M₁~M₄为发射阵元，N₁~N₄为接收阵元，O为反射物，其中反射物的尺寸与各接收阵元和发射阵元相同，其端面与竖直方向逆时针夹角为45°。图中阵元以及反射物只代表其位置，并不代表其几何形状。由于无线信道中信号的传输情况是极其复杂的，为了得到一个比较有用的参数化模型，必需建立关于接收天线阵的假设，即阵元由位于空间已知坐标处的无源阵元按一定的形式排列而成，其接收特性仅与其位置有关而与其尺寸无关(认为其是一个点)，并且阵元都是全向阵元，增益均相等，相互之间的互耦忽略不计，这是阵列信号处理理论的基础^[18]。文献[10, 19]在进行风信号检测时也忽略了阵元结构对风参数估计的影响。因此，本文可以近似认为接收阵元、发射阵元以及反射物为一个点，接收特性只与其位置有关，与尺寸无关。

发射阵元发射的超声波波束经过反射物后，波束方向发生改变，使发射波束被接收传感器阵列所接收，实现多发多收的超声波传感器阵列结构。所有阵元分布在半径为R的圆周上。发射信号M_i发生全反射的位置恰好为N_i。发射阵元位置关于水平线对称，则通过全反射得到4个接收阵元位置关于竖直线对称。由于发射阵元同时发射超声波信号，在无风时刻，接收阵列同时接收信号，不存在时间延迟；在有风时刻，风信号的存在会影响发射信号到达反射物的时间，以及反射信号到达接收阵元的时间，产生时间延迟。根据超声波信号在传播过程中产生的时间延迟构建发射阵列流型矩阵以及接收阵列流型矩阵，通过DOA估计算法进行风参数的测量。虽然超声波在传播过程中存在衰减和遮挡的问题，但对于时间延迟的计算影响不大。因此不需要考虑这类问题。

为了使发射信号M_i被所有接收阵元所接收，以反射物与N_i为水平线其接收阵元N_i的±45°范围内都应存在超声波信号。根据选取的接收阵元个数N以及发射信号角度范围，计算得到接收阵列中相邻两阵元与反射点连线的夹角α=45°/(N-1)。α随着接收阵元个数的增加而减小，考虑到硬件成本和阵列孔径问题，选取四阵元的接收阵列，此时α=15°，M_i发射的信号可以被所有接收阵元N_i所接收。

1.2 阵列流型推导

图 1中V为待测风速的大小，V与竖直顺时针方向的夹角θ为待测风向角，定义风速与水平逆时针方向的夹角θ′，则θ+θ′=π/2。其发射阵列的阵列流型和接收阵列的阵列流型都包含风速、风向信息。由图 1中几何关系，根据矢量分解定理，可以得到风速在反射物与发射传感器M_i连线上的分量：

$ V_{M_{1}}=V \cos \left(\theta^{\prime}-\frac{3}{2} \alpha\right), V_{M_{2}}=V \cos \left(\theta^{\prime}-\frac{1}{2} \alpha\right) $

$ V_{M_{3}}=V \cos \left(\theta^{\prime}+\frac{1}{2} \alpha\right), V_{M_{4}}=V \cos \left(\theta^{\prime}+\frac{3}{2} \alpha\right) $

同理将V_Mi中θ′用θ代替，可以得到风速在反射物与接收传感器N_i连线上的分量V_Ni(i=1, 2, 3, 4)。通过计算得到各发射信号到达反射点的时间为${t_{{M_i}}} = \frac{R}{{c - {V_{{M_i}}}}}$，反射信号到达各接收阵元的时间为${t_{{N_i}}} = \frac{R}{{c + {V_{{N_i}}}}}$，(R为弧形阵列的半径，c为声速)。分别以传感器M₁与传感器N₁为参考阵元，则各发射信号到达反射物的时间相对于到达参考阵元的时间延迟τ_Mi=t_Mi-t_M1，反射信号到达各接收阵元相对于参考阵元的时间延迟τ_Ni=t_Ni-t_N1(i=1, 2, 3, 4)。

根据τ_Mi得到发射阵列的阵列流型为：

$ \boldsymbol{a}(\theta, V)=\left[\mathrm{e}^{-\mathrm{j} 2 {\rm{ \mathsf{ π} }} f \tau_{M_{1}}}\ \mathrm{e}^{-\mathrm{j} 2 {\rm{ \mathsf{ π} }} f \tau_{M_{2}}}\ \mathrm{e}^{-\mathrm{j} 2 {\rm{ \mathsf{ π} }} f \tau_{M_{3}}}\ \mathrm{e}^{-\mathrm{j} 2 {\rm{ \mathsf{ π} }} f \tau_{M_{4}}}\right]^{\mathrm{T}} $

根据τ_Ni得到接收阵列的阵列流型为：

$ \boldsymbol{b}(\theta, V)=\left[\mathrm{e}^{-\mathrm{j} 2 {\rm{ \mathsf{ π} }} f \tau_{N_{1}}}\ \mathrm{e}^{-\mathrm{j} 2 {\rm{ \mathsf{ π} }} f \tau_{N_{2}}}\ \mathrm{e}^{-\mathrm{j} 2 {\rm{ \mathsf{ π} }} f \tau_{N_{3}}}\ \mathrm{e}^{-\mathrm{j} 2 {\rm{ \mathsf{ π} }} f \tau_{N_{4}}}\right]^{\mathrm{T}} $

式中f为超声波信号频率。

1.3 接收信号的匹配滤波处理

发射信号S =[s₁, s₂, s₃, s₄]^T，其中s₁、s₂、s₃、s₄之间相互正交，即相邻两发射信号之间相位差均为90°。目标回波数学模型可表示为：

$ \boldsymbol{x}=\boldsymbol{\eta} \boldsymbol{a}(\theta, V)^{\mathrm{T}} \boldsymbol{S} $

(1)

式中：η ⊂C表示反射信号的复幅度；a (θ, V)为发射信号的导向矢量。接收阵列信号可表示为：

$ \boldsymbol{y}=\boldsymbol{b}(\theta, V) \boldsymbol{\eta} \boldsymbol{a}(\theta, V)^{\mathrm{T}} \boldsymbol{S}+\boldsymbol{E} $

(2)

式中：b (θ, V)为接收信号的导向矢量；y ∈ C^N×L是接收天线的输出; E ∈ C^N×L为零均值，空间白和时间白的复高斯噪声。

对于接收阵元为N的系统，若同时发射M个线性独立正交波形，在接收端通过匹配滤波将来自不同发射天线的信号分离开，因此可以形成MN个相互独立的通道，相当于合成了虚拟阵元，增大了阵列孔径。利用发射信号与阵列接收数据进行匹配滤波，第i个快拍数下得到的匹配滤波输出为：

$ \boldsymbol{X}_{i}=\frac{1}{L} \sum\limits_{n=1}^{L} \boldsymbol{y}_{i}(n) \boldsymbol{S}_{i}(n)^{\mathrm{H}}=\eta_{i} \boldsymbol{b}(\theta, V) \boldsymbol{a}(\theta, V)^{\mathrm{T}}+\boldsymbol{R}_{\mathrm{es}} $

(3)

式中：L为快拍数；η_i为反射系数; R_es= E_i(n)S_i(n)^H; X_i∈ C^N×M，把X_i的列依次排成一列得到：

$ \operatorname{vec}\left(\boldsymbol{X}_{i}\right)=\eta_{i} \boldsymbol{w}+\boldsymbol{v} $

(4)

式中：vec(·)表示列堆栈处理；v ∈ C^MN×1是均值为零; 方差为$\sigma _v^2 = \frac{1}{{\sqrt L }}{\sigma _e}$的复高斯噪声；w =[a (θ, V)⊗ b (θ, V)]为发射阵列与接收阵列的联合导向矢量，w ∈ C^MN×1，⊗表示Kronecker积。当发射信号快拍数为L时，通过匹配滤波器之后的快拍数据表示为：

$ \boldsymbol{Y}=\left[\begin{array}{llll} \operatorname{vec}\left(\boldsymbol{X}_{1}\right) & \operatorname{vec}\left(\boldsymbol{X}_{2}\right) & \cdots & \operatorname{vec}\left(\boldsymbol{X}_{L}\right) \end{array}\right] $

(5)

2 测风算法及性能分析 2.1 二级空间MUSIC算法

多输入多输出系统根据式(5)得到经过匹配滤波处理的阵列输出快拍数据Y。将Y的数据协方差矩阵进行特征值分解，得到信号子空间和噪声子空间，利用2个子空间的正交性估计风速风向参数。由于实际接收数据矩阵Y是有限长的，阵列接收数据的协方差矩阵采用最大似然估计，即：

$ \boldsymbol{R}=\frac{1}{L} \sum\limits_{i=1}^{L} \boldsymbol{Y} \boldsymbol{Y}^{\mathrm{H}} $

(6)

对 R 进行特征值分解得到噪声子空间V₂，进而得到风参数估计的MUSIC谱估计公式为：

$ P(\theta, V)=\frac{1}{\boldsymbol{w}^{\mathrm{H}} \boldsymbol{V}_{2} \boldsymbol{V}_{2}^{\mathrm{H}} \boldsymbol{w}} $

(7)

在理想条件下，发射阵列与接收阵列的联合导向矢量w与噪声子空间V₂正交，P(θ, V)趋于无穷。由于在实际情况下，w与V₂不能完全正交，可以通过寻找空间谱函数P(θ, V)的最大值对应的参数值得到待估计的风速、风向值。

在进行二级谱峰搜索时，先在待估计参数范围内进行一级搜索，即采用大于其估计精度的搜索步长找到空间谱极大值出现的大致位置，粗略估计出风速值θ和风向值V；然后确定θ的邻域[θ-σ, θ+σ]为二级风向角搜索范围，V的邻域[V-ε, V+ε]为二级风速搜索范围，其中，σ和ε均为正数。最后，采用符合其估计精度的搜索步长在二级搜索范围内精确估计出风速风向值。

2.2 算法计算量分析

若发射阵元数为M, 接收阵元数为N，每计算一次空间谱函数时，各步骤需要的计算量如下：

1) a(θ, V)和b(θ, V)的计算：

乘除法次数：9(N-1)+9(M-1)=9(N+M-2)；

加减法次数：3(N-1)+3(M-1)=3(N+M-2)；

指数运算次数：(N-1)+(M-1)=(N+M-2)。

2) Kronecker积的计算：

乘除法次数：NM。

3) 空间谱函数的计算：

乘除法次数：

$ 2 N M(N M-1)+N M+1=M N(2 N M-1)+1 $

加减法次数：

$ \begin{aligned} &(N M-1)^{2}+(N M-2) N M+(N M-1)= \\ N M(2 &N M-3) \end{aligned} $

所以计算一次空间谱函数需要进行的计算量如下：

乘除法次数：D=9(M+N-2)+2(MN)²+1；

加减法次数：A=3(M+N-2)+MN(2MN-3)；

指数运算次数：N+M-2。

风向估计范围为0°~360°，步长为1°，风速估计范围为0~60 m/s，步长为0.1 m/s。在进行计算量分析时，选取一级搜索时风速步长为0.2 m/s，风向步长为2°，二级搜索时σ=10°，ε=1 m/s。MUSIC算法空间谱构建计算量分析如表 1所示。

从表 1可以看出，二级空间谱的构建极大减少了算法计算量，并且随着阵元数的增加、一级搜索步长的增加以及二级搜索参数的减小，该方法对系统实时性的提高效果越明显。

2.3 估计方差与克拉美罗界

在大快拍数的情况下，MUSIC算法的估计误差是均值为零的渐进联合高斯分布, 文献[10]推导了超声波传感器阵列结合MUSIC算法风参数估计方差简化公式，并且证明了估计方差只与阵列流型，噪声功率以及快拍数有关。因此，通过推导得出当风向一定时，风速估计方差为：

$ \operatorname{var}_{V}(\theta, V)=\frac{\sigma^{2} \operatorname{Re}\left\{\boldsymbol{H}_{V} \odot\left(\boldsymbol{w}^{*} \boldsymbol{T} \boldsymbol{w}\right)^{\mathrm{T}}\right\}}{2 L \boldsymbol{H}_{V}^{2}} $

(8)

$ \operatorname{var}_{\theta}(\theta, V)=\frac{\sigma^{2} \operatorname{Re}\left\{\boldsymbol{H}_{\theta} \odot\left(\boldsymbol{w}^{*} \boldsymbol{T} \boldsymbol{w}\right)^{\mathrm{T}}\right\}}{2 L \boldsymbol{H}_{\theta}^{2}} $

(9)

式中：(·)^*表示矩阵的共轭转置运算；⊙为Hadamard积；σ²为噪声功率；L为快拍数；w^*Tw、H_V和H_θ为：

$ \boldsymbol{w}^{*} \boldsymbol{T} \boldsymbol{w}=\boldsymbol{R}_{\mathrm{S}}^{-1}+\sigma^{2} \boldsymbol{R}_{\mathrm{S}}^{-1}\left(\boldsymbol{w}^{*} \boldsymbol{w}\right)^{-1} \boldsymbol{R}_{\mathrm{S}}^{-1} $

(10)

$ \boldsymbol{H}_{V}=\boldsymbol{D}_{V}^{*}\left(\boldsymbol{I}-\boldsymbol{w}\left(\boldsymbol{w}^{*} \boldsymbol{w}\right)^{-1} \boldsymbol{w}^{*}\right) \boldsymbol{D}_{V} $

(11)

$ \boldsymbol{H}_{\theta}=\boldsymbol{D}_{\theta}^{*}\left(\boldsymbol{I}-\boldsymbol{w}\left(\boldsymbol{w}^{*} \boldsymbol{w}\right)^{-1} \boldsymbol{w}^{*}\right) \boldsymbol{D}_{\theta} $

(12)

式中：R_S为信号协方差矩；${\mathit{\boldsymbol{D}}_{\rm{V}}} = \frac{{\partial \mathit{\boldsymbol{w}}}}{{\partial V}}$和${\mathit{\boldsymbol{D}}_{\rm{ \mathsf{ θ} }}} = \frac{{\partial \mathit{\boldsymbol{w}}}}{{\partial \theta }}$分别为联合导向矢量w对风速和风向的偏导数。

对于参数估计问题，方差是针对某种特定的估计量而言的，而克拉美罗界与具体估计方式无关，它反映的是已有信息所能估计参数的最好效果，为任何无偏估计量的方差确定一个下限，估计方差越接近克拉美罗界，估计效果越好。因此，可以用克拉美罗界作为估计性能好坏的标准。

根据文献[20]给出的克拉美罗界定义式，得到风速的克拉美罗界公式以及风向的克拉美罗界公式分别为：

$ \operatorname{CRB}_{V}(\theta, V)=\frac{\sigma^{2}}{2}\left\{\sum\limits_{i=1}^{L} \operatorname{Re}\left[\boldsymbol{S}^{\mathrm{H}}(i) \boldsymbol{H}_{V} \boldsymbol{S}(i)\right]\right\}^{-1} $

(13)

$ \operatorname{CRB}_{\theta}(\theta, V)=\frac{\sigma^{2}}{2}\left\{\sum\limits_{i=1}^{L} \operatorname{Re}\left[\boldsymbol{S}^{\mathrm{H}}(i) \boldsymbol{H}_{\theta} \boldsymbol{S}(i)\right]\right\}^{-1} $

(14)

式中信号源S(i) = diag{s₁(i), s₂(i), s₃(i), s₄(i)}。

3 仿真结果分析

实验条件：超声波传感器阵列结构如图 1所示。选取发射与接收阵列的弧形半径0.1 m，发射阵列中各发射阵元发射的相互正交的超声波频率40 kHz，声速340 m/s；阵元噪声为零均值，复高斯随机噪声；反射系数在采样时刻恒定不变，不同快拍数下独立变化；风向角扫描范围0°~360°，风速扫描范围0~60 m/s，快拍数为5 000。

3.1 可行性实验

为了验证本文所提方法的可行性，设计仿真实验。选取一级搜索风速步长为0.2 m/s，风向步长为2°，二级搜索参数σ=10°，ε=1 m/s。在信噪比为-10 dB时随机选取风速风向参数(45°, 5.4 m/s)，得到如图 2所示的多发多收二级谱峰图和如图 3所示的一发多收全局谱峰图。

由仿真得到的谱峰图可知，对于随机选取的风参数，在低信噪比情况下，多发多收阵列可以无差估计，而一发多收的阵列无法进行风参数估计。

3.2 估计方差与克拉美罗界实验

在信噪比为10 dB，快拍数为5 000的条件下，通过式(8)~(14)计算出测量范围内每个风速和风向值所对应的估计方差与克拉美罗界，并计算二者的差值，如图 4和图 5所示。

由仿真结果可知，所有的风向角估计方差与CRB之差均小于0.04；风速估计方差与CRB之差均小于1×10^-3。风速和风向的估计方差与克拉美罗界较为接近，由此可知，本文所提方法对测量范围内所有风速、风向值的估计具有较小的偏差。

3.3 不同信噪比下风参数估计均方根误差实验

采用蒙特卡洛仿真实验来评估本文所提方法的风速、风向估计性能。均方根误差公式为：

$ \operatorname{RMSE}=\sqrt{\frac{1}{n} \sum\limits_{i=1}^{n}\left(x_{c i}-\bar{x}_{c}\right)^{2}} $

式中：蒙特卡洛实验次数n=100；x_ci为风参数给定值；x_c为第i次仿真结果。应用本文阵列结构在信噪比-40~5 dB内，选取步长为5 dB；并应用目前已经提出的一发四收的弧形超声波传感器阵列结构，弧形半径0.1 m，相邻两接收阵元与发射阵元连线的夹角α=15°，在信噪比-5~15 dB内，选取步长为2 dB。分别随机输入以下3组参数：(45°, -10.5 m/s)、(123°, 23.4 m/s)、(245°, 48.9 m/s)，采用MUSIC算法进行风速、风向估计。估计均方根误差随信噪比变化情况如图 6所示。

由仿真结果可知，随着信噪比的增大，风速和风向角的均方根误差都逐渐减小。本文所提方法当信噪比分别增大到-25 dB和-20 dB时，风向角和风速均方根误差分别减小为零；而一发四收的弧形阵列当信噪比增大到11 dB时，风向角和风速均方根误差才减小为零。因此，在低信噪比情况下本文所提方法有效提高了风参数的估计精度，获得了良好的估计性能。虽然风速、风向参数是随机选取的，但从3.2节估计方差仿真实验中可以看出，在同一信噪比条件下，对整个测量范围内参数的估计误差相差不大。因此，随机选取的风速、风向参数其估计误差可以代表整个测量范围内参数估计误差的趋势。

3.4 不同信噪比下风参数估计成功率实验

实验条件与3.3节相同，风速估计误差小于其步长0.1 m/s时，认为风速估计成功；风向估计误差小于其步长1°时，认为风向估计成功。估计成功率随信噪比变化曲线如图 7所示，图中虚线为采用二级空间MUSIC算法对风参数的估计结果；实线为直接采用符合参数精度的MUSIC算法得到的估计结果。

由仿真结果可以看出，随着信噪比的增加，风参数估计成功率增大。虽然在低信噪比下，二级空间MUSIC算法对参数估计成功率略低于MUSIC算法，当信噪比分别增大到-25 dB和-20 dB时，风向和风速估计成功率分别达到100%。说明二级空间谱的构建方式不仅可以减小算法计算量，而且具有较高的估计精度。通过以上实验和理论分析，验证了本文所提方法具有较好的风参数测量效果。

4 结论

1) 本文针对风速风向测量精度低的问题，基于阵列信号处理理论，提出多输入多输出的超声波传感器阵列结构的风速风向测量方法，提高了低信噪比下的风参数的测量精度。

2) 本文基于MUSIC算法，提出二级空间谱搜索方法，来估计风速风向值，避免了全局谱峰搜索，减小了算法计算量，提高了系统实时性。

3) 为了验证本文所提阵列结构和算法的风速风向测量性能，进行了风参数估计方差和克拉美罗界的理论分析和计算公式的推导，并通过仿真对比试验，验证了所提方法在测量精度以及系统实时性方面的优势。