一种基于惯性系统的探空仪风速测量误差补偿方法

1.本发明属于气象探测技术领域，尤其涉及一种基于惯性系统的探空仪风速测量误差补偿方法。


背景技术：

2.台风是世界上最严重的自然灾害之一。为实现对台风这一自然现象更为本质的科学认知，关键技术瓶颈是对台风进行多要素、长过程、精细化的直接探测，获得台风内核区的温度、湿度、压力、风速、风向等关键气象数据，其中测量台风内核区风速对于预测台风路径与强度尤为重要。目前对于台风风速的测量以风廓线雷达、激光测风雷达、探空仪等为主，而为了满足台风内核区风速的原位探测只能利用探空仪进行测量。但传统的探空仪设备只能依靠对探空仪 gps或北斗系统探测到的位置变化数据解算风速，误差较大。当探空仪处于台风等极端环境中，探空仪本身会在下落过程中绕降落伞急速翻滚、转动，其姿态的剧烈变化会对风速测量造成较大误差影响，此时探空仪的速度不能直接代表实际风场的流速，必须消除由探空仪运动姿态变化所导致的风速测量误差。
3.目前，国内外众多研究高空探测的学者对导航型探空仪风速测量算法做了大量研究。jaakko研究viasala rs92探空仪指出风速测量要进行数据相关性分析和平滑处理。jauhiainen等在总结国际探空仪比对试验中指出，风速测量处理应考虑探空仪摆动因素，但未提及处理的方法。张恩红根据复摆周期的计算方法，研究了探空仪的摆动周期，利用滤波方法去除探空仪摆动对风速测量的影响。然而上述几种方法将探空仪的运动理想化成简单的、具有周期性的锥摆运动，以此来对全部探测数据进行滤波处理，精度较低。当探空仪处在台风这种急剧变化且风速极强的环境中，其运动是复杂且未知的，单纯用周期性的锥摆运动来补偿探测数据将会造成欠补偿或过补偿等问题，无法准确获得探空仪在当前时空的风速风向信息。
4.专利申请cn114019583a，基于惯性补偿的高精度测风系统及方法，公开了一种高精度安全测风系统及方法，通过对惯性模块测得的探空仪加速度积分得到速度，用加权平均的方式对定位模块测得的风速进行补偿；根据惯性模块测得的探空仪姿态信息，解算探空仪摆动周期，并针对不同周期的数据进行多点滤波滤除探空仪摆动的影响。首先，针对第一点，其利用惯性模块测得的探空仪加速度积分得到速度这种做法是不正确且没有意义的。在探空仪系统中，主要存在着降落伞和探空仪两个组件，进行风速补偿的意义主要是补偿探空仪本身围绕降落伞所进行的复杂运动，如果单纯利用惯性模块所得到的加速度信息去解算风速，相当于测量的是整个探空系统的瞬时风速，同gps测量的效果是一样的，并没有解决探空仪绕伞运动所带来的测量误差，无法起到风速补偿的效果。针对第二点，为了解决误差必须利用惯性系统里测得的姿态信息，但该方法仅仅是解算探空仪的摆动周期，这与之前提到的研究文章存在着相同的问题，即利用解算的周期去补偿全部时刻的测风数据会存在过补偿或欠补偿的问题，并不是每一时刻都需要去补偿一致的摆动周期影响，每一时刻都有自己独特的特征，要根据实际情况逐一补偿探空仪测量风速误差。
5.综上，亟需设计一种探空仪风速测量误差补偿方法，补偿因探空仪姿态变化所产生的风速测量误差。


技术实现要素：

6.针对目前探空仪在恶劣工况下因自身姿态变化导致的风速测量误差，本发明提出一种基于惯性系统的探空仪风速测量误差补偿方法。结合探空仪上的惯导信息，解算探空仪真实运动姿态轨迹，补偿因探空仪运动所导致的风速测量误差。本发明的具体技术方案如下：
7.一种基于惯性系统的探空仪风速测量误差补偿方法，风速测量过程中的下投结构包括伞、绳及探空仪，探空仪通过绳与伞底相连，伞为探空仪提供浮空时间，绳在下投过程中保持拉紧状态，所述补偿方法包括以下步骤：
8.s1：在探空仪上集成惯性系统，进行高空下投探测，获取立体空间气象数据，地面接收机实时记录并存储探空仪原位探测数据及三轴惯导信息；
9.s2：利用接收到的探空仪gps/北斗信息，计算探空仪探测到的立体空间风速信息；
10.s3：利用接收到的探空仪三轴惯导信息，计算探空仪下落过程中的姿态信息；
11.s4：将步骤s3的结果通过空间坐标系变换与角度解算，获取探空仪的实时相对空间运动姿态信息；
12.s5：通过步骤s4得到的探空仪的相对空间运动姿态信息，计算探空仪因自身姿态变化导致的风速测量误差，代入步骤s2解算的立体空间风速信息中，完成风速补偿；
13.s6：对经过步骤s5补偿后的风速数据进行滚动时间窗平滑处理，进一步消除环境噪声的影响。
14.进一步地，所述步骤s3中的三轴惯导信息包括三轴加速度、三轴角速度以及三轴角度信息，保留角度范围在(-180
°
,180
°
)的三轴角度信息，通过三次样条插值对三轴角度原始数据插值，再对位移进行差分求解得到运动速度；
15.惯性系统使用欧拉角表示姿态时的坐标系旋转顺序为，先绕z轴旋转，再绕y轴旋转，最后绕x轴旋转，因此，绕y轴旋转的角度范围只有
±
90
°
；
16.以探空仪的质心为原点o、以从探空仪的质心指向绳与探空仪连接点方向为 oxb轴、以从探空仪的质心指向探空仪横向中心线方向为ozb轴、以与oxb轴和ozb轴垂直方向为oyb轴，建立与探空仪固连的坐标系oxbybzb即b系；
17.以伞底为原点o、以从伞底指向地理坐标系东向为oyn轴、以从伞底指向地理坐标系北向为ozn轴、以从伞底指向地理坐标系天向为oxn轴，建立与伞底固连的坐标系ox
nynzn
即n系；在探空仪静止即探空仪与伞不发生相对摆动时，b系能够由n系沿绳平移得到；
18.以下投时刻伞底为原点o、oxi轴指向天向、oyi轴指向东向、ozi轴指向北向，建立空间直角坐标系oxiyizi即i系；
19.定义绳与oxi轴夹角为θ
x
，绳与oyi轴夹角为θy；
20.对原始数据的三次样条插值方法为：
21.s3-1：针对绕oxb轴旋转角度θ，原始数据中有n 1个数据节点：(θ0,t0)， (θ1,t1)，
…
，(θi,ti)，
…
，(θn,tn)，i＝0,1,...,n，其中，ti为当前数据记录时间，θi为 ti时刻对应的oxb轴旋转角度值，在三次样条插值过程中，满足：在每个分段区间上 [ti,t
i 1
]上，f(t)＝fi
(t)都是一个三次方程；插值条件f(ti)＝θi；曲线光滑，即f(t)、 f
′
(t)、f
″
(t)连续；令三次方程各项系数为ai、bi、ci、di，则有：
[0022][0023]
其中，f(t)为经过三次样条插值后的定义域在(t0,tn)上表示t时刻的角度θ＝f(t)的函数，fi(t)为分段区间[ti,t
i 1
]上表示对应时刻角度的函数，fi′
(t)、fi″
(t) 分别为fi(t)的一阶、二阶导函数；
[0024]
s3-2：得到关于未知数的线性方程组：
[0025]
s3-2-1：由fi(ti)＝ai bi(t
i-ti) ci(t
i-ti)2 di(t
i-ti)3＝θi，得到ai＝θi，其中， ai,bi,ci,di均为系数；
[0026]
s3-2-2：令hi＝t
i 1-ti，由边界连续条件fi(t
i 1
)＝θ
i 1
推出：
[0027][0028]
由边界导数连续条件fi′
(t
i 1
)＝f
′
i 1
(t
i 1
)，得到：
[0029][0030]
由fi″
(t
i 1
)＝f
″
i 1
(t
i 1
)，得到：
[0031][0032]
令mi＝fi″
(ti)＝2ci，则式(2)写为mi 6h
idi
＝m
i 1
，则：hi,mi分别为中间变量；
[0033]
s3-2-3：将ai＝ti，代入式(1)中，得到：
[0034][0035]
将ai，bi，ci，di代入式(2)中，得到：
[0036][0037]
s3-2-4：将式(5)拓展至所有时间点，得到以m0,m1,m2,...,mn为未知变量的线性方程组，在插值过程中选择自然边界条件，即m0＝0，mn＝0，则方程组为：
[0038][0039]
至此，在每个分段区间[ti,t
i 1
]上，通过对方程组的求解能够确定ai，bi，ci，di；
[0040]
s3-3：设置新的采样时间间隔则经过插值的角度θ为：
[0041]
θ
i j
＝fi(ti jδti)，j＝1,2,...,n-1
[0042]
通过控制n值即得到相应的时间精度信息；
[0043]
s3-4：对于绕oyb轴的旋转角度γ、绕ozb轴的旋转角度ψ数据，通过相同的三次样条插值方法进行数据插值，插值前的数据共有n 1个，经过插值后有nn 1个。
[0044]
进一步地，所述步骤s4中，在时刻ti，探空仪所在坐标系oxbybzb相对于伞所在的坐标系ox
nynzn
在仅受张紧的绳的限制下发生位移与角度变换关系，通过姿态角解算探空仪相对于伞的运动速度
[0045]
s4-1：按照z-y-x坐标轴旋转顺序，从n系到b系的坐标变换矩阵为：
[0046][0047]
其中，c
θ
是绕oxn旋转θ的旋转矩阵、c
γ
是绕oyn旋转γ的旋转矩阵、c
ψ
是绕 ozn旋转ψ的旋转矩阵，分别表示为：
[0048]
[0049]
s4-2：在只考虑b系相对于n系的旋转关系而不考虑位移关系时， oxn＝[1 0 0]
t
在b系下的表示表示为：
[0050][0051]
θ
x
表示为与oxb之间的夹角：
[0052][0053]
s4-3：在只考虑考虑b系相对于n系的旋转关系而不考虑位移关系时， oxb＝[1 0 0]
t
在b系下的表示表示为：
[0054][0055]
向量oh即在oy
nzn
平面上的投影，n系下表示为：
[0056]
ohn＝[0
ꢀ‑
cosγcosψ sinγ]
t
[0057]
根据几何关系，θy表示为oyn与oh之间的夹角：
[0058][0059]
s4-4：设绳长为l，则tk时刻，k＝0,1,...,nn，经过角度解算得到θ
x
(k)与θy(k)，则探空仪相对于伞底的位置为：
[0060][0061]
通过前向差分，得到tk时刻探空仪相对于伞的东向速度北向速度
[0062][0063]
其中，为tk时刻由经过处理的三轴角速度计算得到的东向速度，为 tk由经过处理的三轴角速度计算得到的北向速度。
[0064]
进一步地，所述步骤s5中，探空仪的gps/北斗信息，包括探空仪相对于地理坐标系的北向速度探空仪相对于地理坐标系东向速度由探空仪的gps/北斗信息得到的风速为：
[0065][0066]
由于探空仪和伞之间依靠绳连接，gps测算得到的风速实际上是探空仪相对于伞
的运动与风速的叠加，因此将gps测算得到的风速分解为：
[0067][0068]
其中，为惯性系统测算出探空仪相对于伞即n系的速度，为在伞下点坐标系即i系下的风速；
[0069]
ti时刻对应的速度与一定能够找到对应时刻的速度与则此时的风速为：
[0070][0071]
对vw进行处理，得到最终经过补偿后的风速v
wind
：
[0072][0073]
其中，floor(
·
)为向下取整。
[0074]
本发明的有益效果在于：
[0075]
1.本发明提出的基于惯性系统的探空仪风速测量误差补偿方法，通过在传统探空仪系统上集成惯性组件，获取探空仪运动的真实姿态信息，进而补偿因探空仪自身运动所导致的风速测量误差。
[0076]
2.本发明利用惯性系统的实测数据进行误差补偿，相较于其他方法估计的姿态误差而言，更具真实性与可靠性，可根据探空仪每次测试的实际情况逐点高精度补偿，补偿探空仪所测的实时风速信息，而不是仅仅通过估计出探空仪大概的摆动周期补偿风速，避免了这种方法过补偿和欠补偿的问题。
附图说明
[0077]
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，通过参考附图会更加清楚的理解本发明的特征和优点，附图是示意性的而不应理解为对本发明进行任何限制，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，可以根据这些附图获得其他的附图。其中：
[0078]
图1是本发明的惯性系统补偿探空仪风速测量误差流程图；
[0079]
图2是探空仪与伞的坐标系定义；
[0080]
图3是空间方位角定义；
[0081]
图4是空间坐标系位置关系；
[0082]
图5是探空仪相对于降落伞的实测运动轨迹；
[0083]
图6是自研搭载惯导系统的探空仪实物图；
[0084]
图7是探空仪上惯导模块安装位置及三轴定义；
[0085]
图8是探空仪东向、北向风速补偿结果，其中，(a)为探空仪东向风速补偿结果，虚线为gps解算得到的东向风速，实线为基于惯性系统补偿后的东向风速，(b)为探空仪北向
风速补偿结果，虚线为gps解算得到的北向风速，实线为基于惯性系统补偿后的北向风速；
[0086]
图9是探空仪东向、北向风速矢量合成后的补偿结果，虚线为gps解算得到的风速，实线为基于惯性系统补偿后的风速。
具体实施方式
[0087]
为了能够更清楚地理解本发明的上述目的、特征和优点，下面结合附图和具体实施方式对本发明进行进一步的详细描述。需要说明的是，在不冲突的情况下，本发明的实施例及实施例中的特征可以相互组合。
[0088]
在下面的描述中阐述了很多具体细节以便于充分理解本发明，但是，本发明还可以采用其他不同于在此描述的其他方式来实施，因此，本发明的保护范围并不受下面公开的具体实施例的限制。
[0089]
本发明在探空仪上集成惯性系统，并分析所捕获的探空仪真实运动姿态信息，以此来补偿探空仪在恶劣工况下因自身运动所导致的风速测量误差。
[0090]
如图1所示，一种基于惯性系统的探空仪风速测量误差补偿方法，风速测量过程中的下投结构包括伞、绳及探空仪，探空仪通过绳与伞底相连，伞为探空仪提供浮空时间，绳在下投过程中保持拉紧状态，补偿方法包括以下步骤：
[0091]
s1：在探空仪上集成惯性系统，进行高空下投探测，获取立体空间气象数据，地面接收机实时记录并存储探空仪原位探测数据及三轴惯导信息；
[0092]
s2：利用接收到的探空仪gps/北斗信息，计算探空仪探测到的立体空间风速信息；
[0093]
s3：利用接收到的探空仪三轴惯导信息，计算探空仪下落过程中的姿态信息；
[0094]
步骤s3中的三轴惯导信息包括三轴加速度、三轴角速度以及三轴角度信息，保留角度范围在(-180
°
,180
°
)的三轴角度信息，通过三次样条插值对三轴角度原始数据插值，再对位移进行差分求解得到运动速度；
[0095]
如图2所示，惯性系统使用欧拉角表示姿态时的坐标系旋转顺序为，先绕z轴旋转，再绕y轴旋转，最后绕x轴旋转，因此，绕y轴旋转的角度范围只有
±
90
°
；
[0096]
以探空仪的质心为原点o、以从探空仪的质心指向绳与探空仪连接点方向为 oxb轴、以从探空仪的质心指向探空仪横向中心线方向为ozb轴、以与oxb轴和ozb轴垂直方向为oyb轴，建立与探空仪固连的坐标系oxbybzb即b系；
[0097]
以伞底为原点o、以从伞底指向地理坐标系东向为oyn轴、以从伞底指向地理坐标系北向为ozn轴、以从伞底指向地理坐标系天向为oxn轴，建立与伞底固连的坐标系ox
nynzn
即n系；在探空仪静止即探空仪与伞不发生相对摆动时，b系能够由n系沿绳平移得到；
[0098]
以下投时刻伞底为原点o、oxi轴指向天向、oyi轴指向东向、ozi轴指向北向，建立空间直角坐标系oxiyizi即i系；
[0099]
定义绳与oxi轴夹角为θ
x
，绳与oyi轴夹角为θy；
[0100]
对原始数据的三次样条插值方法为：
[0101]
s3-1：针对绕oxb轴旋转角度θ，原始数据中有n 1个数据节点：(θ0,t0)， (θ1,t1)，
…
，(θi,ti)，
…
，(θn,tn)，i＝0,1,...,n，其中，ti为当前数据记录时间，θi为 ti时刻对应的oxb轴旋转角度值，在三次样条插值过程中，满足：在每个分段区间上 [ti,t
i 1
]上，f(t)＝fi(t)都是一个三次方程；插值条件f(ti)＝θi；曲线光滑，即f(t)、 f
′
(t)、f
″
(t)连续；令三
次方程各项系数为ai、bi、ci、di，则有：
[0102][0103]
其中，f(t)为经过三次样条插值后的定义域在(t0,tn)上表示t时刻的角度θ＝f(t)的函数，fi(t)为分段区间[ti,t
i 1
]上表示对应时刻角度的函数，fi′
(t)、fi″
(t) 分别为fi(t)的一阶、二阶导函数；
[0104]
s3-2：得到关于未知数的线性方程组：
[0105]
s3-2-1：由fi(ti)＝ai bi(t
i-ti) ci(t
i-ti)2 di(t
i-ti)3＝θi，得到ai＝θi，其中， ai,bi,ci,di均为系数；
[0106]
s3-2-2：令hi＝t
i 1-ti，由边界连续条件fi(t
i 1
)＝θ
i 1
推出：
[0107][0108]
由边界导数连续条件fi′
(t
i 1
)＝f
′
i 1
(t
i 1
)，得到：
[0109][0110]
由fi″
(t
i 1
)＝f
″
i 1
(t
i 1
)，得到：
[0111][0112]
令mi＝fi″
(ti)＝2ci，则式(2)写为mi 6h
idi
＝m
i 1
，则：hi,mi分别为中间变量；
[0113]
s3-2-3：将ai＝ti，代入式(1)中，得到：
[0114][0115]
将ai，bi，ci，di代入式(2)中，得到：
[0116][0117]
s3-2-4：将式(5)拓展至所有时间点，得到以m0,m1,m2,...,mn为未知变量的线性方程组，在插值过程中选择自然边界条件，即m0＝0，mn＝0，则方程组为：
[0118][0119][0120]
至此，在每个分段区间[ti,t
i 1
]上，通过对方程组的求解能够确定ai，bi，ci，di；
[0121]
s3-3：设置新的采样时间间隔则经过插值的角度θ为：
[0122]
θ
i j
＝fi(ti jδti)，j＝1,2,...,n-1
[0123]
通过控制n值即得到相应的时间精度信息；
[0124]
s3-4：对于绕oyb轴的旋转角度γ、绕ozb轴的旋转角度ψ数据，通过相同的三次样条插值方法进行数据插值，插值前的数据共有n 1个，经过插值后有nn 1个。
[0125]
s4：将步骤s3的结果通过空间坐标系变换与角度解算，获取探空仪的实时相对空间运动姿态信息；步骤s4中，在时刻ti，探空仪所在坐标系oxbybzb相对于伞所在的坐标系ox
nynzn
在仅受张紧的绳的限制下发生位移与角度变换关系，如图4 所示，通过姿态角解算探空仪相对于伞的运动速度
[0126]
s4-1：按照z-y-x坐标轴旋转顺序，从n系到b系的坐标变换矩阵为：
[0127][0128]
其中，c
θ
是绕oxn旋转θ的旋转矩阵、c
γ
是绕oyn旋转γ的旋转矩阵、c
ψ
是绕ozn旋转ψ的旋转矩阵，分别表示为：
[0129][0130]
s4-2：在只考虑b系相对于n系的旋转关系而不考虑位移关系时， oxn＝[1 0 0]
t
在b系下的表示表示为：
[0131][0132]
如图3所示，θ
x
表示为与oxb之间的夹角：
[0133][0134]
s4-3：在只考虑考虑b系相对于n系的旋转关系而不考虑位移关系时， oxb＝[1 0 0]
t
在b系下的表示表示为：
[0135][0136]
向量oh即在oy
nzn
平面上的投影，n系下表示为：
[0137]
ohn＝[0
ꢀ‑
cosγcosψ sinγ]
t
[0138]
根据几何关系，θy表示为oyn与oh之间的夹角：
[0139][0140]
s4-4：设绳长为l，则tk时刻，k＝0,1,...,nn，经过角度解算得到θ
x
(k)与θy(k)，则探空仪相对于伞底的位置为：
[0141][0142]
通过前向差分，得到tk时刻探空仪相对于伞的东向速度北向速度
[0143][0144]
其中，为tk时刻由经过处理的三轴角速度计算得到的东向速度，为 tk由经过处理的三轴角速度计算得到的北向速度。
[0145]
s5：通过步骤s4得到的探空仪的相对空间运动姿态信息，计算探空仪因自身姿态变化导致的风速测量误差，代入步骤s2解算的立体空间风速信息中，完成风速补偿；步骤s5
中，探空仪的gps/北斗信息，包括探空仪相对于地理坐标系的北向速度探空仪相对于地理坐标系东向速度由探空仪的gps/北斗信息得到的风速为：
[0146][0147]
由于探空仪和伞之间依靠绳连接，gps测算得到的风速实际上是探空仪相对于伞的运动与风速的叠加，因此将gps测算得到的风速分解为：
[0148][0149]
其中，为惯性系统测算出探空仪相对于伞即n系的速度，为在伞下点坐标系即i系下的风速；
[0150]
ti时刻对应的速度与一定能够找到对应时刻的速度与则此时的风速为：
[0151][0152]
对vw进行处理，得到最终经过补偿后的风速v
wind
：
[0153][0154]
其中，floor(
·
)为向下取整。
[0155]
s6：对经过步骤s5补偿后的风速数据进行滚动时间窗平滑处理，进一步消除环境噪声的影响。
[0156]
为了方便理解本发明的上述技术方案，以下通过具体实施例对本发明的上述技术方案进行详细说明。
[0157]
实施例1
[0158]
如图6-7所示，利用自研搭载惯导系统的探空仪在新疆进行20km下投实验，探测立体空间气象数据。探空仪相对于降落伞的实测运动轨迹如图5所示，解算得到由于探空仪的随机摆动带来的风速测量漂移，分别对gps测算得到的东向速度与北向速度进行补偿。补偿效果如图8所示，其中，(a)为探空仪东向风速补偿结果，虚线为 gps解算得到的东向风速，实线为基于惯性系统补偿后的东向风速，(b)为探空仪北向风速补偿结果，虚线为gps解算得到的北向风速，实线为基于惯性系统补偿后的北向风速。可以看出，经过惯导数据补偿后的东向、北向风速的信噪比有较为明显的提升。将gps测算得到的风速与经过惯导补偿的风速进行矢量合成，得到最终的风速测量结果如图9所示，虚线为gps解算得到的风速，实线为基于惯性系统补偿后的风速，可以看出，经过惯导数据补偿后风速的信噪比有较为明显的提升。
[0159]
以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修
改、等同替换、改进等，均应包含在本发明的保护范围之内。