多经纬仪跟踪探空气球的空中风场反演方法

周 洁,高 航,刘金波,杨伟伟,李健兵

(1.国防科技大学电子科学学院,湖南 长沙 410073；2.海军91236部队,辽宁 葫芦岛 125100)

1 引 言

风是日常生活中最常见的一种自然现象,它无时不在,无处不有,并且对人类生产生活产生巨大的影响,因此对于风场的精确感知意义重大。在航空气象领域,准确掌握精细的三维风场信息,能为飞机安全起降提供重要支撑。目前,用于风场探测技术主要包括探空氢气球、测风仪、风廓线雷达以及激光雷达,其中探空氢气球具有携行方便、操作简单、探测便捷等优点,成为高空气象探测和其他方式真值比对的主要手段之一,在国内外气象行业得到了广泛运用[1-4]。

探空氢气球探测方式是利用经纬仪连续跟踪探空氢气球获得风场信息,目前经纬仪实施空中风探测的主要方法有单经纬仪测风和双经纬仪基线测风方法。其中,单经纬仪测风通过获取特定时间内的探空氢气球的仰角和方位角数据,经过一定算法解算得到相应高度的空中风数据,使用方便灵活、准备时间短、操作人员少,便于在各种地形条件下完成作业任务,在气象野外、伴随保障中仍然发挥着较大作用[5-6]。但单经纬仪测风,需设定恒定气球升速的假设,以此来进行风速的反演,因此测量结果受各项误差源影响较大,准确度稍差[7-8]。双经纬仪基线测风同时使用两台光学测风经纬仪观测氢气球的运动,读出仰角、方位角然后利用三角法或矢量法计算氢气球高度和风向、风速[9]。双经纬仪基线测风方法能够克服恒定升速假设这一问题,有效提高数据的准确度和精度,但其对基线长度和方向的选取直接影响风向风速计算的准确度[10]。

此外,无论是单经纬仪还是双经纬仪在计算空中风场信息时,均预先假定氢气球为风场的良好示踪物,其运动状态与背景风场一致。但氢气球虽然质量较小,但体积较大,在速度变化较大时,较显著的存在氢气球的运动状态与实际风场不一致现象,简单地将二者近似相等与实际不符,因此应当采用一定的反演算法从氢气球运动状态中反演出真实的背景风场数据。本文提出一种基于多经纬仪的风场探测与反演方法(MTLM),通过一定的坐标变换与最小二乘法方法,获取氢气球的精确轨迹信息。随后通过氢气球运动状态方程,基于曳力表达式,实现了由氢气球速度对于背景风场速度的反演。经过仿真与外场实测实验,验证了本方法对于反演背景风场具有较高的性能。

2 基本原理

风传统氢气球观测方法一般基于两个前提假设条件:①气球在整个观测过程中保持恒定升速;②氢气球是背景风场的良好示踪物。但在实际中,这些假设并不能够成立,特别是在复杂风场情况下,可能会导致较大的观测误差。为解决上述问题,本文提出一种基于多经纬仪观测探空氢气球的风场反演方法,该方法主要涉及两个创新点:①多个经纬仪协同观测探空氢气球并通过坐标变化及最小二乘法拟合获得其轨迹和速度;②考虑氢气球的惯性作用,采用运动方程分析准确地反演背景风。

在介绍新方法之前,首先分析传统方法的原理和局限。

2.1 氢气球受力分析

单经纬仪观测探空氢气球方式如图1所示,以经纬仪位置为坐标原点,建立笛卡尔坐标系,X轴以正东为正方向;Y轴以正北为正方向;Z轴以垂直向上为正方向;θ为以正北顺时针方向为正方向的方位角;φ为以垂直向上为正方向的高度角。氢气球在自由大气中做上升运动时,将受到由于空气与氢气间密度差所产生的浮力Ff,氢气球自身及附属物和氢气的总重力G,以及曳力fd。忽略氢气球体积和空气密度的变化,氢气球所受浮力和重力保持不变,曳力大小与氢气球和空气相对速度成正比。因此当氢气球在初始浮力作用下,具有了竖直向上的加速度,曳力随着速度的增加而增大,当浮力、重力和曳力达到平衡状态时,氢气球此时处于匀速上升状态,将保持固定升速上升。

图1 单经纬仪观测氢气球示意图Fig.1 Schematic diagram of hydrogen balloon observation with single theodolite

2.2 空中平均风计算

(1)

(2)

(3)

(4)

2.3 氢气球作为风场示踪物

实际使用单经纬仪进行空中风观测时,使用氢气球作为风场的示踪物,需忽略氢气球本身的惯性和曳力,并且需满足固定的升速vup,但这种假设,实际中并不成立,而惯性则是决定气球的运动状态中的重要参数。惯性是表示流体对物体运动影响的关键参数,通常用斯托克斯数来评估惯性[11]。根据斯托克斯数公式St=t0u0/L0,其中t0=ρdD2/18μ为物体的弛豫时间,u0为流体通过物体时的流速,L0为物体的特征尺寸,ρd为物体密度,D为物体直径,μ为气体的黏度。如果St≪1,惯性则可以忽略；否则,惯性强度很大,物体就不能成为背景风场的良好示踪物。因此在标准气压101.325 kPa、温度为25℃的条件下,空气密度为1.169 kg/m3,粘度为1.8448×10-5Pa·s,氢气密度0.081 kg/m3,直径为0.5 m、运动速度为2 m/s的氢气球斯托克斯数约为243,远大于1,所以实际情况下氢气球的惯性不能再忽略。在氢气球飞行过程中,在背景风场的作用下必然会产生曳力(空气阻力),导致受力的不平衡,无法保持固定的升速vup。若此时仍按单经纬仪的方式进行观测,所得到的空中风场将与真实背景风场具有一定的误差。

3 多经纬仪风场探测

单经纬仪风场探测的局限主要有两点:一是假定氢气球的升速vup固定,但实际升速无法保持固定,因此获得氢气球的高度值存在较大误差；二是相对运动的存在,使得经纬仪观测的氢气球速度不能够直接代表风场的情况。为解决以上问题,本文采取多经纬仪联合观测和背景风场反演方法实现对于风场的探测反演。

3.1 多经纬仪联合观测氢气球速度原理

多经纬仪联合观测的方法放置成如图2所示,经纬仪之间的距离在60～100 m之间均可,夹角控制在60°～90°左右。观测前预先确定相应经纬仪之间的方位角和高度角,以及至少1条边的距离D1,θij表示为第i个经纬仪观测第j个经纬仪的方位角。

图2 三经纬仪观测氢气球示意图Fig.2 Schematic diagram of hydrogen balloon observation with three theodolite

以其中1个经纬仪位置为坐标原点Ο(0,0,0)建立笛卡尔坐标系,(Xm,Ym,Zm)表示为第m个经纬仪的坐标,(xn,yn,zn)表示为第n次观测时氢气球的坐标。根据三角形正弦公式可求得各经纬仪的相对距离,并进一步转化为笛卡尔坐标系下的各经纬仪的坐标值(Xm,Ym,Zm)。

(5)

(6)

当完成多次MTLM观测后,可以得到多个氢气球空中位置,通过多项式拟合方法获得氢气球在观测时间段内的完整运动轨迹S(t),对其关于时间求导即得各个时刻氢气球的速度vb(t)如式(7):

(7)

3.2 背景风场反演原理

在2.1小节中,通过多经纬仪联合观测能够快速准确地获得氢气球的运动参数,根据前文分析,氢气球的运动不能直接代表空中风的情况,因此需采用一定的算法从氢气球的运动信息中提取空中风的信息。如图2所示,根据氢气球受力分析,此时氢气球受重力G、浮力Ff以及曳力fd,根据牛顿第二定律可以得到式(8):

(8)

其中,a为氢气球加速度;m为氢气球自身及附属物和氢气的总质量。空气可看作牛顿流体,因此氢气球的曳力表达式[12]为:

(9)

其中,Cd为氢气球的曳力系数；ρa为空气密度；D为氢气球直径,vw为背景风速；vb为氢气球速度。根据浮力公式Ff=ρagV以及重力公式G=mg,将其代入公式(8)可得:

(10)

其中,V为氢气球体积。

氢气球在低空开始上升时,可以通过控制质量,在竖直方向上使得氢气球处于的平衡状态,此时加速度为a=0,可得到平衡状态下的氢气球上升速度:

(11)

当已知氢气球初始匀速上升的速度时,忽略氢气球体积改变以及空气密度的变化,可以近似认为氢气球的曳力系数在上升过程中保持不变,因此可求得氢气球的曳力系数为:

(12)

将式(9)、式(12)代入式(8)可以得到氢气球运动状态方程:

(13)

根据3.1小节获得的氢气球速度信息,和a=dvb/dt获得加速信息,由上式求出δv的三个分量δvx,δvy,δvz,进一步地,由关系式δv=vw-vb,即vw=vb+δv可以得到背景风速vw。

4 仿真实验

4.1 仿真设置

仿真氢气球在背景风场控制下的移动轨迹,设置仿真背景风场为:

(14)

其中,Uturb代表分布于整个空间中耗散率为5×10-3m2/s3的湍流速度[13]；Ushear代表如图3所示分布于整个空间的非线性风切变的速度。

图3 非线性风切变变化示意图Fig.3 Schematic diagram of nonlinear windshear

仿真氢气球及附属物重70 g,直径为50 cm,假设氢气球为理想,滑球体,曳力系数为0.4[14],标准大气环境(气温0 ℃,1个标准大气压,空气密度1.293 kg/m3),重力加速度为9.8 m/s2。3个经纬仪如图2所示摆放,其中2号、3号经纬仪距1号经纬仪距离均为100 m,相对方位角分别为30°和90°,三个经纬仪坐标分别为(0,0,0),(50,86.6,0),(100,0,0)。根据式(11)获得其初始平衡速度为1.68 m/s。因此氢气球初始速度设为为(2,2,1.68),在背景风场作用下,根据式(13)获得各方向加速度,可通过对加速进行积分获得氢气球的轨迹:

(15)

4.2 误差分析

根据式(6),采用本方法进行观测时需要2个以上经纬仪才能进行观测并求解氢气球位置。考虑方法的适用性、操作性以及经济性因素,通常采用2～4个经纬仪方法进行观测,为进一步验证本方法的鲁棒性和准确性,分别采用2～4个经纬仪进行仿真实验,分别单独引入0.5°的高度角、方位角以及联合观测随机误差进行比对分析,以获得最佳经纬仪数量配置。图4给出了在不同观测误差下,仿真的速度误差比较结果。可以看出,无论是引入哪种误差,增加经纬仪的数量反演的速度误差将会减小,准确度增高。图4(d)、(e)、(f)中可以看出,在只有高度角误差时,增加经纬仪的数量有利于抑制高度角方向上引入的观测误差,但在经纬仪数量大于3个后,速度精度基本保持不变；而图4(g)、(h)、(i)中,只引入了方位角误差时,增加经纬仪的数量,反演得到的速度精度基本保持一致。由此可以看出,增加经纬仪的数量对于观测产生的高度角方向上的误差具有较强的抑制作用,能够有效地消除该方向上的速度波动,但其对方位角上的精度提高作用不明显。经纬仪的数量也不是越多越好,在达到3个经纬仪以后,无论是水平方向上还是垂直方向上的反演精度基本保持不变,因此通常情况下采用3个经纬仪方法即可满足日常观测使用。随后,将以3个经纬仪观测方式进行相应的观测轨迹误差和速度误差进行分析。

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)图4 MTLM方法(2/3/4)观测误差比较,(其中(a)(b)(c) 为存在方位角和高度角观测误差,(d)(e)(f) 为高度角观测误差,(g)(h)(i)为方位角观测误差)Fig.4 Error Comparison of MTPM methods(2/3/4), (where(a)(b)(c)is the error of azimuth and altitude angle, (d)(e)(f)is error of altitude angle,and(g)(h)(i) is the error of azimuth angle)

4.2.1 轨迹误差分析

采用MTLM方法对氢气球进行观测时,在观测方位角和高度角上各加上满足 正态分布的高斯随机误差。在获取3组观测数据后,通过式(6)使用最小二乘法快速计算出氢气球的位置。

如图5所示,其中“观测轨迹”表示加入观测误差后仿真形成的氢气球轨迹,“滤波轨迹”表示采用MTLM方法滤波后获得的氢气球轨迹。根据仿真结果可以看出,采用本方法进行观测时,在复杂风场背景条件下(存在湍流、风切变),定位结果具有较高准确性,各方向均方根误差小于2 m,相对于观测轨迹,各方向精度均提高超过50 %,验证本方法对氢气球轨迹反演具有较高性能。

(a)

(b)

(c)图5 MTLM方法与仿真结果比较Fig.5 Comparison between MTLM method and simulation results

4.2.2 速度误差分析

将MTLM方法观测反演所获得的位置点进行多项式拟合,以获得其空中轨迹的多项式表达形式,再根据式(13)氢气球运动方程分析可分别获得氢气球速度以及背景风场速度。将传统单经纬仪方法((ST),经纬仪观测加多项式拟合方法((TLM),及经MTLM方法并经过运动方程分析(ME)进行对比,具体误差分析如图6所示,其中,图6(a)、(b)、(c)分别为U,V,W风速分量误差结果。由结果可知,通过本方法反演的背景风场与传统单经纬仪相比,水平速度和垂直速度准确度分别提升超过80 %和90 %。此外,经过运动方程分析后,本方法能够有效修正氢气球惯性作用所带来的的偏离背景风场的误差,特别是在垂直方向上,能够更加精准反演出背景风垂直速度分布情况。由于水平方向上,氢气球所受合力仅与曳力大小有关,而曳力大小与相对速度成正比,因此在氢气球与背景风场相对速度较小的情况下,曳力较小,产生的加速度不大,最终氢气球与背景风场速度几乎一致,由图6(a)、(b)可知经反演背景风与氢气球速度相差较小。垂直方向上由于受浮力、重力和曳力的共同作用,其运动状态变化较大,加速度的存在对氢气球运动状态影响较大,因此应通过本方法反演的背景风场代表真实背景风场数据,图5(c))反演背景风精度要远远高于氢气球所代表的风场速度。因此,仿真实验有效验证了本方法在反演背景风场方面的良好性能。此外我们注意到,从图5和图6的速度误差分析中,出现了误差周期性的跳动以及两端的大误差现象,该现象均是由采用的拟合方法所致。本文采用最小二乘法进行拟合,通过最小化总体数据误差的平方和以寻找数据的最佳匹配函数,导致拟合结果与真值会出现“上下波动”现象,这与采用最小二乘法拟合特点吻合。

(a)

(b)

(c)图6 MTLM方法反演氢气球和背景风场误差对比Fig.6 Error comparison of hydrogen balloon and background wind field by MTLM method

5 外场实验验证

为进一步验证本方法观测氢气球风场以及反演背景风场的性能,2021年6月,位于辽宁省葫芦岛市开展了外场多经纬仪联合观测氢气球实验。本次实验采用3部相同的光学经纬仪进行观测,经纬仪布置方式如图2所示,其中2号、3号经纬仪距1号经纬仪距离均为100 m,相对方位角分别为1.7°和62.7°,经纬仪距地面高度均为1.61 m。因此三个经纬仪坐标分别为(0,0,1.61),(2.97,99.96,1.61),(88.86,45.87,1.61)。为验证本方法反演背景风场的准确性,引入激光雷达测量风场数据进行对比。激光雷达与经纬仪观测场距离为50 m,氢气球漂移最大可观测径向距离距1号经纬仪约为800 m,激光雷达最大作用距离为2000 m,其作用距离在可覆盖氢气球观察全过程。

5.1 氢气球速度获取

观测使用的氢气球及附属物重13 g,直径为36 cm,气温28.6 ℃,气压为992.6 hPa,空气密度1.1465 kg/m3,面风速为2.0 m/s,风向270°。采取每30 s同时观测获取1组数据,在经纬仪可见范围内共获得20组数据,如表1所示。使用MTLM方法对氢气球空中轨迹进行定位反演,并经过3阶多项式拟合以获得氢气球空中轨迹表达式,如图7所示。根据多项式拟合获得的曲线,并对获得的曲线通过式(7)计算,以获得各点的速度以及加速度信息。无论是采用单经纬仪还是MTLM方法测量反演所得的风场信息,理论上均为氢气球在背景风场作用下的速度,由于氢气球惯性和曳力的存在,该速度不能代表真实的背景风场速度,因此应当采用3.2小节方法进行反演计算。

(a)X方向位移

(b)Y方向位移

(c)Z方向位移图7 MTLM反演实测数据位移结果Fig.7 MTLM retrieve results of displacement data

5.2 背景风场获取

通过4.1小节获得氢气球各方向的速度信息后,通过代入式(13)求解得到真实背景风场信息,通过与在同一时间段、同一场地进行风场探测的激光雷达的风场信息作对比,相关结果如图7、图8所示。由图8(a)、(b)可以看出,对于水平风速反演来说,采用MTLM方法较单经纬仪性能有较大的提升,风向和风速准确度分别提升86 %和66 %,在高精度风速获取上MTLM方法具有更好的性能。在精度要求不高的条件下,无论是MTLM方法还是单经纬仪方法,都能较准确地得到背景风场信息,两种方法均可以实现对于水平风的快速测量。

(a)水平风向

(b)风速

(c)垂直风速图8 MTLM反演实测数据风速结果Fig.8 MTLM retrieve results of wind speed data

垂直风速的获取结果可由图8(c)表示,MTLM方法不仅可以反演出氢气球垂直速度,也能通过动力学方程反演得出背景真实垂直风速。由图可知,若仍以氢气球速度代表垂直风速,其误差较大幅度增加,通过本方法的到的垂直风速能够较准确地表征真实背景风速信息。

6 结 论

本文提出一种基于多经纬仪的风场探测与反演方法,通过对多经纬仪协同探测获取的氢气球高度角和方位角信息,进行相应的坐标变换与最小二乘法算法,实现了对于氢气球的精确轨迹信息的获取。在此基础上,基于曳力公式对氢气球进行动力学分析,实现了由氢气球速度对于背景风场速度的反演。与传统单经纬仪方法相比,仿真实验验证本方法反演的背景风场在水平速度和垂直速度准确度提升超过80 %和90 %；通过外场实验验证本方法反演的水平风向和风速准确度分别提升86 %和66 %；此外,本方法还能够有效测量垂直风速。经过仿真与外场实测实验,验证了本方法对于反演背景风场具有较好的性能。为实现经纬仪快速、准确、便捷地对高空风场精确测量提供了有效手段。