CN110108902B - 用于三维非正交超声波阵列测风装置的测量误差校正方法 - Google Patents

Abstract

本发明公开了一种用于三维非正交超声波阵列测风装置的测量误差校正方法，先测量单通道风速构造风速矩阵，再利用工具MATLAB计算总风速的平均值及标准差；对于测量结果来说，多次测量结果的平均值应该与实际风场的真实值相等或近似，因此寻找一个合成矩阵，使得利用此矩阵合成的总风速的平均值与风场值相等，但是标准差比原始标准差小很多，从而达到风速测量误差的校正，且测量结果更加稳定。

Description

用于三维非正交超声波阵列测风装置的测量误差校正方法

技术领域

本发明属于光通信技术领域，更为具体地讲，涉及一种用于三维非正交超声波阵列测风装置的测量误差校正方法。

背景技术

目前许多行业都对风速测量有着一定的需求。风的测量方法有许多种，按基本测量原理可以分为：机械式风速仪、超声波风速仪、激光多普勒风速仪等。其中，常用的机械式风速仪如风杯风速计，由于其存在机械部件，各部件易损耗，且存在机械惯性的影响，导致其存在实时精度低、有测量盲区、需要启动风速等缺点，在高精度测量中实用性过低。而激光多普勒风速仪使用复杂，适用场景少，且造价昂贵。超声波风速仪成型后便是一体的工作设备，无机械活动部件，易于安装维护，具有量程广、盲区小、实时精度高、线性度好等优点，成为越来越广泛的风速测量仪器。

利用超声波信号在气体中的传播速度受风的影响来对风速大小进行测量，是超声波应用技术在气体介质中的一种体现。和常用的机械式风速仪不同，超声波测量最大的优点是整个测风系统不需要依靠机械材质的转动，没有惯性的影响，不需要考虑测量器件的磨损，能准确的测量出被测风场的风速信息。与激光多普勒风速仪相比，超声波的测量方法原理简单，且易于制造，在精度方面也丝毫不差于激光多普勒仪。再结合现代数字信号处理与计算机技术，能够精确的得到风速和风向值，可以在更高的层面揭示风矢层的特性，对其他与风速有关的实验研究有着重大意义。

发明内容

本发明的目的在于克服现有技术的不足，提供一种用于三维非正交超声波阵列测风装置的测量误差校正方法，利用测量的单通道风速测量结果，重新寻找一个合成矩阵，使得利用此矩阵合成的新风速的平均值与风场值相等，但是标准差比原始标准差小很多，测量结果更加稳定。

为实现上述发明目的，本发明一种用于三维非正交超声波阵列测风装置的测量误差校正方法，其特征在于，包括以下步骤：

(1)、测量三维非正交超声波阵列测风装置各通道的风速，再利用各通道的风速构成风速矩阵，记为V_阵；

(2)、将各通道的测量时间、各通道测得的风速以及各通道的总风速构成校正前的测量数据矩阵，并导入至工具MATLAB；

(3)、利用工具MATLAB计算总风速的平均值V_ave及标准差V_std；

(4)、结合三维非正交超声波阵列测风装置，利用工具MATLAB产生两个逼近标准值的角度值θ、

(5)、根据θ、

构造合成矩阵A；

(6)、计算A^TA的逆矩阵(A^TA)^-1；

(7)、根据逆矩阵(A^TA)^-1计算V²；

V²＝(V_x V_y V_z)(V_x V_y V_z)^T＝V_阵(A^TA)^-1(V_阵)^T

其中，V_x、V_y、V_z表示在空间直角坐标系下的轴向速度；

(8)、对V²中的每一个元素开方，得到每个时刻校正后的总风速，再对所有校正后的总风速求取平均值，得到拟合后的总风速平均值V_ave1；

(9)、设置阈值V^*，比较V_ave1与V_ave的差值是否小于阈值V^*，如果小于，则进入步骤(10)，否则，返回至步骤(4)；

(10)、计算总风速平均值V_ave1的标准差V_std1，比较V_std1与V_std的大小，如果V_std1小于V_std的λ倍，λ取值为(0,1)，则合成矩阵A为校正后的标准矩阵，返回校正后每个时刻的总风速，误差校正结束；否则，返回步骤(4)。

本发明的发明目的是这样实现的：

本发明一种用于三维非正交超声波阵列测风装置的测量误差校正方法，先测量单通道风速构造风速矩阵，再利用工具MATLAB计算总风速的平均值及标准差；对于测量结果来说，多次测量结果的平均值应该与实际风场的真实值相等或近似，因此寻找一个合成矩阵，使得利用此矩阵合成的总风速的平均值与风场值相等，但是标准差比原始标准差小很多，从而达到风速测量误差的校正，且测量结果更加稳定。

同时，本发明用于三维非正交超声波阵列测风装置的测量误差校正方法还具有以下有益效果：

(1)、结构安装好后仅进行一次校正，只要结构不改变，则以后都可以直接用找到的同一个矩阵来降低误差；

(2)、本次发明不需要用到多次迭代的算法，减少了算法复杂度，实际运用中更加方便，且不占用过多资源。

附图说明

图1是本发明用于三维非正交超声波阵列测风装置的测量误差校正方法流程图；

图2是三维非正交超声波阵列测风装置一种具体实施方式结构图；

图3是换能器安装芯棒的一种具体实施方式结构图；

图4是左右环形支架的一种具体实施方式结构图；

图5是结构安装误差分析图；

图6是拟合前后对比效果图；

图7是验证前后对比效果图一；

图8是验证前后对比效果图二。

具体实施方式

下面结合附图对本发明的具体实施方式进行描述，以便本领域的技术人员更好地理解本发明。需要特别提醒注意的是，在以下的描述中，当已知功能和设计的详细描述也许会淡化本发明的主要内容时，这些描述在这里将被忽略。

实施例

图1是本发明用于三维非正交超声波阵列测风装置的测量误差校正方法流程图。

在本实施例中，如图1所示，本发明一种用于三维非正交超声波阵列测风装置的测量误差校正方法，包括以下步骤：

S1、测量三维非正交超声波阵列测风装置各通道的风速，再利用各通道的风速构成风速矩阵，记为V_阵；

在本实施例中，如图2所示，三维非正交超声波阵列测风装置从上到下包括换能器安装芯棒1、左右环形支架2、主轴3、电路盒4、支撑架5。换能器安装芯棒是中空管状设计，表面根据距离需求设计对应的螺纹，如图3所示。环形支架用于固定换能器芯棒，根据实施方式在环形支架上开螺纹孔，如图4所示。

其中，如图4所示，左右环形支架上分别设置有4个四个螺纹孔，一个螺纹孔位于支架的圆心处，另外三个螺纹孔以中央螺纹孔为几何中心，构成等边三角形，左右环形支架上的孔呈现对应关系。

利用任意三组换能器芯棒上的换能器在空间直角坐标系上三个方向的风速分量，即可求出三维空间中总的风速和风向。

如图2所示，本实施例共利用空间中的8个超声波换能器，两两一组，共构成4条超声波传播通道。通过测量4个通道上的风速后，可通过风速合成获得空间中的真实风速风向。

利用三个非正交的矢量就能合成空间中的一个矢量的原理，我们设换能器1-1和1-2所在通道在空间坐标系中的单位方向向量为i，换能器1-3和1-4的单位方向向量为j，换能器1-5和1-6的单位方向向量为k。由于仅需要三个通道上的速度便可合成最终风速风向，暂不考虑换能器1-7和1-8；

设总风速为V，通道1测得风速为V₁₂，通道2测得风速为V₃₄，通道3测得风速为V₅₆，V_阵＝(V₁₂V₃₄V₅₆)；

S2、将各通道的测量时间、各通道测得的风速以及各通道的总风速构成校正前的测量数据矩阵，并导入至工具MATLAB；

S3、利用工具MATLAB计算总风速的平均值V_ave及标准差V_std；

S4、结合三维非正交超声波阵列测风装置，利用工具MATLAB产生两个逼近标准值的角度值θ、

其取值范围为：θ＝60°±0.5°、

S5、根据θ、

构造合成矩阵A；

S6、计算A^TA的逆矩阵(A^TA)^-1；

S7、根据逆矩阵(A^TA)^-1计算V²；

在本实施例中，风速在空间直角坐标系的轴向速度分别为：V_x、V_y、V_z，则有如下关系式：

(V₁₂ V₃₄ V₅₆)＝(V_x V_y V_z)(i^T j^T k^T)

令矩阵(i^T j^T k^T)为矩阵A，i、j、k为单位方向向量，彼此均不正交，则轴向风速的笛卡尔坐标为下式：

(V_x V_y V_z)＝(V₁₂ V₃₄ V₅₆)A^-1

那么，风速V的平方可表示为式：

V²＝(V_x V_y V_z)(V_x V_y V_z)^T＝(V₁₂ V₃₄ V₅₆)A^-1[(V₁₂ V₃₄ V₅₆)A^-1]^T＝(V₁₂ V₃₄ V₅₆)A^-1(A^-1)^T(V₁₂ V₃₄ V₅₆)^T

根据矩阵运算规则，上式中A^-1(A^-1)^T＝A^-1(A^T)^-1＝(A^TA)^-1；

特殊的，当芯棒安装没有偏差时，根据各换能器的空间位置，可以得到下式：

此时AA^T如下式所示:

在本实施例中，根据结构设计，芯棒的安装螺纹与环形支架的螺纹长度均为12mm，螺纹高度为1.2mm。安装角度误差简易示意图如图5。对于结构加工误差，精度要求安装GB1804-m公差标准，加工误差在±0.1mm。则安装过程中实际安装位置与标准位置最大偏差也为0.1mm，螺纹长度为12mm。

则

arcsin(0.00833)＝Δθ≈0.5°，角度偏差为±0.5°。

求解三维坐标系下三个通道的方向向量，通道1由换能器1-1和1-2的位置决定，其单位方向向量i如式：

其中，θ为芯棒1-1与y轴的竖直夹角，

为芯棒1-1与x轴的水平夹角。

换能器1-3和1-4的单位方向向量j如式：

换能器1-5和1-6的单位方向向量k如式：

由上述三个通道的单位向量i、j、k可以得到在非标准情况下风速风向合成矩阵A如下式：

在标准情况下：θ₀＝60°，

设通道1、2、3的单位方向向量为i₀、j₀、k₀，

设标准矩阵即无安装误差时的风速风向合成矩阵A₀为：

A₀＝(i₀ ^T j₀ ^T k₀ ^T)

设实际矩阵A＝(i^T j^T k^T)，考虑摄动矩阵δA＝A-A₀。

设摄动矩阵δA的元素为a_ij，则δA可表示为式：

当Δθ、

时，sinΔθ、cosΔθ、

的取值如表1。

	Δθ＝-0.5°	Δθ＝0.5°
			sin	-0.00872	0.00872
cos	0.99996	0.99996

表1

已知标准矩阵A₀没有0特征值，所以A₀可逆，||A₀||_a是从属于向量范数||x||_a的算子范数。

对标准矩阵A₀求逆矩阵：

根据上述求逆过程，可以得到(A₀)^-1为下式：

令(A₀)^-1δA＝B，B∈R^3×3，则矩阵(A₀)^-1·δA的谱范数为：

其中，r(B^TB)为矩阵B^TB的谱半径：

将矩阵B^TB改写为：

B^TB＝((A₀)^-1·δA)^T(A₀)^-1·δA＝(δA)^T((A₀)^T)^-1(A₀)^-1δA＝(δA)^T(A₀(A₀)^T)^-1δA

计算A₀(A₀)^T得下式：

求A₀(A₀)^T的逆矩阵并化解得下式：

根据矩阵理论的知识，在一个对角矩阵的左边乘一个矩阵，就是用对角阵的对角元素分别乘这个矩阵的对应各行；在对角阵右边乘一个矩阵，就是用对角阵的对角元素分别乘这个矩阵的对应各列。

令对角矩阵

则矩阵B^TB又可改写为下式：

其中，*表示暂时可忽略。

设B^TB的各元素为b_ij，则主对角元素b_ii的未知项为可通过一定的方式计算出来。b_ii取值的极大值和极小值如表2所示：

表2

矩阵B^TB的特征值λ_Bi与矩阵主对角线元素b_ii的关系如下式：

可以看出B^TB为实对称矩阵，必可相似对角化，即能找到可逆矩阵Q使得：B^TB＝Qλ_BQ^-1，其中，λ_B为B^TB的相似对角阵，Q为3维空间中关于B^TB每一特征值的特征单位向量。在实数空间内有如下定理：

B^TB₂＝Qλ_BQ^-1 ₂＝λ_B2＝max(λ_Bi)

因为B^TB的特征值λ_Bi的和远小于1，可判断其最大特征值小于1，由此可得到矩阵(A₀)^-1δA₂＜1，根据定理，可得A₀+δA可逆。

我们证明了实际合成矩阵A＝A₀+δA是可逆的，则根据风速合成的矩阵表达式，只要能够找到矩阵A，使得由此矩阵合成的风速，其标准差足够小，测量结果比利用标准矩阵合成的更加稳定。

S8、对V²中的每一个元素开方，得到每个时刻校正后的总风速，再对所有校正后的总风速求取平均值，得到拟合后的总风速平均值V_ave1；

S9、设置阈值V^*，比较V_ave1与V_ave的差值是否小于阈值V^*，如果小于，则进入步骤S10，否则，返回至步骤S4；

S10、计算总风速平均值V_ave1的标准差V_std1，比较V_std1与V_std的大小，如果V_std1小于V_std的λ倍，λ取值为(0,1)，则合成矩阵A为校正后的标准矩阵，返回校正后每个时刻的总风速，误差校正结束；否则，返回步骤S4。

实例

在本实施例中，我们用风洞10m/s风速的测量结果进行拟合。将测量数据导入MATLAB，计算得到原始数据的平均值为9.6752m/s，标准差为0.1245。经过上述算法后，找到θ＝59.8368°、

此时拟合结果的平均值为9.6759m/s，与原始平均值大了0.0007m/s，标准差仅为0.0693，小于原始标准差的0.6倍。此时拟合矩阵

拟合效果如图6所示。

将10m/s第二次测量结果和15m/s的测量结果导入MATLAB，用寻找的矩阵A重新计算其合成风速风向，并计算拟合结果平均值和标准差。验证拟合效果如图7和图8所示。拟合前后数据的对比如表3所示。

表3

从拟合效果图6、7、8和表格数据对比可以看出，利用此拟合矩阵对风速进行合成，从对比图上看，拟合后的风速变化曲线比原始数据曲线波动更小。从拟合前后数据对比可以看出，拟合后的数据平均值与原始数据相比变化不大，但是标准差约降至原始标准差的一半，说明其结果比标准矩阵合成的结果更稳定，对测量结果的拟合效果明显。

由以上的寻找和验证的过程，我们就找到了一个合适的合成矩阵A，利用此矩阵对单通道测量结果合成，得到测量风速。用这种方法降低在安装过程中由于角度和距离偏差带来的安装误差，在当前安装条件下的拟合过程完毕。

尽管上面对本发明说明性的具体实施方式进行了描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。

Claims (1)

1.一种用于三维非正交超声波阵列测风装置的测量误差校正方法，其特征在于，包括以下步骤：

(1)、测量三维非正交超声波阵列测风装置各通道的风速，再利用各通道的风速构成风速矩阵，记为V_阵；

三 维非正交超声波阵列测风装置从上到下包括换能器安装芯棒、左右环形支架、主轴、电路盒、支撑架；

利用空间中的8个超声波换能器，两两一组，共构成4条超声波传播通道，设换能器1-1和1-2所在通道在空间坐标系中的单位方向向量为i，换能器1-3和1-4的单位方向向量为j，换能器1-5和1-6的单位方向向量为k；

然后利用任意三组换能器芯棒上的换能器在空间直角坐标系上三个方向的风速分量，即可求出三维空间中总的风速和风向；设总风速为V，通道1测得风速为V₁₂，通道2测得风速为V₃₄，通道3测得风速为V₅₆，V_阵＝(V₁₂ V₃₄ V₅₆)；

(2)、将各通道的测量时间、各通道测得的风速以及各通道的总风速构成校正前的测量数据矩阵，并导入至工具MATLAB；

(3)、利用工具MATLAB计算总风速的平均值V_ave及标准差V_std；

(4)、结合三维非正交超声波阵列测风装置，利用工具MATLAB产生两个逼近标准值的角度值θ、

θ、

分别表示芯棒与y轴的竖直夹角、芯棒与x轴的水平夹角；

(5)、根据θ、

构造合成矩阵A；

(6)、计算A^TA的逆矩阵(A^TA)^-1；

(7)、根据逆矩阵(A^TA)^-1计算V²；

V²＝(V_x V_y V_z)(V_x V_y V_z)^T＝V_阵(A^TA)^-1(V_阵)^T

其中，V_x、V_y、V_z表示在空间直角坐标系下的轴向速度；

(8)、对V²中的每一个元素开方，得到每个时刻校正后的总风速，再对所有校正后的总风速求取平均值，得到拟合后的总风速平均值V_ave1；

(9)、设置阈值V^*，比较V_ave1与V_ave的差值是否小于阈值V^*，如果小于，则进入步骤(10)，否则，返回至步骤(4)；

(10)、计算总风速平均值V_ave1的标准差V_std1，比较V_std1与V_std的大小，如果V_std1小于V_std的λ倍，λ取值为(0,1)，则合成矩阵A为校正后的标准矩阵，返回校正后每个时刻的总风速，误差校正结束；否则，返回步骤(4)。

Images