CN114326814B - 一种无动力飞行器的三维制导系统 - Google Patents

一种无动力飞行器的三维制导系统 Download PDF

Info

Publication number
CN114326814B
CN114326814B CN202111676884.9A CN202111676884A CN114326814B CN 114326814 B CN114326814 B CN 114326814B CN 202111676884 A CN202111676884 A CN 202111676884A CN 114326814 B CN114326814 B CN 114326814B
Authority
CN
China
Prior art keywords
segment
sight
angle
calculating
flight
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
CN202111676884.9A
Other languages
English (en)
Other versions
CN114326814A (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 Aerospace Automatic Control Research Institute
Original Assignee
Beijing Aerospace Automatic Control Research Institute
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 Aerospace Automatic Control Research Institute filed Critical Beijing Aerospace Automatic Control Research Institute
Priority to CN202111676884.9A priority Critical patent/CN114326814B/zh
Publication of CN114326814A publication Critical patent/CN114326814A/zh
Application granted granted Critical
Publication of CN114326814B publication Critical patent/CN114326814B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • 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

  • Traffic Control Systems (AREA)

Abstract

一种无动力飞行器的三维制导系统，包括以下模块：数据获取模块，用于获取当前时刻飞行器的目标距离、视线倾角、视线偏角、速度矢量和速度模值；剩余飞行时间预估模块，用于基于所述视线倾角、视线偏角和速度矢量计算总前置角、前置偏角和前置倾角；基于所述目标距离、视线倾角、视线偏角、速度模值、总前置角、前置偏角和前置倾角，采用基于速度预测的分段迭代法预估剩余飞行时间；制导指令生成模块，用于基于预估剩余飞行时间和期望剩余飞行时间的时间偏差，获取时间偏差控制项，基于当前时刻飞行器的目标距离、视线倾角、视线偏角、速度模值和所述时间偏差控制项生成制导控制指令。

Description

一种无动力飞行器的三维制导系统
技术领域
本发明涉及无动力飞行器制导技术领域，尤其涉及一种无动力飞行器的三维制导系统。
背景技术
飞行器以指定的时间攻击目标，这是实现多飞行器协同饱和攻击的基础。为实现带有时间约束的精确制导，需要较为精确的飞行剩余时间信息以形成控制反馈。现有方法是局域飞行器速度获取变化率已知或可控的前提条件下进行制导的，无动力飞行器在实际飞行中速度变化明显且不受控制，在这种情况下，速度的不确定变化给剩余飞行时间的预测带来了一定的困难，进而使得飞行器的到达时间控制精度低。
发明内容
鉴于上述的分析，本发明实施例旨在提供一种无动力飞行器的三维制导系统，用以解决现有无动力飞行器到达时间控制精度低的问题。
一方面，本发明实施例提供了种无动力飞行器的三维制导系统，包括以下模块：
数据获取模块，用于获取当前时刻飞行器的目标距离、视线倾角、视线偏角、速度矢量和速度模值；
剩余飞行时间预估模块，用于基于所述视线倾角、视线偏角和速度矢量计算总前置角、前置偏角和前置倾角；基于所述目标距离、视线倾角、视线偏角、速度模值、总前置角、前置偏角和前置倾角，采用基于速度预测的分段迭代法预估剩余飞行时间；
制导指令生成模块，用于基于预估剩余飞行时间和期望剩余飞行时间的时间偏差，获取时间偏差控制项，基于当前时刻飞行器的目标距离、视线倾角、视线偏角、速度模值和所述时间偏差控制项生成制导控制指令。
基于上述方案的进一步改进，所述剩余飞行时间预估模块，用于基于所述视线倾角、视线偏角和速度矢量计算总前置角、前置偏角和前置倾角，包括：
根据所述视线倾角和视线偏角计算基准参考坐标系到视线坐标系的方向余弦矩阵；
基于所述方向余弦矩阵计算所述速度矢量在所述视线坐标系中的分量；
根据以下公式计算总前置角：
Figure GDA0004231246460000021
Figure GDA0004231246460000022
η＝arccos(cosηycosηz)
其中，vxL、vyL和vzL分别表示速度矢量在视线坐标系三个坐标轴的分量，ηz表示前置倾角，ηy表示前置偏角，η表示总前置角。
进一步的，所述制导指令生成模块，用于基于当前时刻飞行器的目标距离、视线倾角、视线偏角、速度模值和所述时间偏差控制项生成制导控制指令，包括:
根据公式
Figure GDA0004231246460000023
计算纵向加速度指令；
根据公式
Figure GDA0004231246460000031
计算侧向加速度指令；
其中，Dtc表示目标距离，ηy表示前置偏角，ηz表示前置倾角，KN表示导航比，v表示速度模值，aε表示时间偏差控制项。
进一步地，所述制导指令生成模块，用于基于预估剩余飞行时间和期望剩余飞行时间的时间偏差，获取时间偏差控制项，包括：
根据公式aε＝-KTvsinηycosηyεT计算时间偏差控制项；其中，εT表示时间偏差，ηy表示前置偏角，v表示速度模值，KT表示时间控制系数，aε表示时间偏差控制项。
进一步地，所述剩余飞行时间预估模块基于所述目标距离、视线倾角、视线偏角、速度模值、总前置角、前置偏角和前置倾角，采用基于速度预测的分段迭代法预估剩余飞行时间，包括：
根据所述总前置角判断当前飞行为转弯飞行或直线飞行；
若为转弯飞行，则根据总前置角对转弯段进行分段，采用分段迭代预测转弯段剩余飞行时间；根据直线段起始点的总前置角和目标距离计算剩余直线段航程，根据剩余直线段航程对直线段进行分段，采用分段迭代预测直线段剩余飞行时间；所述转弯段剩余飞行时间与直线段剩余飞行时间的和为预测的总剩余飞行时间；
直线飞行预测模块，用于若为直线飞行，根据所述总前置角和目标距离计算剩余直线段航程，根据剩余直线段航程对直线段进行分段，采用分段迭代预测直线段剩余飞行时间，所述直线段剩余飞行时间为预测的总剩余飞行时间。
进一步地，采用分段迭代预测转弯段剩余飞行时间或采用分段迭代预测直线段剩余飞行时间，包括：
对于每一分段，获取当前分段起始点的状态变量，所述状态变量包括速度模值vi、视线倾角θL,i、视线偏角ψL,i、前置倾角ηz,i、前置偏角ηy,i、地面高度yi；
若当前分段处于转弯段，所述状态变量还包括总前置角ηi和目标距离Dtc,i；若当前分段处于直线段，所述状态变量还包括剩余直线段航程Li；
若当前分段处于转弯段，则根据
Figure GDA0004231246460000041
计算当前分段初始飞行时间预测值
Figure GDA0004231246460000042
若当前分段处于直线段，则根据
Figure GDA0004231246460000043
计算当前分段初始飞行时间预测值
Figure GDA0004231246460000044
其中，Δηi表示当前分段总前置角的变化量，ΔL表示当前分段剩余直线段航程变化量，KN表示导航比；
基于当前分段起始点的状态变量，根据飞行器质心运动的动力学方程计算当前分段的初始速度变化率
Figure GDA0004231246460000045
基于所述初始速度变化率
Figure GDA0004231246460000046
对所述当前分段初始飞行时间预测值
Figure GDA0004231246460000047
进行修正，得到当前分段飞行时间修正值
Figure GDA0004231246460000048
基于当前分段起始点的状态变量和所述当前分段飞行时间修正值
Figure GDA0004231246460000049
预测当前分段结束点的状态变量，当前分段结束点的状态变量为下一分段起始点的状态变量；
所有分段的飞行时间修正值
Figure GDA00042312464600000410
的和作为预测的转弯段或直线段剩余飞行时间。
进一步地，基于当前分段起始点的状态变量，根据飞行器质心运动的动力学方程计算当前分段的初始速度变化率
Figure GDA00042312464600000411
包括：
基于所述视线倾角θL,i、视线偏角ψL,i、前置倾角ηz,i、前置偏角ηy,i和速度模值vi，计算速度矢量在所述基准参考坐标系中各坐标轴的分量；基于速度矢量在所述基准参考坐标系中各坐标轴的分量计算得到当前分段的初始弹道倾角θi；
若当前分段处于转弯段，则根据公式
Figure GDA0004231246460000051
Figure GDA0004231246460000052
计算加速度在弹道坐标系中的分量；若当前分段处于直线段，则根据公式
Figure GDA0004231246460000053
计算加速度在弹道坐标系中的分量；
根据地面高度yi插值得到大气密度ρ和音速aT，根据大气密度ρ和音速aT计算得到马赫数
Figure GDA0004231246460000054
和动压
Figure GDA0004231246460000055
根据以下公式计算升力系数CLd,i：
Fy＝MT·gtccosθi+MT·ay,i
Fz＝MT·az,i
Figure GDA0004231246460000056
Figure GDA0004231246460000057
基于所述马赫数Mach和升力系数CLd,i，采用牛顿迭代法计算攻角指令αc,i；
基于所述马赫数Mach和攻角指令αc,i，根据拟合公式计算阻力系数CD,i；
根据公式
Figure GDA0004231246460000058
计算初始速度变化率，
其中，MT表示飞行器的质量，ST表示飞行器的参考面积，gtc表示引力常数。
进一步地，基于所述初始速度变化率
Figure GDA0004231246460000061
对所述当前分段初始飞行时间预测值
Figure GDA0004231246460000062
进行修正，得到当前分段飞行时间修正值
Figure GDA0004231246460000063
包括：
基于所述初始速度变化率
Figure GDA0004231246460000064
和所述初始飞行时间预测值
Figure GDA0004231246460000065
根据公式
Figure GDA0004231246460000066
得到所述分段结束点的第一速度预测值vp,i；
基于所述分段结束点的第一速度预测值vp,i和所述初始飞行时间预测值
Figure GDA0004231246460000067
根据飞行器质心运动的动力学方程计算得到第一速度变化率
Figure GDA0004231246460000068
根据公式
Figure GDA0004231246460000069
计算得到第一平均修正速度
Figure GDA00042312464600000610
若当前分段处于转弯段，则根据
Figure GDA00042312464600000611
计算当前分段飞行时间修正值
Figure GDA00042312464600000612
若当前分段处于直线段，则根据
Figure GDA00042312464600000613
计算当前分段飞行时间修正值
Figure GDA00042312464600000614
进一步地，基于所述分段结束点的第一速度预测值vp,i和所述初始飞行时间预测值
Figure GDA00042312464600000615
根据飞行器质心运动的动力学方程计算得到第一速度变化率
Figure GDA00042312464600000616
包括：
根据第一速度预测值vp,i和所述初始飞行时间预测值
Figure GDA00042312464600000617
对当前分段的视线倾角和视线偏角进行修正，分别得到第一修正视线倾角θLp,i和第一修正视线偏角ψLp,i；
若当前分段处于转弯段，则根据公式yp,i＝-Dtc,i+1sin(θLp,i)计算得到第一修正地面高度yp,i；若当前分段处于直线段，则根据公式yp,i＝-Li+1sin(θLp,i)计算得到第一修正地面高度yp,i；其中，Dtc,i+1表示i+1段起始点的目标距离，Li+1表示i+1段起始点的剩余直线段航程；
基于所述第一修正视线倾角θLp,i、第一修正视线偏角ψLp,i和第一速度预测值vp,i，计算第一速度矢量在所述基准参考坐标系中各坐标轴的分量；基于第一速度矢量在所述基准参考坐标系中各坐标轴的分量计算得到第一修正弹道倾角θp,i；
基于所述第一速度预测值vp,i、第一修正地面高度yp,i、第一修正弹道倾角θp,i，根据飞行器质心运动的动力学方程计算得到第一速度变化率
Figure GDA0004231246460000071
进一步地，根据第一速度预测值vp,i和所述初始飞行时间预测值
Figure GDA0004231246460000072
对当前分段的视线倾角和视线偏角进行修正，分别得到第一修正视线倾角θLp,i和第一修正视线偏角ψLp,i，包括：
根据以下公式对当前分段的视线倾角进行修正得到第一修正视线倾角θLp,i：
若当前分段处于转弯段，根据
Figure GDA0004231246460000073
计算初始视线倾角变化率
Figure GDA0004231246460000078
和第一修正视线倾角变化率
Figure GDA0004231246460000074
若当前分段处于直线段，则根据
Figure GDA0004231246460000075
计算初始视线倾角变化率
Figure GDA0004231246460000079
和第一修正视线倾角变化率
Figure GDA0004231246460000076
根据公式
Figure GDA0004231246460000077
得到第一修正视线倾角θLp,i；
根据以下公式对当前分段的视线偏角进行修正得到第一修正视线偏角ψLp,i：
若当前分段处于转弯段，则根据
Figure GDA0004231246460000081
Figure GDA0004231246460000082
计算初始视线偏角变化率
Figure GDA0004231246460000088
和第一修正视线偏角变化率
Figure GDA0004231246460000083
若当前分段处于直线段，则根据
Figure GDA0004231246460000084
Figure GDA0004231246460000085
计算初始视线偏角变化率
Figure GDA0004231246460000089
和第一修正视线偏角变化率
Figure GDA0004231246460000086
根据公式
Figure GDA0004231246460000087
得到第一修正视线偏角ψLp,i。
与现有技术相比，本发明通过预测剩余飞行时间，根据时间偏差生成控制项，根据该控制项生成制导控制指令，从而提高制导精确度，使得飞行器能够在约束的时间内到达目标。通过调整侧向加速度控制指令，在侧向增加时间约束控制，使得对原本的纵向运动的影响较小，可以达到良好的控制稳定性，同时形式简单，便于实施。
同时，本发明依据闭环形式的飞行器三维质心运动方程，通过根据总前置角，判断当前飞行为转弯飞行或直线飞行，当为转弯飞行时根据总前置角进行分段，当为直线飞行时根据剩余直线段航程进行分段，采用分段迭代预测剩余飞行时间，通过迭代求解，对无动力飞行器未来速度的大小进行分段预测并迭代修正剩余飞行时间估计公式，使得剩余飞行时间的估计更加精确，从而使得在总前置角大幅变化以及无动力飞行器飞行速度不受控制的情况下仍然可以准确预测剩余飞行时间，从而进一步提高达到时间控制精度。
本发明中，上述各技术方案之间还可以相互组合，以实现更多的优选组合方案。本发明的其他特征和优点将在随后的说明书中阐述，并且，部分优点可从说明书中变得显而易见，或者通过实施本发明而了解。本发明的目的和其他优点可通过说明书以及附图中所特别指出的内容中来实现和获得。
附图说明
附图仅用于示出具体实施例的目的，而并不认为是对本发明的限制，在整个附图中，相同的参考符号表示相同的部件。
图1为本发明的实施例无动力飞行器的三维制导系统的框图；
图2为本发明的实施例的坐标系图；
图3为不同剩余时间估计方法的结果比对图。
具体实施方式
下面结合附图来具体描述本发明的优选实施例，其中，附图构成本申请一部分，并与本发明的实施例一起用于阐释本发明的原理，并非用于限定本发明的范围。
如图2所示，本发明实施例的坐标系包括基准参考坐标系、视线坐标系和弹道坐标系。
基准参考坐标系以飞行器的三维质心为原点O，X轴方向平行于水平面且指向北为正，以XR表示，Y轴垂直地面向上，以YR表示，Z轴平行于水平面按右手定则确定方向，以ZR表示。
视线坐标系以飞行器的三维质心为原点O，X轴方向从飞行器指向目标点，以XL表示，ZL轴方向平行于水平面且指向XL右为正，Y轴按右手定则确定方向，以YL表示。
弹道坐标系以飞行器的三维质心为原点O，X轴方向沿速度矢量的方向，以XV表示，YV轴含速度矢量的铅垂面内与速度矢量垂直并且向上为正，Z轴按右手定则确定方向，以ZV表示。
对本申请实施例中使用的符号进行统一说明：Dtc表示飞行器与目标之间的直线距离；F表示受力；FL表示总升力；CL表示升力系数；CD表示阻力系数；η表示总前置角，即飞行器速度矢量与XL的夹角；ηy表示前置偏角(速度矢量在视线坐标系XL-O-ZL平面内的投影与XL轴之间的夹角)；ηz表示前置倾角(速度矢量与视线坐标系XL-O-ZL平面的夹角)；Δη表示总前置角变化量；Δηy表示前置偏角变化量；Δηz表示前置倾角变化量；Δη为前置角分段量常数；v表示速度模值；
Figure GDA0004231246460000101
表示速度变化率；
Figure GDA0004231246460000102
表示平均速度变化率；y表示距离地面的高度；θL表示视线倾角，即视线矢量XL与XR-O-ZR平面的夹角；ψL表示视线偏角，即视线矢量XL在XR-O-ZR平面内的投影与XR轴的夹角；
Figure GDA0004231246460000103
和
Figure GDA0004231246460000104
分别表示视线倾角变化率和视线偏角变化率；θ表示弹道倾角，即速度矢量与水平面的夹角；ay和az分别表示加速度在弹道坐标系下Yv轴和Zv轴的分量，L表示剩余直线段航程，ΔL表示直线段分段量常数；下标i表示当前分段，下标i+1表示下一个分段。
本发明的一个具体实施例，公开了一种无动力飞行器的三维制导系统，如图1所示，包括以下模块：
数据获取模块，用于获取当前时刻飞行器的目标距离、视线倾角、视线偏角、速度矢量和速度模值；
剩余飞行时间预估模块，用于基于所述视线倾角、视线偏角和速度矢量计算总前置角、前置偏角和前置倾角；基于所述目标距离、视线倾角、视线偏角、速度模值、总前置角、前置偏角和前置倾角，采用基于速度预测的分段迭代法预估剩余飞行时间；
制导指令生成模块，用于基于预估剩余飞行时间和期望剩余飞行时间的时间偏差，获取时间偏差控制项，基于当前时刻飞行器的目标距离、视线倾角、视线偏角、速度模值和所述时间偏差控制项生成制导控制指令。
本发明通过预测剩余飞行时间，根据时间偏差生成控制项，根据该控制项生成制导控制指令，从而提高到达时间精确度，使得飞行器能够在约束的时间内到达目标。通过调整侧向加速度控制指令，在侧向增加时间约束控制，使得对原本的纵向运动的影响较小，可以达到良好的控制稳定性，同时形式简单，便于实施。
其中，获取的速度矢量是在参考基准坐标系中的速度矢量，在基准参考坐标系中三个坐标轴的分量可分别表示为vxR、vyR和vzR。
具体的，剩余飞行时间预估模块，用于基于视线倾角、视线偏角和速度矢量计算总前置角、前置偏角和前置倾角，包括：
S11、根据所述视线倾角和视线偏角计算基准参考坐标系到视线坐标系的方向余弦矩阵；
具体的，根据公式
Figure GDA0004231246460000121
计算基准参考坐标系到视线坐标系的方向余弦矩阵
Figure GDA0004231246460000122
S12、基于所述方向余弦矩阵计算所述速度矢量在所述视线坐标系中的分量；
速度矢量在所述视线坐标系中的分量为
Figure GDA0004231246460000123
其中，vxL、vyL和vzL分别表示速度矢量在视线坐标系中的坐标分量。
S13、根据以下公式计算总前置角：
Figure GDA0004231246460000124
其中，vxL、vyL和vzL分别表示速度矢量在视线坐标系三个坐标轴的分量，ηz表示前置倾角，ηy表示前置偏角，η表示总前置角。
剩余飞行时间预估模块，根据视线倾角、视线偏角和速度矢量计算当前时刻飞行器的总前置角、前置偏角和前置倾角。
剩余飞行时间预估模块基于所述目标距离、视线倾角、视线偏角、速度模值、总前置角、前置偏角和前置倾角，采用基于速度预测的分段迭代法预估剩余飞行时间，包括：
根据总前置角判断当前飞行为转弯飞行或直线飞行。
若为转弯飞行，则根据总前置角对转弯段进行分段，采用基于速度预测的分段迭代预测转弯段剩余飞行时间；根据直线段起始点的总前置角和目标距离计算剩余直线段航程，根据剩余直线段航程对直线段进行分段，采用基于速度预测的分段迭代预测直线段剩余飞行时间；所述转弯段剩余飞行时间与直线段剩余飞行时间的和为预测的总剩余飞行时间；
若为直线飞行，根据所述总前置角和目标距离计算剩余直线段航程，根据剩余直线段航程对直线段进行分段，采用基于速度预测的分段迭代预测直线段剩余飞行时间，所述直线段剩余飞行时间为预测的总剩余飞行时间。
通过根据总前置角，判断当前飞行为转弯飞行或直线飞行，当为转弯飞行时根据总前置角进行分段，当为直线飞行时根据剩余直线段航程进行分段，保证每个分段区间内飞行器的前置倾角的增量均为小角度，采用基于速度预测的分段迭代预测剩余飞行时间，从而使得在总前置角大幅变化以及无动力飞行器飞行速度不受控制的情况下仍然可以准确预测剩余飞行时间。
实施时，前置角分段量常数Δη可取较小的量，例如2°。当总前置角η0＞Δη时，则判断当前飞行为转弯飞行，即当前进行剩余飞行预测时刻飞行器为转弯飞行；否则，判断当前飞行为直线飞行。
若当前飞行为转弯飞行，则根据总前置角对转弯段进行分段，根据公式
Figure GDA0004231246460000131
计算转弯段的分段数为NumP，floor(·)表示向下取整。依次对每一转弯分段预测飞行时间，得到转弯段的剩余飞行时间，然后，根据直线段起始点的总前置角和目标距离计算剩余直线段航程，根据剩余直线段航程对直线段进行分段，依次预测每一直线分段的飞行时间得到直线段的剩余飞行时间，总的剩余飞行时间是转弯段剩余飞行时间与直线段剩余飞行时间的和。其中，直线段的起始点即转弯段的结束点。
若当前飞行为直线飞行，则直接根据剩余直线段航程对直线段进行分段。根据公式
Figure GDA0004231246460000132
计算剩余直线段航程，根据公式
Figure GDA0004231246460000141
计算直线段的分段数NumL+1，依次对每一分段预测飞行时间，得到直线段的剩余飞行时间，即预测的总的剩余飞行时间。
其中，Dtc表示进行剩余飞行预测时刻飞行器与目标之间的直线距离，η表示进行剩余飞行预测时刻的总前置角，KN为采用的比例导引法的导航比，L表示进行剩余飞行预测时刻的剩余直线段航程。
具体的，剩余飞行时间预估模块采用基于速度预测的分段迭代预测转弯段剩余飞行时间或采用基于速度预测的分段迭代预测直线段剩余飞行时间，包括：
S21、对于每一分段，获取当前分段起始点的状态变量，所述状态变量包括速度模值vi、视线倾角θL,i、视线偏角ψL,i、前置倾角ηz,i、前置偏角ηy,i、地面高度yi；
若当前分段处于转弯段，所述状态变量还包括总前置角ηi和目标距离Dtc,i；若当前分段处于直线段，所述状态变量还包括剩余直线段航程Li。
若当前飞行为转弯飞行，对于转弯段的第一分段，其起始点的状态变量即进行剩余飞行预测时刻的状态值，其直线段的第一分段起始点状态变量为转弯段最后一段结束点的状态变量，其中直线段的第一分段起始点的剩余直线段航程
Figure GDA0004231246460000142
Dtc,Nump+1和ηNump+1是转弯段最后一段结束点的目标的距离和总前置角。
若当前飞行为直线飞行，直线段的第一分段起始点的状态变量即进行剩余飞行预测时刻的状态值。
S22、若当前分段处于转弯段，则根据
Figure GDA0004231246460000143
计算当前分段初始飞行时间预测值
Figure GDA0004231246460000144
若当前分段处于直线段，则根据
Figure GDA0004231246460000145
计算当前分段初始飞行时间预测值
Figure GDA0004231246460000146
其中，Δηi表示当前分段总前置角的变化量，ΔL表示当前分段剩余直线段航程变化量，KN表示导航比；
具体的，若当前分段处于转弯段，当前分段总前置角的变化量Δηi根据以下公式计算得到
Figure GDA0004231246460000151
即当i＝NumP时，即第NumP分段取Δηi＝-(ηi-NumP·Δη)，第1～NumP-1段，则取Δηi＝-Δη。
Figure GDA0004231246460000152
根据以下过程得出：
对ηi＝arccos(cosηy,icosηz,i)进行微分，可以得到
Figure GDA0004231246460000153
又因为
Figure GDA0004231246460000154
Figure GDA0004231246460000155
Figure GDA0004231246460000156
有
Figure GDA0004231246460000157
因此，飞行器在三维空间的运动具有如下运动形式，即
Figure GDA0004231246460000158
消去时间变量有
Figure GDA0004231246460000159
然后在区间[ti,t]上积分并化简得
Figure GDA0004231246460000161
因此有
Figure GDA0004231246460000162
然后对sinηt一阶泰勒展开即sinηt＝sinηi+Δηt·cosηt，代入前式有
Figure GDA0004231246460000163
然后在
Figure GDA0004231246460000164
上积分可得
Figure GDA0004231246460000165
Figure GDA0004231246460000166
根据公式
Figure GDA0004231246460000167
对当前分段的飞行时间进行预测，得到当前分段的初始飞行时间预测值
Figure GDA0004231246460000168
若当前分段处于直线段，则根据
Figure GDA0004231246460000169
预测前分段的飞行时间得到当前分段的初始飞行时间预测值
Figure GDA00042312464600001610
S23、基于当前分段起始点的状态变量，根据飞行器质心运动的动力学方程计算当前分段的初始速度变化率
Figure GDA00042312464600001611
S24、基于所述初始速度变化率
Figure GDA00042312464600001612
对所述当前分段初始飞行时间预测值
Figure GDA00042312464600001613
进行修正，得到当前分段飞行时间修正值
Figure GDA00042312464600001614
S25、基于当前分段起始点的状态变量和所述当前分段飞行时间修正值
Figure GDA00042312464600001615
预测当前分段结束点的状态变量，当前分段结束点的状态变量为下一分段起始点的状态变量；
需要说明的是，若当前分段为直线段的最后一段，得到当前分段飞行时间修正值
Figure GDA00042312464600001616
即可得到总的剩余飞行时间，不再预测当前分段结束点的状态变量。
所有分段的飞行时间修正值
Figure GDA00042312464600001617
的和作为预测的转弯段或直线段剩余飞行时间。
在步骤S21之后，步骤S24之前，还包括：
S201、若当前分段处于转弯段，则根据以下公式计算当前分段结束点的前置倾角ηz,i+1、前置偏角ηy,i+1和总前置角ηi+1：
Figure GDA0004231246460000171
Figure GDA0004231246460000172
ηy,i+1＝ηy,i+Δηy,i，ηz,i+1＝ηz,i+Δηz,i，ηi+1＝arccos(cosηy,i+1cosηz,i+1)；
根据公式
Figure GDA0004231246460000173
计算当前分段结束点的目标距离Dtc,i+1；
若当前分段处于直线段，则根据以下公式计算当前分段结束点的剩余直线段航程Li+1、前置倾角ηz,i+1和前置偏角ηy,i+1：
Li+1＝Li-ΔL，
Figure GDA0004231246460000174
其中，KN表示导航比，Δη表示前置角分段量常数。
步骤S23中，基于当前分段起始点的状态变量，根据飞行器质心运动的动力学方程计算当前分段的初始速度变化率，包括步骤S231-S237：
S231、基于所述视线倾角θL,i、视线偏角ψL,i、前置倾角ηz,i、前置偏角ηy,i和速度模值vi，计算速度矢量在所述基准参考坐标系中各坐标轴的分量；基于速度矢量在所述基准参考坐标系中各坐标轴的分量计算得到当前分段的弹道倾角θi；
具体的，根据以下公式计算弹道倾角θ：
Figure GDA0004231246460000175
Figure GDA0004231246460000181
其中，
Figure GDA0004231246460000182
表示基准参考坐标系到视线坐标系的方向余弦矩阵，v表示速度模值，ηz表示前置倾角，ηz表示前置偏角，θL表示视线倾角，ψL表示视线偏角，vpx,vpy,vpz分别表示速度矢量在视线坐标系各轴的分量和vqx,vqy,vqz分别表示速度矢量在参考坐标系各轴的分量。
将上述计算弹道倾角θ的公式表示为θ＝fθ(θL,ψL,ηy,ηz,v)。
将根据当前分段起始点的状态变量视线倾角θL,i、视线偏角ψL,i、速度模值vi、前置倾角ηz,i和前置偏角ηy,i带入θ＝fθ(θL,ψL,ηy,ηz,v)可计算出当前分段的初始弹道倾角θi。
S232、若当前分段处于转弯段，则根据公式
Figure GDA0004231246460000183
Figure GDA0004231246460000184
计算当前加速度在弹道坐标系中的分量；若当前分段处于直线段，则根据公式
Figure GDA0004231246460000185
计算当前加速度在弹道坐标系中的分量；
ay表示加速度在弹道坐标系的Yv轴的分量，az表示加速度在弹道坐标系的Zv轴的分量。
S233、根据地面高度yi插值得到大气密度ρ和音速aT，根据大气密度ρ和音速aT计算得到马赫数
Figure GDA0004231246460000186
和动压
Figure GDA0004231246460000187
具体的，根据美国标准大气模型，可根据当前分段距离地面的高度插值得到大气密度ρ和音速aT，根据
Figure GDA0004231246460000191
计算马赫数Mach，根据
Figure GDA0004231246460000192
计算动压。
S234、根据以下公式计算升力系数CLd,i：
Fy＝MT·gtccosθi+MT·ay,i，Fz＝MT·az,i
Figure GDA0004231246460000193
其中，MT表示飞行器质量，ST表示飞行器参考面积，gtc表示引力常数，sign表示符号函数，FL表示总升力，Fy表示升力在弹道坐标系的Yv轴的分量，Fz表示升力在弹道坐标系的Zv轴的分量；根据升力的计算原理：升力＝升力系数*动压*飞行器参考面积，计算升力系数CLd,i。
S235、基于所述马赫数Mach和升力系数CLd,i，采用牛顿迭代法计算攻角指令αc,i。
具体的，牛顿迭代法的迭代公式为：
Figure GDA0004231246460000194
迭代初值为5°即α0为5°，迭代精度满足精度要求，例如1e-5，或达到最大迭代次数，例如100次，则迭代结束。
具体的，牛顿迭代法用到的函数f(α)和导数f'(α)如下：
f(α)＝kCl_00+kCl_10α+kCl_01Mach+kCl_20α2+kCl_11α·Mach+kCl_02Mach2+kCl_30α3+kCl_21α2·Mach+kCl_12α·Mach2+kCl_03Mach3+kCl_31α3·Mach+kCl_22α2·Mach2+kCl_13α·Mach3+kCl_ 04Mach4+kCl_32α3·Mach2+kCl_23α2·Mach3+kCl_14α·Mach4+kCl_05Mach5-CLd
f'(α)＝kCf_00+kCf_10α+kCf_01Mach+kCf_20α2+kCf_11α·Mach+kCf_02Mach2+kCf_21α2·Mach+kCf_12α·Mach2+kCf_03Mach3+kCf_22α2·Mach2+kCf_13α·Mach3+kCf_04Mach4
其中，kCl_00、kCl_10、kCl_01、kCl_20、kCl_11、kCl_02、kCl_30、kCl_21、kCl_12、kCl_03、kCl_31、kCl_22、kCl_13、kCl_04、kCl_32、kCl_23、kCl_14、kCl_05，kCf_00、kCf_10、kCf_01、kCf_20、kCf_11、kCf_02、kCf_21、kCf_12、kCf_03、kCf_22、kCf_13、kCf_04为事先根据飞行器的气动参数拟合的系数。
将步骤S232计算的马赫数Mach和步骤S234计算的升力系数CLd,i带入上式迭代计算得到攻角指令αc,i。牛顿迭代法的收敛速度快，通过设置收敛条件可以得到较为精确的数值解，从而快速准确的计算公攻角指令。
S236、基于所述马赫数Mach和攻角指令αc,i，根据拟合公式计算阻力系数CD,i；
具体的，按下面的公式计算阻力系数：
Figure GDA0004231246460000201
其中，ktc_00、ktc_10、ktc_01、ktc_20、ktc_11、ktc_02、ktc_30、ktc_21、ktc_12、ktc_03、ktc_40、ktc_31、ktc_22、ktc_13、ktc_04、ktc_50、ktc_41、ktc_32、ktc_23、ktc_14、ktc_05为事先根据飞行器的气动参数拟合的系数。飞行器的气动特性一般是用大量的离散数据描述的，通过对数据的拟合，可以得到阻力系数和攻角与马赫数之间的近似的拟合关系式，方法简单，计算速度快。
将步骤S232计算的马赫数Mach和步骤S236计算的攻角指令αc,i带入上式迭代计算得到阻力系数CD,i。
S237、根据公式
Figure GDA0004231246460000202
计算初始速度变化率。
具体的，根据飞行器质心运动的动力学方程
Figure GDA0004231246460000211
计算得到初始速度变化率。无动力飞行器受力主要来自阻力和重力，因此其动力学方程表示为
Figure GDA0004231246460000212
其中
Figure GDA0004231246460000213
表示阻力对速度变化率的影响，gtc sinθi表示重力对速度变化率的影响。
其中，MT表示飞行器的质量，ST表示飞行器的参考面积，gtc表示引力常数。
将步骤233-步骤237计算速度变化率的方法表示为
Figure GDA0004231246460000214
v表示基准参考坐标系中的速度模量，ay表示加速度在弹道坐标系的Yv轴的分量，az表示加速度在弹道坐标系的Zv轴的分量，θ表示弹道倾角，y表示地面高度。
具体的，步骤S24基于所述初始速度变化率
Figure GDA0004231246460000215
对所述当前分段初始飞行时间预测值
Figure GDA0004231246460000216
进行修正，得到当前分段飞行时间修正值
Figure GDA0004231246460000217
包括步骤S241-S244。
S241、基于所述初始速度变化率
Figure GDA00042312464600002114
和所述初始飞行时间预测值
Figure GDA0004231246460000218
根据公式
Figure GDA0004231246460000219
得到所述分段结束点的第一速度预测值vp,i；
S242、基于所述分段结束点的第一速度预测值vp,i和所述初始飞行时间预测值
Figure GDA00042312464600002110
根据飞行器质心运动的动力学方程计算得到第一速度变化率
Figure GDA00042312464600002111
具体的，步骤S242中，基于所述分段结束点的第一速度预测值vp,i和所述初始飞行时间预测值
Figure GDA00042312464600002112
根据飞行器质心运动的动力学方程计算得到第一速度变化率
Figure GDA00042312464600002113
包括步骤S2421-S2424：
S2421、根据第一速度预测值vp,i和所述初始飞行时间预测值
Figure GDA0004231246460000221
对当前分段的视线倾角和视线偏角进行修正，分别得到第一修正视线倾角θLp,i和第一修正视线偏角ψLp,i；
具体的，根据第一速度预测值vp,i和所述初始飞行时间预测值
Figure GDA0004231246460000222
对当前分段的视线倾角和视线偏角进行修正，分别得到第一修正视线倾角θLp,i和第一修正视线偏角ψLp,i，包括：
根据以下公式对当前分段的视线倾角进行修正得到第一修正视线倾角θLp,i：
若当前分段处于转弯段，则根据
Figure GDA0004231246460000223
计算初始视线倾角变化率
Figure GDA0004231246460000224
和第一修正视线倾角变化率
Figure GDA0004231246460000225
若当前分段处于直线段，则根据
Figure GDA0004231246460000226
计算初始视线倾角变化率
Figure GDA0004231246460000227
和第一修正视线倾角变化率
Figure GDA0004231246460000228
根据公式
Figure GDA0004231246460000229
计算第一修正视线倾角θLp,i。
需要说明的是，若当前分段为直线段的最后一段，此时Li+1为0，因此第一修正视线倾角变化率
Figure GDA00042312464600002210
直接取0。
对当前分段起始点的视线倾角变化率
Figure GDA00042312464600002211
和第一修正视线倾角变化率
Figure GDA00042312464600002214
取平均，得到平均视线倾角变化率
Figure GDA00042312464600002212
从而根据初始飞行时间预测值
Figure GDA00042312464600002213
对当前分段的视线倾角进行修正，得到第一修正视线倾角θLp,i。
根据以下公式对当前分段的视线偏角进行修正得到第一修正视线偏角ψLp,i：
若当前分段处于转弯段，则根据
Figure GDA0004231246460000231
Figure GDA0004231246460000232
计算初始视线偏角变化率
Figure GDA0004231246460000233
和第一修正视线偏角变化率
Figure GDA0004231246460000234
若当前分段处于直线段，则根据
Figure GDA0004231246460000235
Figure GDA0004231246460000236
计算初始视线偏角变化率
Figure GDA0004231246460000239
和第一修正视线偏角变化率
Figure GDA0004231246460000237
根据公式
Figure GDA0004231246460000238
得到第一修正视线偏角ψLp,i。
需要说明的是，若当前分段为直线段的最后一段，此时Li+1为0，因此第一修正视线偏角变化率
Figure GDA00042312464600002310
直接取0。
S2422、若当前分段处于转弯段，则根据公式yp,i＝-Dtc,i+1sin(θLp,i)计算得到第一修正地面高度yp,i；若当前分段处于直线段，则根据公式yp,i＝-Li+1sin(θLp,i)计算得到第一修正地面高度yp,i。
S2423、基于所述第一修正视线倾角θLp,i、第一修正视线偏角ψLp,i和第一速度预测值vp,i，计算第一速度矢量在所述基准参考坐标系中各坐标轴的分量；基于第一速度矢量在所述基准参考坐标系中各坐标轴的分量计算得到第一修正弹道倾角θp,i；
具体的，根据步骤S231中计算弹道倾角θ的公式表示为θ＝fθ(θL,ψL,ηy,ηz,v)，将第一修正视线倾角θLp,i、第一修正视线偏角ψLp,i、第一速度预测值vp,i、当前分段结束点的前置倾角ηz,i+1、和当前分段结束点的前置偏角ηy,i+1带入fθ，计算得到第一修正弹道倾角θp,i。
其中当前分段结束点的前置倾角ηz,i+1和当前分段结束点的前置偏角ηy,i+1可根据步骤S201计算得到。
S2424、基于所述第一速度预测值vp,i、第一修正地面高度yp,i、第一修正弹道倾角θp,i，根据飞行器质心运动的动力学方程计算得到第一速度变化率
Figure GDA0004231246460000241
具体的，若当前分段处于转弯段，则根据公式
Figure GDA0004231246460000242
计算加速度在弹道坐标系中的分量；若当前分段处于直线段，则根据公式
Figure GDA0004231246460000243
计算加速度在弹道坐标系中的分量。
需要说明的是，若当前分段为直线段最后一段，则由于Li+1为0，因此不采用上述公式计算ayp,i和azp,i，而直接将ayp,i和azp,i置为0。
其中当前分段结束点的前置倾角ηz,i+1、当前分段结束点的前置偏角ηy,i+1、当前分段结束点的目标距离Dtc,i+1、当前分段结束点的剩余直线段航程Li，可根据步骤可根据步骤S201计算得到。
根据步骤S233-S237中计算速度变化率的公式
Figure GDA0004231246460000244
将第一修正弹道倾角θp,i、第一速度预测值vp,i、ayp,i、azp,i和第一修正地面高度yp,i带入
Figure GDA0004231246460000251
计算得到第一速度变化率
Figure GDA0004231246460000252
S243、根据公式
Figure GDA0004231246460000253
计算得到第一平均修正速度
Figure GDA0004231246460000254
即，对初始速度变化率
Figure GDA0004231246460000255
和第一速度变化率
Figure GDA0004231246460000256
取平均，得到平均速度变化率
Figure GDA0004231246460000257
根据初始速度vi、平均速度变化率
Figure GDA0004231246460000258
和初始飞行时间预测值
Figure GDA0004231246460000259
对速度进行修正，得到第一平均修正速度
Figure GDA00042312464600002510
通过取平均值提高估计的精度。
S244、若当前分段处于转弯段，则根据
Figure GDA00042312464600002511
计算当前分段飞行时间修正值
Figure GDA00042312464600002512
若当前分段处于直线段，则根据
Figure GDA00042312464600002513
计算当前分段飞行时间修正值
Figure GDA00042312464600002514
根据得到第一平均修正速度
Figure GDA00042312464600002520
重新计算当前分段的飞行时间，从而得到当前飞行时间修正值
Figure GDA00042312464600002515
通过根据初始速度变化率对视线倾角、视线偏角、弹道倾角进行修正，进一步修正飞行时间预测值，从而是预测更加准确。
步骤S25，基于当前分段起始点的状态变量和所述当前分段飞行时间修正值
Figure GDA00042312464600002516
预测当前分段结束点的状态变量，包括：
S251、根据当前分段的初始速度变化率
Figure GDA00042312464600002519
和所述当前分段飞行时间修正值
Figure GDA00042312464600002517
根据公式
Figure GDA00042312464600002518
得到当前分段结束点的第二速度预测值vq,i；
以修正后的飞行时间估计值再次计算速度预测值，提高估计的精度。
S252、基于所述分段结束点的第二速度预测值vq,i和当前分段飞行时间修正值
Figure GDA0004231246460000261
根据飞行器质心运动的动力学方程计算得到第二速度变化率
Figure GDA0004231246460000262
具体的，根据
Figure GDA0004231246460000263
计算第二修正视线倾角θLq,i，根据
Figure GDA0004231246460000264
计算得到第二视线偏角ψLq,i；根据yq,i＝-Dtc,i+1sinθLq,i计算第二修正地面高度yq,i；
根据步骤S231中计算弹道倾角θ的公式表示为θ＝fθ(θL,ψL,ηy,ηz,v)，将第二速度预测值vq,i、第二修正视线倾角θLq,i、第二修正视线偏角ψLq,i、当前分段结束点的前置倾角ηz,i+1、和当前分段结束点的前置偏角ηy,i+1带入fθ，计算得到第二修正弹道倾角θq,i。
根据步骤S2424中计算ayp,i和azp,i的公式，将其中vp,i替换为vq,i，计算得到ayq,i和azq,i
根据步骤S233-S237中计算速度变化率的公式
Figure GDA0004231246460000265
将第二修正弹道倾角θq,i、第二速度预测值vq,i、ayq,i、azq,i和第二修正地面高度yq,i带入
Figure GDA0004231246460000266
计算得到第二速度变化率
Figure GDA0004231246460000267
S253、对初始速度变化率
Figure GDA0004231246460000268
和第二速度变化率
Figure GDA0004231246460000269
取平均，得到第二平均修正速度
Figure GDA00042312464600002610
Figure GDA00042312464600002611
S254、基于所述第二平均修正速度
Figure GDA00042312464600002612
计算当前分段的结束点的第三速度预测值
Figure GDA00042312464600002613
所述第三速度预测值为预测得到的当前分段的结束点的速度vi+1；
S255、根据当前分段起始点的视线倾角θL,i、当前分段的结束点的速度vi+1和所述当前分段飞行时间修正值
Figure GDA0004231246460000271
计算得到当前分段的结束点的视线倾角θL,i+1；根据当前分段起始点的视线偏角ψL,i、当前分段的结束点的速度vi+1和所述当前分段飞行时间修正值
Figure GDA0004231246460000272
计算得到当前分段的结束点的视线偏角ψL,i+1。
具体的，根据以下公式计算得到当前分段的结束点的视线倾角θL,i+1；
若当前分段处于转弯段，根据公式
Figure GDA0004231246460000273
计算第二修正视线倾角变化率
Figure GDA0004231246460000274
若当前分段处于直线段，根据公式
Figure GDA0004231246460000275
计算第二修正视线倾角变化率
Figure GDA0004231246460000276
根据公式
Figure GDA0004231246460000277
计算当前分段的结束点的视线倾角θL,i+1。
若当前分段处于转弯段，根据公式
Figure GDA0004231246460000278
计算第二修正视线篇角变化率
Figure GDA0004231246460000279
若当前分段处于直线段，根据公式
Figure GDA00042312464600002710
计算第二修正视线篇角变化率
Figure GDA00042312464600002712
根据公式
Figure GDA00042312464600002711
计算当前分段的结束点的视线偏角ψL,i+1。
S256、若当前分段处于转弯段，则根据公式yi+1＝-Dtc,i+1sinθL,i+1计算当前分段的结束点的地面高度；若当前分段处于直线段，则根据yi+1＝-Li+1sinθL,i+1计算当前分段的结束点的地面高度yi+1。
通过采用分段的方式，对每个分段预测飞行时间，所有分段的时间总和即为总的剩余飞行时间，从而使得无动力飞行器速度显著下降以及总前置角大幅变化情况下剩余飞行时间的估算更准确。
使用本申请的预测方法和其他方法的剩余飞行时间估计结果如附图3所示。可以看出，实际的剩余飞行时间是斜率为-1的直线，传统估计方法和常速小角度估计方法由于不考虑速度变化和大前置角，很明显地估计偏差较大。而采用本申请提高的剩余飞行时间预测方法，其剩余飞行时间估计曲线大致和实际的剩余飞行时间曲线重合，随着时间推移接近程度越来越高，显示出较高的估计精度。
得到预估剩余飞行时间后，制导指令生成模块根据公式
Figure GDA0004231246460000281
计算预估剩余飞行时间和期望剩余飞行时间的时间偏差。基于预估剩余飞行时间和期望剩余飞行时间的时间偏差，获取时间偏差控制项，包括：
根据公式aε＝-KTvsinηycosηyεT计算时间偏差控制项；其中，εT表示时间偏差，ηy表示前置偏角，v表示速度模值，KT表示时间控制系数，aε表示时间偏差控制项。
将时间偏差控制项引入制导律，得到制导控制指令。
具体的，基于当前时刻飞行器的目标距离、视线倾角、视线偏角、速度模值和所述时间偏差控制项生成制导控制指令，包括:
根据公式
Figure GDA0004231246460000282
计算纵向加速度指令；
根据公式
Figure GDA0004231246460000283
计算侧向加速度指令；
其中，Dtc表示目标距离，ηy表示前置偏角，ηz表示前置倾角，KN表示导航比，v表示速度模值，aε表示时间偏差控制项。
得到纵向加速度指令和侧向加速度指令后，将控制指令发送至飞行器姿态控制系统，通过姿态控制系统调整飞行器的飞行姿态，从而实现飞行器按照预期的飞行时间飞向目标。
本发明中，上述各技术方案之间还可以相互组合，以实现更多的优选组合方案。本发明的其他特征和优点将在随后的说明书中阐述，并且，部分优点可从说明书中变得显而易见，或者通过实施本发明而了解。本发明的目的和其他优点可通过说明书以及附图中所特别指出的内容中来实现和获得。
本领域技术人员可以理解，实现上述实施例方法的全部或部分流程，可以通过计算机程序来指令相关的硬件来完成，所述的程序可存储于计算机可读存储介质中。其中，所述计算机可读存储介质为磁盘、光盘、只读存储记忆体或随机存储记忆体等。
以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。

Claims (1)

1.一种无动力飞行器的三维制导系统，其特征在于，包括以下模块：
数据获取模块，用于获取当前时刻飞行器的目标距离、视线倾角、视线偏角、速度矢量和速度模值；
剩余飞行时间预估模块，用于基于所述视线倾角、视线偏角和速度矢量计算总前置角、前置偏角和前置倾角；基于所述目标距离、视线倾角、视线偏角、速度模值、总前置角、前置偏角和前置倾角，采用基于速度预测的分段迭代法预估剩余飞行时间；
制导指令生成模块，用于基于预估剩余飞行时间和期望剩余飞行时间的时间偏差，获取时间偏差控制项，基于当前时刻飞行器的目标距离、视线倾角、视线偏角、速度模值和所述时间偏差控制项生成制导控制指令；
所述剩余飞行时间预估模块基于所述视线倾角、视线偏角和速度矢量计算总前置角、前置偏角和前置倾角，包括：
根据所述视线倾角和视线偏角计算基准参考坐标系到视线坐标系的方向余弦矩阵；
基于所述方向余弦矩阵计算所述速度矢量在所述视线坐标系中的分量；
根据以下公式计算总前置角：
Figure FDA0004231246440000011
Figure FDA0004231246440000012
η＝arccos(cosηycosηz)
其中，vxL、vyL和vzL分别表示速度矢量在视线坐标系三个坐标轴的分量，ηz表示前置倾角，ηy表示前置偏角，η表示总前置角；
所述制导指令生成模块，用于基于当前时刻飞行器的目标距离、视线倾角、视线偏角、速度模值和所述时间偏差控制项生成制导控制指令，包括:
根据公式
Figure FDA0004231246440000021
计算纵向加速度指令；
根据公式
Figure FDA0004231246440000022
计算侧向加速度指令；
其中，Dtc表示目标距离，ηy表示前置偏角，ηz表示前置倾角，KN表示导航比，v表示速度模值，aε表示时间偏差控制项；
所述制导指令生成模块，用于基于预估剩余飞行时间和期望剩余飞行时间的时间偏差，获取时间偏差控制项，包括：
根据公式aε＝-KTvsinηycosηyεT计算时间偏差控制项；其中，εT表示时间偏差，ηy表示前置偏角，v表示速度模值，KT表示时间控制系数，aε表示时间偏差控制项；
所述剩余飞行时间预估模块基于所述目标距离、视线倾角、视线偏角、速度模值、总前置角、前置偏角和前置倾角，采用基于速度预测的分段迭代法预估剩余飞行时间，包括：
根据所述总前置角判断当前飞行为转弯飞行或直线飞行；
若为转弯飞行，则根据总前置角对转弯段进行分段，采用分段迭代预测转弯段剩余飞行时间；根据直线段起始点的总前置角和目标距离计算剩余直线段航程，根据剩余直线段航程对直线段进行分段，采用分段迭代预测直线段剩余飞行时间；所述转弯段剩余飞行时间与直线段剩余飞行时间的和为预测的总剩余飞行时间；
直线飞行预测模块，用于若为直线飞行，根据所述总前置角和目标距离计算剩余直线段航程，根据剩余直线段航程对直线段进行分段，采用分段迭代预测直线段剩余飞行时间，所述直线段剩余飞行时间为预测的总剩余飞行时间；
采用分段迭代预测转弯段剩余飞行时间或采用分段迭代预测直线段剩余飞行时间，包括：
对于每一分段，获取当前分段起始点的状态变量，所述状态变量包括速度模值vi、视线倾角θL,i、视线偏角ψL,i、前置倾角ηz,i、前置偏角ηy,i、地面高度yi；
若当前分段处于转弯段，所述状态变量还包括总前置角ηi和目标距离Dtc,i；若当前分段处于直线段，所述状态变量还包括剩余直线段航程Li；
若当前分段处于转弯段，则根据
Figure FDA0004231246440000031
计算当前分段初始飞行时间预测值
Figure FDA0004231246440000032
若当前分段处于直线段，则根据
Figure FDA0004231246440000033
计算当前分段初始飞行时间预测值
Figure FDA0004231246440000034
其中，Δηi表示当前分段总前置角的变化量，ΔL表示当前分段剩余直线段航程变化量，KN表示导航比；
基于当前分段起始点的状态变量，根据飞行器质心运动的动力学方程计算当前分段的初始速度变化率
Figure FDA0004231246440000035
基于所述初始速度变化率
Figure FDA0004231246440000036
对所述当前分段初始飞行时间预测值
Figure FDA0004231246440000037
进行修正，得到当前分段飞行时间修正值
Figure FDA0004231246440000038
基于当前分段起始点的状态变量和所述当前分段飞行时间修正值
Figure FDA0004231246440000039
预测当前分段结束点的状态变量，当前分段结束点的状态变量为下一分段起始点的状态变量；
所有分段的飞行时间修正值
Figure FDA00042312464400000310
的和作为预测的转弯段或直线段剩余飞行时间；
基于当前分段起始点的状态变量，根据飞行器质心运动的动力学方程计算当前分段的初始速度变化率
Figure FDA0004231246440000041
包括：
基于所述视线倾角θL,i、视线偏角ψL,i、前置倾角ηz,i、前置偏角ηy,i和速度模值vi，计算速度矢量在基准参考坐标系中各坐标轴的分量；基于速度矢量在所述基准参考坐标系中各坐标轴的分量计算得到当前分段的初始弹道倾角θi；
若当前分段处于转弯段，则根据公式
Figure FDA0004231246440000042
Figure FDA0004231246440000043
计算加速度在弹道坐标系中的分量；若当前分段处于直线段，则根据公式
Figure FDA0004231246440000044
计算加速度在弹道坐标系中的分量；
根据地面高度yi插值得到大气密度ρ和音速aT，根据大气密度ρ和音速aT计算得到马赫数
Figure FDA0004231246440000045
和动压
Figure FDA0004231246440000046
根据以下公式计算升力系数CLd,i：
Fy＝MT·gtccosθi+MT·ay,i
Fz＝MT·az,i
Figure FDA0004231246440000047
Figure FDA0004231246440000048
基于所述马赫数Mach和升力系数CLd,i，采用牛顿迭代法计算攻角指令αc,i；
基于所述马赫数Mach和攻角指令αc,i，根据拟合公式计算阻力系数CD,i；
根据公式
Figure FDA0004231246440000051
计算初始速度变化率，
其中，MT表示飞行器的质量，ST表示飞行器的参考面积，gtc表示引力常数；
基于所述初始速度变化率
Figure FDA0004231246440000052
对所述当前分段初始飞行时间预测值
Figure FDA0004231246440000053
进行修正，得到当前分段飞行时间修正值
Figure FDA0004231246440000054
包括：
基于所述初始速度变化率
Figure FDA0004231246440000055
和所述初始飞行时间预测值
Figure FDA0004231246440000056
根据公式
Figure FDA0004231246440000057
得到所述分段结束点的第一速度预测值vp,i；
基于所述分段结束点的第一速度预测值vp,i和所述初始飞行时间预测值
Figure FDA0004231246440000058
根据飞行器质心运动的动力学方程计算得到第一速度变化率
Figure FDA0004231246440000059
根据公式
Figure FDA00042312464400000510
计算得到第一平均修正速度
Figure FDA00042312464400000511
若当前分段处于转弯段，则根据
Figure FDA00042312464400000512
计算当前分段飞行时间修正值
Figure FDA00042312464400000513
若当前分段处于直线段，则根据
Figure FDA00042312464400000514
计算当前分段飞行时间修正值
Figure FDA00042312464400000515
基于所述分段结束点的第一速度预测值vp,i和所述初始飞行时间预测值
Figure FDA00042312464400000516
根据飞行器质心运动的动力学方程计算得到第一速度变化率
Figure FDA00042312464400000517
包括：
根据第一速度预测值vp,i和所述初始飞行时间预测值
Figure FDA00042312464400000518
对当前分段的视线倾角和视线偏角进行修正，分别得到第一修正视线倾角θLp,i和第一修正视线偏角ψLp,i；
若当前分段处于转弯段，则根据公式yp,i＝-Dtc,i+1sin(θLp,i)计算得到第一修正地面高度yp,i；若当前分段处于直线段，则根据公式yp,i＝-Li+1sin(θLp,i)计算得到第一修正地面高度yp,i；其中，Dtc,i+1表示i+1段起始点的目标距离，Li+1表示i+1段起始点的剩余直线段航程；
基于所述第一修正视线倾角θLp,i、第一修正视线偏角ψLp,i和第一速度预测值vp,i，计算第一速度矢量在基准参考坐标系中各坐标轴的分量；基于第一速度矢量在所述基准参考坐标系中各坐标轴的分量计算得到第一修正弹道倾角θp,i；
基于所述第一速度预测值vp,i、第一修正地面高度yp,i、第一修正弹道倾角θp,i，根据飞行器质心运动的动力学方程计算得到第一速度变化率
Figure FDA0004231246440000061
根据第一速度预测值vp,i和所述初始飞行时间预测值
Figure FDA0004231246440000062
对当前分段的视线倾角和视线偏角进行修正，分别得到第一修正视线倾角θLp,i和第一修正视线偏角ψLp,i，包括：
根据以下公式对当前分段的视线倾角进行修正得到第一修正视线倾角θLp,i：
若当前分段处于转弯段，根据
Figure FDA0004231246440000063
计算初始视线倾角变化率
Figure FDA0004231246440000064
和第一修正视线倾角变化率
Figure FDA0004231246440000065
若当前分段处于直线段，则根据
Figure FDA0004231246440000066
计算初始视线倾角变化率
Figure FDA0004231246440000067
和第一修正视线倾角变化率
Figure FDA0004231246440000068
根据公式
Figure FDA0004231246440000071
得到第一修正视线倾角θLp,i；
根据以下公式对当前分段的视线偏角进行修正得到第一修正视线偏角ψLp,i：
若当前分段处于转弯段，则根据
Figure FDA0004231246440000072
Figure FDA0004231246440000073
计算初始视线偏角变化率
Figure FDA0004231246440000074
和第一修正视线偏角变化率
Figure FDA0004231246440000075
若当前分段处于直线段，则根据
Figure FDA0004231246440000076
Figure FDA0004231246440000077
计算初始视线偏角变化率
Figure FDA0004231246440000078
和第一修正视线偏角变化率
Figure FDA0004231246440000079
根据公式
Figure FDA00042312464400000710
得到第一修正视线偏角ψLp,i。
CN202111676884.9A 2021-12-31 2021-12-31 一种无动力飞行器的三维制导系统 Active CN114326814B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111676884.9A CN114326814B (zh) 2021-12-31 2021-12-31 一种无动力飞行器的三维制导系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111676884.9A CN114326814B (zh) 2021-12-31 2021-12-31 一种无动力飞行器的三维制导系统

Publications (2)

Publication Number Publication Date
CN114326814A CN114326814A (zh) 2022-04-12
CN114326814B true CN114326814B (zh) 2023-06-16

Family

ID=81022823

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111676884.9A Active CN114326814B (zh) 2021-12-31 2021-12-31 一种无动力飞行器的三维制导系统

Country Status (1)

Country Link
CN (1) CN114326814B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115079723B (zh) * 2022-07-20 2024-06-14 中国人民解放军国防科技大学 一种任意时间到达的固定翼无人机制导方法
CN115560638B (zh) * 2022-09-02 2023-05-09 北京理工大学 基于修正量与控制时间关系的抗卫星失效弹道控制方法

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104077469A (zh) * 2014-05-28 2014-10-01 中国人民解放军海军航空工程学院 基于速度预测的分段迭代剩余时间估计方法
CN107422744A (zh) * 2017-05-02 2017-12-01 中国科学院声学研究所 一种基于径向速度控制的交会时间控制方法
CN107943079A (zh) * 2017-11-27 2018-04-20 西安交通大学 一种剩余飞行时间在线估计方法
CN111474948A (zh) * 2019-12-25 2020-07-31 中国人民解放军海军潜艇学院 一种带时间控制的前置导引与姿态控制制导的方法
CN111506101A (zh) * 2019-10-21 2020-08-07 北京理工大学 基于通信网络拓扑结构的飞行器协同制导控制方法及系统

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104077469A (zh) * 2014-05-28 2014-10-01 中国人民解放军海军航空工程学院 基于速度预测的分段迭代剩余时间估计方法
CN107422744A (zh) * 2017-05-02 2017-12-01 中国科学院声学研究所 一种基于径向速度控制的交会时间控制方法
CN107943079A (zh) * 2017-11-27 2018-04-20 西安交通大学 一种剩余飞行时间在线估计方法
CN111506101A (zh) * 2019-10-21 2020-08-07 北京理工大学 基于通信网络拓扑结构的飞行器协同制导控制方法及系统
CN111474948A (zh) * 2019-12-25 2020-07-31 中国人民解放军海军潜艇学院 一种带时间控制的前置导引与姿态控制制导的方法

Also Published As

Publication number Publication date
CN114326814A (zh) 2022-04-12

Similar Documents

Publication Publication Date Title
CN108036676B (zh) 一种基于三维再入弹道解析解的全射向自主再入制导方法
CN114326814B (zh) 一种无动力飞行器的三维制导系统
CN104778376B (zh) 一种临近空间高超声速滑翔弹头跳跃弹道预测方法
CN105573340B (zh) 一种固定翼无人机抗侧风的飞行控制方法
CN105486308A (zh) 估计弹目视线角速率的快收敛Kalman滤波器的设计方法
CN114020019A (zh) 飞行器的制导方法与装置
CA2699137A1 (en) Hybrid inertial system with non-linear behaviour and associated method of hybridization by multi-hypothesis filtering
CN104121930B (zh) 一种基于加表耦合的mems陀螺漂移误差的补偿方法
CN106406333A (zh) 一种基于积分型终端滑模的平流层飞艇俯仰角跟踪方法
CN115248038B (zh) 一种发射系下的sins/bds组合导航工程算法
CN114964226B (zh) 噪声自适应强跟踪扩展卡尔曼滤波器四旋翼姿态解算方法
CN114326813B (zh) 一种无动力飞行器的剩余飞行时间预测方法和系统
CN112034879B (zh) 一种基于高度-射程比的标准轨迹跟踪制导方法
CN110647161B (zh) 基于状态预测补偿的欠驱动uuv水平面轨迹跟踪控制方法
CN113110527A (zh) 一种自主水下航行器有限时间路径跟踪的级联控制方法
CN109445283B (zh) 一种用于欠驱动浮空器在平面上定点跟踪的控制方法
CN108398883B (zh) 一种rlv进场着陆轨迹快速推演及确定方法
CN115268475B (zh) 基于有限时间扰动观测器的机器鱼精确地形跟踪控制方法
CN106483967B (zh) 一种基于角速度信息测量与滑模的飞艇俯仰角稳定方法
CN115342815A (zh) 反大气层内或临近空间机动目标视线角速率估计方法
CN109634301B (zh) 一种考虑视场角限制以及运动不确定性的结合记忆的旋翼飞行器高速飞行避障方法
CN102508819B (zh) 基于角速度的飞行器极限飞行时四元数勒让德近似输出方法
CN102495830B (zh) 基于角速度的飞行器极限飞行时四元数Hartley近似输出方法
KR102114051B1 (ko) 무게중심의 이동을 고려한 항공기의 비선형 제어방법
CN102495831B (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