CN101750514A - 高空驻空飞行器风速风向实时测量方法与装置 - Google Patents

Abstract

本发明是一种高空驻空飞行器风速风向实时测量方法与装置，该方法由测量装置和与风速风向提取算法构成。测量装置由无刷电机驱动、以恒定角速率转动的刚性旋杆，刚性旋杆两端装有测压探头，测压探头连接在一个压差传感器上。当风吹向恒定转动的刚性旋杆时，两测压探头间的压差为一周期性余弦信号。该信号的幅值与风速大小成正比，信号的相位包含了风向角信息。通过本发明的提取算法，可从压差信号中提取出风速风向信息。本发明的特点是在既定风速下，压差信号的幅值与大气密度和旋臂线速度之积成正比。因此，通过合理地选择刚性旋杆转速和旋臂长度，在大气密度很低的情况下也能得到较大的信号量，从而能可靠地提取出风速风向信息。

Description

高空驻空飞行器风速风向实时测量方法与装置

技术领域

本发明主要用于高空驻空飞行器的飞行控制。其作用是提供飞行控制所需的实时相对风速/风向信息(空速/侧滑角)和真风速/风向信息。本发明的领域属于航空测控技术。

背景技术

高空驻空飞行器包括平流层飞艇、临近空间浮升一体飞行器、高空气球等飞行器。近年来在区域大气环境监测、防灾减灾、高分辨率实时监视、预警和导弹防御、反恐、区域通信等需求的驱动下，以平流层飞艇为代表的高空驻空型飞行器引起了各国的普遍重视，美国、日本、俄罗斯、韩国、欧盟、中国等主要国家陆续启动了相关的研究计划，开始了较深入的研究开发工作。

高空驻空飞行器体积庞大，运动速度低，其动力学特性受风的影响很大。此外，驻空要求，使得这类飞行器经常要迎风飞行或悬停。因此对风速、风向进行实时、准确的测量是实施长期驻空以及提高飞行控制品质的前提。

当风速传感器安装在飞行器上时，所测量的风速为相对风速，即大气相对于飞行器的速度。通常更习惯用飞行器相对于大气的速度来描述相对风速，称之为空速。空速(或相对风速)与真风速之间有以下关系：

真风速矢量＝地速矢量-空速矢量

用数学公式表述为

${\stackrel{\→}{V}}_w={\stackrel{\→}{V}}_G-{\stackrel{\→}{V}}_a---\left(1\right)$

其中，

分别为真风速、地速和空速矢量。

地速

可通过飞行器配备的GPS/INS组合导航系统精确测量，所以获取了相对风速(空速)

也就意味着获取了真风速

对大气飞行器来说，准确测量空速的意义不仅限于获取真风速，空速/侧滑角本身就是飞行控制中一个重要的参考量。

在航空领域空速测量的标准方法是采用压差空速计，通过测量动压得到空速。低速时动压与空速之间满足伯努利方程

$\mathrm{\Δp}=\frac{1}{2}\mathrm{\ρ}{V}_a^2---\left(2\right)$

可见，动压与大气密度ρ成正比，与空速的平方成正比。

在高空，由于大气非常稀薄，加上高空驻空飞行器航速较低，导致动压非常小，精确测量空速非常困难。另一方面，由于动压与空速的平方成正比，为了应付极端大风情况，压差传感器的量程还不能太小，这进一步加剧了常规飞行时空速提取的难度。

表1给出了20km高度不同空速所对应的动压值。当速度小于4m/s时，动压＜1Pa，速度为8m/s时，动压也仅有2.8Pa。从如此微弱的信号中，准确提取空速非常困难，要求压差传感器有极高的精度和分辨率。

表1  20km高度不同空速所对应的动压

V(m/s)	1	4	8	12	20	40
							Δp(Pa)	0.0445	0.711	2.85	6.40	17.8	71.1

与空速大小的测量相比，确定风向更为困难，传统的风标式测向装置以及压差测风向装置，在平流层高度均因动压过低无法使用。

航空仪表之外的其他测风装置例如超声波风速计、热线\热球风速仪也因平流层环境过于苛刻而无法使用。以超声波风速计为例，环模实验表明现有产品在海拔7km以上的高度无法正常使用。热线\热球风速仪对环境的要求更为苛刻，一般只用于室内洁净气流的测量。

目前在平流层及以上高度，还没有现成的产品能直接用于为低速飞行器提供飞行控制所需要的实时风速、风向信息。

发明内容

(一)要解决的技术问题

有鉴于此，本发明的主要目的在于提供一种能直接用于为低速飞行器提供飞行控制所需要的实时风速、风向信息的高空驻空飞行器风速风向实时测量方法与装置，以解决现有技术不能对低速飞行器提供飞行控制所需要的实时风速、风向信息的问题。

(二)技术方案

为达到上述目的，本发明第一方面是提供高空驻空飞行器风速风向实时测量装置，在高空稀薄大气环境中，为低速、驻空型飞行器提供风速风向测量信息，该测量装置包括：压差传感器，第一测压探头、第二测压探头，刚性旋杆、传感器安装盘，无刷电机及其控制器、微处理器和附件组成，其中：刚性旋杆对称安装在无刷电机的转动轴上；第一测压探头和第二测压探头分别安装在刚性旋杆长边的两端面上，且第一测压探头和第二测压探头与刚性旋杆旋转平面平行，且第一测压探头、第二测压探头分别与刚性旋杆长边面垂直，第一测压探头和第二测压探头到无刷电机转动轴的距离相等；传感器安装盘安装在刚性旋杆上；压差传感器安装在传感器安装盘上，压差传感器与第一测压探头和第二测压探头连接；无刷电机带动刚性旋杆以恒定角速率ω转动；无刷电机上设有绝对式旋转编码器，绝对式旋转编码器用于输出刚性旋杆的各个时刻转动角；当风吹向转动的刚性旋杆时，第一测压探头和第二测压探头之间产生并输出压差；附件包含气压计和温度计，用于测量大气压强和大气温度；微处理器根据附件给出的气压强和大气温度信号计算出大气密度；微处理器与压差传感器连接，压差传感器输出的压差信号给入微处理器，无刷电机与微处理器连接，微处理器对无刷电机输出的刚性旋杆转动角信号进行高速采样；微处理器根据大气密度、压差信号和转动角信号并采用风速风向提取算法得到飞行器相对大气的相对风速、风向信息，该相对风速、风向信息与导航系统给出的地速信息融合，得到真风速风向信息。

为达到上述目的，本发明第二方面是提供高空驻空飞行器风速风向实时测量方法，利用高空驻空飞行器风速风向实时测量装置方法，在高空稀薄大气环境中，为低速、驻空型飞行器提供风速风向测量信息，该方法的步骤包括：

步骤S1：启动高空驻空飞行器风速风向实时测量装置，等待无刷电机驱动转速恒定；

步骤S2：取测量周期为整数个转动周期，在测量周期内微处理器U6对压差传感器的压差Δp和无刷电机转角θ进行m次采样；

步骤S3：微处理器对附件进行采样，得到大气压强P和大气温度T，确定大气密度为 $\mathrm{\ρ}=\frac{P}{\mathrm{RT}};$

步骤S4：微处理器采用风速风向提取算法对采样点(Δp_i，θ_i)进行处理，i＝1，2，…，m，提取出相对风速大小V_w和相对风向角ψ_w信息，或者体坐标系分量V_wx，V_wy；

步骤S5：微处理器将相对风速V_w、ψ_w或者V_wx、V_wy送至飞行器飞控计算机，由飞控计算机将相对风速和导航系统给出的地速信息融合，提取出真风速；

步骤S6：返回步骤S1，解算下一输出周期的风速、风向信息。

(三)有益效果

本发明给出的方法和装置可在高空稀薄大气环境下，实现相对风速风向(空速/侧滑角)的实时、准确测量，和导航系统给出的地速信息融合后，可推算出准确的真风速风向信息。该发明具有以下优点：

1、能同时给出风速、风向信息；

2、测量精度高，可实现对1m/s及以下风速、风向的准确测量；

3、测量信号与风速大小呈线性关系，传感器量程不会因最大可测风速的增大而急剧增大；

4、相对风向的测量精度不受大气密度的影响；

5、易于实现、成本低廉。

附图说明

图1为本发明装置Q示意图，图中的椭圆代表测压探头的运动轨迹；

图2为本发明装置Q在高空驻空飞行器Q’上的安装示意图；

图3为本发明装置工作原理说明示意图；

图4为本发明装置本体坐标系、刚性旋杆转角和体系风向角的定义图；

图5a为本发明装置工作时压差传感器和转角编码器输出的压差-时间曲线；

图5b为本发明装置工作时压差传感器和转角编码器输出转角-时间曲线；

图6为风速/风向输出周期、转动周期和压差采样周期示意图。

图7a为1m/s风速和图7b12m/s风速时，4种风速风向提取算法给出的风速测量结果的比较的仿真曲线；

图8为图7风速测量结果对应的95％置信区间半长；

图9为1m/s风速和12m/s风速时，4种风速风向提取算法给出的风向测量结果；

图10为图9风向测量结果对应的95％置信区间半长。

具体实施方式

下面结合附图详细说明本发明技术方案中所涉及的各个细节问题。应指出的是，所描述的实施例仅旨在便于对本发明的理解，而对其不起任何限定作用。

本发明测量装置Q如附图1所示，图2为本发明测量装置Q在高空驻空飞行器Q’上的安装示意图；

本发明装置Q由主装置和附件构成。主装置包括：压差传感器U1、测压探头U2包括：第一测压探头U2A和第二测压探头U2B、刚性旋杆U3、传感器安装盘U4、无刷电机U5及其控制器和微处理器U6组成。附件U7包括气压计和温度传感器，用于推算大气密度。若飞行器已经有大气数据传感器，则附件U7可不要，直接由大气数据传感器提供大气密度信息。刚性旋杆U3对称安装在无刷电机U5的转动轴上；第一测压探头U2A和第二测压探头U2B分别安装在刚性旋杆U3的长边窄面的两端，且第一测压探头U2A和第二测压探头U2B与刚性旋杆U3旋转平面平行，且第一测压探头U2A、第二测压探头U2B分别与刚性旋杆U3的长边窄面垂直，第一测压探头U2A和第二测压探头U2B到无刷电机U5转动轴的距离相等；传感器安装盘U4安装在刚性旋杆U3上；压差传感器U1安装在传感器安装盘U4上，压差传感器U1与第一测压探头U2A和第二测压探头U2B连接；无刷电机U5带动刚性旋杆U3以恒定角速率ω转动；无刷电机U5上设有绝对式旋转编码器，绝对式旋转编码器用于输出刚性旋杆U3的各个时刻转动角；当风吹向转动的刚性旋杆时，第一测压探头U2A和第二测压探头U2B之间产生并输出压差；附件U7包含气压计和温度计，用于测量大气压强和大气温度；微处理器U6根据附件U7给出的气压强和大气温度信号计算出大气密度；微处理器U6与压差传感器U1连接，压差传感器U1输出的压差信号给入微处理器U6，无刷电机U5与微处理器U6连接，微处理器U6对无刷电机U5输出的刚性旋杆U3转动角信号进行高速采样；微处理器U6根据大气密度、压差信号和转动角信号并采用风速风向提取算法得到飞行器相对大气的相对风速、风向信息，该相对风速、风向信息与导航系统给出的地速信息融合，得到真风速风向信息。

压差传感器U1为双极型压差传感器或两只单极型压差传感器组成，第一测压探头U2A和第二测压探头U2B是皮托管，所述的这些部件仅仅是实施例的一种，其他实施例的具体形式在此不再赘述。

第一测压探头U2A、第二测压探头U2B分别与刚性旋杆U3垂直，第一测压探头U2A、第二测压探头U2B到无刷电机U5转动轴的距离相等，称该距离为刚性旋杆U1的旋臂半径，用符号r表示。装置工作时，无刷电机U5带动刚性旋杆U3以恒定角速率ω转动。无刷电机U5配有绝对式旋转编码器，既用于无刷电机U5的转速控制，也用于输出各个时刻，刚性旋杆U3的转动角。微处理器U6对双极型压差传感器U1进行高速采样，利用实施方法一节中的给出的算法，即可提取出相对风速、风向信息。该信息和导航系统给出的地速信息融合后，又可得到真风速风向信息。

方法原理说明

下面结合附图3、4来阐述该发明的工作原理，并借以说明该发明为何能在极低密度的大气环境下，实现对风速风向的准确测量。

当刚性旋杆U3以恒定旋转角速率ω转动时，第一测压探头U2A和第二测压探头U2B相对转动轴的线速度V为：

V＝ωr                                        (3)

其中，r为刚性旋杆U1的旋臂半径。假定第一测压探头U2A、第二测压探头U2B正好位于附图3所示的与风向平行的位置，此时第一测压探头U2A和第二测压探头U2B相对于大气的速度V_A和V_B分别为：

$\left\{\begin{array}{c}{V}_A=V+V_w=\mathrm{\ωr}+V_w\\ V_B=V-V_w=\mathrm{\ωr}-V_w\end{array}\right.---\left(4\right)$

其中，风速大小为V_w；根据伯努利方程，第一测压探头U2A、第二测压探头U2B在A、B处的动压p_A和p_B分别为：

$\left\{\begin{array}{c}{p}_A=\frac{1}{2}\mathrm{\ρ}{V}_A^2=\frac{1}{2}\mathrm{\ρ}[V^2+V_w^2+2{\mathrm{VV}}_w]\\ p_B=\frac{1}{2}\mathrm{\ρ}{V}_B^2=\frac{1}{2}\mathrm{\ρ}[V^2+V_w^2-2{\mathrm{VV}}_w]\end{array}\right.---\left(5\right)$

其中，ρ是大气密度，所以对压差传感器U1所测的压差为：

$\mathrm{\Δp}=p_A-p_B=\frac{1}{2}\mathrm{\ρ}(V_A^2-V_B^2)=2\mathrm{\ρ}{\mathrm{VV}}_w=k\left(\mathrm{\ρ}\right)V_w---\left(6\right)$

其中，k(ρ)＝2ρV                            (7)

式(6)表明由风wind引起的压差测量值与风速大小V_w成正比，比例系数为大气密度ρ和第一测压探头U2A、第二测压探头U2B线速度V乘积的两倍。

若取无刷电机U5的转速为n＝2400rpm＝40rps，并取刚性旋杆U1的旋臂半径r＝30cm，则线速度V为：

V＝ωr＝2πnr＝24π＝75.4m/s                 (8)

在20km处，大气密度ρ＝0.08891kg/m³，由此可计算出

k(ρ)＝2ρV＝13.4kg/[m²s]

也即，1m/s的风速也会引起13.4Pa的压差。这种压差量级，常规级别的压差传感器U1即可准确测量，无需采用高精度、高分辨率的压差传感器U1。

以上就是该方法测速的核心原理，虽然平流层高度大气密度很低，但由于第一测压探头U2A和第二测压探头U2B相对转动轴的线速度V＝ωr的补偿作用，即使是很低的空速，也可产生足够大的压差信号。

关于旋臂半径r和旋转角速度ω的选择，关键要考虑以下因素：

1)若第一测压探头U2A和第二测压探头U2B处线速度V过大，和当地音速之比超过0.3，则伯努利方程无法使用，需要考虑气体的压缩效应，这必将造成数据分析和处理方面的困难；

2)若旋臂半径r太大，则转动引起的动载和气动弹性效应可能会导致刚性旋杆U3、第一测压探头U2A和第二测压探头U2B激烈抖颤，影响测量精度。

所以旋臂半径r和旋转角速度ω均不宜过大，我们建议r的选择在30～60cm之间，旋转角速度ω在600～3000rpm之间，且约束ω·r＜0.3×295≈88m/s。其中，295为平流层音速。

下面再结合附图4，推导发明装置Q工作时，压差传感器输出信号的具体形式。

首先定义测量装置的本体坐标系oxyz。该坐标系是测量基准，所测风速、风向均表述在该坐标系。该坐标系同时也是所发明装置在飞行器上进行安装的基准。安装时，要求该坐标系各轴与飞行器的本体系各对应轴平行。

本体系的xoy平面与刚性旋杆U3的旋转平面重合，z轴与无刷电机U5轴重合。图中，-y轴方向为描述刚性旋杆U3转动角度θ的基准边。-y轴位置对应无刷电机U5旋转编码器的零角度位置。一旦标记了测角基准边，则x轴、y轴和z轴的指向就确定了。

设风速wind与x轴的夹角为ψ_w，其风速大小为Vw。在某个瞬时t，刚性旋杆U3转动到θ转角的位置。

由速度合成定理，第一测压探头U2A相对于大气的速度

为：测压探头U2A相对于转动轴o点的速度加上o点相对大气的速度

后者即为空速

同样，第二测压探头U2B相对于大气的速度

等于第二测压探头U2B相对于转动轴o点的速度

加上o点相对大气的速度

写成矢量的形式，即

${\stackrel{\→}{V}}_A={\stackrel{\→}{V}}_{\mathrm{Ao}}+{\stackrel{\→}{V}}_o,{\stackrel{\→}{V}}_B={\stackrel{\→}{V}}_{\mathrm{Bo}}+{\stackrel{\→}{V}}_o---\left(9\right)$

将

分解成与刚性旋杆U3相切的分量和沿着刚性旋杆U3的分量。由于测压探头U2与刚性旋杆U3垂直，因此仅与刚性旋杆U3相切的分量可被测压探头U2所敏感，沿刚性旋杆U3的分量不能被敏感。

记与刚性旋杆U3相切的分量为V_AT、V_BT，根据附图4，很容易导出：

$\left\{\begin{array}{c}{V}_{\mathrm{AT}}=\mathrm{\ωr}+V_w\cos ({\mathrm{\ψ}}_w-\mathrm{\θ})\\ V_{\mathrm{BT}}=\mathrm{\ωr}-V_w\cos ({\mathrm{\ψ}}_w-\mathrm{\θ})\end{array}\right.---\left(10\right)$

从而，第一测压探头U2A、第二测压探头U2B所敏感处的动压P_A、P_B分别为：

$\left\{\begin{array}{c}{p}_A=\frac{1}{2}\mathrm{\ρ}{V}_{\mathrm{AT}}^2=\frac{1}{2}\mathrm{\ρ}[V^2+V_w^2{\cos }^2({\mathrm{\ψ}}_w-\mathrm{\θ})+2{\mathrm{VV}}_w\cos ({\mathrm{\ψ}}_w-\mathrm{\θ})]\\ p_B=\frac{1}{2}\mathrm{\ρ}{V}_{\mathrm{BT}}^2=\frac{1}{2}\mathrm{\ρ}[V^2+V_w^2{\cos }^2({\mathrm{\ψ}}_w-\mathrm{\θ})-2{\mathrm{VV}}_w\cos ({\mathrm{\ψ}}_w-\mathrm{\θ})]\end{array}\right.---\left(11\right)$

其中，ω为测量装置刚性旋杆U1恒定转动的角速率；V＝ωr为第一测压探头U2A和第二测压探头U2B相对于转动轴的线速度；r为刚性旋杆U1的旋臂半径；

压差传感器U1最终所测压差Δp的输出为：

Δp＝p_A-p_B＝2ρωrV_wcos(ψ_w-θ)＝k(ρ)V_wcos(ψ_w-θ)    (12)

其中k(ρ)＝2ρV，ρ大气密度；

Δp式中表明由风引起的压差测量值为一周期性余弦信号，该余弦信号的幅值与风速大小成正比，该余弦信号的相位包含了风向角ψ_w的信息；由于幅值比例系数k(ρ)为大气密度ρ、第一测压探头U2A和第二测压探头U2B线速度V之积的2倍，通过合理地选择刚性旋杆转速和旋臂长度，在大气密度很低的情况下也能得到较大的信号量，从而能可靠地提取出风速风向信息。

可见，压差传感器U1的输出为一周期性的正弦曲线，附图5a为发明装置工作时压差传感器和转角编码器输出的压差-时间曲线；图5b为发明装置工作时压差传感器和转角编码器输出转角-时间曲线；压差-时间曲线中，幅值绝对值最大点对应风速大小，该点处的转角对应风向角(或反风向角)。

从式(12)可以看出，该信号的幅值包含了风速大小信息，信号的相位包含了风向角信息。通过专门设计的算法即可从输出信号中提取出风速风向信息，具体将在实施方法一节中介绍。

总实施流程：

将本发明测量装置Q按照附图2所示的方式安装在飞行器Q’上，要求测量装置前方和左右两侧无遮挡。

在高空稀薄大气环境中，为低速、驻空型飞行器提供风速风向测量信息，该方法的步骤包括：

步骤(1)：当发明装置启动，并等待无刷电机U5转速稳定后，按照下述步骤提取出相对风速风向信息；

步骤(2)：按照附图6给出的方式，所述周期性余弦信号是刚性旋杆U1转动的整数个转动周期，以采样周期h对压差传感器U1和无刷电机U5转角进行采样。采样n个转动周期，在测量周期内对压差Δp和无刷电机U5转角θ进行m次采样点，微处理器U6对采样点(Δp_i，θ_i)进行处理，并由风速风向提取算法，得到风速风向信息，i＝1，2，…，m。

步骤(3)：微处理器U6对附件的压强计和温度计进行采样，得到大气压强P和大气温度T，由下式确定出大气密度ρ：

$\mathrm{\ρ}=\frac{P}{\mathrm{RT}}=\frac{P}{287\·T}---\left(13\right)$

其中，R＝287，为空气气体常数。

步骤(4)：微处理器U6采用风速风向提取算法对采样点(Δp_i，θ_i)进行处理，i＝1，2，…，m，利用随后介绍的4种提取算法之一(推荐采用算法2或算法3)，提取出相对风速大小V_w和相对风向角ψ_w信息，或者体坐标系分量V_wx，V_wy。

步骤(5)：微处理器U6将相对风速V_w、ψ_w或者V_wx、V_wy送至飞行器飞控计算机(FCC)，由FCC将相对风速和导航系统给出的地速信息融合，按公式(1)提取出真风速。

步骤(6)：返回步骤(1)，解算下一输出周期的风速、风向信息。

算法1：

若在刚性旋杆U3转动周期T内可对压差进行足够密集的采样，则可采用下述方法确定风速、风向：

如附图5所示，找到周期T内压差的最大值Δp_max和最小值Δp_min及其对应的刚性旋杆转动角来提取风速大小和方向信息，利用下式即可计算出风速V_w：

$V_w=\frac{\mathrm{\Δ}{p}_{\max }-\mathrm{\Δ}{p}_{\min }}{4\mathrm{\ρ\ωr}}=\frac{\left|\mathrm{\Δ}{p}_{\max }\right|+\left|\mathrm{\Δ}{p}_{\min }\right|}{4\mathrm{\ρ\ωr}}---\left(14\right)$

又，当θ＝ψ_w时Δp＝Δp_max，当θ＝ψ_w+π时Δp＝Δp_min，所以风向角ψ_w可用下式确定：

ψ_w＝(θ_p++θ_p--π)/2                             (15)

式中，θ_p+和θ_p-分别对应压差值为Δp_max和Δp_min时的旋杆转动角(具体见附图5a和附图5b)，当风速风向输出周期T_wm内包含n个转动周期时，用式(14)、(15)分别计算出每个转动周期T内的(V_wi，ψ_wi)，i＝1，2，…，n。用下式取平均后，作为风速和风向的最终输出值，即：

$V_w=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{V}_{\mathrm{wi}},{\mathrm{\ψ}}_w=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\ψ}}_{\mathrm{wi}}---\left(16\right)$

式中，V_wi为第i个转动周期T内的风速大小，ψ_wi为第i个转动周期T内的风向角；

算法1原理简明、计算简单，并且式(14)巧妙地消除了压差传感器U1可能存在的直流漂移，但要达到一定的精度需要在转动周期T内进行足够次数的采样。以风向的测量为例，要达到2°以下的测向精度，T内采样数至少要大于360/4＝90次。

该方法另一个潜在的问题是所利用的信息较少，这使其容易受到测量野值的影响。虽然有这些不足，但由于简明扼要，且便于人工判读，该方法适合在地面测试和飞行数据分析等后处理等环节中应用。

算法2：

算法2采用傅立叶级数逼近法计算风速、风向，是将压差Δp信号按傅里叶级数展开至一阶项，一阶傅里叶系数与风速在体坐标系上的分量成正比，利用这一特征计算风速风向信息；

首先将式(12)改写成下述形式

Δp＝2ρωrV_wcosψ_wcosθ+2ρωrV_wsinψ_wsinθ      (17)

上式是理想的情况，实际压差信号

必然包含直流漂移、高频毛刺等测量噪声。但

仍为周期信号，所以可采用傅立叶级数展开

$\mathrm{\Δ}\tilde{p}{=A}_0+\underset{k=1}{\overset{\∞}{\mathrm{\Σ}}}(A_k\cos \mathrm{k\ωt}+B_k\sin \mathrm{k\ωt}),t\∈(0,T]---\left(18\right)$

其中，A_k，B_k为傅立叶级数系数。

对比式(17)和(18)，显然只需计算A₀，A₁，B₁即可。由傅立叶系数计算公式：

$A_0=\frac{1}{2\mathrm{\π}}{\∫}_0^{2\mathrm{\π}}\mathrm{\Δ}\tilde{p}\left(t\right)\mathrm{d\θ}=\frac{1}{T}{\∫}_0^T\mathrm{\Δ}\tilde{p}\left(t\right)\mathrm{dt}\≈\frac{1}{m}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\mathrm{\Δ}{\tilde{p}}_i---\left(19\right)$

$A_1=\frac{1}{\mathrm{\π}}{\∫}_0^{2\mathrm{\π}}\mathrm{\Δ}\tilde{p}\left(t\right)\cos \mathrm{\θd\θ}=\frac{2}{T}{\∫}_0^T\mathrm{\Δ}\tilde{p}\left(t\right)\cos \mathrm{\θdt}\≈\frac{2}{m}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\mathrm{\Δ}{\tilde{p}}_i\cos {\mathrm{\θ}}_i---\left(20\right)$

$B_1=\frac{1}{\mathrm{\π}}{\∫}_0^{2\mathrm{\π}}\mathrm{\Δ}\tilde{p}\left(t\right)\sin \mathrm{\θd\θ}=\frac{2}{T}{\∫}_0^T\mathrm{\Δ}\tilde{p}\left(t\right)\sin \mathrm{\θdt}\≈\frac{2}{m}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\mathrm{\Δ}{\tilde{p}}_i\sin {\mathrm{\θ}}_i---\left(21\right)$

其中，m为风速风向输出周期T_wm内的总采样数。T_wm必须是转动周期T的整数倍，即附图6所示的形式：T_wm＝nT＝mh。但m不必被n整除。

求出傅立叶系数A₀、A₁和B₁后，风速、风向可采用下式计算

$V_w=\frac{\sqrt{A_1^2+B_1^2}}{2\mathrm{\ρ\ωr}},{\mathrm{\ψ}}_w={\tan }^{-1}(B_1/A_1)---\left(22\right)$

风向角ψ_w的具体象限由傅立叶系数A₁，B₁的符号确定。风速V_w在体坐标系上的分量V_wx为：

$V_{\mathrm{wx}}=\frac{A_1}{2\mathrm{\ρ\ωr}},V_{\mathrm{wy}}=\frac{B_1}{2\mathrm{\ρ\ωr}}---\left(23\right)$

对于飞行控制来说，采用体系分量实施控制更为方便。

计算表明，只要转动周期T有10个以上的采样点，即m/n≥10，就可取得满意的估计精度。

算法3：

算法3采用广义最小二乘拟合计算风速、风向。

假定在风速风向输出周期T_wm内对压差传感器U1进行了m次采样。取{1，sinθ，cosθ}为基函数，由m个采样点(Δp_i，θ_i)数据可构造如下的广义最小二乘方程：

$\left[\begin{array}{ccc}1& \cos {\mathrm{\θ}}_1& \sin {\mathrm{\θ}}_1\\ 1& \cos {\mathrm{\θ}}_2& \sin {\mathrm{\θ}}_2\\ \·& \·& \·\\ \·& \·& \·\\ \·& \·& \·\\ 1& \cos {\mathrm{\θ}}_m& \sin {\mathrm{\θ}}_m\end{array}\right]\left[\begin{array}{c}{a}_0\\ a_1\\ b_1\end{array}\right]=\left[\begin{array}{c}\mathrm{\Δ}{\tilde{p}}_1\\ \mathrm{\Δ}{\tilde{p}}_2\\ \·\\ \·\\ \·\\ \mathrm{\Δ}{\tilde{p}}_m\end{array}\right]---\left(24\right)$

式中，θm是第m个采样点的转角，a₀，a₁，b₁为广义最小二乘拟合系数，

是第m个采样点的实际压差信号。

该方程的最小二乘解为：

$\left[\begin{array}{c}{a}_0\\ a_1\\ b_1\end{array}\right]={\left[\begin{array}{ccc}m& \mathrm{\Σ}\cos {\mathrm{\θ}}_i& \mathrm{\Σ}\sin {\mathrm{\θ}}_i\\ \mathrm{\Σ}\cos {\mathrm{\θ}}_i& \mathrm{\Σ}\cos {\mathrm{\θ}}_i^2& \mathrm{\Σ}\sin {\mathrm{\θ}}_i\cos {\mathrm{\θ}}_i\\ \mathrm{\Σ}\sin {\mathrm{\θ}}_i& \mathrm{\Σ}\sin {\mathrm{\θ}}_i\cos {\mathrm{\θ}}_i& \mathrm{\Σ}\sin {\mathrm{\θ}}_i^2\end{array}\right]}^{-1}\left[\begin{array}{c}\mathrm{\Σ\Δ}{\tilde{p}}_i\\ \mathrm{\Σ\Δ}{\tilde{p}}_i\cos {\mathrm{\θ}}_i\\ \mathrm{\Σ\Δ}{\tilde{p}}_i\sin {\mathrm{\θ}}_i\end{array}\right]---\left(25\right)$

求出a₀，a₁，b₁后，风速、风向可采用下式计算

$V_w=\frac{\sqrt{a_1^2+b_1^2}}{2\mathrm{\ρ\ωr}},{\mathrm{\ψ}}_w={\tan }^{-1}(b_1/a_1)---\left(26\right)$

ψ_w的具体象限由a₁，b₁的符号确定。风速在体坐标系上的分量V_wx、V_wx如下表示为：

$V_{\mathrm{wx}}=\frac{a_1}{2\mathrm{\ρ\ωr}},V_{\mathrm{wx}}=\frac{b_1}{2\mathrm{\ρ\ωr}}---\left(27\right)$

通常只要转动周期T内的采样点数m/n＞10，即可获得较好的估计精度。

和算法2不同，算法3并不要求采样时间窗T_wm为T的整数倍，m次采样可分布在任意的[0，t]时间段内，甚至采样间隔也可以是非均匀的。所以算法3不但精度高，还具有很大的灵活性，但计算量较算法2大。

算法4：

算法4采用快速傅立叶变换(FFT)提取风速、风向。利用快速傅立叶变换变换来估计风速大小和方向；从快速傅立叶变换谱线中主频谱所对应的频谱参数提取风速风向信息。取m为2的幂次方，利用FFT程序计算：

$P_{k+1}=\underset{i=0}{\overset{m-1}{\mathrm{\Σ}}}\mathrm{\Δ}{\tilde{p}}_{i+1}\exp (-j\frac{2\mathrm{\πk}}{m}i),k=\mathrm{0,1},...,m-1---\left(28\right)$

记P_max＝|P_J|＝max{|P_k+1|}，k＝0，1，…，m-1     (29)

其中，P_k+1为FFT频谱参数，P_J为主频谱参数，P_max为主频谱参数的幅值。

风速、风向采用下式计算：

$V_w=\frac{P_{\max }}{\mathrm{m\ρ\ωr}}---\left(30\right)$

${\mathrm{\ψ}}_w={\tan }^{-1}\left[\frac{-\mathrm{imag}\left(P_J\right)}{\mathrm{real}\left(P_J\right)}\right]---\left(31\right)$

风速在体坐标系上的分量为：

$V_{\mathrm{wx}}=\frac{\mathrm{real}\left(P_J\right)}{\mathrm{m\ρ\ωr}},V_{\mathrm{wy}}=\frac{-\mathrm{imag}\left(P_J\right)}{\mathrm{m\ρ\ωr}}---\left(32\right)$

其中，imag(·)和real(·)分别表示复数P_J的虚部和实部。

采用FFT方法存在的一个问题是，由于

为周期信号，因此采样时间窗必须是转动周期T的整数倍，即T_wm满足附图6所示形式：

T_wm＝nT＝mh                                       (33)

否则对应的FFT频谱参数会发生畸变，无法准确复现原信号。不过这样一来，确保m为2的幂次方可能会比较棘手。

FFT所计算的离散频谱参数中相邻谱线的频率间隔为：

$\mathrm{\Δf}=\frac{1}{T_{\mathrm{wm}}}=\frac{1}{\mathrm{mh}}=\frac{f_p}{m}---\left(34\right)$

其中，f_p为对压差传感器U1的采样频率。故由式(29)所计算的主频f_m为：

$f_m=\mathrm{J\Δf}=J\frac{f_p}{m}---\left(35\right)$

理论上2πf_m应等于无刷电机U5转速ω，故通过计算二者间的误差ε，可评估转速ω的稳定性，前提是相邻谱线的频率间隔Δf足够小。

$\mathrm{\ϵ}=|2\mathrm{\π}{f}_m-\mathrm{\ω}|=|2\mathrm{\πJ}\frac{f_p}{m}-\mathrm{\ω}|---\left(36\right)$

FFT方法并未显式利用压差信号主频ω已知这一特点，这使得其对风速、风向的估计精度不如算法2和算法3。但不依赖转速ω某些时候也是一个优点，尤其在系统测试和性能评估时非常有用。FFT方法的另一个优点是有非常成熟的模块可用。

图7a～图10b为实施方法中给出的4种风速风向提取算法比较的仿真曲线。仿真条件为：刚性旋杆U3转速2400rpm，旋臂半径r的长30cm，压差传感器U1量程±500Pa，压差传感器U1精度为满量程的2.5％(3σ)。待测风速2组：风速1m/s，风向角30°；风速12m/s，风向角30°，其中：

图7a为1m/s风速和图7b12m/s风速时，4种风速风向提取算法给出的风速测量结果；其中圆点标记曲线为真实值，三角标记曲线为算法1测量仿真结果，圆形标记曲线为算法2测量仿真结果，矩形标记为算法3测量仿真结果，菱形标记曲线为算法4测量仿真结果。

图8a和图8b为图7a和图7b风速测量结果对应的95％置信区间半长；其中圆点标记曲线为零值，三角标记曲线为算法1测量仿真结果，圆形标记曲线为算法2测量仿真结果，矩形标记为算法3测量仿真结果，菱形标记曲线为算法4测量仿真结果。

图9a为1m/s风速和图9b12m/s风速时，4种风速风向提取算法给出的风向测量结果；其中圆点标记曲线为真实值，三角标记曲线为算法1测量仿真结果，圆形标记曲线为算法2测量仿真结果，矩形标记为算法3测量仿真结果，菱形标记曲线为算法4测量仿真结果。

图10a和图10b为图9a和图9b风向测量结果对应的95％置信区间半长。其中圆点标记曲线为零值，三角标记曲线为算法1测量仿真结果，圆形标记曲线为算法2测量仿真结果，矩形标记为算法3测量仿真结果，菱形标记曲线为算法4测量仿真结果。

从上述仿真曲线中可以看出算法2和算法3效果最好，算法1较差，算法4精度介于算法1和算法2、3之间，但算法4在若干频率点上精度存在跳跃，这主要是采样时间窗不满足式(33)所造成。当采样频率＞400Hz后，采样频率对算法2和算法3的精度影响较小，但对算法1和算法4影响较大。算法1的精度随采样频率的提高趋于稳定。采样频率对算法4精度的影响主要体现在采样时间窗上，当采样时间窗不满足式(33)时，精度呈跳跃状。

以上所述，仅为本发明中的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉该技术的人在本发明所揭露的技术范围内，可理解想到的变换或替换，都应涵盖在本发明的包含范围之内，因此，本发明的保护范围应该以权利要求书的保护范围为准。

Claims (8)

1.一种高空驻空飞行器风速风向实时测量装置，其特征在于：在高空稀薄大气环境中，为低速、驻空型飞行器提供风速风向测量信息，该测量装置包括：压差传感器(U1)，第一测压探头(U2A)、第二测压探头(U2B)，刚性旋杆(U3)、传感器安装盘(U4)，无刷电机及其控制器(U5)、微处理器(U6)和附件(U7)组成；

刚性旋杆(U3)对称安装在无刷电机(U5)的转动轴上；第一测压探头(U2A)和第二测压探头(U2B)分别安装在刚性旋杆(U3)的长边窄面两端，且第一测压探头(U2A)和第二测压探头(U2B)与刚性旋杆(U3)旋转平面平行，且第一测压探头(U2A)、第二测压探头(U2B)分别与刚性旋杆(U3)的长边窄面垂直，第一测压探头(U2A)和第二测压探头(U2B)到无刷电机(U5)转动轴的距离相等；传感器安装盘(U4)安装在刚性旋杆(U3)上；压差传感器(U1)安装在传感器安装盘(U4)上，压差传感器(U1)与第一测压探头(U2A)和第二测压探头(U2B)连接；无刷电机(U5)带动刚性旋杆(U3)以恒定角速率ω转动；无刷电机(U5)上设有绝对式旋转编码器，绝对式旋转编码器用于输出刚性旋杆(U3)的各个时刻转动角；当风吹向转动的刚性旋杆时，第一测压探头(U2A)和第二测压探头(U2B)之间产生并输出压差；附件(U7)包含气压计和温度计，用于测量大气压强和大气温度；微处理器(U6)根据附件(U7)给出的气压强和大气温度信号计算出大气密度；微处理器(U6)与压差传感器(U1)连接，压差传感器(U1)输出的压差信号给入微处理器(U6)，无刷电机(U5)与微处理器(U6)连接，微处理器(U6)对无刷电机(U5)输出的刚性旋杆(U3)转动角信号进行高速采样；微处理器(U6)根据大气密度、压差信号和转动角信号并采用风速风向提取算法得到飞行器相对大气的相对风速、风向信息，该相对风速、风向信息与导航系统给出的地速信息融合，得到真风速风向信息。

2.如权利要求1所述高空驻空飞行器风速风向实时测量装置，其特征在于，当风相对于飞行器的以速度V_w、方向ψ_w吹向测量装置时，根据伯努利方程，第一测压探头(U2A)、第二通测压探头(U2B)位于处的动压P_A、P_B分别为：

$\left\{\begin{array}{c}{p}_A=\frac{1}{2}\mathrm{\ρ}[V^2+V_w^2{\cos }^2({\mathrm{\ψ}}_w-\mathrm{\θ})+2{\mathrm{VV}}_w\cos ({\mathrm{\ψ}}_w-\mathrm{\θ})]\\ p_B=\frac{1}{2}\mathrm{\ρ}[V^2+V_w^2{\cos }^2({\mathrm{\ψ}}_w-\mathrm{\θ})-{2\mathrm{VV}}_w\cos ({\mathrm{\ψ}}_w-\mathrm{\θ})]\end{array}\right.$

其中，ω为测量装置刚性旋杆(U1)恒定转动的角速率；V＝ωr为测压探头(U2)相对于转动轴的线速度；r为刚性旋杆(U1)的旋臂半径；

所以压差传感器(U1)所测压差Δp为：

Δp＝p_A-p_B＝2ρωrV_wcos(ψ_w-θ)＝k(ρ)V_wcos(ψ_w-θ)

其中k(ρ)＝2ρV，ρ大气密度；

Δp式中表明由风引起的压差测量值为一周期性余弦信号，该余弦信号的幅值与风速大小成正比，该余弦信号的相位包含了风向角ψ_w的信息；由于幅值比例系数k(ρ)为大气密度ρ、第一测压探头(U2A)和第二测压探头(U2B)线速度V之积的2倍，通过合理地选择刚性旋杆(U1)转速和旋臂长度，在大气密度很低的情况下也能得到较大的信号量，从而能可靠地提取出风速风向信息。

3.如权利要求2所述高空驻空飞行器风速风向实时测量装置，其特征在于，所述旋臂半径r选择在30～60cm之间，旋转角速度ω在600～3000rpm之间，且约束ω·r＜0.3×295≈88m/s，其中，295为平流层音速。

4.如权利要求2所述高空驻空飞行器风速风向实时测量装置，其特征在于，所述周期性余弦信号是刚性旋杆(U1)转动的整数个转动周期，在测量周期内对压差Δp和无刷电机(U5)转角θ进行m次采样；微处理器(U6)对采样点(Δp_i，θ_i)进行处理，并由风速风向提取算法，得到风速风向信息，i＝1，2，…，m。

5.如权利要求1所述高空驻空飞行器风速风向实时测量装置，其特征在于，所述测量装置的本体坐标系oxyz是测量基准，所测风速、风向均表述在该坐标系；该坐标系也是所述测量装置在飞行器上进行安装的基准，安装时，要求该坐标系各轴与飞行器的本体系各对应轴平行。

6.如权利要求1所述高空驻空飞行器风速风向实时测量装置，其特征在于，所述本体坐标系oxyz的本体系的xoy平面与刚性旋杆(U3)的旋转平面重合，z轴坐标与无刷电机(U5)轴重合，-y轴坐标方向为描述刚性旋杆(U3)转动角度θ的基准边；-y轴坐标位置对应无刷电机(U5)绝对式旋转编码器的零角度位置；标记了测角基准边，则确定x轴坐标、y轴坐标和z轴坐标的指向。

7.一种如权利要求1所述装置的高空驻空飞行器风速风向实时测量方法，其特征在于：在高空稀薄大气环境中，为低速、驻空型飞行器提供风速风向测量信息，该方法的步骤包括：

步骤S1：启动高空驻空飞行器风速风向实时测量装置，等待无刷电机驱动转速恒定；

步骤S2：取测量周期为整数个转动周期，在测量周期内微处理器U6对压差传感器U1的压差Δp和无刷电机U5转角θ进行m次采样；

步骤S3：微处理器对附件进行采样，得到大气压强P和大气温度T，确定大气密度为 $\mathrm{\ρ}=\frac{P}{\mathrm{RT}};$

步骤S4：微处理器采用风速风向提取算法对采样点(Δp_i，θ_i)进行处理，i＝1，2，…，m，提取出相对风速大小V_w和相对风向角ψ_w信息，或者体坐标系分量V_wx，V_wy；

步骤S5：微处理器将相对风速V_w、ψ_w或者V_wx、V_wy送至飞行器飞控计算机，由飞控计算机将相对风速和导航系统给出的地速信息融合，提取出真风速；

步骤S6：返回步骤S1，解算下一输出周期的风速、风向信息。

8.如权利要求7高空驻空飞行器风速风向实时测量方法，其特征在于：所述微处理器采用风速风向提取算法包括四种算法，该四种算法的特征分别是：

(1)算法1利用转动周期内的压差的最大值和最小值及其对应的刚性旋杆转动角来提取风速大小和方向信息；

(2)算法2将压差信号按傅里叶级数展开至一阶项，一阶傅里叶系数与风速在体坐标系上的分量成正比，利用这一特征计算风速风向信息；

(3)算法3取{1，sinθ，cosθ)为基函数，由m个采样点(Δp_i，θ_i)构造广义最小二乘方程；从解算该方程所得到的基函数拟合系数提取风速风向信息；

(4)算法4利用快速傅立叶变换变换来估计风速大小和方向；从快速傅立叶变换谱线中主频谱所对应的频谱参数提取风速风向信息。

Images