CN112257186A - 针对小型四旋翼飞行器气动参数的时域辨识方法 - Google Patents

针对小型四旋翼飞行器气动参数的时域辨识方法 Download PDF

Info

Publication number
CN112257186A
CN112257186A CN202011255237.6A CN202011255237A CN112257186A CN 112257186 A CN112257186 A CN 112257186A CN 202011255237 A CN202011255237 A CN 202011255237A CN 112257186 A CN112257186 A CN 112257186A
Authority
CN
China
Prior art keywords
small
rotor aircraft
time domain
identification method
data
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.)
Granted
Application number
CN202011255237.6A
Other languages
English (en)
Other versions
CN112257186B (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.)
Beijing Institute of Technology BIT
Original Assignee
Beijing Institute of Technology BIT
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 Beijing Institute of Technology BIT filed Critical Beijing Institute of Technology BIT
Priority to CN202011255237.6A priority Critical patent/CN112257186B/zh
Publication of CN112257186A publication Critical patent/CN112257186A/zh
Application granted granted Critical
Publication of CN112257186B publication Critical patent/CN112257186B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/10—Geometric CAD
    • G06F30/15—Vehicle, aircraft or watercraft design
    • 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
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Geometry (AREA)
  • General Physics & Mathematics (AREA)
  • Evolutionary Computation (AREA)
  • General Engineering & Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Automation & Control Theory (AREA)
  • Aviation & Aerospace Engineering (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Navigation (AREA)

Abstract

本发明公开了一种针对小型四旋翼飞行器气动参数的时域辨识方法，包括以下步骤：模型结构确定、输入信号设计、通过飞行试验完成数据采集、数据预处理、模型气动参数估计。本发明所述的针对小型四旋翼飞行器气动参数的时域辨识方法，将适用于多采样频率的RTS平滑扩展卡尔曼滤波算法与基于最大似然的输出误差法相结合，解决了现有时域方法无法有效应用于小型四旋翼气动参数辨识的问题具有成本低、周期短、精度高等诸多优点。

Description

针对小型四旋翼飞行器气动参数的时域辨识方法
技术领域
本发明涉及一种气动参数的时域辨识方法，具体涉及一种针对小型四旋翼飞行器气动参数的时域辨识方法，属于无人机控制领域。
背景技术
近年来，小型四旋翼飞行器因具有目标小、造价低、重量轻、隐蔽性好等优势，能更好的适应复杂多变的环境，广泛应用于军民产业之中。
随着微小型无人飞行器的广泛应用，对其飞行性能的要求越来越高，对其控制算法的研究也越来越多。控制算法的研究过程中，无论是否是基于模型的控制方法，其控制律的开发、仿真与验证都离不开准确的飞机本体动力学模型。
然而，现有的小型四旋翼飞行器动力学模型中气动参数大都采用机理建模方法和风洞建模方法得到。风动建模方法需要进行风洞试验，成本非常昂贵。机理建模过程中对无人机各个部分的组件进行基于物理原的建模，需要进行大量的台架实验来调节各部件的建模增益系数以匹配飞行数据，方法耗时很长，且由于部分过于理想化的建模假设导致飞行器动力学模型精确性和实用性降低。
相比之下，时域辨识方法具有成本低、周期短、准确性较高等优点，十分契合小型四旋翼飞行器研究迭代快速的特点。但是，时域气动参数辨识方法应用于动力学特性不稳定且机载传感器精度精度有限的小型四旋翼飞行器时，易发生数值积分发散，极易受到量测噪声和系统误差影响，导致辨识精度下降，无法达到有效的气动参数辨识。
因此，亟待设计一种针对小型四旋翼飞行器的气动参数辨识方法，解决时域气动参数辨识法无法有效应用于小型四旋翼飞行器的问题，从而获得一种成本低、周期短、准确率高的小型四旋翼飞行器动力学模型。
发明内容
为了克服上述问题，本发明人进行了锐意研究，设计出一种针对小型四旋翼飞行器时域气动参数辨识方法，通过将时域方法结合RTS平滑扩展卡尔曼滤波，能够有效解决上文提到的现有时域方法存在的问题。
具体地，所述方法，包括以下步骤：
S1、模型结构确定；
S2、输入信号设计；
S3、通过飞行试验完成数据采集；
S4、数据预处理；
S5、模型气动参数估计。
在步骤S1中，所述模型为小型四旋翼飞行器状态空间模型。
进一步地，所述模型通过下述子步骤获得：
S11、获得飞行器动力学模型；
S12、获得小型四旋翼飞行器状态空间模型。
在步骤S2中，所述输入信号采用3211多阶跃信号的方式输入，在3211多阶跃信号的信号中，比例为1的脉宽对应于占主导地位的系统预期固有频率满足以下公式:
Figure BDA0002772909460000021
其中，fn为占主导地位的系统期望固有频率固。
在步骤S3中，所述飞行试验是在小型四旋翼的悬停状态下进行的，小型四旋翼悬停状态下的性能对于四旋翼的飞行能力来说是最为关键的，飞行试验采用悬停状态下进行，能够获得最具代表性的试验数据。
在步骤S4中，所述数据预处理包括飞行数据兼容性分析和传感器位置校正。
进一步地，飞行数据兼容性分析为通过对飞行试验所得输出信号进行RTS平滑-扩展卡尔曼滤波实现的，RTS平滑-扩展卡尔曼滤波的状态方程和量测方程可进行如下表示：
Figure BDA0002772909460000031
其中w,v为不相关的零均值白噪声，从测量数据中获得，两者的方差分别为Q,R。
优选地，所述输出信号进行RTS平滑-扩展卡尔曼滤波包括前向卡尔曼滤波和后向卡尔曼滤波过程。
根据本发明，在步骤S5中，所述模型气动参数估计通过对步骤S4中预处理后的数据进行输出误差法辨识获得。
在一个优选的实施方式中，所述输出误差法辨识包括如下子步骤：
S51、求解负对数似然函数获得，获得优化目标函数；
S52、求解优化步长；
S53、根据步长进行优化；
S54、以松弛优化的思想交替调整直至达到收敛。
本发明所具有的有益效果包括：
(1)根据本发明提供的针对小型四旋翼飞行器气动参数的时域辨识方法，将适用于多采样频率的RTS平滑扩展卡尔曼滤波算法与基于最大似然的输出误差法相结合，实现了时域气动参数辨识在小型四旋翼飞行器上的高精度应用；
(2)根据本发明提供的针对小型四旋翼飞行器气动参数的时域辨识方法，通过量测数据对输入信号进行了设计，可以得到重点频率附近的输入信号并在短时间内获取包含丰富信息的数据
(3)根据本发明提供的针对小型四旋翼飞行器气动参数的时域辨识方法，无需高精度传感器及昂贵复杂的建模手段即可得到四旋翼的气动模型，成本低、周期短，适应小型四旋翼飞行器的快速迭代研究。
附图说明
图1示出根据本发明一种优选实施方式的针对小型四旋翼飞行器气动参数的时域辨识方法流程示意图；
图2示出一种优选实施方式中针对小型四旋翼飞行器气动参数的时域辨识方法示意图；
图3示出RTS平滑扩展卡尔曼滤波流程示意图；
图4示出实验例1中姿态角拟合效果图。
具体实施方式
下面通过附图和实施例对本发明进一步详细说明。通过这些说明，本发明的特点和优点将变得更为清楚明确。
在这里专用的词“示例性”意为“用作例子、实施例或说明性”。这里作为“示例性”所说明的任何实施例不必解释为优于或好于其它实施例。尽管在附图中示出了实施例的各种方面，但是除非特别指出，不必按比例绘制附图。
本发明提供了一种针对小型四旋翼飞行器气动参数的时域辨识方法，如图1所示，包括以下步骤：
S1、模型结构确定；
S2、输入信号设计；
S3、通过飞行试验完成数据采集；S4、数据预处理；
S5、模型气动参数估计。
在步骤S1中，所述模型为小型四旋翼飞行器状态空间模型。具体地，通过下述子步骤获得：
S11、获得飞行器动力学模型；
S12、获得小型四旋翼飞行器状态空间模型；
在步骤S11中，所述飞行器动力学模型，可以通过如下形式的飞行器六自由度运动方程组来表示，分别为：
平行动力学方程组
Figure BDA0002772909460000051
旋转动力学方程组
Figure BDA0002772909460000052
平移运动学方程组
Figure BDA0002772909460000053
旋转运动学方程组
Figure BDA0002772909460000054
其中，u,v,w分别机体系下飞行器质心速度在机体轴方向的速度分量；
Figure BDA0002772909460000061
机体系下分别为飞行器质心加速度在机体轴方向的分量；
p,q,r代表机体坐标系相对地面坐标系转动的角速度在机体坐标系各轴上的分量；
φ,θ,ψ分别表示滚转角、俯仰角和偏航角三个欧拉角；
Figure BDA0002772909460000062
为机体系下飞行器滚转、俯仰、偏航方向的角速度；
L,M,N分别表示俯仰力矩、滚转力矩、偏航力矩；
Ix,Iy,Iz为飞行器对机体坐标系中各坐标轴的转动惯量；
ax,ay,az为机体系相对于地面坐标系加速度沿地面坐标系XYZ轴的分量；
Figure BDA0002772909460000063
为地面坐标系下飞行器质心速度沿地面坐标系坐标轴的分量。
在步骤S12中，所述小型四旋翼飞行器状态空间模型为通过在飞行器动力学模型中增加输入信号影响后获得，具体的，将飞行器六自由度运动方程组中增加输入信号影响，并进行小角度线性化，获得能够用于气动参数辨识的小型四旋翼飞行器状态空间模型，所述状态空间模型可以通过下式表示：
Figure BDA0002772909460000071
其中，Xu,Yv,Zw,Lv,Lp,Mu,Mq,Nr,
Figure BDA0002772909460000072
代表动力学特性相关的力或力矩导数，均为待辨识参数；
MCthr,MClat,MClon,MCdir分别为垂向通道，偏航通道，俯仰通道和滚转通道的输入信号，其具体数值在步骤S3中数据采集时获取；下标thr代表垂向、lat代表偏航向、lon代表俯仰向、dir代表滚转向。
为了获得精确的可靠的结果,输入信号要按照被研究系统的特性在某些限制范围内调整。在有非线性反应的情况下,输入信号设计特别重要，在这种情况下,输入信号的幅值要适应特定的非线性,也就是说输入信号要有几种幅值,且具有不同的比率。
在步骤S2中，所述输入信号采用3211多阶跃信号的方式输入，所述3211多阶跃信号是一种已被证实用于识别研究的较好的单输入信号，例如其在[聂雪媛,刘中玉,杨国伟.基于CFD气动力辨识模型的气动弹性数值计算[J].振动与冲击,2014,33(020):20-25.]中已被成熟应用。
这种输入信号一方面比较简单，另一方面，其又能激发相当宽的频率范围，无需按照被研究系统的特征频率精确调整转换间隔时间。在本发明中，通过3211多阶跃信号的方式，不仅能够在短时间获得较为丰富的数据信息，还起到辨识激励的作用，产生频率在感兴趣频率附近的信号，在满足实际约束前提下充分地激励系统的要求。
进一步地，所述3211多阶跃信号为信号带宽为3:2:2:1的交替脉冲组成；在3211多阶跃信号中，比例为1的脉宽对应于占主导地位的系统预期固有频率满足以下公式:
Figure BDA0002772909460000081
其中，fn为占主导地位的系统期望固有频率固。
优选地，发明人根据经过多次试验确定，fn为1Hz，以1Hz的主频进行各测量数据的3211多阶跃输入信号设计，能够获得更好的激励辨识效果。
在步骤S3中，所述飞行试验是在小型四旋翼的悬停状态下进行的，小型四旋翼悬停状态下的性能对于四旋翼的飞行能力来说是最为关键的，飞行试验采用悬停状态下进行，能够获得最具代表性的试验数据。
进一步地，所述数据的采集，如图2所示，包括采集输入信号和输出信号，所述输入信号包括飞手输入信号和步骤S2中涉及的3211多阶跃信号，所述输出信号包括无人机旋翼电机输入信号、欧拉角、角速率、线性加速度、航向、气压高度、高度和地速。
其中，所述飞手输入信号、电机输入信号和欧拉角通过飞控进行采集，所述角速率和线性加速度通过IMU进行采集，所述航向和气压高度分别通过地磁极和气压计进行采集，所述高度和地速通过GPS进行采集。在一个优选的实施方式中，对飞手输入信号、电机输入信号、欧拉角、角速率、线性加速度、航向和气压高度的采集采用较高的频率进行，对高度和地速的采集采用较低的频率进行。
进一步地，所述较高的频率是指大于20Hz的频率，所述较低频率是指小于20Hz的频率，优选地，采集频率选取被采样传感器主频率的5倍。由于不同的传感器的主频不同，采用不同的采样频率能够获得更为准确的数据。
在步骤S4中，所述数据预处理包括飞行数据兼容性分析和传感器位置校正。
其中，飞行数据兼容性分析为通过对飞行试验所得输出信号进行RTS平滑-扩展卡尔曼滤波实现的。
根据本发明，RTS平滑-扩展卡尔曼滤波的状态方程和量测方程可进行如下表示：
Figure BDA0002772909460000091
其中w,v为不相关的零均值白噪声，与传感器硬件相关，可从测量数据中获得，两者的方差分别为Q,R，Q,R根据w,v值通过最优滤波方法获得，具体获得方法可参考E.A.Morelli,Estimating noise characteristics from flight test data usingopti-mal Fourier smoothing,J.Aircr.32(1995)689–695，在本发明中不再赘述；
状态变量xa为：
xa＝[u v w φ θ ψ h]T (八)量测变量za为：
za＝[Vxm Vym Vzm hm φm θm ψm]T (九)
Vxm,Vym,Vzm为地面坐标系下飞行器沿坐标轴的速度分量，h为传感器测量得到的高度，hm为相对地面坐标系的绝对高度；
状态方程中的函数f(xa)为：
Figure BDA0002772909460000101
量测方程函数h(xa)为：
Figure BDA0002772909460000102
为求解扩展卡尔曼滤波的增益K及RTS平滑中的增益Ks，对f(xa)和h(xa)按照以下形式求偏导得到线性系统矩阵F(xa)，线性量测矩阵H(xa)：
Figure BDA0002772909460000111
其中，F(xa),H(xa)矩阵中的
Figure BDA0002772909460000112
分别代表f(xa)矩阵，h(xa)矩阵元素对xa矩阵元素求偏导，下标w代表矩阵的第w行，w∈[1，7]。
由线性系统矩阵F(xa)，线性量测矩阵H(xa)离散后可得状态转移矩阵Φ(xa)和量测转移矩阵ξ(xa)：
Figure BDA0002772909460000113
其中，Ts为离散步长，取Ts＝0.002s，I表示单位矩阵。
根据状态转移矩阵和量测转移矩阵可以得到与噪声协方差矩阵相关的常值矩阵R(k),Q(k)：
Figure BDA0002772909460000114
其中，t为积分变量，无实际意义。
所述输出信号进行RTS平滑-扩展卡尔曼滤波包括前向卡尔曼滤波和后向卡尔曼滤波过程，如图3所示。
具体地，所述前向卡尔曼滤波通过下述步骤进行：
S411、前向滤波初值选取，根据采集到数据的前10组取平均值确定状态变量的初值xa0，状态变量协方差矩阵的初值P0本领域技术人员可根据经验选择，例如状态变量协方差矩阵的初值P0可以设置为：
Figure BDA0002772909460000121
S412、通过下式进行状态量的更新：
Figure BDA0002772909460000122
其中，
Figure BDA0002772909460000123
为第k-1组数据对应的状态量估计值；
Figure BDA0002772909460000124
为第k组数据对应的状态量预测值，且
Figure BDA0002772909460000125
S413、通过下式对状态变量协方差矩阵更新：
Figure BDA0002772909460000126
其中，
Figure BDA0002772909460000127
表示根据第k-1组数据的估计值
Figure BDA0002772909460000128
进一步地，通过式十四计算得到的用于第k步的离散后的状态转移矩阵。
S414、通过下式计算扩展卡尔曼滤波增益K：
Figure BDA0002772909460000129
其中
Figure BDA00027729094600001210
表示用于第k步的线性状态转移矩阵，其是根据第k-1组数据的估计值
Figure BDA00027729094600001211
通过式十二计算得到的。
S415、根据测量到的量测变量数值，与根据状态量预测值
Figure BDA00027729094600001212
计算得到的量测变量值做差乘以卡尔曼增益作为修正量，对
Figure BDA00027729094600001213
进行修正，即可得到了状态变量的估计值
Figure BDA00027729094600001214
对状态变量进行修正可通过下式表示：
Figure BDA00027729094600001215
S416、通过下式对状态变量协方差矩阵修正：
Figure BDA00027729094600001216
在步骤S412～S416中，k代表对第k组数据，带下标k或(k)的变量代表与第k组数据相关的变量，N1为总的样本数据组数；将步骤S415和S416中获得的
Figure BDA0002772909460000131
和
Figure BDA0002772909460000132
再应用到S412中，从k＝1开始循环计算直至k＝N1，即可得到最终的前向滤波估计值
Figure BDA0002772909460000133
所述后向卡尔曼滤波通过下述步骤进行：
S421、后向滤波初值选取，后向滤波状态变量的估计初值与状态变量协方差矩阵初值为：
Figure BDA0002772909460000134
S422、通过下式获得平滑增益Ks矩阵
Figure BDA0002772909460000135
S423、获得平滑后的滤波结果：
具体地，第k组数据对应的RTS平滑扩展卡尔曼滤波后状态量的估计值
Figure BDA0002772909460000136
为：
Figure BDA0002772909460000137
第k个数据点处RTS平滑扩展卡尔曼滤波后的状态量的协方差矩阵
Figure BDA0002772909460000138
为：
Figure BDA0002772909460000139
其中，变量
Figure BDA00027729094600001310
分别为前向扩展卡尔曼滤波中得到的状态估计量，一步状态估计量，状态量协方差矩阵估计值，一步状态量协方差矩阵估计值，离散的一步状态转移矩阵。
在本发明中，RTS平滑扩展卡尔曼滤波后状态量的估计值
Figure BDA00027729094600001311
即为更为精确的飞行试验输出数据。在本发明中，通过RTS平滑-扩展卡尔曼滤波，解决了采样频率不同导致的同一时间各数据可用性不一致的问题，从而减少了扩展卡尔曼滤波的内在误差，从而得到了更为准确的数据。
所述传感器位置校正是为了解决特定传感器，如IMU等，因测量需求安装在机身特定位置的情况，或其它原因导致原点无法对齐的问题，通过传感器位置误差校正进行纠正。
具体校正公式本领域技术人员可根据实际需要自由设置，例如通过下式进行加速度位置误差校正：
Figure BDA0002772909460000141
其中，[xa ya za]T表示传感器相对质心的位置矢量，g0为重力加速度取9.8m/s2
在步骤S5中，所述模型气动参数估计通过对步骤S4中预处理后的数据进行输出误差法辨识获得。
具体地，将步骤S1中获得的状态空间模型拆分为四个通道，各通道状态空间模型如下:
俯仰通道
Figure BDA0002772909460000142
偏航通道
Figure BDA0002772909460000143
滚转通道
Figure BDA0002772909460000144
垂向通道
Figure BDA0002772909460000151
在本发明中，通过对四个通道分别进行辨识完成对整个状态空间模型的辨识。
下面以俯仰通道为例，详述具体的辨识方法，其它三个通道的辨识与俯仰通道辨识方法相同。
俯仰通道状态方程、输出方程和量测方程如下：
Figure BDA0002772909460000152
其中，vb为零均值白噪声，其方差为Rb，xb为状态量，yb为输出量，zb为量测量。
进一步地，根据式五可得：
Figure BDA0002772909460000153
其中，下标g代表量测得到的带噪声数据。
更进一步地，状态方程和输出方程的系数矩阵可以表示为：
Figure BDA0002772909460000154
待辨识参数向量可以表示为：
Figure BDA0002772909460000161
残差η表达式为：
η(ib)＝zb(ib)-yb(ib)b＝1,2,...,N2 (三十二)
其中，ib代表不同组的辨识数据，N2代表总的辨识数据组数。
根据本发明，所述输出误差法辨识包括如下子步骤：
S51、求解负对数似然函数，获得优化目标函数；
通过下式求解负对数似然函数方程：
Figure BDA0002772909460000162
其中，
Figure BDA0002772909460000163
表示辨识参数为μ,量测方程为z的负对数似然函数，
Figure BDA0002772909460000164
代表取括号内量的似然函数，n0代表输出变量的个数。σ为残差的协方差矩阵，估计值可以表示为
Figure BDA0002772909460000165
|σ|表示矩阵的对角线元素平方的加和。
由于负对数似然函数方程中只有
Figure BDA0002772909460000166
这部分为可变量，故优化目标函数可表示为：
Figure BDA0002772909460000167
S52、求解优化步长；
依据修正牛顿-拉夫逊算法，忽略二阶梯度矩阵中的二阶偏导数项，得到对应第ik个待辨识参数向量的步长表达式为：
Figure BDA0002772909460000171
其中
Figure BDA0002772909460000172
代表第ik个待辨识参数向量，N3＝N2+1。
S53、根据步长进行优化；
所述优化过程如下：
Figure BDA0002772909460000173
进一步地，优化过程的初值通过最小二乘方法根据输出数据获得状态变量和状态变量协方差矩阵获得的。
在一个优选的实施方式中，优化过程的初值为：
μ0＝[0 0 0 0]T
Figure BDA0002772909460000174
S54、以松弛优化的思想交替调整直至达到收敛。
在本发明中，以松弛优化的思想最小化负对数似然代价函数，将μ和σ交替调整最终得到效果最好的未知参数集，具体地，先将μ固定，求解得到当前步的
Figure BDA0002772909460000175
再将σ值固定为估计得到的
Figure BDA0002772909460000176
依据代价函数求解当前步骤的
Figure BDA0002772909460000177
最小值，如此重复直到满足以下收敛条件中的一条或多条:
收敛条件一：μ变化量中的元素绝对值小于1.0×10-5：
Figure BDA0002772909460000178
其中，jn代表待辨识参数向量中的第jn个元素，np代表带辨识参数向量中元素的个数。
收敛条件二：代价函数
Figure BDA0002772909460000181
的变化对于连续迭代小于0.001：
Figure BDA0002772909460000182
收敛条件三：代价函数梯度中元素的绝对值小于0.05：
Figure BDA0002772909460000183
收敛条件四：
Figure BDA0002772909460000184
中对角线元素相对变化小于0.05：
Figure BDA0002772909460000185
其中
Figure BDA0002772909460000186
代表
Figure BDA0002772909460000187
第l个对角线元素。
当收敛条件满足后，所得到的θ即为最终辨识出的模型气动参数，以此气动参数作为动力学模型的参数，即可获得一种准确率高的小型四旋翼飞行器动力学模型。
实施例
实施例1
对一个小型四旋翼飞行器动力学模型进行辨识，小型四旋翼状态空间模型表示为：
Figure BDA0002772909460000191
输入信号设计为3211多阶跃信号的方式输入，3211多阶跃信号中比例为1的脉宽对应于占主导地位的系统预期固有频率满足以下公式:
Figure BDA0002772909460000192
其中fn为1Hz。
在小型四旋翼悬停状态下进行飞行试验完成数据采集，获取输入信号和输出信号，其中对飞手输入信号、电机输入信号、欧拉角、角速率、线性加速度、航向和气压高度采用100Hz的频率采集，对高度和地速采用5Hz的频率采集。。
完成数据采集后，对传感器位置误差进行校正和对输出信号进行RTS平滑扩展卡尔曼滤波，其中RTS平滑扩展卡尔曼滤波前向卡尔曼滤波和后向卡尔曼滤波过程，在前向卡尔曼滤波过程中，根据采集到输出信号中数据的前10组取平均值确定状态变量的初值xa0，状态变量协方差矩阵的初值P0设置为：
Figure BDA0002772909460000201
通过下式进行状态量的更新：
Figure BDA0002772909460000202
通过下式对状态变量协方差矩阵更新：
Figure BDA0002772909460000203
通过下式计算扩展卡尔曼滤波增益K：
Figure BDA0002772909460000204
根据测量到的量测变量数值，与根据状态量预测值
Figure BDA0002772909460000205
计算得到的量测变量值做差乘以卡尔曼增益作为修正量，对
Figure BDA0002772909460000206
进行修正，得到状态变量的估计值
Figure BDA0002772909460000207
通过下式对状态变量进行修正：
Figure BDA0002772909460000208
将获得的
Figure BDA0002772909460000209
和
Figure BDA00027729094600002010
再次代入式十五～式十五，从k＝1开始循环计算直至k＝N1，得到最终的前向滤波估计值
Figure BDA00027729094600002011
在后向卡尔曼滤波过程中，后向滤波状态变量的估计初值与状态变量协方差矩阵初值选为：
Figure BDA00027729094600002012
通过下式获得平滑增益Ks矩阵
Figure BDA00027729094600002013
通过下式获得平滑后的滤波结果：
Figure BDA00027729094600002014
Figure BDA0002772909460000211
进而获得精确的飞行试验输出数据
Figure BDA0002772909460000212
对精确的飞行试验输出数据进行输出误差法辨识，将状态空间模型拆分为四个通道，对四个通道分别进行辨识，以俯仰通道为例，俯仰通道状态方程、输出方程和量测方程如下：
Figure BDA0002772909460000213
通过下式求解负对数似然函数方程：
Figure BDA0002772909460000214
获得优化目标函数：
Figure BDA0002772909460000215
通过下式求解优化目标函数的优化步长：
Figure BDA0002772909460000216
根据步长按照下式进行优化：
Figure BDA0002772909460000217
优化过程的初值为：
μ0＝[0 0 0 0]T
Figure BDA0002772909460000218
以松弛优化的思想交替调整直至达到收敛，当收敛条件满足后，所得到的θ即为最终辨识出的模型气动参数，以此气动参数作为动力学模型的参数，获得小型四旋翼飞行器动力学模型。
对比例1
进行与实施例1相同的辨识，区别在于，在数据的预处理过程中，不进行RTS平滑扩展卡尔曼滤波，直接将获得的数据采用松弛原理的输出误差法处理。
实验例1
将四旋翼飞行器的实际运动情况与实施例1和对比例1中获得的四旋翼飞行器动力学模型模拟值进行比对，其姿态角拟合效果图结果如图4所示。
从图中可以看出,实施例1中获得的四旋翼飞行器动力学模型与四旋翼飞行器的实际运动情况吻合度高，通过四旋翼飞行器动力学模型可以准确的反应小型四旋翼飞行器的实际状况。
在本发明的描述中，需要说明的是，术语“上”、“下”、“内”、“外”、“前”、“后”等指示的方位或位置关系为基于本发明工作状态下的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。此外，术语“第一”、“第二”、“第三”、“第四”仅用于描述目的，而不能理解为指示或暗示相对重要性。
在本发明的描述中，需要说明的是，除非另有明确的规定和限定，术语“安装”“相连”“连接”应作广义理解，例如，可以是固定连接，也可以是可拆卸连接，或一体的连接普通；可以是机械连接，也可以是电连接；可以是直接连接，也可以通过中间媒介间接连接，可以是两个元件内部的连通。对于本领域的普通技术人员而言，可以具体情况理解上述术语在本发明中的具体含义。
以上结合了优选的实施方式对本发明进行了说明，不过这些实施方式仅是范例性的，仅起到说明性的作用。在此基础上，可以对本发明进行多种替换和改进，这些均落入本发明的保护范围内。

Claims (10)

1.一种针对小型四旋翼飞行器气动参数的时域辨识方法，包括以下步骤：
S1、模型结构确定；
S2、输入信号设计；
S3、通过飞行试验完成数据采集；
S4、数据预处理；
S5、模型气动参数估计。
2.根据权利要求1所述的针对小型四旋翼飞行器气动参数的时域辨识方法，其特征在于，
在步骤S1中，所述模型为小型四旋翼飞行器状态空间模型。
3.根据权利要求2所述的针对小型四旋翼飞行器气动参数的时域辨识方法，其特征在于，
所述模型通过下述子步骤获得：
S11、获得飞行器动力学模型；
S12、获得小型四旋翼飞行器状态空间模型。
4.根据权利要求1所述的针对小型四旋翼飞行器气动参数的时域辨识方法，其特征在于，
在步骤S2中，所述输入信号采用3211多阶跃信号的方式输入，在3211多阶跃信号的信号中，比例为1的脉宽对应于占主导地位的系统预期固有频率满足以下公式:
Figure FDA0002772909450000011
其中，fn为占主导地位的系统期望固有频率固。
5.根据权利要求1所述的针对小型四旋翼飞行器气动参数的时域辨识方法，其特征在于，
在步骤S3中，所述飞行试验是在小型四旋翼的悬停状态下进行的，小型四旋翼悬停状态下的性能对于四旋翼的飞行能力来说是最为关键的，飞行试验采用悬停状态下进行，能够获得最具代表性的试验数据。
6.根据权利要求1所述的针对小型四旋翼飞行器气动参数的时域辨识方法，其特征在于，
在步骤S4中，所述数据预处理包括飞行数据兼容性分析和传感器位置校正。
7.根据权利要求1所述的针对小型四旋翼飞行器气动参数的时域辨识方法，其特征在于，
飞行数据兼容性分析为通过对飞行试验所得输出信号进行RTS平滑-扩展卡尔曼滤波实现的，RTS平滑-扩展卡尔曼滤波的状态方程和量测方程可进行如下表示：
Figure FDA0002772909450000021
其中w,v为不相关的零均值白噪声，从测量数据中获得，通过w,v可解算出两者的方差为Q,R。
8.根据权利要求1所述的针对小型四旋翼飞行器气动参数的时域辨识方法，其特征在于，
所述输出信号进行RTS平滑-扩展卡尔曼滤波包括前向卡尔曼滤波和后向卡尔曼滤波过程。
9.根据权利要求1所述的针对小型四旋翼飞行器气动参数的时域辨识方法，其特征在于，
在步骤S5中，所述模型气动参数估计通过对步骤S4中预处理后的数据进行输出误差法辨识获得。
10.根据权利要求9所述的针对小型四旋翼飞行器气动参数的时域辨识方法，其特征在于，
所述输出误差法辨识包括如下子步骤：
S51、求解负对数似然函数获得，获得优化目标函数；
S52、求解优化步长；
S53、根据步长进行优化；
S54、以松弛优化的思想交替调整直至达到收敛。
CN202011255237.6A 2020-11-11 2020-11-11 针对小型四旋翼飞行器气动参数的时域辨识方法 Active CN112257186B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011255237.6A CN112257186B (zh) 2020-11-11 2020-11-11 针对小型四旋翼飞行器气动参数的时域辨识方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011255237.6A CN112257186B (zh) 2020-11-11 2020-11-11 针对小型四旋翼飞行器气动参数的时域辨识方法

Publications (2)

Publication Number Publication Date
CN112257186A true CN112257186A (zh) 2021-01-22
CN112257186B CN112257186B (zh) 2023-05-02

Family

ID=74265460

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011255237.6A Active CN112257186B (zh) 2020-11-11 2020-11-11 针对小型四旋翼飞行器气动参数的时域辨识方法

Country Status (1)

Country Link
CN (1) CN112257186B (zh)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112949216A (zh) * 2021-02-03 2021-06-11 中国空气动力研究与发展中心高速空气动力研究所 一种基于混合性能函数的在线寻峰数据处理方法
CN114266103A (zh) * 2021-09-16 2022-04-01 中国空气动力研究与发展中心计算空气动力研究所 一种飞行器参数和噪声特性在线估计方法及存储介质
CN114527796A (zh) * 2022-02-25 2022-05-24 国网山东省电力公司临沂供电公司 一种无人机仿输电线飞行的方法和系统
CN115309178A (zh) * 2021-05-06 2022-11-08 北京理工大学 一种小型四旋翼无人机前飞实验系统设计及模型辨识方法

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107101636A (zh) * 2017-05-23 2017-08-29 南京航空航天大学 一种使用卡尔曼滤波器辨识多旋翼动力学模型参数的方法
CN110187713A (zh) * 2019-04-12 2019-08-30 浙江大学 一种基于气动参数在线辨识的高超声速飞行器纵向控制方法

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107101636A (zh) * 2017-05-23 2017-08-29 南京航空航天大学 一种使用卡尔曼滤波器辨识多旋翼动力学模型参数的方法
CN110187713A (zh) * 2019-04-12 2019-08-30 浙江大学 一种基于气动参数在线辨识的高超声速飞行器纵向控制方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
李正楠 等: "基于卡尔曼滤波最大似然参数估计的气动参数辨识", 《四川兵工学报》 *
李超: "《中国优秀硕士学位论文全文数据库 (工程科技II辑)》", 15 August 2014 *

Cited By (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112949216A (zh) * 2021-02-03 2021-06-11 中国空气动力研究与发展中心高速空气动力研究所 一种基于混合性能函数的在线寻峰数据处理方法
CN112949216B (zh) * 2021-02-03 2023-07-14 中国空气动力研究与发展中心高速空气动力研究所 一种基于混合性能函数的在线寻峰数据处理方法
CN115309178A (zh) * 2021-05-06 2022-11-08 北京理工大学 一种小型四旋翼无人机前飞实验系统设计及模型辨识方法
CN114266103A (zh) * 2021-09-16 2022-04-01 中国空气动力研究与发展中心计算空气动力研究所 一种飞行器参数和噪声特性在线估计方法及存储介质
CN114266103B (zh) * 2021-09-16 2023-05-19 中国空气动力研究与发展中心计算空气动力研究所 一种飞行器参数和噪声特性在线估计方法及存储介质
CN114527796A (zh) * 2022-02-25 2022-05-24 国网山东省电力公司临沂供电公司 一种无人机仿输电线飞行的方法和系统
CN114527796B (zh) * 2022-02-25 2024-05-28 国网山东省电力公司临沂供电公司 一种无人机仿输电线飞行的方法和系统

Also Published As

Publication number Publication date
CN112257186B (zh) 2023-05-02

Similar Documents

Publication Publication Date Title
CN112257186B (zh) 针对小型四旋翼飞行器气动参数的时域辨识方法
CN108827299B (zh) 一种基于改进四元数二阶互补滤波的飞行器姿态解算方法
Hamel et al. Attitude estimation on SO [3] based on direct inertial measurements
CN111024064B (zh) 一种改进Sage-Husa自适应滤波的SINS/DVL组合导航方法
CN103630137B (zh) 一种用于导航系统的姿态及航向角的校正方法
US6498996B1 (en) Vibration compensation for sensors
CN110081878B (zh) 一种多旋翼无人机的姿态及位置确定方法
CN110398257A (zh) Gps辅助的sins系统快速动基座初始对准方法
CN104764467B (zh) 空天飞行器惯性传感器误差在线自适应标定方法
CN111141313B (zh) 一种提高机载局部相对姿态匹配传递对准精度的方法
CN110553642B (zh) 一种提高惯性制导精度的方法
CN111156987A (zh) 基于残差补偿多速率ckf的惯性/天文组合导航方法
CN109443342A (zh) 新型自适应卡尔曼无人机姿态解算方法
CN112945225A (zh) 基于扩展卡尔曼滤波的姿态解算系统及解算方法
CN110702113B (zh) 基于mems传感器的捷联惯导系统数据预处理和姿态解算的方法
CN102680004A (zh) 一种挠性陀螺位置姿态测量系统pos的标度因数误差标定与补偿方法
CN113175926B (zh) 一种基于运动状态监测的自适应水平姿态测量方法
CN108827288A (zh) 一种基于对偶四元数的降维捷联惯性导航系统初始对准方法及系统
CN108458709A (zh) 基于视觉辅助测量的机载分布式pos数据融合方法和装置
CN112857398A (zh) 一种系泊状态下舰船的快速初始对准方法和装置
CN116007620A (zh) 一种组合导航滤波方法、系统、电子设备及存储介质
CN116182871B (zh) 一种基于二阶混合滤波的海缆探测机器人姿态估计方法
Koehl et al. Wind-disturbance and aerodynamic parameter estimation of an experimental launched micro air vehicle using an EKF-like observer
CN111412887B (zh) 一种基于卡尔曼滤波的攻角、侧滑角辨识方法
Jategaonkar et al. Identification of moderately nonlinear flight mechanics systems withadditive process and measurement noise

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