CN112785052B - 基于粒子滤波算法的风速风向预测方法 - Google Patents

基于粒子滤波算法的风速风向预测方法 Download PDF

Info

Publication number
CN112785052B
CN112785052B CN202110047680.2A CN202110047680A CN112785052B CN 112785052 B CN112785052 B CN 112785052B CN 202110047680 A CN202110047680 A CN 202110047680A CN 112785052 B CN112785052 B CN 112785052B
Authority
CN
China
Prior art keywords
wind speed
wind
state
wind direction
particle
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202110047680.2A
Other languages
English (en)
Other versions
CN112785052A (zh
Inventor
李新波
王晓玉
姜良旭
左昕雨
范书源
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Jilin University
Original Assignee
Jilin University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Jilin University filed Critical Jilin University
Priority to CN202110047680.2A priority Critical patent/CN112785052B/zh
Publication of CN112785052A publication Critical patent/CN112785052A/zh
Application granted granted Critical
Publication of CN112785052B publication Critical patent/CN112785052B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
    • G06Q10/00—Administration; Management
    • G06Q10/04—Forecasting or optimisation specially adapted for administrative or management purposes, e.g. linear programming or "cutting stock problem"
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01P—MEASURING LINEAR OR ANGULAR SPEED, ACCELERATION, DECELERATION, OR SHOCK; INDICATING PRESENCE, ABSENCE, OR DIRECTION, OF MOVEMENT
    • G01P13/00—Indicating or recording presence, absence, or direction, of movement
    • G01P13/02—Indicating direction only, e.g. by weather vane
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01P—MEASURING LINEAR OR ANGULAR SPEED, ACCELERATION, DECELERATION, OR SHOCK; INDICATING PRESENCE, ABSENCE, OR DIRECTION, OF MOVEMENT
    • G01P5/00—Measuring speed of fluids, e.g. of air stream; Measuring speed of bodies relative to fluids, e.g. of ship, of aircraft
    • G01P5/24—Measuring speed of fluids, e.g. of air stream; Measuring speed of bodies relative to fluids, e.g. of ship, of aircraft by measuring the direct influence of the streaming fluid on the properties of a detecting acoustical wave
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/004—Artificial life, i.e. computing arrangements simulating life
    • G06N3/006—Artificial life, i.e. computing arrangements simulating life based on simulated virtual individual or collective life forms, e.g. social simulations or particle swarm optimisation [PSO]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
    • G06F2218/12—Classification; Matching

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Data Mining & Analysis (AREA)
  • Mathematical Optimization (AREA)
  • General Engineering & Computer Science (AREA)
  • Pure & Applied Mathematics (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Software Systems (AREA)
  • Operations Research (AREA)
  • Business, Economics & Management (AREA)
  • Databases & Information Systems (AREA)
  • Algebra (AREA)
  • Evolutionary Computation (AREA)
  • Strategic Management (AREA)
  • Human Resources & Organizations (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Economics (AREA)
  • Probability & Statistics with Applications (AREA)
  • Computer Hardware Design (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Development Economics (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Game Theory and Decision Science (AREA)
  • Aviation & Aerospace Engineering (AREA)
  • Multimedia (AREA)
  • Entrepreneurship & Innovation (AREA)
  • Marketing (AREA)
  • Quality & Reliability (AREA)
  • Tourism & Hospitality (AREA)
  • General Business, Economics & Management (AREA)
  • Evolutionary Biology (AREA)
  • Acoustics & Sound (AREA)
  • Geometry (AREA)
  • Health & Medical Sciences (AREA)
  • Artificial Intelligence (AREA)
  • Biomedical Technology (AREA)

Abstract

本发明公开了一种基于粒子滤波算法的风速风向预测方法，属于风速风向测量技术领域，包括构建弧形阵列测风结构，得到阵列流型矢量；以风速风向参数作为待估计量，建立用于预测风速风向的状态空间模型，所述状态空间模型由状态方程和观测方程组成；根据状态方程和观测方程，利用粒子滤波算法估计风速风向的状态值，实现风速风向的动态预测；所述粒子滤波算法的似然函数采用多重信号分类MUSIC算法空间谱函数，本发明实现了风速风向值的动态预测，并且具有较高的风速风向预测精度。

Description

基于粒子滤波算法的风速风向预测方法
技术领域
本发明属于风速风向测量技术领域，具体涉及了一种基于粒子滤波算法的风速风向预测方法。
背景技术
风是由空气流动而产生的一种自然现象，在现今生活中有着广泛的应用，风速风向的准确测量及动态预测在天气预报、环境监控、风力发电、航空航天等领域具有非常重要的意义。目前，常用的测风仪有机械式风速仪、超声波式风速仪、皮托管式风速仪、热线式风速仪、激光多普勒风速仪、雷达多普勒风速仪等。相比于其他种类的测风仪，超声波式风速仪具有无机械磨损、测量速度快、维护成本低、性价比高的优势，但是其测风精度完全取决于超声波信号在顺风和逆风时传播时间的测量精度。
为了应用多维空域信息，提高风参数的测量精度，基于阵列信号处理思想，现有技术中提出了一种风速风向测量方法，详见专利文献号CN104897925B，该方法主要采用一发多收的均匀线性阵列结构，结合多重信号分类MUSIC算法估计风速风向值。并且通过仿真实验验证了所述方法的可行性，但是测量精度受信噪比影响较大。为了提高系统抗干扰能力，现有技术中还公开了一种风向风速测量装置，该装置采用一发多收的弧形超声波传感器阵列结构，采用波束形成算法估计风速和风向，并且该方案通过仿真实验证明了在信噪比大于10db时，可以实现测量范围内所有风速和风向均方根误差为零的参数估计。但是该方案是基于静态的风速风向估计，即假设风速风向值为常数，当风速风向值改变时，需采取批处理算法，不断地进行二维谱峰搜索，计算量较大，不便于实现风速和风向的连续测量。由以上内容可知，虽然现有阵列测风结构和方法利用了多维空域信息，具有提高测量精度的优势，但是只能对静态风速风向值进行估计，不能对变化的风速风向值进行动态预测。如何提高风速风向的测量精度以及实现动态预测是测风领域急需解决的问题。
发明内容
本发明的目的是针对现有的风速风向测量只能对静态风速风向值进行估计，不能对变化的风速风向值进行动态预测的问题，而提出了一种基于粒子滤波算法的风速风向预测方法，该方法在阵列信号数学模型的基础上建立了状态方程和观测方程，然后将多重信号分类MUSIC算法的空间谱函数作为粒子滤波算法的似然函数，实现风速风向值的动态预测，具有较高的风速风向预测精度。
为实现上述目的，本发明采用如下技术方案：基于粒子滤波算法的风速风向预测方法，其特征在于，包括如下步骤：
步骤1：构建弧形阵列测风结构，得到阵列流型矢量；
弧形阵列测风结构由五个超声波传感器组成，五个超声波传感器中的一个为发射阵元，另外四个为接收阵元，以发射阵元为中心，四个接收阵元处于同一弧线上，α为相邻两接收阵元与发射阵元连线夹角，V为待测风速的大小，θ为待测风向角，风的来向为北偏东θ°方向，根据矢量的分解得到待测风速V在发射阵元与四个接收阵元连线上的分量，分别为V1、V2、V3及V4，其中
Figure BDA0002897882250000021
Figure BDA0002897882250000022
得到发射阵元所发射信号到第i个接收阵元的时间为：
Figure BDA0002897882250000023
R为弧形半径，c为超声波传播速度，Vi为待测风速V在发射阵元与第i个接收阵元连线上的分量，以第1个接收阵元为参考阵元，得到发射信号到达第i个接收阵元相对于参考阵元的时延为：
τ1＝0
Figure BDA0002897882250000024
Figure BDA0002897882250000031
Figure BDA0002897882250000032
因此，阵列流型矢量为：
Figure BDA0002897882250000033
其中f为发射阵元所发射的超声波频率，j表示复数域中的虚数单位；
步骤2：建立预测风速风向的状态空间模型
以风速风向参数作为待估计量，建立用于预测风速风向的状态空间模型，所述状态空间模型由状态方程和观测方程组成；
所述状态方程表示待估计量的变化规律，具体表达形式如下：
X(n+1)＝FcvX(n)+Gcvv(n)
其中X(n+1)为(n+1)时刻风速风向的状态向量；X(n)为n时刻风速风向的状态向量，
Figure BDA0002897882250000034
V(n)为n时刻风速，
Figure BDA0002897882250000035
为n时刻风速变化率，θ(n)为n时刻风向，
Figure BDA0002897882250000036
为n时刻风向变化率，T表示转置，
Figure BDA0002897882250000037
ΔT表示相邻两次风速风向预测的时间间隔；v(n)＝[v1(n),v1(n),v2(n),v2(n)]T为状态噪声矢量，其中v1(n)为状态空间模型迭代风速时的噪声，v2(n)为状态空间模型迭代风向时的噪声；
所述观测方程表示待估计量与观测量之间的函数关系，具体表达形式如下：
Z(n)＝a(θ,V)s(n)+W(n)
其中n＝1,2,…,L，L为快拍数，Z(n)为阵列接收信号矢量，a(θ,V)为阵列流型矢量，s(n)为发射源的发射信号；W(n)为接收阵元接收到的噪声矢量，W(n)＝[w1(n),w2(n),w3(n),w4(n)]T，T为矩阵转置，w1(n)、w2(n)、w3(n)及w4(n)分别为四个接收阵元接收到的均值为零、方差为δ2的高斯白噪声；
步骤3：根据状态方程和观测方程，利用粒子滤波算法估计风速风向的状态值，实现风速风向的动态预测；
所述粒子滤波算法的似然函数采用多重信号分类MUSIC算法空间谱函数；
具体构建多重信号分类MUSIC算法空间谱函数如下：
Figure BDA0002897882250000041
其中U为噪声子空间，a(θ,V)为阵列流型矢量，H为共轭转置，aH(θ,V)为a(θ,V)的共轭转置，UH为U的共轭转置；
因此，粒子滤波算法的似然函数表达形式如下：
Figure BDA0002897882250000042
式中i为粒子序号，
Figure BDA0002897882250000043
为k时刻，第i个粒子的风速风向状态向量；
Figure BDA0002897882250000044
为k时刻，第i个粒子的风速风向观测数据；a(θi,Vi)为第i个粒子的阵列流型矢量，aH(θi,Vi)为a(θi,Vi)的共轭转置，Uk为k时刻根据风速风向的状态向量得到的噪声子空间，Uk H为Uk的共轭转置。
进一步，步骤3中利用粒子滤波算法估计风速风向的状态值的过程如下：
①初始化：k时刻，采样NS个粒子，第i个粒子的风速风向状态向量为
Figure BDA0002897882250000051
得到集合
Figure BDA0002897882250000052
各粒子权重赋值为1/NS，各粒子初始状态风速风向值采用MUSIC算法进行估计，风速变化率
Figure BDA0002897882250000053
为风速测量精度0.1m/s，风向变化率
Figure BDA0002897882250000054
为风向测量精度1°；
②更新权值：k时刻，按照弧形阵列的协方差矩
Figure BDA0002897882250000055
E[g]为计算数学期望的表达式，本发明中
Figure BDA0002897882250000056
Zk(n)为k时刻阵列接收信号矢量，
Figure BDA0002897882250000057
为Zk(n)的共轭转置；将其进行特征值分解得到噪声子空间，根据在k时刻获得先验样本，得到第i个粒子的风速风向状态为
Figure BDA0002897882250000058
表示状态方程迭代方式，即先验概率密度函数，
Figure BDA0002897882250000059
为经过迭代得到的第i个粒子前k-1时刻的风速风向的状态向量，抽取粒子得到粒子状态，得到第i个粒子的阵列流型矢量a(θ(i),V(i))，利用似然函数
Figure BDA00028978822500000510
根据粒子权重更新方程
Figure BDA00028978822500000511
计算粒子权重；
Figure BDA00028978822500000512
为k时刻第i个粒子的权重，
Figure BDA00028978822500000513
为k-1时刻第i个粒子的权重；
③重采样：按照公式
Figure BDA0002897882250000061
计算归一化权重，重新采样NS个粒子并分配权值为1/NS；
④状态估计：按照后验概率密度函数
Figure BDA0002897882250000062
计算风速风向状态的估计值，即风速风向的状态值；
Xk为k时刻风速风向状态向量，
Figure BDA0002897882250000063
为第i个粒子通过迭代过程经过前k个时刻的观测得到的观测向量；
Figure BDA0002897882250000064
为Xk的狄拉克函数；
⑤预测：根据X(n+1)＝FcvX(n)+Gcvv(n)预测下一时刻的状态；
⑥k＝k+1，重复步骤②至步骤⑥，当k大于观测长度，终止预测。
通过上述设计方案，本发明可以带来如下有益效果：本发明提出了一种基于粒子滤波算法的风速风向预测方法，该方法构建的弧形阵列测风结构采用超声波传感器发射和接收超声波信号，由于风信号的存在，使得发射信号到达各接收阵元存在时间延迟，得到阵列接收信号矢量并将其作为观测方程，由于弧形阵列测风结构利用了多维空域信息，避免了基于时差法原理的超声测风仪由于测量超声波传播时间的准确性而影响测风准确性的问题，提高了风参数的测量精度。并根据风速风向匀速变化规律构建二阶匀速模型，将其作为状态方程。结合状态方程和观测方程得到状态空间模型。将粒子滤波算法与阵列信号处理中多重信号分类MUSIC算法相结合实现了风速风向值的动态预测。通过仿真实验验证了本发明所提方法可以对风速风向值进行动态预测，并且在信噪比大于5dB时，风速和风向预测的均方根误差均小于其变化率，具有较高的风速风向预测精度。
附图说明
图1为基于粒子滤波算法的风速风向预测方法的流程图；
图2为弧形阵列测风结构示意图；
图3为风向预测结果图；
图4为风速预测结果图；
图5为不同信噪比下风速和风向预测的均方根误差关系图。
具体实施方式
为了使本发明的目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。
本发明提出了一种基于粒子滤波算法的风速风向预测方法，该方法包括如下过程：
一、建立预测风速风向的状态空间模型：
以风速风向参数作为待估计量，建立状态方程和观测方程，组成用于预测风速风向的状态空间模型。其中状态方程表示待估计量的变化规律，观测方程表示待估计量与观测量之间的函数关系。
1、状态方程
通常风速风向的变化是相对缓慢的，因此，可以采用CV模型，即二阶匀速模型，模拟风速风向的变化规律，在所述二阶匀速模型中风速风向值的变化速度为常数，变化的加速度为高斯白噪声，其中加速度反应了因扰动而产生的无法预测的建模误差。
设n时刻的风速风向状态向量为：
Figure BDA0002897882250000071
其中，T为矩阵转置，V(n)和
Figure BDA0002897882250000072
分别为n时刻风速和风速变化率，风速加速度
Figure BDA0002897882250000081
其中v1(n)为风速加速度，服从均值为零，v1(n)的方差为δv1的高斯分布；θ(n)和
Figure BDA0002897882250000082
分别为n时刻风向和风向变化率，风向加速度
Figure BDA0002897882250000083
其中v2(n)为风向加速度，服从均值为零，v2(n)的方差为δv2的高斯分布。由CV模型可以得出：
Figure BDA0002897882250000084
Figure BDA0002897882250000085
Figure BDA0002897882250000086
Figure BDA0002897882250000087
其中，ΔT表示前一个时间步长和当前时间步长之间的时间周期，即相邻两次风速风向预测的时间间隔(以秒为单位)。
在进行风速风向预测时，令风速风向状态向量为
Figure BDA0002897882250000088
则相应的状态方程为：
X(n+1)＝FcvX(n)+Gcvv(n) (5)
其中，
Figure BDA0002897882250000091
v(n)＝[v1(n),v1(n),v2(n),v2(n)]T为状态噪声矢量，其中v1(n)为状态空间模型迭代风速时的噪声，v2(n)为状态空间模型迭代风向时的噪声。
2、观测方程
如图2所示，弧形阵列测风结构由五个超声波传感器组成，五个超声波传感器中的一个为发射阵元，另外四个为接收阵元，α为相邻两接收阵元与发射阵元连线夹角，V为待测风速的大小，θ为待测风向角(风的来向为以竖直正北方向顺时针偏θ°方向)，由于发射的超声波频率大于20KHz，进行离散化的采样超声波频率大于20KHz，因此，在一个极小的采样周期内，即采样周期小于50微秒，风速和风向的幅值不变。根据矢量的分解得到待测风速V在发射阵元与四个接收阵元连线上的分量，分别为V1、V2、V3及V4，其中
Figure BDA0002897882250000092
Figure BDA0002897882250000093
得到发射信号到第i个接收阵元的时间为：
Figure BDA0002897882250000094
(i＝1,2,3,4)，R为弧形半径，c为超声波传播速度，Vi为待测风速V在发射阵元与第i个接收阵元连线上的分量，以第1个接收阵元为参考阵元，得到发射信号到达第i个接收阵元相对于参考阵元的时延为：
τ1＝0 (6)
Figure BDA0002897882250000101
Figure BDA0002897882250000102
Figure BDA0002897882250000103
因此，弧形阵列测风结构的阵列流型矢量为：
Figure BDA0002897882250000104
其中：f为发射的超声波频率，j表示复数域中的虚数单位。
由于只有一个发射信号，阵列流型矩阵A(θ,V)＝a(θ,V)，进而得到阵列接收信号矢量：
Z(n)＝a(θ,V)s(n)+W(n) (11)
其中n＝1,2,…,L，L为快拍数，Z(n)为阵列接收信号矢量，a(θ,V)为阵列流型矢量，s(n)为发射源的发射信号；W(n)为接收阵元接收到的噪声矢量，W(n)＝[w1(n),w2(n),w3(n),w4(n)]T，T为矩阵转置，w1(n)、w2(n)、w3(n)及w4(n)分别为四个接收阵元接收到的均值为零、方差为δ2的高斯白噪声；(5)式表明了待估计风速和风向的变化规律，其作为状态方程，(11)式表明了待估计量与观测量之间的关系，将其作为观测方程，联合(5)式即为状态空间模型，此状态空间模型可以对风速风向的变化进行预测。
二、粒子滤波算法
粒子滤波是一种基于蒙特卡洛仿真的近似贝叶斯滤波算法，其主要思想是用一群随机样本(粒子)来描述后验概率分布，通过测量模型来调节各样本的权值和位置，以逼近真实的概率分布，以样本均值代替积分运算，从而获得状态的最小方差，得到系统的估计值，原则上适用于任意非线性、非高斯随机系统的状态估计。
1、序贯蒙特卡罗估计
根据阵列接收信号矢量，得到k时刻，L个快拍数下的观测数据为：
Zk＝{Z(1+kL),Z(2+kL),...,Z(L+kL)} (12)
公式(12)中Z(1+kL)表示第1+kL快拍接收阵列接收到的数据；Z(2+kL)第2+kL快拍接收阵列接收到的数据，Z(L+kL)表示第L+kL快拍接收阵列接收到的数据；L为快拍数，k为第k时刻；
用Z1:K＝{Z1,Z2,...,ZK}表示离散时间的测量序列，K为采样周期，所提出的粒子滤波MUSIC算法基于状态方程和观测方程递归估计风速风向值，实现风参数的预测。整个过程分为预测和更新两个阶段，预测阶段根据先验概率密度函数生成一组采样粒子
Figure BDA0002897882250000111
每个样本通过状态方程，在k时刻获得先验样本，得到第i个粒子的风速风向状态为：
Figure BDA0002897882250000112
Figure BDA0002897882250000113
表示状态方程迭代方式，即先验概率密度函数，
Figure BDA0002897882250000114
表示第i个粒子k时刻的风速风向状态向量；
Figure BDA0002897882250000115
表示第i个粒子前k-1时刻迭代得到的风速风向状态向量，更新阶段，在收到测量结果ZK后，评估每个先验样本的可能性，并获得每个样本的权重：
Figure BDA0002897882250000121
其中，
Figure BDA0002897882250000122
为粒子滤波算法的似然函数，重要性密度函数
Figure BDA0002897882250000123
表示第i个粒子k时刻的风速风向观测数据；
Figure BDA0002897882250000124
表示第i个粒子k时刻的风速风向状态向量；
Figure BDA0002897882250000125
表示第i个粒子前k-1时刻迭代得到的风速风向状态向量；Z1:k-1为经过迭代得到的第i个粒子前k-1时刻的风速风向的观测向量；
Figure BDA0002897882250000126
为k时刻第i个粒子的权重，
Figure BDA0002897882250000127
为k-1时刻第i个粒子的权重；
因此，粒子权重更新方程为：
Figure BDA0002897882250000128
得到归一化权重为：
Figure BDA0002897882250000129
经过重采样后，合成粒子，得到后验概率密度函数
Figure BDA00028978822500001210
Xk为k时刻风速风向状态向量，
Figure BDA0002897882250000131
为第i个粒子通过迭代过程经过前k个时刻的观测得到的观测向量；
Figure BDA0002897882250000132
为Xk的狄拉克函数；
2、粒子滤波算法的似然函数
MUSIC算法作为高分辨算法，在进行静态风速风向估计时，其主要思想是利用协方差矩阵特征值分解，得到噪声子空间，构建空间谱函数
Figure BDA0002897882250000133
其中U为噪声子空间，a(θ,V)为阵列流型矢量，H为共轭转置，aH(θ,V)为a(θ,V)的共轭转置，UH为噪声子空间U的共轭转置；，根据PMUSIC(θ,V)通过谱峰搜索得到风速风向值。
k时刻，根据离散时间的测量序列Z1:K＝{Z1,Z2,...,ZK}得到L快拍数下，弧形阵列的协方差矩阵为：
Figure BDA0002897882250000134
请给出E[g]为计算数学期望的表达式，本发明
Figure BDA0002897882250000135
将测风阵列的协方差矩阵进行特征值分解得到噪声子空间Uk，根据此时粒子的风速风向状态
Figure BDA0002897882250000136
结合状态方程得到阵列流型矢量a(θ(i),V(i))。在遍历测量范围内所有风速风向值进行谱峰搜索时，它反映的规律为：当待估计风速风向值与此时的遍历值接近程度越高时，MUSIC算法空间谱函数越大，反之越小。权重的计算是粒子滤波算法的核心，通常用于计算粒子权重的似然函数为高斯函数，它反映了滤波器的一步预测值与超声波传感器测量值的接近程度越高，其权重越大，反之越小的规律。这一规律恰好与MUSIC算法风速风向测量规律一致。因此，可以将空间谱函数PMUSIC(θ,V)代替粒子滤波算法中的似然函数，得到粒子滤波算法的似然函数为：
Figure BDA0002897882250000141
式中i为粒子序号，a(θ(i),V(i))为粒子i的阵列流型矢量，aH(θi,Vi)为粒子i的阵列流型矢量共轭转置；Uk为k时刻根据风速风向状态得到的噪声子空间，Uk H噪声子空间Uk的共轭转置；当粒子的θ(i)、V(i)与真实值越接近，似然函数越大，则权值越大。根据权值的大小能实现“优质”粒子的大量复制，并对“劣质”粒子实现淘汰制。经过权值的计算，重新指导粒子空间分布。
粒子滤波算法(PF)步骤为：
(1)初始化：采样NS个粒子，各粒子权重赋值为1/NS，各粒子初始状态风速风向值采用MUSIC算法进行估计，风速变化率
Figure BDA0002897882250000142
为风速测量精度0.1m/s，风向变化率
Figure BDA0002897882250000143
为风向测量精度1°；
(2)更新权值：k时刻，按照式(19)计算阵列协方差矩阵，将其进行特征值分解得到噪声子空间，根据式(13)抽取粒子得到粒子状态，得到阵列流型矢量a(θ(i),V(i))，根据式(20)计算似然函数，根据式(15)计算粒子权重；
(3)重采样：按照式(16)计算归一化权重，重新采样NS个粒子并分配权值为1/NS；
(4)状态估计：按照式(17)计算状态的估计值；
(5)预测：根据状态方程式(5)预测下一时刻的状态；
(6)k＝k+1，重复步骤(2)至步骤(6)，当k大于观测长度，终止预测。
三、仿真实验
1、可行性实验
超声波传感器发射的超声波信号频率为40KHz；弧形半径为0.1m；各阵元噪声为高斯白噪声；快拍数为2500；信噪比SNR＝10dB；初始风速为5m/s，风向为10°；无噪声情况下，风速变化率为0.1m/s，风向变化率为1°；粒子数为500；预测时间长度为50s。应用弧形阵列测风结构结合以多重信号分类MUSIC算法空间谱函数作为似然函数的粒子滤波算法，得到的风速风向预测仿真结果如图3和图4所示。
由仿真结果可以看出，风向预测值与真实值重合，风速预测值与真实值重合，验证了本发明所提方法可以进行风速风向值的预测。
2、不同信噪比下均方根误差实验
均方根误差公式：
Figure BDA0002897882250000151
其中n为蒙特卡洛实验次数，n＝100；在进行风参数预测时，
Figure BDA0002897882250000152
为每次蒙特卡洛实验在预测时间长度内风参数估计误差的均值。
从仿真结果可以看出，随着信噪比的增加，风速和风向预测的均方根误差随着信噪比的增加而减小，并且在信噪比大于5dB时，风速和风向预测的均方根误差均小于其变化率，验证了本发明具有很好的风速风向预测精度。
综上，本发明将阵列信号处理中的多重信号分类MUSIC算法结合粒子滤波算法用到测风领域，实现了风速风向的动态预测；构建风速风向匀速变化模型，作为粒子滤波算法的状态方程；采用弧形超声波传感器阵列结构作为测风阵列结构，将其阵列接收信号数学模型作为粒子滤波算法的观测方程；由于多重信号分类MUSIC算法的空间谱函数与粒子滤波算法中的似然函数具有相似的特性，因此，将多重信号分类MUSIC算法的空间谱函数PMUSIC(θ,V)代替粒子滤波算法的似然函数，用于更新粒子权重，进行粒子滤波算法重采样。

Claims (2)

1.基于粒子滤波算法的风速风向预测方法，其特征在于，包括如下步骤：
步骤1：构建弧形阵列测风结构，得到阵列流型矢量；
弧形阵列测风结构由五个超声波传感器组成，五个超声波传感器中的一个为发射阵元，另外四个为接收阵元，以发射阵元为中心，四个接收阵元处于同一弧线上，α为相邻两接收阵元与发射阵元连线夹角，V为待测风速的大小，θ为待测风向角，风的来向为北偏东θ°方向，根据矢量的分解得到待测风速V在发射阵元与四个接收阵元连线上的分量，分别为V1、V2、V3及V4，其中
Figure FDA0003534084160000011
Figure FDA0003534084160000012
得到发射阵元所发射信号到第i个接收阵元的时间为：
Figure FDA0003534084160000013
其中i的取值为1、2、3或4，R为弧形半径，c为超声波传播速度，Vi为待测风速V在发射阵元与第i个接收阵元连线上的分量，以第1个接收阵元为参考阵元，得到发射信号到达第i个接收阵元相对于参考阵元的时延为：
τ1＝0
Figure FDA0003534084160000014
Figure FDA0003534084160000015
Figure FDA0003534084160000021
因此，阵列流型矢量为：
Figure FDA0003534084160000022
其中f为发射阵元所发射的超声波频率，j表示复数域中的虚数单位；
步骤2：建立预测风速风向的状态空间模型
以风速风向参数作为待估计量，建立用于预测风速风向的状态空间模型，所述状态空间模型由状态方程和观测方程组成；
所述状态方程表示待估计量的变化规律，具体表达形式如下：
X(n+1)＝FcvX(n)+Gcvv(n)
其中X(n+1)为n+1时刻风速风向的状态向量；X(n)为n时刻风速风向的状态向量，
Figure FDA0003534084160000023
V(n)为n时刻风速，
Figure FDA0003534084160000024
为n时刻风速变化率，θ(n)为n时刻风向，
Figure FDA0003534084160000025
为n时刻风向变化率，T表示转置，
Figure FDA0003534084160000026
ΔT表示相邻两次风速风向预测的时间间隔；v(n)＝[v1(n),v1(n),v2(n),v2(n)]T为状态噪声矢量，其中v1(n)为状态空间模型迭代风速时的噪声，v2(n)为状态空间模型迭代风向时的噪声；
所述观测方程表示待估计量与观测量之间的函数关系，具体表达形式如下：
Z(n)＝a(θ,V)s(n)+W(n)
其中n＝1,2,…,L，L为快拍数，Z(n)为阵列接收信号矢量，a(θ,V)为阵列流型矢量，s(n)为发射源的发射信号；W(n)为接收阵元接收到的噪声矢量，W(n)＝[w1(n),w2(n),w3(n),w4(n)]T，T为矩阵转置，w1(n)、w2(n)、w3(n)及w4(n)分别为四个接收阵元接收到的均值为零、方差为δ2的高斯白噪声；
步骤3：根据状态方程和观测方程，利用粒子滤波算法估计风速风向的状态值，实现风速风向的动态预测；
所述粒子滤波算法的似然函数采用多重信号分类MUSIC算法空间谱函数；
具体构建多重信号分类MUSIC算法空间谱函数如下：
Figure FDA0003534084160000031
其中U为噪声子空间，a(θ,V)为阵列流型矢量，H为共轭转置，aH(θ,V)为a(θ,V)的共轭转置，UH为U的共轭转置；
因此，粒子滤波算法的似然函数表达形式如下：
Figure FDA0003534084160000032
式中i为粒子序号，
Figure FDA0003534084160000033
为k时刻，第i个粒子的风速风向状态向量；
Figure FDA0003534084160000034
为k时刻，第i个粒子的风速风向观测数据；a(θi,Vi)为第i个粒子的阵列流型矢量，aH(θi,Vi)为a(θi,Vi)的共轭转置，Uk为k时刻根据风速风向的状态向量得到的噪声子空间，Uk H为Uk的共轭转置。
2.根据权利要求1所述的基于粒子滤波算法的风速风向预测方法，其特征在于：步骤3中利用粒子滤波算法估计风速风向的状态值的过程如下：
①初始化：k时刻，采样NS个粒子，第i个粒子的风速风向状态向量为
Figure FDA0003534084160000041
得到集合
Figure FDA0003534084160000042
各粒子权重赋值为1/NS，各粒子初始状态风速风向值采用MUSIC算法进行估计，风速变化率
Figure FDA0003534084160000043
为风速测量精度0.1m/s，风向变化率
Figure FDA0003534084160000044
为风向测量精度1°；
②更新权值：k时刻，按照弧形阵列的协方差矩
Figure FDA0003534084160000045
，Zk(n)为k时刻阵列接收信号矢量，
Figure FDA0003534084160000046
为Zk(n)的共轭转置；将其进行特征值分解得到噪声子空间，根据在k时刻获得先验样本，得到第i个粒子的风速风向状态为
Figure FDA0003534084160000047
表示状态方程迭代方式，即先验概率密度函数，
Figure FDA0003534084160000048
为经过迭代得到的第i个粒子前k-1时刻的风速风向的状态向量，抽取粒子得到粒子状态，得到第i个粒子的阵列流型矢量a(θ(i),V(i))，利用似然函数
Figure FDA0003534084160000049
根据粒子权重更新方程
Figure FDA00035340841600000410
计算粒子权重；
Figure FDA00035340841600000411
为k时刻第i个粒子的权重，
Figure FDA00035340841600000412
为k-1时刻第i个粒子的权重；
③重采样：按照公式
Figure FDA0003534084160000051
计算归一化权重，重新采样NS个粒子并分配权值为1/NS；
④状态估计：按照后验概率密度函数
Figure FDA0003534084160000052
计算风速风向状态的估计值，即风速风向的状态值；
Xk为k时刻风速风向状态向量，
Figure FDA0003534084160000053
为第i个粒子通过迭代过程经过前k个时刻的观测得到的观测向量；
Figure FDA0003534084160000054
为Xk的狄拉克函数；
⑤预测：根据X(n+1)＝FcvX(n)+Gcvv(n)预测下一时刻的风速风向状态；
⑥k＝k+1，重复步骤②至步骤⑥，当k大于观测长度，终止预测。
CN202110047680.2A 2021-01-14 2021-01-14 基于粒子滤波算法的风速风向预测方法 Active CN112785052B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110047680.2A CN112785052B (zh) 2021-01-14 2021-01-14 基于粒子滤波算法的风速风向预测方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110047680.2A CN112785052B (zh) 2021-01-14 2021-01-14 基于粒子滤波算法的风速风向预测方法

Publications (2)

Publication Number Publication Date
CN112785052A CN112785052A (zh) 2021-05-11
CN112785052B true CN112785052B (zh) 2022-04-19

Family

ID=75755931

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110047680.2A Active CN112785052B (zh) 2021-01-14 2021-01-14 基于粒子滤波算法的风速风向预测方法

Country Status (1)

Country Link
CN (1) CN112785052B (zh)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113901701A (zh) * 2021-09-30 2022-01-07 南通职业大学 一种基于改进粒子滤波的非线性结构动载荷识别方法
CN114397474B (zh) * 2022-01-17 2022-11-08 吉林大学 基于fcn-mlp的弧形超声传感阵列风参数测量方法
CN116991841B (zh) * 2023-09-25 2023-12-19 温州市工业与信息技术发展有限公司 一种用于混合风数据模型的数据智能清洗方法

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104897925A (zh) * 2015-06-24 2015-09-09 吉林大学 超声波风速风向测量装置及测量方法
CN106621614A (zh) * 2016-12-30 2017-05-10 上海高盾科技发展有限公司 一种双向低风阻颗粒物滤网及其制备方法
CN106778846A (zh) * 2016-12-02 2017-05-31 华北电力大学 一种基于支持向量机的风电场风速预测方法
CN107576819A (zh) * 2017-08-29 2018-01-12 吉林大学 一种测量风速和风向的方法及系统
CN107765347A (zh) * 2017-06-29 2018-03-06 河海大学 一种高斯过程回归和粒子滤波的短期风速预测方法
CN108169511A (zh) * 2018-01-11 2018-06-15 吉林大学 三维空间来风的风速测量系统及方法
CN110543929A (zh) * 2019-08-29 2019-12-06 华北电力大学（保定） 一种基于Lorenz系统的风速区间预测方法及系统

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104897925A (zh) * 2015-06-24 2015-09-09 吉林大学 超声波风速风向测量装置及测量方法
CN106778846A (zh) * 2016-12-02 2017-05-31 华北电力大学 一种基于支持向量机的风电场风速预测方法
CN106621614A (zh) * 2016-12-30 2017-05-10 上海高盾科技发展有限公司 一种双向低风阻颗粒物滤网及其制备方法
CN107765347A (zh) * 2017-06-29 2018-03-06 河海大学 一种高斯过程回归和粒子滤波的短期风速预测方法
CN107576819A (zh) * 2017-08-29 2018-01-12 吉林大学 一种测量风速和风向的方法及系统
CN108169511A (zh) * 2018-01-11 2018-06-15 吉林大学 三维空间来风的风速测量系统及方法
CN110543929A (zh) * 2019-08-29 2019-12-06 华北电力大学（保定） 一种基于Lorenz系统的风速区间预测方法及系统

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
Modified Particle Filtering Algorithm for Single Acoustic Vector Sensor DOA Tracking;Xinbo Li等;《Sensors》;20151016;第15卷(第10期);26198-26211 *
基于波束形成算法风速风向的测量研究;西继东;《中国优秀硕士学位论文全文数据库-基础科学辑》;20180115(第1期);A009-43 *

Also Published As

Publication number Publication date
CN112785052A (zh) 2021-05-11

Similar Documents

Publication Publication Date Title
CN112785052B (zh) 基于粒子滤波算法的风速风向预测方法
CN106842128B (zh) 运动目标的声学跟踪方法及装置
CN108594190B (zh) 一种高分辨海杂波的仿真方法
CN106646417A (zh) 广义帕累托分布参数的迭代最大似然估计方法
CN111638521A (zh) 一种成像高度计遥感数据的海面风速反演方法
CN116087908B (zh) 一种基于协同作业的雷达高精度物位计测量方法
CN114330163B (zh) 高频地波超视距雷达台风-电离层扰动动力学模型建模方法
CN107064893B (zh) 基于对数矩的广义帕累托分布参数估计方法
CN115356702A (zh) 一种基于海上多源雷达回波的蒸发波导深度学习联合反演方法
CN117309079B (zh) 基于时差法的超声飞渡时间测量方法、装置、设备及介质
CN113359130B (zh) 一种低散射高速运动目标的探测方法
CN116794649A (zh) 基于波形选择的杂波中机动目标跟踪方法
CN106443623B (zh) 一种天波超视距雷达目标与电离层参数联合估计方法
Khurjekar et al. Uncertainty aware deep neural network for multistatic localization with application to ultrasonic structural health monitoring
CN110851788B (zh) 基于神经网络的超声背散射零差k模型参数估算方法
CN113514809B (zh) 矢量脱靶量参数的测量方法、装置、电子设备和存储介质
Li et al. Ocean Surface Wind Direction Inversion from HFSWR Based on BP Neural Network
CN114035172B (zh) 基于粒子滤波的雷达回波资料反演大气波导的方法
CN116911082B (zh) 基于降水雷达和同化资料的降水粒子质量和数量估算方法
Bychkova et al. Ultrasonic multipath monitoring of turbulent gas flows: pulse signals correlation processing
CN116577734B (zh) 基于先验知识的机载雷达精细化杂波仿真方法与装置
CN117907997A (zh) 一种基于改进粒子滤波的声源跟踪与环境参数反演方法
Vasyukov et al. Investigation of principles of simulation of space-time processing of wideband signals
CN118050740B (zh) 一种调频连续波激光雷达信号分析方法及系统
CN114580615B (zh) 一种基于神经网络的分布式小平台水下纯方位定位方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant