CN110243359B - 基于低空风预测模型的安全航迹规划方法 - Google Patents
基于低空风预测模型的安全航迹规划方法 Download PDFInfo
- Publication number
- CN110243359B CN110243359B CN201910467594.XA CN201910467594A CN110243359B CN 110243359 B CN110243359 B CN 110243359B CN 201910467594 A CN201910467594 A CN 201910467594A CN 110243359 B CN110243359 B CN 110243359B
- Authority
- CN
- China
- Prior art keywords
- wind
- distance
- point
- aircraft
- node
- 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
Links
- 238000000034 method Methods 0.000 title claims abstract description 149
- 239000013598 vector Substances 0.000 claims abstract description 63
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 43
- 239000011159 matrix material Substances 0.000 claims abstract description 30
- 238000012545 processing Methods 0.000 claims abstract description 6
- 238000012216 screening Methods 0.000 claims abstract description 4
- 230000008569 process Effects 0.000 claims description 29
- 230000000694 effects Effects 0.000 claims description 13
- 230000006870 function Effects 0.000 description 36
- 238000012937 correction Methods 0.000 description 26
- 230000009194 climbing Effects 0.000 description 25
- 230000008859 change Effects 0.000 description 19
- 238000004458 analytical method Methods 0.000 description 16
- 238000004088 simulation Methods 0.000 description 13
- 238000004364 calculation method Methods 0.000 description 10
- 238000011160 research Methods 0.000 description 6
- 230000001174 ascending effect Effects 0.000 description 5
- 238000006073 displacement reaction Methods 0.000 description 5
- 238000004891 communication Methods 0.000 description 4
- 230000003247 decreasing effect Effects 0.000 description 4
- 101100098799 Magnaporthe oryzae (strain 70-15 / ATCC MYA-4617 / FGSC 8958) TAS2 gene Proteins 0.000 description 3
- 101150002629 TAS1 gene Proteins 0.000 description 3
- 238000010586 diagram Methods 0.000 description 3
- VYQNWZOUAUKGHI-UHFFFAOYSA-N monobenzone Chemical compound C1=CC(O)=CC=C1OCC1=CC=CC=C1 VYQNWZOUAUKGHI-UHFFFAOYSA-N 0.000 description 3
- 238000005457 optimization Methods 0.000 description 3
- 238000012876 topography Methods 0.000 description 3
- 241000255925 Diptera Species 0.000 description 2
- 238000013473 artificial intelligence Methods 0.000 description 2
- 230000015572 biosynthetic process Effects 0.000 description 2
- 230000007423 decrease Effects 0.000 description 2
- 230000007547 defect Effects 0.000 description 2
- 238000001514 detection method Methods 0.000 description 2
- 238000011161 development Methods 0.000 description 2
- 239000000446 fuel Substances 0.000 description 2
- 230000002068 genetic effect Effects 0.000 description 2
- 238000005259 measurement Methods 0.000 description 2
- 239000002245 particle Substances 0.000 description 2
- 230000002194 synthesizing effect Effects 0.000 description 2
- 230000000007 visual effect Effects 0.000 description 2
- LFQSCWFLJHTTHZ-UHFFFAOYSA-N Ethanol Chemical compound CCO LFQSCWFLJHTTHZ-UHFFFAOYSA-N 0.000 description 1
- 206010034719 Personality change Diseases 0.000 description 1
- 230000001133 acceleration Effects 0.000 description 1
- 230000003044 adaptive effect Effects 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 238000013528 artificial neural network Methods 0.000 description 1
- 238000007664 blowing Methods 0.000 description 1
- 239000003086 colorant Substances 0.000 description 1
- 238000010835 comparative analysis Methods 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 239000005043 ethylene-methyl acrylate Substances 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 238000010801 machine learning Methods 0.000 description 1
- 239000000463 material Substances 0.000 description 1
- 238000002715 modification method Methods 0.000 description 1
- 238000012544 monitoring process Methods 0.000 description 1
- 230000035515 penetration Effects 0.000 description 1
- 230000002093 peripheral effect Effects 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 230000000630 rising effect Effects 0.000 description 1
- 238000000926 separation method Methods 0.000 description 1
- 230000003746 surface roughness Effects 0.000 description 1
- 230000002459 sustained effect Effects 0.000 description 1
- 238000001308 synthesis method Methods 0.000 description 1
- 101150052846 tas3 gene Proteins 0.000 description 1
- 238000012795 verification Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
- G01C21/00—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00
- G01C21/005—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 with correlation of navigation data from several sources, e.g. map or contour matching
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
- G01C21/00—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00
- G01C21/20—Instruments for performing navigational calculations
- G01C21/203—Specially adapted for sailing ships
-
- 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"
- G06Q10/047—Optimisation of routes or paths, e.g. travelling salesman problem
-
- 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
- G06Q50/00—Information and communication technology [ICT] specially adapted for implementation of business processes of specific business sectors, e.g. utilities or tourism
- G06Q50/10—Services
- G06Q50/26—Government or public services
- G06Q50/265—Personal security, identity or safety
-
- 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
- Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
- Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
- Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation
Landscapes
- Engineering & Computer Science (AREA)
- Business, Economics & Management (AREA)
- Remote Sensing (AREA)
- Radar, Positioning & Navigation (AREA)
- Physics & Mathematics (AREA)
- Human Resources & Organizations (AREA)
- General Physics & Mathematics (AREA)
- Strategic Management (AREA)
- Tourism & Hospitality (AREA)
- Economics (AREA)
- Theoretical Computer Science (AREA)
- Automation & Control Theory (AREA)
- Development Economics (AREA)
- General Business, Economics & Management (AREA)
- Marketing (AREA)
- Computer Security & Cryptography (AREA)
- Game Theory and Decision Science (AREA)
- Educational Administration (AREA)
- Health & Medical Sciences (AREA)
- General Health & Medical Sciences (AREA)
- Primary Health Care (AREA)
- Entrepreneurship & Innovation (AREA)
- Quality & Reliability (AREA)
- Operations Research (AREA)
- Traffic Control Systems (AREA)
Abstract
本发明公开一种基于低空风预测模型的安全航迹规划方法,步骤是:将救援环境的三维网格的网格中心点作为航空器的飞行路径节点,以航空器的垂直安全距离和侧向安全距离为准,对飞行路径节点进行筛选,得到可通行矩阵A;引入距离权重,得到带距离权重的邻接矩阵C;采用Dijkstra算法对邻接矩阵C求解最短路径,途经的节点即为最短航迹节点;利用基于地形坡度的低空风场预测方法得到每个最短航迹节点的风速和风向,结合航空器的航速和航向,得到实际航迹;判断实际航迹是否满足航空器飞行的垂直安全距离和侧向安全距离,若满足则不作处理,若不满足则进行航迹修正,最终得到低空救援安全航迹。此种方法能够准确体现地形坡度对风矢量的影响,保证航迹安全。
Description
技术领域
本发明涉及一种基于低空风预测模型的安全航迹规划方法。
背景技术
目前,在发达国家当中,航空应急救援已经形成比较完善的救援体系,而我国的航空救援体系的发展起步较晚,但随着通用航空产业的发展以及目前救援方式弊端的日趋凸显,以直升机、无人机为主的空中救援模式将成为未来抗震救灾的主要方式。由于救援活动通常包括目视搜救、伤员运输和物资投放等行动,通常需要用于空中救援的通用航空器保持较低高度进行大范围飞行,需要实时的飞行冲突解脱[1],风险系数高;另一方面,发生灾害的地区又通常为地形地貌较为复杂的山地区域或偏远地区,地面通信导航服务难以对救援航空器提供有效支持和保障,甚至无法覆盖低空救援区域。此时,低空风场以及复杂地形将严重威胁救援航空器的飞行安全,救援效率也严重受限。因此,为了降低飞行触地或与障碍物相撞的可能性,提高救援效率,考虑低空风与地形的安全航迹规划对保证救援航空器安全运行有着重要意义。
低空救援航迹规划的研究,主要围绕如下三类方法开展:
(1)基于栅格的航迹规划技术。基于栅格的算法主要包括RRT算法、A*算法、D*算法、Dijkstra算法、动态规划法。该类方法通过对三维空间栅格化基础上,借助上述算法实现三维航迹规划方法,使生成的航迹可以满足包括最低飞行高度、最大拐爬升/下滑角等在内的各种约束条件,实现飞行冲突回避。如Bazhenov[2]提出了一种从当前飞机位置到预定轨迹变化点的无冲突轨迹形成的算法,以解决航空器和地形冲突问题。算法提供了预定轨迹变化点到达时间同步与飞机前方和安全的高空区域的地形,并提出了一种沿这些轨迹计算航空器速度的算法。Yin[3]研究了在跨大西洋飞行中避免形成持续的航迹时,横向和纵向变化对飞行轨迹的影响。利用复杂的地球系统模型(EMAC),对飞行时间和飞行距离进行航迹优化。但是这种航迹规划方法更适用于民用航空的航班轨迹规划,且运用场景中没有考虑地形阻碍,与低空救援的情景特点有较大区别。Franco[4]研究了受风不确定性影响的飞机轨迹预测问题,提出飞行时间和燃料消耗的概率分析。该方法侧重于对燃油消耗的分析,这在低空救援任务中不是首要考虑因素。Li[5]提出了一种规划无人机在已知室内环境下的通用路径模型。该方法的特点在于得出的航迹规划路径与障碍保持一定距离,但是这样的航迹路径放在低空救援环境中,却不一定是航程距离最短的路径,忽视了救援效率。
(2)基于人工智能的航迹规划技术。该类方法将现代人工智能算法用于航迹规划的避障研究中,主要算法涉及启发式搜索[6]、神经网络、遗传算法[7]、机器学习[8,9]、粒子群算法[10]等。如Wang[6]提出了一种具有威胁和地形障碍规避的远程飞行实时混合路径规划方法。该方法针对一系列成本因素和约束条件,寻找最优路径。该算法仅在需要时调用,从而节省计算资源,降低了时间成本,但也因此只探测航线前方有限距离的威胁规避,所以是一种实时动态规划方法,而不能作为预先规划方法。Kulida[7]针对机载航迹安全系统速度控制与其他飞机无法达到分离标准的情况,提出了一种考虑复杂地形的、长度一定的低空飞行航迹规划的遗传算法,最后给出了该算法在耶利佐沃机场附近山区的仿真轨迹。Chen[10]提出了一种基于dubins曲线和粒子群算法的路径规划方法。但是,该方法运用在低空救援的三维空间航迹规划中,降维处理会损失大量算法求解空间,降低获取最优路径的可能性。Li[11]提出了一种改进的无人机低空突防弹道规划概率路线图(PRM)方法,该方法的威胁对象是地对空导弹,判断威胁依据一种概率计算方法,由于航空器性能受限,以及出于安全考虑,近地贴地的航迹规划结果不适用于大部分救援航空器。
(3)基于图形的航迹规划技术。基于图形的算法主要包括Voronoi图法、通视图法、子目标网络法、随机路线图法等。该方法可通过上述算法实现不同飞行策略约束下的飞行航迹规划,确定全局最优航迹。Aleshin[12]提出了一种飞机飞行预测、潜在冲突态势检测和冲突解决路径可实施性评估的方法。然而这种方法侧重于原定航迹的局部调整,对于其他潜在更优航迹缺少了全局搜索的考虑。Takeichi[13]提出了一种基于飞行和气象条件的水平飞行时间不确定性自适应预测模型。利用二次监视雷达S模式采集的实际飞行数据和数值天气预报进行处理,获得大量飞行时间误差和飞行气象条件。根据不确定度传播规律,建立了飞行时间不确定度随马赫数、飞行距离、风速和温度变化的自适应预测模型。但是,该方法需要依赖二次雷达对实时数据的采集才能实现[14],而在低空救援环境中,气象雷达、通信支持保障服务难以覆盖到地形复杂、高度较低的低空空域。
低空救援航空活动具有飞行高度低、受到低空天气因素的影响比较大、活动空间和范围受限于高度层限制以及地形复杂等特点。以往的研究主要侧重于飞行中的威胁回避、地形跟踪方面的航迹规划和实时冲突探测与解脱,所采用的方法及考虑的情景因素未能突显低空救援活动的特征。在真实的低空救援环境当中,地形环境对救援航空器的飞行路线影响非常大,低空风场对救援航空器的航迹也会产生诸如偏航等安全影响。
涉及的参考文献:
[1]Zhang M,Yu J,Zhang Y,Wang S,Yu H,Flight conflict resolution duringlow-altitude rescue operation based on ensemble conflict models,Advances inMechanical Engineering,2017,9(4),1–16.
[2]Bazhenov S G,Egorov N A,Kulida E L,et al.Control of aircrafttrajectory and speed to avoid terrain and traffic conflicts during approachmaneuvering.Automation&Remote Control,2016,77(10),1827-1837.
[3]Yin F,Grewe V,Christine,et al.Impact on flighttrajectory characteristics when avoiding the formation of persistentcontrails for transatlantic flights.Transportation Research Part D:Transportand Environment,2018,65,466-484.
[4]Franco A,Rivas D,Valenzuela A,et al.Probabilistic aircrafttrajectory prediction in cruise flight considering ensemble windforecasts.Aerospace Science and Technology,2018,350-362.
[5]Fangyu Li,Sisi Zlatanova,Martijn Koopman,Xueying Bai,AbdoulayeDiakité.Universal path planning for an indoor drone.Automation inConstruction,2018,98,275-283.
[6]Wang H,Li Q,Cheng N.Real-time path planning for low altitudeflight based on A*algorithm and TF/TA algorithm,in:2012IEEE InternationalConference on Automation Science and Engineering(CASE),Seoul,South Korea,20-24Aug.,2012.
[7]Kulida E L,Lebedev V G.Genetic Algorithm for GeneratingTrajectories of Specified Length for the Onboard Flight Path Safety System,Procedia Computer Science,2017,112,1015-1022.
[8]Qi P,Zhao X,Wang Y,et al.Aeroelastic and trajectory control ofhigh altitude long endurance aircraft,IEEE Transactions on Aerospace andElectronic Systems,2018,54,6,2992-3003.
[9]Zhou LH,Zhang JQ,et al.A fuzzy adaptive backstepping control basedon mass observer for trajectory tracking of a quadrotor UAV,InternationalJournal of Adaptive Control and Signal Processing,2018,32(12),1675-1693.
[10]Chen QY,Lu YF,et al.Path planning for UAVs formationreconfiguration based on Dubins trajectory,Journal of Central SouthUniversity,2018,25(11),2664-2676.
[11]Li Q,Wei C,Wu J,Zhu X.Improved PRM Method of Low AltitudePenetration Trajectory Planning for UAVs,In:Proceedings of 2014 IEEE ChineseGuidance,Navigation and Control Conference,Yantai,China,8-10 Aug.,2014.
[12]Aleshin B.S.,Bazhenov S.G.,Lebedev V.G.,and Kulida E.L.Assessmentof Implementability and Safety of Aircraft Flight Paths by an OnboardMathematical Model,Automation and Remote Control,2014,75,4,745–754.
[13]Takeichi N.Adaptive prediction of flight time uncertainty forground-based4D trajectory management,Transportation Research Part C-emergingTechnologies,2018,95,335-345.
[14]Delahaye D,Rabut C,Stéphane Puechmorel.Wind Field Evaluation byUsing Radar Data and Vector Spline Interpolation,in:2011 9th IEEEInternational Conference on Control and Automation(ICCA),Santiago,Chile,19-21Dec.,2011.
[15]Cao S,Wang T,Ge Y,et al.Numerical study on turbulent boundarylayers over two-dimensional hills-effects of surface roughness and slope,Journal of Wind Engineer and Industrial Aerodynamics,2012,104-106,342-349.
[16]Cao S,Tamura T.Experimental study on roughness effects onturbulent boundary layer flow over a two-dimensional steep hill,Journal ofWind Engineering and Industrial Aerodynamics,2006,94,1,1-19.
[17]Li W,Ding P,Duan C,Qiu R,Lin J,Shi X.Comparison of spatialinterpolation approaches for in-core power distribution reconstruction,Nuclear Engineering and Design,2018,337,66-73.
[18]Lavrinenko A.V.,Moldovanova E.A.,Mymrina D.F.,Popova A.I.,PopovaK.Y.,Popov Y.B.,Spatial interpolation of meteorological fields using amultilevel parametric dynamic stochastic low-order model,Journal ofAtmospheric and Solar-Terrestrial Physics,2018,181,Part A,38-43.
[19]Foehn A,Hernández J G,Schaefli B,Cesare G D,Spatial interpolationof precipitation from multiple rain gauge networks and weather radar data foroperational applications in Alpine catchments,Journal of Hydrology,2018,563,1092-1110.
[20]Francisco González-Longatt,Humberto Medina,Javier Serrano González,Spatial interpolation and orographic correction to estimate wind energyresource in Venezuela,Renewable and Sustainable Energy Reviews,2015,48,1-16.
[21]EbubekirGüler,Seyit Ahmet/>Investigation of windshear coefficients and their effect on electrical energy generation,AppliedEnergy,2011,88,11,4097-4105.
[22]Maduako E.Okorie,Freddie Inambao,Zivayi Chiguvare,Evaluation ofWind Shear Coefficients,Surface Roughness and Energy Yields over InlandLocations in Namibia,Procedia Manufacturing,2017,7,630-638.
[23]Bazhenov S G,Egorov N A,Kulida E L,et al.Generation of Traffic/Terrain Conflict Free Trajectories and Their Analysis by Means of the On-board Model of the Airplane,IFAC-PapersOnLine,2016,49(12),1430-1435.
发明内容
本发明的目的,在于提供一种基于低空风预测模型的安全航迹规划方法,其同时考虑距离、地形坡度、高程对风矢量的影响,能够准确体现地形坡度对风矢量的影响,考虑不同风向的风场对航空器航行轨迹的影响,进行航迹修正,保证航迹安全。
为了达成上述目的,本发明的解决方案是:
一种基于低空风预测模型的安全航迹规划方法,包括如下步骤:
步骤1,建立救援环境的三维网格,将网格中心点作为救援航空器的飞行路径节点,以航空器的垂直安全距离和侧向安全距离为准,对飞行路径节点进行筛选,得到可通行矩阵A;
步骤2,引入距离权重,得到带距离权重的邻接矩阵C;
步骤3,采用Dijkstra算法对带距离权重的邻接矩阵C求解最短路径,途经的节点即为最短航迹节点;
步骤4,利用基于地形坡度的低空风场预测方法得到每个最短航迹节点对应的风速和风向,结合航空器的航速和航向,从而得到航空器在各航迹节点的实际航迹;
步骤5,判断步骤4得到的实际航迹是否满足航空器飞行的垂直安全距离和侧向安全距离,若满足则不作处理,若不满足则进行航迹修正,最终得到低空救援安全航迹。
采用上述方案后,本发明在低空救援安全航迹规划中主要关注低空风和地形地貌两个方面。在低空风预测模型构建方面中,分析处于不同坡度的山体风的规律得出拟合函数式,并与反距离插值方法结合,得到的风场内插方法可同时考虑距离、山体坡度、高程对风矢量的影响。与以往的低空航迹规划不同,本发明考虑山地地形坡度预测低空风,提出构建保证航空器安全的规划航迹方法,具有如下创新:
(1)针对传统反距离风场内插方法忽略了山体坡度对风矢量的影响,分析风矢量在不同坡度的山体上的变化规律,提出基于坡度的风场内插方法,该方法的插值结果能够更好的体现地形坡度对风矢量的影响,使得风场插值结果更贴近实际,实验结果表明,本发明方法对于风速的预测误差最大降低了13.37%,平均降低2.74%;
(2)针对低空风场影响和地形阻碍限制等问题,分析不同风向的风场对航空器航行轨迹的影响,并得出相应的航路节点修正方法以保证航迹安全;将Dijkstra算法与三维网格空域划分方法相结合构成改进的Dijkstra算法,对该算法的联通矩阵求解做优化,解决了二维平面算法难以求解三维地形阻碍空间的最优路径的问题;结合路径求解算法和航路节点修正方法,得出低空救援安全航迹规划方法,并通过算例对比验证本发明方法的安全性。
附图说明
图1是风速在缓坡上随高度变化情况示意图;
图2是风速在陡坡上随高度变化情况示意图;
图3是竖直风矢量对飞行状态的影响示意图;
图4是航空器在顺风风场下的起飞轨迹;
图5是航空器在顺风风场下的降落轨迹;
图6是水平侧风对航空器的影响示意图;
图7是基于Dijkstra算法的轨迹修正方法流程图;
图8是山体轮廓剖面图;
图9是阿坝汶川周边地形仿真图;
图10是风矢量待测点位置图;
图11是救援区域三维网格;
图12是无风情形下的救援航迹;
其中,(a)是主视图,(b)是左视图;
图13是低空风影响下的救援轨迹;
图14是冲突隐患位置示意图;
图15是修正后的安全规划航迹示意图;
图16是航迹规划对比图;
图17是两处冲突隐患点示意图。
具体实施方式
1低空风预测模型的构建
1.1现有风场内插方法的原理分析
Cao[15,16]提出了考虑地形因素的风场内插方法,此方法是在反距离加权法的基础上引入地形高程变化度的因子,从而构造新的权重函数。权重函数形式如下:
公式中的w表示权重函数,r表示两点之间的距离,h表示两点之间的高程差,a和b分别表示r和h的指数。两点之间距离较近,根据风场的待求点与观测点的相对位置,该公式表示的情况可以分为以下几种:
第一种:两点之间距离较近,且地势变化不大;
第二种:两点之间距离较远,且地势变化不大;
第三种:两点之间距离较近,且地势变化较大;
第四种:两点之间距离较远,且地势变化较大。
通过分析可知,两点之间地势变化不大时,随着距离的贴近,两者的风速及风向发生变化的可能性越小,相关性越大;而当两点之间的距离一定时,地势起伏变化越大,两者的风速及风向不相同的可能性越高,相关性越小。因此在这四种情况中,观测点与待求点的相关性应该是依次递减的,为了方便比较第二种和第三种情况的权重系数,需要对权重函数的距离因子和地势变化因子去量纲化。
上式中,rmax和hmax为所有r和h的最大值。为了保证第三种情况的权重值小于第二种,需要满足a<b。在反距离加权法中,a通常取2,故此处只需要b取3即可满足要求[17,18]。
通过分析可以发现,该权重函数主要考虑了两点之间的距离及绝对高度差对风矢量的影响,而在实际情况中,山地起伏的坡度也是影响风速变化的一个重要因素,而式(2)中并没有体现这一点。本文把坡面的垂直高度和水平方向的距离比称作坡度。通过阅读大量相关文献[15,18,19],发现山体坡度对风速有如下影响:当山坡坡度在0.15-0.65时,视作缓坡考虑,此时山脚风速沿山坡逐渐增长,在山顶处风速达到最大值,风速最大值约为山脚风速的两倍;当山坡坡度在0.5-0.75时,视作陡坡考虑,山脚风速沿山坡逐渐增长,在山顶处风速达到最大值,风速最大值约为山脚风速的1.8倍;但是要想将坡体对风速的影响放入插值中进行改进,还需要进一步寻找风速、山体坡度以及坡体高度的关系。
1.2基于山体坡度的风矢量规律分析
参考文献[15,16]的研究方法,通过收集汶川附近的风矢量数据,分析不同山体坡度与风矢量的关系。现收集汶川气象测风塔中的风矢量数据,记录的风矢量分布于缓坡以及陡坡,其中缓坡的坡度分别有0.3/0.4/0.5/0.65,而陡坡的坡度分别有0.8/1/1.2/1.5。分布于缓坡的风矢量,位于山坡坡底处的风速基本在14m/s至17m/s的范围内,而分布于缓坡的风矢量,位于山坡坡底处的风速基本在26m/s至30m/s的范围内;为了方便讨论,我们把位于山坡坡底处的风速称作入口风速,通过观察可以发现,缓坡的入口风速远小于陡坡的入口风速,这是由于缓坡坡体附近地形较平缓,而陡峭坡体所处的地形较前者更起伏且复杂,根据伯努利原理,狭窄陡峭地带更容易形成“穿堂风”,这是造成两者入口风速差异的主要原因。但是入口风速的不同不会对接下来的分析结论造成太大影响,因为我们主要关注风速在坡体上的变化趋势。为了便于分析,下面给出风矢量在不同坡度的坡体随高度上升的变化过程图。
图1给出的风速变化过程中,可以看出入口风在缓坡上升至150m高度前,风速一直在呈线性递增;在高度150m至550m范围内,风速基本维持在一个稳定水平;当高度大于550m时,风速呈线性递减。其中代表坡度0.3的蓝色曲线由于坡体较平缓,在最高处附近的风矢量没有形成涡流或紊乱现象,故风速继续保持在20.3m/s左右。图中不同颜色的实线对应不同坡度下风速随高度变化的趋势,黑色虚线是根据采集的风数据拟合得到的曲线,拟合函数为:
f(x)=a1×sin(b1×x×c1)+a2×sin(b2×x×c2)+a3×sin(b3×x×c3) (3)
其中a1=47.53,b1=0.007575,c1=-1.14,a2=49.3,b2=0.01263,c2=0.1726,a3=22.4,b3=0.01541,c3=2.298。该拟合函数的决定系数R2=0.9849,决定系数越接近1,则拟合效果越好,和方差SSE=0.2809,均方根误差RMSE=0.3748,总体误差较低。
图2给出的风速变化过程中,可以发现坡度为0.8情况下的入口风速与其他三条虚线的入口风速差别较大,这是由于周边地形的差异造成入口风速的不同,但不影响对比分析风速变化趋势。当坡度为0.8时,在风矢量上升到高度150m之前,呈线性递增趋势,而在高度150m至550m范围内,风速缓慢下降并维持在20m/s附近;当高度大于550m时,风速呈线性递减趋势。另外三种坡度下的风矢量则在高度100m至200m范围内呈线性递减,其余变化趋势与坡度为0.8的情况基本一致。图中唯一黑色实线是根据采集的风数据拟合得到的曲线,拟合函数为:
f(x)=d1×sin(e1×x×f1)+d2×sin(e2×x×f2)+d3×sin(e3×x×f3) (4)
其中d1=139.8,e1=0.00383,f1=1.448,d2=125.7,e2=0.006422,f2=3.735,d3=28.76,e3=0.009827,f3=5.668。该拟合函数的决定系数R2=0.9975,决定系数越接近1,则拟合效果越好,和方差SSE=0.1287,均方根误差RMSE=0.2537,总体误差较低。
对于坡度与风速之间关系的讨论,得到了两种拟合函数式(3)和式(4),两种函数可以高度近似地表达风矢量在平缓和陡峭地形上的变化趋势,该函数可以运用到风矢量插值模型当中,从而使得风矢量插值模型可以同时体现地形坡度、地形高程变化、两点之间的距离三个主要因素对风速的影响。
1.3考虑地形因素的风场内插方法改进
综上,得到的风矢量插值方法为:
(1)确定风矢量参考点j,确定风矢量待求点i,确定地形基准面及其高程;
(2)计算每个待求点所在坡体的坡度,坡度由待求点i所在高度Hi以及该点与基准点在主要风向上的水平距离Li共同确定,其中基准点指该点所在山坡的坡脚外延轮廓与主风向延长线的交点;
(3)若待求点所在坡体坡度小于0.2,则视为平原地区考虑[17],直接利用反距离权重插值法计算风矢量,权重函数是在式(3)的基础上,去除高程因子得出,并且距离因子的幂指数取a=2。确定权重函数后,代入风场内插公式,具体如下:
(4)若待求点所在坡体坡度大于0.2且小于0.7,则视为缓坡地区考虑,若坡度大于0.7,则视为陡坡地区考虑[17]。计算缓坡或陡坡地区的风矢量时,先用反距离权重插值法计算坡体山脚处风速,再利用相应的拟合函数计算坡体上的风速,拟合函数见式(3)和式(4),实际待测点在山坡上的高度位置对应拟合函数中相对高度,函数值的结果需要进行对应插值修正,修正值为UX-|Ucal|或VX-|Vcal|,其中UX、VX表示在山坡相对位置的南北风分量或东西风分量,最终风矢量正负号与参考点风矢量的正负号一致;
(5)近地风速计算完毕后,需要计算低空风速,而风随海拔高度变化一般遵循指数律公式或对数律公式[20],此处我们采用指数律公式估算不同高度的低空风速:
式中:Ui、Uj分别为处于高度Zi、Zj的风速,m为风随高度的切变系数,此处的切变系数由地表情况决定,详情如表1[21-22]:
表1不同地表情况对应的切变系数
至此,改进的风场内插方法通过将传统的反距离内插方法以及基于坡度的拟合函数相结合,可以同时考虑到地形高程差、坡度、距离以及地面粗糙度对风矢量的影响。
2低空风影响下的航迹修正
救援区域的风场数据可以通过预测模型得出,而初步安全航迹可以基于三维网格的建立以及Dijkstra算法求得,此外还需要将地形、风场两个方面结合考虑,判断航空器在飞行时受低空风影响下所产生的航迹是否能够规避地形,确保航迹的安全。依据风向,此处把低空风场分为水平风场和竖直风场,在以往研究中,主要考虑水平侧风对航空器的影响,故接下来主要分析竖直风场以及水平顺逆风对航空器航迹的影响,从而得出为确保航迹安全所需的飞行角度、航空器应采取的空速以及计划航路节点的修正距离,此结果将作为后续调整航迹节点的计算依据。
2.1竖直风场下的航迹修正
竖直方向的风矢量可分为上升气流和下沉气流,它们不改变航空器的平飞速度,但是会影响航空器的爬升率以及下降率。航空器在上升气流中爬升,其上升角和爬升率都会增大,反之,在下沉气流中上升,会使上升角和爬升率减少;在其他飞行阶段,竖直方向的气流均会对航空器的飞行状态产生类似影响,具体可见图3。竖直风场下航迹修正中使用的相关参数见表2。
表2竖直风场下航迹修正相关参数表
在上升阶段,若航空器遭遇下沉气流,则需要以空速VTAS1和爬升角θ1飞行,才可以达到实际爬升空速VGSC和实际爬升角度θ的航行轨迹;若航空器遭遇上升气流,则需要以空速VTAS2和爬升角θ2飞行,才可以达到实际爬升空速VGSC和实际爬升角度θ的航行轨迹。在平飞阶段和下降阶段,均可按照上述逻辑修正航空器的航行轨迹。然而要想达到这样的修正目的,需要确定空速VTAS1或VTAS2、以及爬升角θ1或θ2,具体计算过程如下:
以航空器在下沉气流中的爬升阶段为例。设航空器预计到达的位置点所要求的时间为T,预计到达的位置点距离起始点的高度为H,距离起始点的水平距离为L,那么航空器从起始点至预计位置点的实际爬升距离为S,实际爬升速度为VGSC,实际爬升角度为θ。现在需要求:应采取的空速VTAS1,应采取的爬升角度θ1。
对于航速矢量和空速矢量组成的三角形,运用余弦定理可得下列关系:
类似的,航空器在上升气流中的爬升阶段,要使得航空器以实际爬升空速VGSC和实际爬升角度θ到达位置点,那么应采取的空速VTAS2和应采取的爬升角度θ2可以通过下列公式求出:
航空器在下沉气流中的平飞阶段,已知VGSL和Vsink,那么要使得航空器以平飞空速VGSL并且水平飞行到达预计位置点,那么应采取的空速VTAS3和应采取的爬升角度θ3可以通过下列公式求出:
航空器在上升气流中的平飞阶段,已知VGSL和Vrise,那么要使得航空器以平飞空速VGSL并且水平飞行到达预计位置点,那么应采取的空速VTAS4和应采取的下滑角度θ4可以通过下列公式求出:
假设航空器从开始下降点到预计位置点的高度差为ΔH,两点之间的水平距离差为ΔL,时间差为ΔT,要使得航空器以实际下降速度VGSF和实际下降角度θ0到达预计位置点,那么当航空器处于下沉气流中时,应采取的空速VTAS5和应采取的下滑角度θ5可以通过下列公式求出:
处于上升气流中时,要使得航空器以实际下降速度VGSF和实际下降角度θ0到达预计位置点,应采取的空速VTAS6和应采取的下滑角度θ6可以通过下列公式求出:
2.2水平顺逆风下的航迹修正
此处我们主要讨论水平风场内的顺风和逆风。水平顺风和水平逆风主要影响航空器的平飞速率,在航空器起降过程中,水平风场不改变其爬升率或下滑率。以往主要研究水平风场对航空器到达预计位置点的时间差并对其进行修正,而航空器在低空救援环境中,面临低空风以及起伏地形,仅仅修正航空器到达预计位置点的最终时间,是不足以保证其全过程的航迹安全的。因此,需要进一步分析航空器在水平顺风以及水平逆风的航迹。水平风场下航迹修正中使用的相关参数见表3。
表3水平风场下航迹修正相关参数表
在航空器的巡航阶段,航空器以巡航速度VTAS水平飞行,若进入顺风风场,则航空器的地速为VGS=VTAS+Vwind,起始点至预计位置点之间的水平距离为L,则该航段所需的飞行时间为T=L/(VTAS+Vwind);若进入逆风风场,则航空器的地速为VGS=VTAS-Vwind,起始点至预计位置点之间的水平距离为L,则该航段所需的飞行时间为T=L/(VTAS-Vwind)。
对于起飞阶段,在无风情形下,固定翼航空器的起飞过程为:由A点处开始滑跑加速,以速度VTAS爬升至B点。直升机则没有滑跑过程,直接爬升至C点。但若被顺风风场笼罩,则起飞的滑跑距离增加,地速增大,爬升距离增大,爬升角变小,此时起飞轨迹虽然由A点开始,但是当航空器飞至与B点相同的高度时,其真实位置点应在C点,这是因为在地面滑跑过程中受到顺风影响,航空器因升力损失需要加速,增加的速度大小为顺风风速大小,如此才可以满足升力损失,而在爬升阶段,航空器同样受顺风影响,水平位移距离会增加。实质上,不管是从升力补偿的角度来理解,还是从速度矢量和的角度来理解,航空器的水平位移都会增加。L2指的是航空器在地面滑跑阶段,风速Vwind在该时段的位移量;L3指的是航空器在爬升阶段,风速Vwind在该时段的位移量;C点与B点之间的距离为L2+L3。
显然,若在顺风风场的笼罩下,航空器仍然以无风情形下的起飞策略执行,不能继续认为航空器到达位置为B点,若前方存在类似山体等障碍物,如图4所示,则很可能导致航空器触地坠毁。因此需要对航空器的预先规划航迹进行修正,提升爬升率,或者将滑跑起始点提前,根据上述分析可知滑跑起始点提前量为L2+L3。现假设起飞滑跑所需时间为t1,爬升阶段所需时间为t2,则开始起飞滑跑位置点应提前的距离为:L2+L3=(t1+t2)Vwind。
对于降落阶段,在无风情形下,固定翼航空器的降落过程为:由D点处开始下降,以速度VTAS下滑至跑道,减速滑跑经过S的距离后停至E点。直升机则没有减速滑跑过程,直接降落至地面。但若被顺风风场笼罩,则在降落在跑道上的地速增大,下滑距离增大,下滑角变小,减速距离增加。此时起飞轨迹虽然由D点开始,但是航空器降落在跑道上的位置点不同,两个降落点之间距离S2;航空器在跑道上减速直至停止所需距离也不同,距离差为S1。同样的,在顺风风场的笼罩下,航空器不能以无风情形下的降落策略执行,否则很可能冲出跑道或者撞向前方障碍物,如图5所示。因此需要对航空器的预先规划航迹进行修正,加大下滑率,或者将降落起始点提前,根据上述分析可知提前量为S1+S2。现假设下滑至跑道所需时间为t3,滑跑减速阶段所需时间为t4,则下降起始点应提前的距离为:S1+S2=(t3+t4)Vwind。
与顺风起降的情况相反的是,逆风有利于提高航空器在起降过程中的安全性。与无风情形相比,逆风起降的爬升角或下滑角增大,起飞助跑距离或减速滑跑距离减少,航空器能够在更短的水平距离中完成起飞或降落过程。在降落过程中,由于航空器在逆风下降的轨迹有更高的安全性,故不需要改变飞行策略;而在起飞过程中,虽然航空器在逆风起飞的安全性比无风情形下更高,但是完成起飞进入平飞巡航阶段的起始点位置会提前,这对后续的整体航迹有影响,故需要对其修正,即航空器进入巡航阶段的起始点位置将会提前:(t1+t2)Vwind。
2.3水平侧风下的航迹修正
本文的航迹修正同样考虑了水平侧风对航空器的影响。由于以往课题组已经对水平侧风有充分研究,故此处只作简单分析。
图6中,AB为预定航迹,若在飞行过程中,航空器遭遇水平侧风Vwind,则会发生航迹偏移,为使航空器按照原定航迹飞行,应该改变航向,使航向偏向上风面。
此时航空器的地速VGS=VTAS×cosθTAS+Vwind×cosθwind,式中,VGS表示地速,VTAS表示真空速,Vwind表示风速,θTAS表示航向与航迹的夹角,θwind表示风向与航迹的夹角。水平侧风使得航空器造成的航迹偏移量为L=Vwind×(S/VTAS),其中S表示航空器从第一个航迹节点至第二个航迹节点之间的距离。关于水平风场的航迹修正,只需将航迹节点朝着上风向修正与偏移量相当的距离即可。
2.4基于Dijkstra算法的航迹修正方法
基于上述讨论分析,以下给出完整的低空安全航迹求解方法:
1、根据地形数据生成模拟地形,并在救援环境中建立三维网格,网格中心点将作为救援航空器的飞行路径节点;
2、以航空器的垂直安全距离和侧向安全距离为准,从所有的路径节点中筛选出可以通行的节点,并作为算法求解的空间。判断依据:每个节点到该节点坐标所在地面的距离应大于航空器的垂直安全距离,同时每个节点到所在高度层最近起伏地形的直线距离应大于航空器的侧向安全距离。保留所有满足判断依据的网格节点,记为可通行矩阵A;
3、由可通行矩阵A得到带距离权重的邻接矩阵C。矩阵元素表示对应的两个航路节点之间的距离,航路节点之间的距离由坐标点求得。此时得到的带有距离权重的矩阵记为B,判断矩阵B中的相邻航路节点是否联通。判断依据:相邻航路节点之间的x坐标之差不可大于三维网格的长;相邻航路节点之间的y坐标之差不可大于三维网格的宽;相邻航路节点之间的z坐标之差不可大于三维网格的高;三个条件同时满足则视为两节点联通。若两个航路节点之间不可直接联通,则距离无穷大;若两个航路节点之间可直接联通,则保留原距离。最终得到带距离权重的邻接矩阵C。
4、通过步骤1至步骤3,已经为Dijkstra算法运用到三维空间的求解做好铺垫。接下来运用Dijkstra算法对邻接矩阵C求解最短路径。初始化航路的起点s0和终点t,设已求出最短路径的节点集合为S,其余未确定最短路径的节点集合为U;
5、从U中选取一个距离s0最小的节点sk,把sk加入S中;
6、以sk为新考虑的中间点,修改U中各节点的距离;若从起点s0到终点t的距离(经过点sk)比原来距离(不经过顶点sk)短,则修改最短距离值;
7、重复步骤5和6直到所有节点都包含在S中,此时的最短距离值即为起点s0到终点t的最短距离,途径的航路节点即为最短航迹节点;
8、在每个航路节点上引入低空风,通过已知的气象站坐标和航路节点坐标,利用基于地形坡度的低空风场预测方法可以得到每个航路节点上的风速和风向信息;
9、根据风速、风向和航空器的航速和航向,可以得到航迹节点在低空风影响下的实际航迹;
10、判断受低空风影响后的实际航迹是否满足航空器飞行的垂直安全距离和侧向安全距离;
11、对于不满足要求的航迹节点,根据不同风向,使用前文提到的不同修正方法。若受到竖直风场影响导致航迹节点存在安全隐患,则将原航迹节点的前序航迹节点的航空器航向修正,具体策略为航空器航向朝上风向调整θ(θ根据风向和航空器状态取θ1-θ6,对应公式(10)、(12)、(14)、(16)、(18)、(20))。但调整的航向不得超过航空器的最大俯仰角,具体数值由航空器性能数据决定,一般情况下不超过±15度。(可参考民航局第195号令CCAR-121-R4第121.177条);若受到水平风场影响导致航迹节点存在安全隐患,则将原航迹节点朝着偏移相反方向调整相等的偏移距离,具体偏移距离已经在前文给出;
12、修正完毕,最终得到低空救援安全航迹,总体流程见图7。
3低空救援航迹规划仿真验证
3.1数字方程模拟地形
起伏的地形不仅会对航空器的飞行构成一定的阻碍,而且会极大地影响近地风速和风向,发生变化的风速及风向会改变航空器的飞行姿态,若飞行姿态变化剧烈,航空器的飞行安全将难以保证。因此,地形地貌特征在低空救援安全航迹规划中是不可回避的重要方面。
首先需要建立数字方程来模拟二维情形下的山体轮廓,常用的函数模型有高斯模型、正弦模型、余弦模型和贝尔模型[18]。与其他模型相比,贝尔模型在坡度的变化上比较缓慢,不便于较小的仿真地图上表示更广阔的地形地势环境,而正弦模型和余弦模型本质上可以通过平移横坐标得到彼此的图像,故此处选择高斯模型(21)和余弦模型(22),以下给出它们的函数方程:
f(x)=cos(x) (22)
当方程(21)取均值为2,方差为0的时候,可以使得两个图像的对称轴均为Y轴,其曲线与方程(22)的曲线对比如图8。图中点画线表示高斯模型,实线表示余弦模型,图像表示的区域是截取靠近波峰的位置,观察可以发现,高斯模型模拟的山体比余弦模型模拟的山体在坡度上更加平缓。而对于余弦模型,其山脚处位置在地势上出现急剧变化,因此在接下来的研究中,采用高斯模型模拟平缓山坡,采用余弦模型模拟陡峭山峰。
本次地形仿真选取阿坝汶川与成都、绵阳交接处的地形,选取范围在北纬30°45′至31°43′、东经102°51′至103°44′,地形数据来源于全球人造卫星大地测绘网站,通过MATLAB精细仿真可以得到成都、绵阳、汶川交界处的地形地貌。以东西走向建立x坐标,x坐标由西至东逐渐增加,以南北走向建立y坐标,y由北至南逐渐增加,坐标原点选取地点为:东距绵阳130公里、南距成都210公里。图9中经纬坐标通过插值调整改为横纵坐标表示,插值方法为三次样条插值,即给定的数据点区间范围内,每个数据节点之间的曲线可以用一个三次多项式函数表示,每个节点可以达到二阶连续,且要求正数第一个与第二个的多项式函数的三阶导数相当,根据此要求可求出每段三次多项式函数的表达式,从而得到函数曲线上任意位置的插值,曲线相连最终形成地形曲面。
阿坝汶川附近地形的最高峰海拔达到4200m左右,地形走势上,由西北至东南逐步降低,至成都和绵阳的交接开阔地带,地势海拔降低至2000m左右。此次仿真过程中,救援航空器起降基准面选取成都和绵阳的开阔平原地带,该地带海拔为2000m。
3.2基于地形坡度的低空风预测
从中央气象台获取阿坝汶川、成都、绵阳地区5月及9月风力及风向数据,部分数据来源于汶川雁门乡气象监测站,统计发现在汶川地区的最大风力达10级,最小风力为3级,高频次风力为5-6级,风势走向多为东南风和西北风,且每天中的白昼多为东南风向,时间段在上午10点至下午6点;夜间多为西北风向,时间段为晚上9点至凌晨2点。成都和绵阳地区的夏季最大风力达到7级,秋季最大风力达到9级,高频风力为3-5级。
此次仿真背景模拟5.12汶川大地震灾区救援状况,假设阿坝汶川灾区的风力风向不明,已知成都、绵阳附近参考点J1、J2、J3的风力风向,预测汶川山地低空风A至E点。已知参考点J1至A点距离130km,参考点J3至E点距离120km。通过坐标位置可以确定每个待测点与参考点之间的距离。图10中的红色字母表示风矢量待测点。表4给出风矢量参考点的坐标,表5给出风矢量待测点的坐标。
表4风矢量参考点信息
待测点与J1的水平距离为R1,高度差为H1;待测点与J2的水平距离为R2,高度差为H2;待测点与J3的水平距离为R3,高度差为H3。
表5风矢量待测点信息(单位:m)
参考点J2的风矢量可以分解为东西风向的风矢量U2和南北方向的风矢量V2,大小均为7.71m/s。运用本文的风场内插方法计算风矢量过程如下:
缓坡上的风矢量计算以A点为例,首先计算山脚处风矢量:
最终得到的风矢量数值正负号与参考点风分量正负号保持一致,原因是保持风分量方向不变,故UA=8.49m/s,VA=-7.30m/s。
陡坡上的风矢量计算以D点为例,首先计算山脚处风矢量:
最终得到的风矢量数值正负号与参考点风分量正负号保持一致,原因是保持风分量方向不变,故UD=8.1m/s,VD=-7.1m/s。
其余待测点的计算过程不再赘述,计算结果由表6和表7给出,其中表6是使用文献[15]的方法得出的结果,表7是使用本文方法计算得出的结果。
表6传统风场内插方法计算结果
表7基于坡度的内插方法计算结果
风分量需要进一步合成为风矢量,合成方法依据航行速度三角形原理,最终得到每个点的风速和风向。此次选取点的风速实测值为气象记录站对该点十分钟内的风速取平均值,风向实测值为该点风向角度在十分钟内的平均值,表8中的风向以角度表示,风矢量指向正北方向的角度为0°,指向正西方向的角度为90°,并且逆时针依次递增。风速1和风向1表示使用文献[15]的方法计算得到的结果,风速2和风向2表示使用本文方法得到的计算结果。
表9中的误差由该式算出:
其中,x表示计算值,μ表示实测值。
表8两种方法计算结果对比
表9两种方法计算误差对比
通过分析比较发现,运用本文方法对风场进行内插计算得到的结果,整体上进一步贴近实测值,尤其对于缓坡的部分的风矢量预测更精确,但是可以发现D点的风速误差较大,这是由于该点处于峡谷地带,气流变化大,因此预测难度也较大。从风向方面看,两种方法对于山谷中的风向预测精度相差不大,其中D点的风向误差较大,这说明现实中风向的变化受地形影响非常大,狭窄的山谷、陡峭的地势等都会对风向产生很大影响,因此更难预测。总的来说,本文方法相对于文献[15]方法,对于风速的预测误差最小降低0.34%,最大降低了13.37%,平均降低2.74%;对于风向的预测误差最小降低了0.22%,最大降低了2.38%,平均降低0.52%。
近地风速计算完毕后,需要计算低空风速,而运用指数律公式需要确定切变系数,即下垫面情况,此次仿真救援环境的地貌多为林地,故取切变系数m=0.3。首先需要将三个参考点的风分量插值到每个待测点所在的高度层:
以参考点J1插值到A点所在高度层为例:
J1的风分量只有南北方向V1,J1所在高度ZJ1=2255m,A所在高度为ZA=2618m,则坐标A(113,71,2618)的风分量为V1×(ZA/ZJ1)m=6.7×(2618/2255)0.3=7.01m/s。
将低空风分量合成为最终风矢量,即可得到各个航路节点上的风速和风向。
3.3低空救援安全航迹规划仿真
由于低空空域中没有设置航线,故此次仿真实验中根据航空器性能和目视飞行安全间隔设置空间立体三维网格,每个网格设置中心点作为航路节点,航路节点作为算法求解航行路径的预选点,后续可根据山体地形、低空风速等因素对航路节点位置进行调整。网格的水平长和宽为35000m,高度440m,地图高度与现实高度比例一致。三维网格效果如图11所示。
图11中的点为网格中心点,边框为网格边缘,三维网格建立后可覆盖整个救援区域。现假定有救援任务,派遣直升机从绵阳附近出发,前往灾区附近搜救伤员(在后续的路径求解过程中,直升机到达灾区上空并保持悬停状态表示在灾区附近搜救的过程),随后将伤员送往成都,直升机型号为直9,巡航速度为250km/h,爬升率为7.7m/s,且救援路程不超过其最大航程,求直升机救援航迹。
判断直升机运行安全的两个衡量标准可分为飞行超障高度和侧向安全裕度。查找民航局飞行标准司发布AC-91-FS-2014-22的咨询通告《直升机安全运行指南》发现,对于救援作业中的超越障碍高度和侧向安全裕度,文件没有明确给出,故参考民航总局批准巡线作业的最低要求:巡线飞行真高不低于50m,侧向距离不小于旋翼直径的1.5倍。直9旋翼直径为11.93m,根据文件要求,侧向距离不小于18m。
运用Dijkstra算法求解路径,得到的航路节点坐标依次为A-B-C-D-E-F-G-H-I-J-K-L。航程为88.98km,所用时间为21.3min,航迹见图12。
该航迹为无风状态下的安全航迹,直升机始终保持足够的安全距离和超障距离飞行,但这种情况为理想状况,接下来运用低空风预测模型计算航路节点上的风速和风向,计算结果如表10所示。
表10航路节点的风速和风向
将航路节点上的风速与航空器在该节点的空速合成,并按照无风情形下的计划航线飞行,最终得到低空风影响下的航迹,见图13。
图13中的红色实现为无风情形下的航迹,蓝色虚线为低空风影响下的航迹,蓝色航迹的航路节点依次为A'-B'-C'-D'-E'-F'-G'-H'-I'-J'-K'-L'。表11给出低空风影响下的航路节点,每个航路节点与无风情形下的航路节点相比较的偏移量也在表中给出。
表11低空风影响下航路节点的偏移量
蓝色航迹有两个航段存在严重冲突隐患,第一个冲突隐患位于节点B'至节点D'之间的航段上,该航段的冲突位置如图14(a)所示。由于附近存在较强烈的水平风场,使得直升机在水平方向上的位移量增加,导致在爬升过程中,过早接近山顶,且通过测算发现,最近处的垂直距离不足42m,一旦驾驶员操作失误,或出现下击气流,则很可能导致航空器触地发生意外。
第二个冲突隐患位于节点H'至节点I'之间的航段上,该航段的冲突位置如图14(b)所示。按照原计划直升机选择绕开高海拔山峰飞行,但通过低空风预测模型发现,附近存在较强烈的吹往西北方向的风场,且随着海拔的上升,风速越来越大,使得直升机距离原规划航迹的偏离量逐步增加,使得航迹与山峰外延最短垂直距离不足30m,存在极大的危险。
针对两个严重潜在威胁点,需要在相应的航段位置调整飞行策略。为了调整后的航迹不再发生类似冲突,需要计算低空风对航迹造成的偏移量。在此次飞行过程中,有10个航路节点发生偏移,其中有3个航路节点产生的偏移量对安全隐患产生的贡献较大,依据航空器性能、航路节点的位置和节点上的风矢量信息,可以求得相应的偏移量,见表12。
表12关键航路节点的偏移量
除了这三个航路节点坐标,其余的航路节点由于受到侧风影响,最终到达的航路终点位置与原计划有所偏离,这些偏离量一并放入新的航迹规划中修正,根据偏离量和航空器在该航路节点空速,可以计算出航空器的飞行角度,最终求得航空器到达的位置点,得到的安全航迹规划如图15所示,航迹节点的具体坐标见表13。
表13低空风影响下的安全航迹节点及偏移量
图15中的黑色虚线航迹为修正后的规划航迹,航路节点依次为A”-B”-C”-D”-E”-F”-G”-H”-I”-J”-K”-L”,航程比原计划增加8233m,预计时间增加127s。修正过程中,对于不影响航空器安全的航路节点予以保留,对于影响航空器安全、或者可能造成与终点距离偏差较大的航路节点,予以修正,沿着预先规划的安全航迹飞行,可确保航空器在近地飞行过程中,即使遭受低空风的影响,仍然留有足够的安全距离,同时,所增加的航行路程和时间成本分别为9.3%和9.9%,增加幅度不大。
3.4航迹规划对比分析
为评估本章航迹规划方法的有效性,与Bazhenov S G[23]所提的方法进行对比分析。仿真环境依旧使用上述的汶川附近地形环境,救援行动使用的航空器依旧为直9型号直升机,通过本文方法与文献[23]方法求解得到的低空救援航迹如图16所示,具体坐标值见表14。
表14文献方法求解得到的低空救援航迹
图16中的黑色虚线航迹为本文方法得到的低空救援航迹规划结果,其航路节点坐标和其他相关信息不再赘述。洋红色实线为文献[23]方法得到的航迹规划结果,其航路节点坐标依次为M-N-O-P-R-S-T-U-T'-V-W-X-Y-Z-M'-N',航程为110.82km,所用时间为26.3min。由于该文献的路径求解算法的本质是基于二维平面的路径寻优,故航迹在经过灾区附近时,其航路节点的顺序为S-T-U-T'-V。
图17为对比文献的航迹规划结果存在的两处冲突隐患点。图17(a)为第一个冲突隐患点,航迹与山体最近处的垂直距离仅余15m,图17(b)为第二个冲突隐患点,航迹与山体最近处的垂直距离仅余23m,两处冲突隐患点均存在极大的航空器撞地威胁,且该情况为未引入低空风场的理想情形。显然,若引入低空风场,该航迹结果的风险只增不减。
表15为两种航迹规划方法的性能比较。其中,未修正航迹表示航迹规划结果在低空风场影响下的航迹偏移情况,而安全航迹表示航迹规划结果已经考虑低空风场对航线轨迹的影响并提前做出相应修正。水平距离表示在航迹上,航空器与地形障碍的水平距离最近处,该距离若小于18m则视为航空器存在安全隐患;垂直距离表示在航迹上,航空器与地形障碍的垂直距离最近处,该距离若小于50m则视为航空器存在安全隐患;水平距离和垂直距离若不满足安全距离条件,则显示具体的数值。本文方法未修正航迹与文献[23]方法的航迹均存在安全隐患。文献[23]的航迹规划方法主要由9个步骤的算法实现,核心的最优路径求解算法同样为Dijkstra算法,但由于路径求解过程局限于航空器当前所处高度层,且对于垂直越障高度的距离判定存在滞后性,故造成规划结果存在安全风险。
表15两种航迹规划方法性能比较
与文献[23]所提方法不同,本文方法的核心算法虽然也基于Dijkstra算法,但是本文规划方法在处理联通矩阵过程中,将可通行的三维网格航路节点提前筛选出来,形成坐标矩阵并在该矩阵中求解最优路径,克服了二维路径求解算法运用在三维坐标矩阵中的困难,不仅充分利用了可行解空间,且避免了路径求解结果出现地形碰撞的可能,保证了航迹安全。从结果上看也印证了这一点,与文献[23]方法相比,本文方法的安全航迹在航程上缩短了12.3%,在耗时上减少了11%,同时还保证了航空器在救援过程的安全,体现了本文所提规划方法的安全性和优越性。
4结论
本发明提出了一种基于坡度山风规律分析的风场内插方法,该方法可以将地形高程、地形坡度、下垫面粗糙度和距离等因素对风矢量的影响综合考虑得出风场内插值,相比传统的反距离风场内插方法,本章提出的方法对于缓坡地形的风场数值预测,其准确度有明显提高;提出了一种低空风影响下的安全航迹规划方法,该方法克服以往将低空风场和复杂地势相互割裂的弊端,将两者结合考虑进行航迹的安全规划,分析不同风向的风场对航行轨迹的影响,并根据航路节点的偏移状况,提出相应的修正策略,从而确保航迹安全。
本发明提出的风场内插方法对于缓坡地形的风矢量预测,其准确度有明显提高,而处于更复杂的狭窄山谷地形的风场,以及其他陡峭地势的风场,预测精度还有待提高。在低空航迹规划方面,后续可以加入低空风切变等其他气象因素的考虑,或者进一步分析低空风对航空器飞行状态的影响,基于空气动力与航空器性能,构建更全面、更精确的低空救援航迹规划方法。
Claims (6)
1.一种基于低空风预测模型的安全航迹规划方法,其特征在于包括如下步骤:
步骤1,建立救援环境的三维网格,将网格中心点作为救援航空器的飞行路径节点,以航空器的垂直安全距离和侧向安全距离为准,对飞行路径节点进行筛选,得到可通行矩阵A;
步骤2,引入距离权重,得到带距离权重的邻接矩阵C;
步骤3,采用Dijkstra算法对带距离权重的邻接矩阵C求解最短路径,途经的节点即为最短航迹节点;
步骤4,利用基于地形坡度的低空风场预测方法得到每个最短航迹节点对应的风速和风向,结合航空器的航速和航向,从而得到航空器在各航迹节点的实际航迹;
步骤5,判断步骤4得到的实际航迹是否满足航空器飞行的垂直安全距离和侧向安全距离,若满足则不作处理,若不满足则进行航迹修正,最终得到低空救援安全航迹;
在步骤4中,利用基于地形坡度的低空风场预测方法得到每个最短航迹节点对应的风速和风向的过程是:
步骤41,确定风矢量参考点j及待求点i,确定地形基准面及其高程;
步骤42,计算每个待求点所在坡体的坡度,坡度由待求点i所在高度Hi以及该点与基准点在主要风向上的水平距离Li共同确定,其中基准点指该点所在山坡的坡脚外延轮廓与主风向延长线的交点;
步骤43,根据待求点所在坡体坡度,采用不同方法计算风矢量:
若待求点所在坡体坡度小于0.2,则视为平原地区考虑,利用反距离权重插值法计算风速;具体如下:
若待求点所在坡体坡度大于0.2且小于0.7,则视为缓坡地区考虑,若坡度大于0.7,则视为陡坡地区考虑;计算缓坡或陡坡地区的风矢量时,先用反距离权重插值法计算坡体山脚处风速,再利用相应的拟合函数计算坡体上的风速;
风向与参考点风向一致;
第一种拟合函数为:
f(x)=a1×sin(b1×x×c1)+a2×sin(b2×x×c2)+a3×sin(b3×x×c3) (3)
其中a1=47.53,b1=0.007575,c1=-1.14,a2=49.3,b2=0.01263,c2=0.1726,a3=22.4,b3=0.01541,c3=2.298,该拟合函数的决定系数R2=0.9849,决定系数越接近1,则拟合效果越好,和方差SSE=0.2809,均方根误差RMSE=0.3748;
第二种拟合函数为:
f(x)=d1×sin(e1×x×f1)+d2×sin(e2×x×f2)+d3×sin(e3×x×f3) (4)
其中d1=139.8,e1=0.00383,f1=1.448,d2=125.7,e2=0.006422,f2=3.735,d3=28.76,e3=0.009827,f3=5.668,该拟合函数的决定系数R2=0.9975,决定系数越接近1,则拟合效果越好,和方差SSE=0.1287,均方根误差RMSE=0.2537;
步骤44,基于步骤43求得的近地风速,根据航迹节点的海拔高度,采用指数律公式或对数律公式估算航迹节点所在高度的低空风速。
2.如权利要求1所述的基于低空风预测模型的安全航迹规划方法,其特征在于:所述步骤1中,可通行矩阵A包含满足判断依据的网格节点,所述判断依据的内容是:每个节点到该节点坐标所在地面的距离大于航空器的垂直安全距离,且每个节点到所在高度层最近起伏地形的直线距离大于航空器的侧向安全距离。
3.如权利要求1所述的基于低空风预测模型的安全航迹规划方法,其特征在于:所述步骤2的具体过程是:
步骤21,根据可通行矩阵A中节点的坐标点求得两个节点之间的距离,将此作为带距离权重的矩阵B的矩阵元素;
步骤22,判断矩阵B中相邻节点是否联通,若两个节点不能够直接联通,设定距离为无穷大;若两个节点能够直接联通,则保留原距离,从而得到带距离权重的邻接矩阵C。
4.如权利要求3所述的基于低空风预测模型的安全航迹规划方法,其特征在于:所述步骤22中,判断矩阵B中相邻节点是否联通的方法是:判断如下3个条件,相邻节点之间的x坐标之差不大于三维网格的长,相邻节点之间的y坐标之差不大于三维网格的宽,相邻节点之间的z坐标之差不大于三维网格的高;若3个条件同时满足则判断为相邻节点联通,否则判断为不联通。
5.如权利要求1所述的基于低空风预测模型的安全航迹规划方法,其特征在于:所述步骤3中,采用Dijkstra算法对带距离权重的邻接矩阵C求解最短路径的过程是:
步骤31,初始化航路的起点s0和终点t,设已求出最短路径的节点集合为S,其余未确定最短路径的节点集合为U;
步骤32,从U中选取一个距离起点s0最小的节点sk,把sk加入S中;
步骤33,以sk为新考虑的中间点,若从起点s0到终点t、且经过sk的距离比不经过sk的距离短,那么将最短距离值更新为从起点s0到终点t、且经过sk的距离:
步骤34,重复步骤32-33,直至所有节点都包含在S中,此时的最短距离值即为起点s0到终点t的最短距离,途径的航路节点即为最短航迹节点。
6.如权利要求1所述的基于低空风预测模型的安全航迹规划方法,其特征在于:所述步骤5中,进行航迹修正的方法是:若受到竖直风场影响导致航迹节点存在安全隐患,则将原航迹节点的前序航迹节点的航空器航向修正,具体是将航空器航向朝上风向调整θ,θ=θ1-θ6,θ1、θ6分别表示风向和航空器状态,且调整的航向不得超过航空器的最大俯仰角;若受到水平风场影响导致航迹节点存在安全隐患,则将原航迹节点朝着偏移相反方向调整相等的偏移距离。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201910467594.XA CN110243359B (zh) | 2019-05-31 | 2019-05-31 | 基于低空风预测模型的安全航迹规划方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201910467594.XA CN110243359B (zh) | 2019-05-31 | 2019-05-31 | 基于低空风预测模型的安全航迹规划方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN110243359A CN110243359A (zh) | 2019-09-17 |
CN110243359B true CN110243359B (zh) | 2023-03-24 |
Family
ID=67885596
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201910467594.XA Active CN110243359B (zh) | 2019-05-31 | 2019-05-31 | 基于低空风预测模型的安全航迹规划方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN110243359B (zh) |
Families Citing this family (18)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN114521248A (zh) * | 2019-09-30 | 2022-05-20 | 索尼集团公司 | 信息处理设备、信息处理方法和程序 |
CN110794874B (zh) * | 2019-10-11 | 2022-12-30 | 东南大学 | 一种定位误差约束下的飞行器航迹快速规划方法 |
CN110906936B (zh) * | 2019-12-18 | 2022-11-18 | 哈尔滨工程大学 | 一种水下机器人路径规划方法 |
CN111428916B (zh) * | 2020-03-12 | 2023-04-07 | 南京大学 | 一种海上救援船舶巡航路径规划方法 |
CN111580552B (zh) * | 2020-05-09 | 2023-08-04 | 陕西飞机工业(集团)有限公司 | 一种飞机圆航迹自动飞行控制方法 |
CN111895998B (zh) * | 2020-06-17 | 2022-07-15 | 成都飞机工业(集团)有限责任公司 | 一种大型固定翼无人飞行器分段堆栈式航路规划方法 |
CN111854754B (zh) * | 2020-06-19 | 2023-01-24 | 北京三快在线科技有限公司 | 无人机航线规划方法、装置、无人机及存储介质 |
CN112362060B (zh) * | 2020-08-28 | 2022-08-26 | 中国南方航空股份有限公司 | 一种民航飞行航路规划方法 |
CN112506226B (zh) * | 2020-12-24 | 2022-02-01 | 中国人民解放军军事科学院国防科技创新研究院 | 一种基于温度约束条件的长航时无人机飞行航线规划方法 |
CN113435010B (zh) * | 2021-05-26 | 2023-09-12 | 中国再保险(集团)股份有限公司 | 大尺度精细化地形数字化模拟方法和装置 |
CN113624237B (zh) * | 2021-08-10 | 2024-01-02 | 北京理工大学 | 一种基于Dubins路径的航程调整无人机航迹规划方法 |
CN114545931B (zh) * | 2022-01-27 | 2024-05-24 | 大连海事大学 | 一种基于改进人工势场法引导的Bi-RRT算法的水面无人艇路径规划方法 |
CN114625170B (zh) * | 2022-03-24 | 2023-05-12 | 中国民用航空飞行学院 | 一种山区火灾直升机救援飞行路径动态规划方法 |
CN114460972B (zh) * | 2022-04-13 | 2022-06-07 | 中国民航大学 | 一种无人机城区运行管控方法 |
CN115290927A (zh) * | 2022-08-02 | 2022-11-04 | 厦门航空有限公司 | 跑道风分量的修正方法、装置、计算机设备和存储介质 |
CN116481545B (zh) * | 2023-06-21 | 2023-08-25 | 广州中海电信有限公司 | 一种基于卫星通信的船舶通导规划方法、存储介质及系统 |
CN117649785B (zh) * | 2023-11-28 | 2024-06-07 | 中国民航管理干部学院 | 一种无人机多运行人分布式协同冲突化解方法及系统 |
CN117634719B (zh) * | 2024-01-25 | 2024-05-24 | 中国电子科技集团有限公司电子科学研究院 | 一种基于航行概率与空间拓扑约束的航线预测方法及系统 |
Family Cites Families (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN102929285A (zh) * | 2012-11-16 | 2013-02-13 | 中国民用航空飞行学院 | 多救援直升机多目标分配与航迹规划方法 |
CN103577702A (zh) * | 2013-11-13 | 2014-02-12 | 中国航空工业集团公司西安飞机设计研究所 | 一种在低空风切变中飞机临界规避参数的确定方法 |
CN107357999A (zh) * | 2017-07-19 | 2017-11-17 | 云南电网有限责任公司电力科学研究院 | 一种风场的数值模拟方法及系统 |
-
2019
- 2019-05-31 CN CN201910467594.XA patent/CN110243359B/zh active Active
Also Published As
Publication number | Publication date |
---|---|
CN110243359A (zh) | 2019-09-17 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN110243359B (zh) | 基于低空风预测模型的安全航迹规划方法 | |
CN110930770B (zh) | 一种基于管制意图和飞机性能模型的四维航迹预测方法 | |
EP2631732B1 (en) | Method for flying an aircraft along a flight path | |
US11619953B2 (en) | Three dimensional aircraft autonomous navigation under constraints | |
US20170076614A1 (en) | System and method for optimizing an aircraft trajectory | |
CN100535684C (zh) | 用于防止航空器进入涡流发生器尾涡流危险区的方法和系统 | |
Murrieta-Mendoza et al. | Performance DataBase creation using a level D simulator for Cessna Citation X aircraft in cruise regime | |
Schilke et al. | Dynamic route optimization based on adverse weather data | |
Coombes et al. | Landing site reachability in a forced landing of unmanned aircraft in wind | |
Hauf et al. | Adverse weather diversion model DIVMET | |
CN111192481A (zh) | 一种基于碰撞风险的进离场程序无人机管控区边界确定方法 | |
Xue et al. | Intent Modeling and Conflict Probability Calculation for Operations in Upper Class E Airspace | |
Krozel et al. | Analysis of clear-air turbulence avoidance maneuvers | |
CN111142555A (zh) | 一种基于碰撞风险的机场无人机管控区域划设方法 | |
Su et al. | A comprehensive flight plan risk assessment and optimization method considering air and ground risk of UAM | |
Ramée et al. | Aircraft flight plan optimization with dynamic weather and airspace constraints | |
Menon et al. | A modeling environment for assessing aviation safety | |
KR20230078097A (ko) | 도심항공교통 실증 관리를 위한 디지털 트윈 기술에 기반한 3d 가시화 방법 | |
Pleter et al. | A review of flight trajectory optimisations | |
Zhang | Safe path correction method for ambulance aircraft based on low-altitude wind prediction | |
Wells et al. | Minimising emissions from flights through realistic wind fields with varying aircraft weights | |
Thipphavong | Reducing aircraft climb trajectory prediction errors with top-of-climb data | |
Carmona et al. | Fuel savings through missed approach maneuvers based on aircraft reinjection | |
Dimitriou et al. | Preliminary design of a BWB UAV for highway traffic monitoring | |
Visser et al. | Optimization of rotorcraft simultaneous noninterfering noise abatement approach procedures |
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 |