CN108663536A - 一种测量三维风速和超声温度仪器的风速数据修正方法 - Google Patents

Abstract

本发明提供一种测量三维风速和超声温度仪器的风速数据修正方法，所述方法包括：S1、获取观测系统中测量三维风速和超声温度的仪器测量的第一风速数据和第一温度数据，以及获取观测系统中另一温度传感器测量的第二温度数据；S2、将第一温度数据和第二温度数据进行比较，判断所述仪器是否处于数据异常状态；S3、在所述仪器发生异常时，根据所述第一风速数据、和在数据异常状态下所述仪器的性能参数，获取真实风速数据。使用本发明的方法恢复得到真实的观测数据，经济、高效、快速、精度高。

Description

一种测量三维风速和超声温度仪器的风速数据修正方法

技术领域

本发明涉及一种测量三维风速和超声温度仪器的风速数据修正方法。

背景技术

三维超声风速测量仪，能准确测出大气中的立体风速风向，进而能真实反映气流的流动状况，及时捕捉空气中的湍流、瞬间的阵风。

如，现有技术中提供一种可以测量三维风速和超声温度的仪器，该仪器对超声风速和超声温度测量的基本原理是测量超声脉冲在一对相互分开但两者间距离为常值的传感器/探头之间的传输时间差。超声波在两个相对方向上的传输时间均依赖于声速和测量路径中的空气速度，可以采用下述的公式(1)和公式(2)测量传感器/探头两端往返的时间，如图1所示。

上述公式中，d₃是路径长度,t_u3和t_d3分别是从传感器/探头两端分别往返的时间。公式(1)为风速与超声波方向一致时，公式(2)为风速与超声波速度方向相反时。

为保证观测结果的精确性，在一次测量过程中需要让超声波在一对超声传感器/探头之间来回往返一次；并同时由3对这样的传感器/探头同时测量。其中，公式(1)和(2)表示的是第3对传感器/探头的测量原理。第1对和第2对的传感器/探头的测量原理与第3对一致，但每对传感器/探头之间的距离可能有所变化。

由图1可知，直接测量的风速方向分别与一对探头的设置方向一致；而这三对传感器/探头的设置方向在空间并不是正交的。

为便于展示3维风场，如图6所示，以仪器上固定的一处为原点，建立仪器坐标系，该仪器坐标系为右手正交立体坐标系；其X方向与该仪器的第一对探头的设置方向一致。

因此，基于每对探头两端的距离以及三对探头相对于该仪器坐标系的方位角和顶角存在唯一对应的一组坐标转换矩阵A，用于将基于这三对探头直接获取的风速数据转化为该仪器坐标系下的3维风速数据。作为这类三维超声风速测量仪(即测量三维风速和超声温度的仪器，简称测量仪器)的关键参数，这组坐标转换矩阵作为关键的出厂参数在该仪器出厂或每次校准标定时会被精确测量并固化在其专用的测量软件中。

当该仪器由于各种意外因素导致某一对或多对探头两端的距离发生改变，或至少一个或一对探头的方位角或顶角发生改变。仪器使用人员没有及时发现仪器异常；或发现了异常，但没有机会及时返厂进行校准标定。则利用原有的出厂参数获得的观测数据将较大地偏离真实数据，导致测量的数据不可用。

发明内容

针对现有技术中的问题，本发明提供一种测量三维风速和超声温度仪器的风速数据修正方法。

第一方面，本发明提供一种测量三维风速和超声温度仪器的风速数据修正方法，包括：

S1、获取观测系统中测量三维风速和超声温度的仪器测量的第一风速数据和第一温度数据，以及获取观测系统中另一温度传感器测量的第二温度数据；

S2、将第一温度数据和第二温度数据进行比较，判断所述仪器是否处于数据异常状态；

S3、在所述仪器发生异常时，根据所述第一风速数据、和在数据异常状态下所述仪器的性能参数，获取真实风速数据。

可选地，所述步骤S3包括：

所述第一风速数据为所述仪器处于数据异常状态时获取的错误的风速数据u_x,u_y.u_z；

S31、根据所述第一风速数据，获取原始测量的用于计算风速的第一测量数据u₁,u₂,u₃；

S32、根据数据异常状态下所述仪器的性能参数，获取所述仪器对应的数据异常状态的形状系数矩阵A_T；

S33、根据所述第一测量数据和所述形状系数矩阵A_T,获得修复后的作为真实风速数据的第二风速数据u_cx,u_cy,u_cz。

可选地，所述步骤S3包括：

所述第一风速数据为所述仪器处于数据异常状态时获取的错误的风速数据u_x，u_y，u_z；

S031、根据所述第一风速数据，获取原始测量的用于计算风速的第二测量数据u_{t1_1},u_{t2_1},u_{t3_1}；

S032、根据数据异常状态下所述仪器的性能参数，获取所述仪器对应的数据异常状态的形状系数矩阵A_T；

S033、根据修正后的第一测量数据和所述形状系数矩阵A_T,获得修复后的作为真实风速数据的第二风速数据u_cx,u_cy,u_cz。

可选地，所述性能参数包括：处于数据异常状态下所述仪器中每一对探头之间的距离、相邻探头之间的夹角。

可选地，所述步骤S31包括：

所述第一风速数据为根据下述公式获取的；

其中，A为嵌入在所述仪器中的该仪器处于正常状态时使用的形状系数矩阵；

所述u_x，u_y，u_z为第一风速数据，所述u₁,u₂,u₃为所述仪器的每对探头的测量数据；

相应地，根据公式：获取真实测量数据u₁,u₂,u₃；

所述步骤S33包括：

根据公式获取第二风速数据u_cx,u_cy,u_cz；

其中，A_T为所述仪器对应的数据异常状态的形状系数矩阵。

可选地，所述步骤S031包括：

所述第一风速数据为根据下述公式获取的；

其中，A为嵌入在所述仪器中的该仪器处于正常状态时使用的形状系数矩阵；

所述u_x,u_y，u_z为第一风速数据，所述u_{t1_1},u_{t2_1},u_{t3_1}为所述仪器的每对探头的测量数据通过阴影校正后得到的真实测量数据；

相应地，根据公式获取真实测量数据u_{t1_1},u_{t2_1},u_{t3_1}；

步骤S033包括：

根据公式获取第二风速数据u_cx,u_cy,u_cz；

其中，A_T为所述仪器对应的数据异常状态的形状系数矩阵。

可选地，

其中，θ₁为第一条声波路径相对于Z轴的顶点角，为第一条声波路径相对于Z轴的方位角；

θ₂为第二条声波路径相对于Z轴的顶点角，为第二条声波路径相对于Z轴的方位角；

θ₃为第三条声波路径相对于Z轴的顶点角，为第三条声波路径相对于Z轴的方位角；

X、Y、Z为预先定义的三维正交仪器坐标系；

根据获取形状系数矩阵A的方式，获取处于数据异常状态的所述仪器的各个声波路径相对于Z轴的顶点角和方位角，根据所述获取形状系数矩阵A的方式和获取的处于数据异常状态的所述仪器的顶点角和方位角，获取数据异常状态的形状系数矩阵A_T。

第二方面，本发明还提供一种测量三维风速的超声温度仪器的风速数据修正方法，包括：

M1、若确定测量三维风速的超声温度仪器的风速发生异常时，则根据所述仪器当前状态下获取的第一风速数据和异常状态下该仪器的性能参数，获取真实风速数据。

可选地，所述M1包括：

所述第一风速数据为所述仪器处于数据异常状态时获取的错误的风速数据u_x，u_y，u_z；

M01、根据所述第一风速数据，获取原始测量的用于计算风速的第一测量数据u₁,u₂,u₃；

M02、根据数据异常状态下所述仪器的性能参数，获取所述仪器对应的数据异常状态的形状系数矩阵A_T；

M03、根据所述第一测量数据和所述形状系数矩阵A_T,获得修复后的作为真实风速数据的第二风速数据u_cx,u_cy,u_cz。

或者，

M11、根据所述第一风速数据，获取原始测量的用于计算风速的第二测量数据u_{t1_1},u_{t2_1},u_{t3_1}；

M12、根据数据异常状态下所述仪器的性能参数，获取所述仪器对应的数据异常状态的形状系数矩阵A_T；

M13、根据修正后的第一测量数据和所述形状系数矩阵A_T,获得修复后的作为真实风速数据的第二风速数据u_cx,u_cy,u_cz。

可选地，所述性能参数包括：处于数据异常状态下所述仪器中每一对探头之间的距离、相邻探头之间的夹角。

本发明具有的有益效果：

本发明提供的数据修正方法适用于该仪器设置的超声探头发生几何参数变化、但超声测量保持正常的使用场景。在该仪器轻微损坏后且不能及时返厂时，仍然可以使用该仪器进行现场测量。并在返厂校正后，使用本发明提供的方法基于测量数据快速地恢复得到真实的观测数据，经济、高效、快速、精度高。

附图说明

图1为现有技术中测量仪器的测量过程示意图；

图2为测量仪器的一个超声探头的支撑支架因外力撞击而脱漆的示意图；

图3为本发明实施例中风速测量及风影校正的示意图；

图4为本发明实施例中声速测量及侧向风校正的示意图；

图5为根据本发明实施例中说明的方法对超声温度进行修正前、修正后以及另一温度测试设备同时测量的低频温度数据对比图；

图6为测量仪器正常时的坐标系的示意图。

具体实施方式

为了更好的解释本发明，以便于理解，下面结合附图，通过具体实施方式，对本发明作详细描述。

本发明的数据修正方法是针对Campbell Scientific Inc公司的三维风速和超声温度测量仪而提出的。IRGASON风速仪是一台精密仪器，出厂或校准均需采用坐标测量机(CMM，Coordinate Measurement Machine)进行出厂参数的校准及标定。

以下以IRGASON风速仪在某野外观测基地进行风速、升温及热辐射观测时获取的数据为例，对本发明提出的数据校正方法进行说明。

在执行野外现场观测任务时，经过数据比对，发现该IRGASON风速仪数据异常。根据下述的公式17至公式23，利用超声风速数据间接地精确计算出高频的气温波动数据(这里的高频是指10HZ、20HZ甚至更高)。

在进行数据验证工作时，分别从IRGASON测量的高频率气温数据中选取了三个典型观测日的每半小时间隔的平均温度数据，如图5(a)中的蓝色线所示；相应地，选取该普通探头同时期测量的半小时平均温度值，如图5(a)中的绿色线所示。鉴于普通探头没有任何外力损坏，可以判断为数据是正常的，那么通过比较蓝色线的温度数据和绿色线的温度数据，可以直观地判断出IRGASON观测的气温数据是异常的，进而可以推断得出IRGASON观测的风速数据也是异常的。

在图5中，选取了三个典型观测日的半小时的平均数据，图中蓝色线为IRGASON测量的高频率气温数据的半小时平均值，绿色线为普通探头测量的半小时平均值，因为普通探头没有任何外力损坏，可以判断为数据是正常的，那么比较蓝线的数据和绿线的数据，可以直观判断出IRGASON观测的气温数据是异常的，进而可以判断出测量的风速数据也是异常的。因为超声风速温度计本身是有运行状态监测系统的，会通过一些运行指标判断数据是否正常，并通过诊断信息提示。但这是这种轻微的磕碰，并没有实际影响设备运行状态(设备本身的监测并没有觉得异常)，但通过数据比较发现，观测结果已经异常了，所以这种条件下的异常判断和数据修复是比较有难度也很有意义的。

因此一旦由于外部力量导致一对探头两端的距离发生改变，如果没有及时调整测量系统中固化的参数，距离变化之后的所有测量数据不再准确，而作为一种高频、精确的测量设备，对关键计算参数细微的改变都会导致最后计算结果无法使用。

如果没有及时调整测量系统中固化的参数，距离变化之后的所有测量数据不再准确，而作为一种高频、精确的测量设备，对关键计算参数细微的改变都会导致最后计算结果无法使用。

在野外观测中认为一对探头两端距离发生改变，主要是由于探头在运输、安装或者观测期间发生了外力磕碰。如果是客观确认距离变化，必须要回到实验室或者工厂，用专门的设备进行重新测量，并比较测量前后两次的数据变化进行判断。

特别地，以目前的经验判断，如果仅靠目视就能明显发现探头距离发生改变，并发生形变，那么该设备已不太可能进行正常观测计算，已无数据输出，也就没有数据恢复的意义了。

如果由于外力导致探头两端距离发生变化时，如果该种情况较轻，如图2所示，没有影响超声波信号的强度，那么设备运行状态不会发出警告，但是最明显的变化是三个方向的风速数据已经偏离了正常值，但一般观测站点上不会再同时配备同类型的超声仪器，所以超声数据的异常不太容易发现。

不过一般观测站点会同时配备空气温度观测设备，通过该类型设备与三维超声风速/温度仪测量的温度数据相比较，可以比较直观地发现由于距离变化而导致的数据异常。

在野外观测中认为一对探头两端距离发生改变，主要是由于探头在运输、安装或者观测期间发生了外力磕碰，并且数据异常而进行的自我判断和分析。如果是客观确认距离变化，必须要回到实验室或者工厂，用专门的设备进行重新测量，并比较测量前后两次的数据变化进行判断。以我们目前的经验判断，如果仅靠目视就能明显发现探头距离发生改变，并发生形变，那么该设备已不太可能进行正常观测计算，已无数据输出，也就没有数据恢复的意义了。

另一方面，因为该仪器是有运行状态监测软硬件系统的，可通过预先设定的运行指标判断该仪器是否正常工作，及初步判断数据是否正常，并提示诊断信息。

通过检查该IRGASON仪器，如图2所示，发现该IRGASON仪器在一处的探头被轻微磕碰。尽管这种磕碰并没有实际影响设备的运行状态，因为设备本身的运行状态监测系统并没有诊断出设备工作异常。但通过观测的温度数据比较可以发现，该IRGASON仪器的观测结果是异常的。

在野外现场观测时，大气观测适用的测量时间窗口一旦错过将不可再现。在确认该仪器的超声测量结果无异常的前提下，仍旧利用该IRGASON仪器进行了执行了既定的现场观测计划，并在随后将其返厂重新进行了校准和标定。

该IRGASON(序列号:1131)几何损坏前之前的出厂参数(2014.9.9)及几何损坏之后重新校准后的出厂参数(2016.3.6)如下表1所示。

表1

本次实验发明人发现数据异常的原因是：理论上根据下述的公式17、18、19、20、21、22、23，超声风速的数据可以精确恢复出高频的气温波动数据(10HZ、20HZ甚至更高)。

具体地，在得到错误的三维超声的数据，即公式5中的u_x,u_y,u_z都是错误的，其错误来源其实就是公式5中的矩阵A(即表征三维超声探头三对探头之间距离构建的坐标系，因为工艺的问题，每个产品本身都有其各自的矩阵A，而磕碰正是造成了矩阵A的变化，所以再按照以前的系数计算，得到的u_x,u_y,u_z必然是错误的)，因此基于公式5将u_x,u_y,u_z反算出u₁,u₂.u₃，再将设备返回工厂后测量计算出新的矩阵A_T，再将原始数据u₁,u₂.u₃，与新的矩阵A_T计算出新的三维风速u_xT,u_yT,u_zT，进而计算出新的空气温度。

进一步地，在于在原始风速观测中，设备探头本身会对u₁,u₂.u₃，测量本身造成影响(类似于树荫对阳光照射的影响，称为阴影订正)，如图3所示，那么需要安装公式11订正，因此在数据修复中，同样需要进行订正，但如果按照公式11进行计算的话，由于未知数个数大于方程个数，就需要进行迭代，但是这样的迭代计算会非常繁琐和复杂，因为本身的数据是高频数据，10HZ，这样进行长时间迭代的话，对于数据恢复订正工作来说无疑是非常大的影响，于是采用了一个简单的假设，例如可以比较公式B1与公式11，大大简化了迭代的步骤，同时根据实际工作经验和理论，设置了符合实际的判断阈值，极大地减少了迭代步骤，但没有影响最后的校正结果。同样如图5(c)所示，灰色线就是订正后的温度数据，与绿色线已经基本吻合，在观测误差标准内浮动。

同样的道理，温度数据修复中，在原始风速观测中有个侧向风的影响，如图4所示，所以进行侧向风的订正，如公式33所示。

本实施例中修复风速的步骤如下：

步骤M1、基于公式10，基于错误的风速数据，和由于磕碰造成形变产生的错误的形状系数矩阵A，得到原始的超声风速仪坐标系下的风速数据u_{t1_1},u_{t2_1},u_{t3_1},

步骤M2、基于公式11的迭代过程，重新再现阴影效应对原始风速数据u_{t1_1},u_{t2_1},u_{t3_1}的影响，并进行订正；如图3所示；

步骤M3、基于公式12和公式13，基于三维超声观测的基于原理与公式3、4、5、6类似，得出订正后的新的原始风速u_c1,u_c2.u_c3；

步骤M4、重新测量磕碰后的超声的新的几何形状，得出新的形状系数矩阵A’,基于公式14得到新的订正后的三维超声数据u_cx,u_cy,u_cz.上述的下标C为修复的意思。

本实施例中温度正常测量的步骤如下：

步骤N1、基于公式17、18表现了超声温度的基本测量原理，通过记录信号的反射时间计算原始超声信号(该信号为测量风速和温度的原始信号，风速和温度都是基于该信号的来回反射时间计算得到的)的速度；

步骤N2、基于公式19和公式20，修正了侧向风对于原始超声信号速度的影响；

步骤N3、基于公式21、22和23，计算了修正后的超声信号速度计算大气温度的过程。

以上是超声温度正常测量原理。

温度数据修复的步骤是：

步骤N4、在已知修复后的超声原始风速数据的基础上，基于上述基本原理，在公式24的基础上，进行推导，得到温度修正的公式27.

步骤N5、公式28到公式31为公式27中一个不确定项的近似化简推导过程

步骤N6、公式33为最终修正后的超声温度的计算公式。

具体地，以下详细说明上述的风速修复的过程。

如图1所示，声速风速计会测量在每条声波路径中超声信号上行飞行的时间t_ui和下行飞行的时间(t_di)，下标i表示声波路径的顺序，可以是1,2或3。超声信号在第三声波路径上的上行飞行时间如前面的公式(1)，公式(1)中沿着第三声波路径的d₃是在出厂时或重新校准过程中使用坐标测量机(CMM)精确测量的长度，c₃是声速，u₃是沿着第三声波路径的气流速度。

超声信号在第三声波路径上的下行飞行时间如公式(2)所示：

进一步地，声波风速计输出三维风数据的算法如下：

联立方程(1)和(2)，可以得到u₃，如式(3)所示：

同样地，可以得到沿着第一声波路径的气流速度u₁，沿着第二声波路径的气流速度u₂；也即沿着第i声波路径的气流速度u_i可以如式(4)所示：

与d₃类似，沿着第一声波路径的d₁、沿着第二声波路径的d₂也是在出厂时或重新校准过程中使用坐标测量机(CMM)精确测量的长度。

从图1和图6中可知，第一声波路径、第二声波路径、第三声波路径在空间不满足正交关系。

因此，如式(5)所示，将从这三个非正交声波路径测得的风速u_i(i＝1,2或3)通过变换矩阵A变换到一个右手正交坐标系中，从而将这3个声波路径方向上的风速表达为更为常用的3维正交风速。

具体地，该右手正交坐标系为仪器坐标系。该仪器坐标系的定义如下：该坐标系的原点位于声源测量体积中心的原点，其中，ux-uy平面与风速计结构中内部装置所构建的图像平面相平行，并将风速计结构中的内置灯泡调平，该第一声波路径位于uy-uz平面内。

在本发明实施例中，变换矩阵A是一个3×3的矩阵，通过精确测量得到的这3条声波路径相对与该仪器坐标系的角度参数构建得到。因此，每台风速仪所适用的变换矩阵A都是其独有的，并嵌入在其软件系统的固件中。

经过如式(5)的处理后，最终该风速计输出的3维正交风速既是其在该仪器坐标系中的u三个分量u_x，u_y和u_z。

源自超声探头的阴影效应，测试得到的每一声波路径上的风速u_i在幅值上是小于其真实值u_Ti的。

由于来自声波换能器的阴影效应，通常认为测试得到的每一声波路径上的风速u_i在幅值上是小于其真实值u_Ti的。下述公式中符号u_{Ti_n}的下标T表示“真”，下标_n表示u_{Ti_n}是在阴影校正步骤中的第n次迭代中估计得到的u_{i_n}，真值u_{Ti_n}是在阴影校正步骤中利用u_i来逼近的。

目前，在IRGASON OS 5或更新的版本中，阴影校正步骤已经是一个常规的可选项。

若在该仪器的固件中已经默认增加了阴影校正步骤时，在式(5)中，参与坐标系转换的风速数据就是经过n次迭代后逼近的真实风速u_{Ti_n},

在IRGASON声波风速计中，与声波换能器的尺寸、每一声波路径的几何参数有关的校正方程由式(7)给出：

其中，α_i是总真实风向量U_T与沿着声波路径i的真实风向量u_Ti两者之间的夹角，该夹角的值在这两个向量被精确估计之前是未知的。参考图1和图3，式(7)可以改写为：

其中，U_T为总真实风向量的幅值，由式(9)给出

式(8)应用于阴影修正有两个问题：1)尽管u_Ti在方程(8)的左边以解析解的形式给出，但是u_Ti在方程(8)右边的平方根下存在二阶项，因此u_Ti的解析解并不容易得到。2)U_T是是未知的，因为u_x，u_y和u_z均是在在阴影校正之前从u₁，u₂和u₃中计算得到的。

从式(8)可以看出，u_i的数值与u_Ti的数值比较接近；因此，可以认为，根据式(5)得到的u_x,u_y和u_z与根据式(6)得到的u_x,u_y和u_z是非常接近的。因此，利用式(8)采取迭代算法来估计u_Ti从而进行阴影校正是一个可行的途径。

在第一次迭代时，作为u_Ti的估计，方程(8)右边的u_Ti可以用u_i来替代。考虑到应该使用方程式中的u_x,u_y和u_z来计算U_T。U_T可以根据式(5)中的u_x，u_y和u_z来计算。

多次迭代，可以确保上次迭代和先前迭代之间的u_x，u_y和u_z的偏差小于1mm s^-1，也即小于1.96σ，其中σ是u_x，u_y和u_z之间的最大精度(即恒定风速时的标准偏差)。

最后一次迭代中确定的u_{T1_n}，u_{T2_n}和u_{T3_n}最终用于式(6)来计算三维风速u_x，u_y和u_z以作为声波风速计的输出。

此外，如式(4)至(6)所示，声波风速计测量通过t_ui和t_di来计算u_x,u_y,和u_z这三维风速，同时计算声波温度(Ts)。因此，式(1)中的声波路径长度(di)及在方程(5)和(6)中的转换矩阵A在出厂或校准过程中是嵌入在该声波风速仪的固件中的。

仪器属于数据异常状态时，则音速路径长度和音速路径角度必然已经不同于在d_i和A的数值被嵌入到固件中时的情况。此时，公式(4)中的di和(5)、(6)对应的矩阵A对于此风速计不再有效。所以u_x,u_y,和u_z的输出仍然基于生产或校准过程中的嵌入式di和A是错误的。纠正错误的输出，u_x,u_y,和u_z需要转换回t_ui和t_di并重新计算。

为纠正错误的输出，u_x,u_y,和u_z需要被转换回t_ui和t_di，并根据现场测试t_ui和t_di时的真实声波路径长度和真实声波路径角度使用t_ui和t_di重新计算u_x,u_y,和u_z；该仪器在现场部署前确定的该仪器的几何参数与经过物理变形后使用该仪器在现场实际测试时，该仪器真实的几何参数是不同的。

在现场试验结束后，IRGASON(SN：1131)返回到厂家。与制造过程相同，使用CMM重新测量了IRGASON(SN：1131)的声波路径长度和声波路径长度和角度。

重新测量的声波路径长度由d_Ti(i＝1,2或3)表示，并且重新测量的声波路径角度用于将转换矩阵A重建为A_T。d_Ti和A_T已经更新到IRGASON的固件内以供将来现场使用。d_Ti和A_T也用于修正在返厂重新测量几何参数之前在测试现场获取的u_x,u_y和u_z和Ts的数据。针对输出的u_x,u_y,和u_z，是否包括换能器阴影校正，是不同的修正步骤。

1)包括换能器阴影校正

将三维坐标系统中的u_x,u_y和u_z转换为探头阴影校正之后沿声波路径的风速。

利用式(B5)，ui可以表示为

U_T可以用公式(9)来计算；u_{Ti_m}可以使用u_{Ti_n}进行合理近似，因为u_{Ti_m}和u_{Ti_n}足够接近，可以确保u_x、u_y、u_z收敛于它们的测量精度。

使用ui和di，可以得到方程(4)中方括号内的时间项；

同样地，根据公式(4)，使用u_Ti，可计算得到沿第i个声道的气风速度u_ci：

在附录B中的声波换能器阴影修正迭代算法中，用u_ci替换u_i，且u_ci被进一步地修正为u_{cTi_n}。

根据公式(6)，修正后的三维坐标系中的3维风速矢量可以由下式得到：

2)不包括换能器阴影校正

将三维坐标系统中的ux，uy和uz转换为沿每一声波路径的风速：

根据公式(12)和公式(13)，重新计算沿每一声波路径的风速。

根据公式(5)，修正后的三维坐标系中的3维风速矢量可以由下式得到：

另外，以下详细说明上述的温度数据修复的过程。

声波风速计输出声波温度的算法

根据公式(1)和公式(2)，还可以得到：

相似地，可以得到c₁和c₂(如图1和图4所示)。

对应于ci的式(17)，i＝1,2,3,如下：

如图4所示，c_i是沿着声道i测得的声速。

当侧风(u⊥i)或与声道i相垂直的风分量是零时，c_i就是真实声速(c_Ti)。但是，侧风很少是零，因此，需要将c_i修正为c_Ti。

如图1和图5所示，可以由下式计算真实声速(c_Ti)：

再次参见图5，有：

根据声波温度的定义，单位为K，沿第i声道的声速(Tsi)应该表示为：

其中，其中γ_d(值为1.4003)是恒定压力(1,004J K^-1kg^-1)下的干燥空气比热与恒定体积(717J K^-1kg^-1)下干燥空气比热的比值；Rd(值为287.04J K^-1kg^-1)是干燥空气的气体常数。

最终，声波风速计输出的声波温度(Ts，以℃为单位)是对三个声道中的声波温度数值求平均值，如下式：

利用式(20)替换c_Ti，利用式(18)替换c_i，则有：

恢复声波温度数据的过程说明

方程(23)表明，给定di，声波风速计使用其测量的飞行时间t_ui和t_di估计声波温度；参见式(4)，声波风速计也使用t_ui和t_di计算沿声波路径的风速(u_i或u_Ti，如果选择了换能器阴影校正选项)；利用式(9)计算总风速(U_T)，其中，沿三维坐标系方向的各分量数据是从u_i或u_Ti根据式(5)或式(6)，式(16)或式(14)，转换得到的。

在声波风速计意外地发生几何变形后，其声速路径长度和声速路径角度可能会与di和A在嵌入固件时的数值不同，因此，方程(23)中使用的di和在方程式(5)或式(6)的A的数值并不再适于此风速计在测量时使用。因此，该仪器仍然基于之前嵌入到固件中的di和A而输出的u_x，u_y，u_z和T_s必然不能代表现场测试时的风速。

前述已经公开了使用重新测量的声波路径长度(d_Ti)和使用重新测量的声波路径角度而构建的A_T来恢复三维风数据的步骤。以下提供使用d_Ti和从上述步骤中修正的三维风速数据来修正Ts的步骤。

基于公式(20)，声波路径i上经过侧风校正后的声音速度可表示为：

其中，c_ci是声波路径i上未知的修正声音速度，UcT为

在式(22)中，利用c_cTi ²替换c_Ti ²,修正声音温度(Tcs in℃)表示为：

现在，需要推导c_cTi ²。

式(24)-式(20)，可以得到：

用式(26)去替代式(25)中的c_cTi ²，将记为 记为可以得到：

式(27)中，只有是未知的。

根据式(18)，c_ci ²可以由下式得到：

从而，该未知项可以由下式(29)得到：

在式(29)中，唯一的一个未知项为c_i ²。

根据式(20)，式(29)可以被表示为：

在该公式(30)的等号右边，c_Ti ²是未知的。

但是，式(29)的等式右边从整体上是一个微分项，其中，c_Ti ²可以由尽可能接近c_Ti ²的临近值合理地逼近，该值可以由式(22)中对c_T1 ²,c_T2 ²,c_T3 ²的平均而确定。

在没有测量误差和随机误差时，这三个c_Ti均与风速(流体速度)无关；因为，这三个数值均是声音速度的真值而并非是沿着各自的声音路径测量得到的声音速度。

因此，c_Ti ²可以合理地由这三个值c_Ti ²(c_T ²)的平均值来逼近，如式(31)：

其中，c_T ²可以由式(22)计算得到：

由于需要用c_T ²来替换c_Ti ²，式(31)等式右侧作为整体，其相对误差将小于4％；即使在最糟糕的情形下，如，由于c_Ti ²值之间的差异导致超声温度的变化量在空气温度-30℃时达到了10℃(即U_T＝0，u_Ti＝0)。

将式(31)代入式(27)，可以得到：

在式(33)等式的右边，Ts之后的全部项即为修正声波温度项。

5、应用

在不考虑阴影修正时，依次利用式(15)、(12)、(13)、(16)即可获得修正风速数据。

在考虑阴影修正时，利用式(10)、(11)、(12)、(13)、(14)(15)(16)即可获得修正风速数据。

基于以上修正风速数据，利用公式(9)、(32)和(33)来获得修正超声温度数据。

数据验证

如下附录C所示，利用空气温度、相对湿度和大气压力获取声波温度。

使用HMP155A和IRGASON内置气压计的研究级传感器测量了空气温度，相对湿度和大气压力，并且这些变量的数据也以10Hz(每秒10条记录)存储。这些数据可用于估计声波温度(参见附录C：从空气温度，相对湿度和大气压力得到的声波温度)。使用公式(33)修正的声波温度数据与测量声速范围内的计算声波温度的三个代表值进行比较：图5(a)中的-20.01±0.14℃，图5(b)中的-9.06±0.13℃和图5(c)中的-1.90±0.22℃。测得的(即未恢复的)和计算出的声波温度值之间的差异，图5(a)中的9.60±0.14℃，图5(b)中的9.53±0.17℃和图5(c)中的8.93±0.24℃被缩小到0.99±0.14℃，0.57±0.17℃和-0.25±0.24℃，分别作为回收和计算声波温度值之间的差值。声波风速计IRGASON的准确度为±0.5℃(与声波风速计设计师的个人通信)以及来自HMP155A(Campbell Scientific Inc.，1990)的空气温度低于0℃和相对湿度1.2％的精度为±0.2～0.3℃，从中得出计算出的声波温度(附录C)，如果恢复值和计算值之间的平均声波温度差异在±0.80℃范围内或甚至可能被认为是可能的更宽范围内，则可以合理地判断回收的声波温度数据是令人满意的正确测量和计算声波温度之间的可能差异范围。如图6所示，公式(33)在采用声速测量仪测量声波温度数据的几何形变方面显然做了很好的工作，但在图5(a)的情况下几乎不能令人满意(即0.99±0.14℃，这是回收和计算之间声波温度的差异)，虽然0.99±0.14℃的范围与±0.80℃没有显着差异。如果所有的c_Ti在其几何变形中都不是由声速风速仪测量的话，c_T的c_Ti近似值可能是由于c_Ti的近似值所致，但在本研究中并非如此。根据公式(22)中，不可能从Ts得到单独的c_Ti，这是来自任何声速风速计的声波温度的唯一输出。现在c_T1 ²、c_T2 ²、c_T3 ²，的平均值是已知的，并且声道长度的变化是已知的。可以估计三种声速之间的差异，并将它们的平均值(c_T2)调整到c_T1 ²、c_T2 ²、c_T3 ²近似值，尽管精确值是不可能的。调整值至少可以在一定程度上反映c_Ti ²之间的变异性，并可合理预期改善数据修正。

调整

侧风校正后的声速(c_Ti)与风速无关。考虑到空气密度和大气压力，没有风，c_Ti等于从声道i测得的声速(c_i)(见公式(19))。在这种情况下，没有风，t_ui和t_di在公式(18)是相同的并且可以用t_i表示。因此，公式(18)在这种情况下相当于

在公式(33)，c_T ²是三平方c_Ti的平均值(见公式(22)和(32))，但是个体c_Ti是未知的；因此，为了恢复改善，必须通过合理的调整从c_T ²估计。由于几何变形，c_T ²和c_Ti ²之间的幅度差异必须与c_Ti ²误差相关。平方公式(34)的两边得到

cTi2的总分化由下式给出

考虑到飞行时间即使在其几何变形中也可以通过声速风速计(Δt_i＝0)正确测量，该等式变为

在数学上的分化，c_Ti ²可以合理地用c_T近似，由下式给出

这是c_Ti ²偏移c_T ²的误差。这个误差可以合理地用来表示c_Ti ²远离c_T ²的偏差。三个c_Ti ²值远离c_T ²的偏差是三个远离c_T ²的c_Ti ²之间变异性的量度。

尽管个体c_Ti ²未知，但三个c_Ti ²的平均值被称为c_T ²。这个平均值在调整之后应该保持不变，因为c_Ti ²与c_T ²之间的变化在不同的范围内调整。如果调整后的c_Ti ²的平均值不等于c_T ²，则所有调整后的c_Ti ²都应该加上或减去相同的不断地将三个调整后的c_Ti ²值的平均值作为c_T ²，但c_Ti ²值之间的变化保持不变。这个常数必须是三个c_Ti ²值的平均值。基于这些分析，将c_T ²调整为c_Ti ²可以被构建为

用这个公式代替公式(30)中的c_Ti ²，然后用公式30的结果替换公式(27)中的公式(27)可以写成

在方程的右边，Ts之后的整个项是调整后的声波温度恢复项。这个数据用公式33计算过，又用公式40计算过。显然，这个公式40比公式33计算的好。如图5中的图a至c所示，修正和计算之后的声波温度差异分别降低到0.81±0.14℃，0.38±0.17℃和-0.45±0.24℃。统计意义上，这些差异值落入±0.80℃的范围内。公式(40)被认为比公式33更好。这是令人满意的。最终，公式(40)用于数据恢复，并且如图5所示被合并到软件中。

验证三维风速修正

虽然没有明确验证，但是在恢复的声波温度数据的验证中隐含地验证了恢复的三维风速数据，因为恢复的三维风速数据用于恢复声波温度数据(见公式(33)和(40))。声波温度的横风校正需要3D风资料。如果3D风没有得到很好的恢复，声波温度数据就不能得到满意的恢复。本研究中声波温度数据的满意恢复隐含验证了三维风资料的满意恢复。

恢复到计算声波温度的可比性

恢复的声波温度来自快速响应声速风速仪的测量结果，计算出的声速温度来自慢响应空气温度和相对湿度探头以及IRGASON内置气压计的测量结果；因此前者反映了高频声波温度的波动，后者反映的波动并不像前者那么高。这样，来自同时测量的一对恢复和计算的声波温度值(即，时间序列数据文件中的相同记录)不具有可比性。两组之间的差异是没有意义的。因此，数据比较中使用了半小时内恢复的和计算的声波温度值之间的平均差值。

为了更好地比较，通过滞后时间对准的恢复和计算的声波温度值计算差异。根据空气温度，相对湿度和大气压计算出的声波温度；使用慢响应传感器进行测量；被认为是滞后于回复之一，以响应声波温度的波动。大约10秒的滞后时间在经恢复的声波温度的时间序列与计算的声波温度的不同时间延迟(即，时移)序列的交叉相关的最大值的经验最大值上发现。恢复的和计算的声波温度值之间的差异使用恢复的声波温度而没有时间滞后并且计算声波温度并且时间滞后10秒来计算。使用滞后时间可能是不必要的，但可能会使比较尽可能合理。

在较低温度下恢复高于计算的声波温度

参见图5。与计算的声波温度相比，由公式(40)在-20.01℃时高出0.81±0.14℃(图5a)，在-9.06℃时高出0.38±0.17℃(图5b)，然而在-1.90℃时甚至是0.45±0.24℃(图5c)。这种与温度差异的趋势可能与声速测量仪在不同温度下的性能和温度精度较低有关和湿度探头在较低的温度范围内。

在空气温度约20℃(即制造温度)的制造环境下测量声波风速计的声波路径长度和几何形状，并将其嵌入到其固件中以用于现场应用。然而，在制造温度以上或以下，由于声波风速计结构的热膨胀或收缩，声波路径长度必须在制造温度下长于或短于声波路径长度，而声波路径的长度值在固件内部不变。因此，声速测量仪可能会低估或高估声速，从而导致声速过低。当温度不高于或低于制造温度时，低估或高估可能不明显，而风速仪必须工作得最好围绕制造温度。在这项研究中，声波风速计的工作空气温度低于-20℃内，声道在一定程度上变短，从而可能影响其测量性能。虽然在文献中未发现低气温或高气温的情况下对声速测量仪的测量性能进行评估，在低于制造温度下估计的声波温度如图6a至图6c所示。

尽管在不同的空气温度下，温度和相对湿度探头以及IRGASON内置气压计的性能(其测量值用于计算声波温度)的性能比声波风速计更稳定，但它们的准确度是最佳温度为20℃，温度远离20℃时会降低。例如，HMP155A的空气温度准确度在20℃时为±0.1℃，在-20℃时为±0.25℃，相对湿度(RH)的准确度为±(1.0+0.008RH)，在-20℃时为±(1.2+0.012RH)％。图5a中较低温度下回收声波温度值和计算声波温度值之间的较大不一致可能是由于空气温度越低，HMP155A和IRGASON内置气压计的精度越低所致。

计算声波温度时的辐射

看图5c。与使用公式40恢复的声波温度相比时，计算出的声波温度在12:00至12:30期间升高0.45±0.24℃，在12:15至12:27期间甚至升高0.65±0.19℃，这可能部分原因是在HMP155A的辐射屏蔽层上有较高的750W m-2短波太阳辐射进入。计算得出声波温度来源于HMP155A和IRGASON内置气压表的气温，相对湿度和大气压力的测量结果。安放在防辐射罩内的HMP155A会受到太阳辐射的污染。即使是使用防辐射罩遮挡HMP155A的阳光，当使用这种防护罩时，在太阳光下和电子电源下的传感器产生的任何热量都被低效地消散。因此，屏蔽罩内的空气和HMP155A传感元件比周围的环境空气温度高。辐射屏蔽内部的空气温度取决于屏蔽结构，环境风速和其他环境条件。在图5c的情况下，输入的短波辐射量为750Wm-2时，辐射屏蔽层内的空气含量较高的情况并不少见。在我们的研究中，这种较高的气温可能直接导致计算声波温度的高估(见公式C1)

前述的3D风的阴影校正是用户的可编程选项，并且声波温度的侧风校正被应用。在选项阴影修正的情况下，三维风数据可以使用公式10至16恢复。如果不是，按照公式(15)，(12)，(13)和(16)的顺序。根据三维风恢复过程的数据，声波温度数据可以用公式(9)，(32)和(38)至(40)恢复；因此，三维风资料和声波温度的满意恢复情况能通过声波温度数据恢复的满意度反映出来。

本实施例可以在声波风速计内部进行或不进行阴影校正的情况下恢复三维风速数据，也可以在声速测量仪内部恢复侧风校正的声波温度数据。

附录A转换矩阵

在微气象应用中，风速用仪器或自然风的三维(3D)直角坐标系统表示，但声波风速计测量沿其三条非正交声波路径的流速(见图1和图6)。因此，对于应用，需要将沿着三条声波路径的流速转换为3D正交坐标系，参考如图6所示的声速测量仪的几何结构(即3D正交仪器坐标系)。如果u_x和u_y分别是x和y方向上的两个水平速度，u_z是z方向的垂直速度(图6)，该系统的定义是xy平面平行于风速计球泡平坦的仪器平面，yz平面上的第一条声道和源于中心的测量量。由于声波路径很容易在系统中显示，所以沿着声波路径的流向量可以在三维中正交表示。在定义的三维正交仪器坐标系中(见图6)，假设第一条声波路径相对于z轴的顶点角为θ1，方位角φ1等于90°，固定在xy平面上(见图6)，沿第一路径(u1)的流速可以分解为u_x，u_y和u_z三个正交分量：

u₁[0 sinθ₁ cosθ₁]'＝[u_x u_y u_z]' (A1)

类似地，沿着第二和第三路径的风矢量可以分别分解为三个正交分量，如下所示：

其中θ2和φ2是第二声道的顶点角和方位角，并且

其中θ₃和是第三声波路径的顶点角和方位角。结合公式(A1)至(A3)得到：

名义上为IRGASON的声道，θ1，θ2和θ3全部为30°，和为330°和210°(见图6)。给定这些角度使用坐标测量机的测量数据进行计算，并与声波路径长度一起在几何变形前后的IRGASON序列号1131中列出，如上表1。

使用表1中的数据，公式A4中的矩阵A和它的转置在几何变形之前(即，在IRGASON固件中使用，尽管在变形之后在现场无效)得到

和

IRGASON几何变形后的两个矩阵变为：

其中下标T表示“真”，在IRGASON变形之后，它应该在现场使用，尽管它没有被使用，并且

Matrix A-1，A_T和A^-1 _T用于我们的数据恢复，A也用于固件内部以及阴影校正是可选的情况。

附录B声波换能器阴影修正的迭代算法

给定变换矩阵A，使用公式(5)，测得的矢量[u₁ u₂ u₃]′在三维正交仪器坐标系中转化为风向量[u_x u_y u_z]′，随后，使用公式(9)中U_T。在公式(8)右边的平方根下用u_i代替u_Ti。给出了第一次迭代的近似方程：

其中下标i是1,2或3，并且u_Ti的下标_1表示它是从第一次迭代计算的。

第一次迭代

方程(B1)用于第一次迭代中的换算器阴影修正。

第二次迭代

使用公式(9)，U_T被重新计算。在公式(B1)右边的平方根下用uTi_1替换ui。给出了第二次迭代的近似方程：

第三次迭代……第n次迭代

其中下标m＝n-1，使用公式(9)，U_T再次被重新计算。,类似于u_{Ti_2}的计算，u_{Ti_n}使用方程计算

以确保u_x，u_y或u_z在上一次迭代和上上一次迭代之间的差值为小于1mm s-1≈1.96σ,其中σ是u_x，u_y和u_z中最大精度(即恒定风的标准偏差)。在u_x，u_y和u_z测量范围内的数字测试结果表明，迭代大多收敛于n＝2，全部收敛于在n≤3。

附录C空气温度，相对湿度和大气压力下的声波温度

如果在现场测量空气温度(以℃为单位)，相对湿度(RH以百分比表示)和大气压力(P以千帕表示)，则可以使用众所周知的方法计算声波温度(以℃为单位)方程(Kaimal和Gaynor，1991)：

其中e是空气水蒸汽压力(kPa)，可由T，RH和P计算如下。给定T和P，可用Buck(1981)计算饱和水蒸汽压力(单位kPa)：

其中fw(Td，P)是增强因子：

f_w(T,P)＝1.00041+P[3.48×10^-5+7.4×10^-9(T+30.6-0.38P)²] (C3)

使用空气相对湿度的定义，空气水蒸汽压力由下式给出：

将测得的T和P以及计算的e提交到公式(C1)中，可以计算声波温度。

上述各个实施例可以相互参照，本实施例不对各个实施例进行限定。

最后应说明的是：以上所述的各实施例仅用于说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述实施例所记载的技术方案进行修改，或者对其中部分或全部技术特征进行等同替换；而这些修改或替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的范围。

Claims (10)

1.一种测量三维风速和超声温度仪器的风速数据修正方法，其特征在于，包括：

S1、获取观测系统中测量三维风速和超声温度的仪器测量的第一风速数据和第一温度数据，以及获取观测系统中另一温度传感器测量的第二温度数据；

S2、将第一温度数据和第二温度数据进行比较，判断所述仪器是否处于数据异常状态；

S3、在所述仪器发生异常时，根据所述第一风速数据、和在数据异常状态下所述仪器的性能参数，获取真实风速数据。

2.根据权利要求1所述的方法，其特征在于，所述步骤S3包括：

所述第一风速数据为所述仪器处于数据异常状态时获取的错误的风速数据u_x,u_y.u_z；

S31、根据所述第一风速数据，获取原始测量的用于计算风速的第一测量数据u₁,u₂,u₃；

S32、根据数据异常状态下所述仪器的性能参数，获取所述仪器对应的数据异常状态的形状系数矩阵A_T；

S33、根据所述第一测量数据和所述形状系数矩阵A_T,获得修复后的作为真实风速数据的第二风速数据u_cx,u_cy,u_cz。

3.根据权利要求1所述的方法，其特征在于，所述步骤S3包括：

所述第一风速数据为所述仪器处于数据异常状态时获取的错误的风速数据u_x，u_y，u_z；

S031、根据所述第一风速数据，获取原始测量的用于计算风速的第二测量数据u_{t1_1},u_{t2_1},u_{t3_1}；

S032、根据数据异常状态下所述仪器的性能参数，获取所述仪器对应的数据异常状态的形状系数矩阵A_T；

S033、根据修正后的第一测量数据和所述形状系数矩阵A_T,获得修复后的作为真实风速数据的第二风速数据u_cx,u_cy,u_cz。

4.根据权利要求2或3所述的方法，其特征在于，所述性能参数包括：处于数据异常状态下所述仪器中每一对探头之间的距离、相邻探头之间的夹角。

5.根据权利要求2所述的方法，其特征在于，

所述步骤S31包括：

所述第一风速数据为根据下述公式获取的；

其中，A为嵌入在所述仪器中的该仪器处于正常状态时使用的形状系数矩阵；

所述u_x，u_y，u_z为第一风速数据，所述u₁,u₂,u₃为所述仪器的每对探头的测量数据；

相应地，根据公式：获取真实测量数据u₁,u₂,u₃；

所述步骤S33包括：

根据公式获取第二风速数据u_cx,u_cy,u_cz；

其中，A_T为所述仪器对应的数据异常状态的形状系数矩阵。

6.根据权利要求2所述的方法，其特征在于，

所述步骤S031包括：

所述第一风速数据为根据下述公式获取的；

其中，A为嵌入在所述仪器中的该仪器处于正常状态时使用的形状系数矩阵；

所述u_x,u_y，u_z为第一风速数据，所述u_{t1_1},u_{t2_1},u_{t3_1}为所述仪器的每对探头的测量数据通过阴影校正后得到的真实测量数据；

相应地，根据公式获取真实测量数据u_{t1_1},u_{t2_1},u_{t3_1}；

步骤S033包括：

根据公式获取第二风速数据u_cx,u_cy,u_cz；

其中，A_T为所述仪器对应的数据异常状态的形状系数矩阵。

7.根据权利要求5或6所述的方法，其特征在于，

其中，θ₁为第一条声波路径相对于Z轴的顶点角，为第一条声波路径相对于Z轴的方位角；

θ₂为第二条声波路径相对于Z轴的顶点角，为第二条声波路径相对于Z轴的方位角；

θ₃为第三条声波路径相对于Z轴的顶点角，为第三条声波路径相对于Z轴的方位角；

X、Y、Z为预先定义的三维正交仪器坐标系；

根据获取形状系数矩阵A的方式，获取处于数据异常状态的所述仪器的各个声波路径相对于Z轴的顶点角和方位角，根据所述获取形状系数矩阵A的方式和获取的处于数据异常状态的所述仪器的顶点角和方位角，获取数据异常状态的形状系数矩阵A_T。

8.一种测量三维风速的超声温度仪器的风速数据修正方法，其特征在于，包括：

M1、若确定测量三维风速的超声温度仪器的风速发生异常时，则根据所述仪器当前状态下获取的第一风速数据和异常状态下该仪器的性能参数，获取真实风速数据。

9.根据权利要求8所述的方法，其特征在于，所述M1包括：

所述第一风速数据为所述仪器处于数据异常状态时获取的错误的风速数据u_x，u_y，u_z；

M01、根据所述第一风速数据，获取原始测量的用于计算风速的第一测量数据u₁,u₂,u₃；

M02、根据数据异常状态下所述仪器的性能参数，获取所述仪器对应的数据异常状态的形状系数矩阵A_T；

M03、根据所述第一测量数据和所述形状系数矩阵A_T,获得修复后的作为真实风速数据的第二风速数据u_cx,u_cy,u_cz。

或者，

M11、根据所述第一风速数据，获取原始测量的用于计算风速的第二测量数据u_{t1_1},u_{t2_1},u_{t3_1}；

M12、根据数据异常状态下所述仪器的性能参数，获取所述仪器对应的数据异常状态的形状系数矩阵A_T；

M13、根据修正后的第一测量数据和所述形状系数矩阵A_T,获得修复后的作为真实风速数据的第二风速数据u_cx,u_cy,u_cz。

10.根据权利要求9所述的方法，其特征在于，所述性能参数包括：处于数据异常状态下所述仪器中每一对探头之间的距离、相邻探头之间的夹角。