CN115016547A - 基于剖分网格的飞行器航迹规划方法 - Google Patents
基于剖分网格的飞行器航迹规划方法 Download PDFInfo
- Publication number
- CN115016547A CN115016547A CN202210767445.7A CN202210767445A CN115016547A CN 115016547 A CN115016547 A CN 115016547A CN 202210767445 A CN202210767445 A CN 202210767445A CN 115016547 A CN115016547 A CN 115016547A
- Authority
- CN
- China
- Prior art keywords
- subdivision
- aircraft
- radar
- grid
- detection probability
- 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.)
- Pending
Links
- 238000000034 method Methods 0.000 title claims abstract description 51
- 238000001514 detection method Methods 0.000 claims abstract description 126
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 80
- 230000008859 change Effects 0.000 claims abstract description 9
- 238000004364 calculation method Methods 0.000 claims description 46
- 230000008520 organization Effects 0.000 claims description 20
- 150000001875 compounds Chemical class 0.000 claims description 7
- 230000000007 visual effect Effects 0.000 claims description 5
- 230000011218 segmentation Effects 0.000 claims description 4
- 238000006073 displacement reaction Methods 0.000 claims description 3
- 230000035515 penetration Effects 0.000 abstract description 6
- 230000000875 corresponding effect Effects 0.000 description 52
- 235000019580 granularity Nutrition 0.000 description 30
- 238000009826 distribution Methods 0.000 description 23
- 238000004088 simulation Methods 0.000 description 18
- 238000013461 design Methods 0.000 description 9
- 238000010586 diagram Methods 0.000 description 8
- 230000000694 effects Effects 0.000 description 8
- 238000004458 analytical method Methods 0.000 description 6
- 229910052792 caesium Inorganic materials 0.000 description 4
- TVFDJXOCXUVLDH-UHFFFAOYSA-N caesium atom Chemical compound [Cs] TVFDJXOCXUVLDH-UHFFFAOYSA-N 0.000 description 4
- 230000007123 defense Effects 0.000 description 4
- 230000008569 process Effects 0.000 description 4
- 238000007796 conventional method Methods 0.000 description 3
- 238000002474 experimental method Methods 0.000 description 3
- 230000014509 gene expression Effects 0.000 description 3
- 230000006872 improvement Effects 0.000 description 3
- 238000011161 development Methods 0.000 description 2
- 230000018109 developmental process Effects 0.000 description 2
- 238000011156 evaluation Methods 0.000 description 2
- 238000012954 risk control Methods 0.000 description 2
- 230000009471 action Effects 0.000 description 1
- 230000003044 adaptive effect Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000012512 characterization method Methods 0.000 description 1
- 238000006243 chemical reaction Methods 0.000 description 1
- 230000009194 climbing Effects 0.000 description 1
- 239000003086 colorant Substances 0.000 description 1
- 230000000052 comparative effect Effects 0.000 description 1
- 230000002596 correlated effect Effects 0.000 description 1
- 230000009189 diving Effects 0.000 description 1
- 230000007613 environmental effect Effects 0.000 description 1
- 239000004744 fabric Substances 0.000 description 1
- 230000003993 interaction Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 239000002245 particle Substances 0.000 description 1
- 230000035479 physiological effects, processes and functions Effects 0.000 description 1
- 238000009877 rendering Methods 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
- 238000012795 verification Methods 0.000 description 1
- 238000012800 visualization Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G05—CONTROLLING; REGULATING
- G05D—SYSTEMS FOR CONTROLLING OR REGULATING NON-ELECTRIC VARIABLES
- G05D1/00—Control of position, course, altitude or attitude of land, water, air or space vehicles, e.g. using automatic pilots
- G05D1/12—Target-seeking control
Landscapes
- Engineering & Computer Science (AREA)
- Aviation & Aerospace Engineering (AREA)
- Radar, Positioning & Navigation (AREA)
- Remote Sensing (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Automation & Control Theory (AREA)
- Radar Systems Or Details Thereof (AREA)
Abstract
本发明公开了基于剖分网格的飞行器航迹规划方法,包括两大步骤,先基于剖分网格,进行多雷达探测区环境建模,对飞行器的飞行物理空间进行表征;然后在基于剖分网格组织的物理空间中,利用改进的A*算法进行飞行器航迹规划。本发明中采用立体剖分网格理论这一新的底层表征框架以有效解决同一空间不同飞行器带来的空间基准距离变化需求导致的同一空间重复建模的问题,并基于该理论对传统A*算法进行改进,形成了更加适应于飞行器突防背景、更加通用、高效自主的航迹规划方法,辅助飞行员进行航迹决策。
Description
技术领域
本发明涉及航迹规划技术领域,尤其涉及基于剖分网格的飞行器航迹规划方法。
背景技术
随着信息化战争的不断演进,空天战场中针对飞行器的各类威胁飞速发展。这就需要飞行员在面临威胁时不仅要进行复杂的手工操作,还要做出航迹选择等判断决策,但受限于生理和心理等限制,单纯依靠飞行员完成任务变得越来越困难,迫切需要更加通用、高效自主的航迹规划能力辅助决策。军事应用中突防背景下的飞行器航迹规划需求就是在综合考虑算法规划时间、最终航迹长度以及雷达威胁区分布等因素后,在出发前或任务过程中迅速为飞行器规划出一条符合作战意图的有效航迹,并且在该过程中实现风险可控的突防。
飞行器种类繁复多样,特性不一。在传统航迹规划方法中,数据组织结构通常采取基于经纬网或局部空域直角坐标系网格划分的方式,在目标物理空间中依据一定数据组织方法离散取点计算威胁区并支撑后续航迹规划应用。假设在某一个固定的时间基准下,当飞行器的飞行速率不同导致空间基准距离必须改变时,传统方法只能重新对同一环境中的威胁以新的空间尺度计算威胁区,无法将不同速率的飞行器有效地整合在同一个方法和框架中关联计算,不能层次关联不同的空间组织粒度以适应多型飞行器航迹规划的需要。
发明内容
针对上述存在的问题,本发明旨在提供一种基于剖分网格的飞行器航迹规划方法,利用剖分网格解决空域威胁建模问题,并在此基础之上利用网格特点改进A*算法进行航迹规划,从而得到更加通用、高效自主的航迹规划方法辅助飞行员进行航迹决策。
为了实现上述目的,本发明所采用的技术方案如下:
基于剖分网格的飞行器航迹规划方法,其特征在于,包括以下步骤,
S1:基于剖分网格,进行多雷达探测区环境建模,对飞行器的飞行物理空间进行表征;
S2:在步骤S1中基于剖分网格组织的物理空间中,利用改进的A*算法进行飞行器航迹规划。
进一步的,步骤S1中多雷达探测区环境建模的具体操作包括以下步骤,
S101:计算单雷达探测概率;
S102:根据步骤S101中计算出来的单雷达探测概率,计算多雷达探测区内多雷达联合探测概率;
S103:基于剖分网格组织,依据多雷达联合探测概率,计算任一粒度层级雷达探测概率;
S104:利用剖分网格组织下不同粒度层级之间的关联关系,计算新的粒度层级下剖分体元对应的探测概率大小;
S105:依据相应剖分网格体元内探测概率值大小,对基于剖分网格的雷达探测区进行可视化建模。
进一步的,步骤S101的具体操作包括以下步骤,
S1012:将标准飞行器反射截面积与实际目标反射截面积关联计算,得到式中,假设目标i的截面积为标准截面积,目标m的截面积为实际截面积,分别表示i和m对应的单雷达平均信号功率,和分别表示i和m对应的探测概率;
进一步的,步骤S102中,假设待计算物理空间内有n台相同体制的雷达,它们之间相互独立,那么其联合探测概率为在直角坐标系中,假设雷达的坐标为(xi,yi,0),目标的坐标为(x,y,z),则多雷达联合探测概率可表示为
进一步的,步骤S103的具体操作包括以下步骤,
S1031:利用剖分网格对多雷达探测区进行区域网格划分,视标准体元为飞行器当前参数下所需网格粒度相对应的剖分网格赤道标准体元大小;
S1032:利用剖分编码的二进制三维标识数据和空间体元关系计算方法中的体元位移运算方法,分别求取当前位置坐标编码与雷达所在位置编码的三维剖分体元数量差,再分别乘以基础体元的长、宽、高,得到步骤S102中多雷达联合探测概率计算公式中的坐标差值;
S1033:利用步骤S1032中求得的坐标差值,结合步骤S102中的多雷达联合探测概率计算公式,计算得到该粒度层级的雷达探测概率。
进一步的,步骤S104中利用公式计算新的粒度层级下剖分体元对应的探测概率大小,式中,code_L代表当前所处粒度层级为L时剖分体元块所对应的编码标识;code_M代表以M层级对空域环境进行划分时每一个该层级剖分体元对应的编码标识;P代表该剖分编码对应的联合雷达探测概率值;Pcode_L代表飞行器在L层级标准体元活动时至少被发现一次的概率,也即飞行器在L层级标准体元活动时的雷达探测概率;8L-M代表三维八叉树组织下,L层级下的基础剖分体元由多少个M层级剖分体元组成;j代表相应的体元序列;T代表当前L层级基础剖分体元结构下包含的M层级剖分体元中含有探测概率信息的体元个数。
进一步的,步骤S2的具体操作包括以下步骤,
S201:设置飞行器的约束条件,根据需要粒度确定相应网格层级;
S202:确定起点和终点,建立开集和闭集数据表;
S203:将飞行器的起点作为父节点,而后依据分段变步长思想通过剖分编码位变化确定子节点,并分别计算每个子节点的f(n)值后放入开集中;
S204:在开集中找出最小f(n)对应的所有子节点;
S205:若在开集最小f(n)对应的子节点有多个,继续进行子节点前进方向二次选择,找出更优前进方向及其对应的子节点。
进一步的,步骤S203分段变步长方法搜索所有的子节点的具体操作包括以下步骤,
S2031:计算父节点的启发式代价因子h(n’),并结合任务背景设定判断值D;
S2032:当h(n’)值小于D值时,采用小步长高层级网格搜索所有的子节点;
S2033:当h(n’)值大于D值时,采用大步长低层级网格搜索所有的子节点。
进一步的,步骤S203中子节点的f(n)的计算方法为
式中,Pcode_L代表飞行器在L层级标准体元活动时的雷达探测概率,P(n)为相应子节点安全通行的概率。
进一步的,步骤S205中子节点前进方向二次选择的具体操作包括以下步骤,
S2051:在基于空间八叉树组织结构的空间划分中,找到包含目标和父节点的最小八叉树结构,采用公式对八个子块的风险进行评估,式中,q代表包含当前位置和目标位置的最小八叉树网格的八个子块,q=1,2,…,8;Nobstacle、Ntotal代表在最小粒度层级网格体元建模时,包含目标和父节点的最小八叉树结构下,当前层级八个剖分体元块内分别包含的最小粒度体元中威胁体元的总数量和应包含的所有最小粒度体元的总数量;
S2052:在包含目标和父节点的最小八叉树结构中,选择Psparsity(q)值最大的子块,作为可行方向,也即飞行器航迹的目标点。
本发明的有益效果是:
1、本发明中基于剖分网格的飞行器航迹规划方法采用立体剖分网格理论这一新的底层表征框架以有效解决同一空间不同飞行器带来的空间基准距离变化需求导致的同一空间重复建模的问题,并基于该理论对传统A*算法进行改进,形成了更加适应于飞行器突防背景、更加通用、高效自主的航迹规划方法,辅助飞行员进行航迹决策。
2、本发明中基于剖分网格,进行多雷达探测区环境建模,对飞行器的飞行物理空间进行表征,能够形成多级关联,层级可调,重复可用的雷达威胁区表征能力,可将不同速率的飞行器有效地整合在同一个方法和框架中关联计算,有效解决涉及不同飞行器或飞行器不同速率时威胁区重复建模的问题。
3、本发明中将传统的A*算法进行改进,一是在节点搜索效率提升上,提出分段变步长搜索设计,减少节点计算数量;二是对节点f(n)的计算方法进行改进,通过风险可控设计使子节点探测概率与航迹搜索相关联,使整个算法搜索可以考虑穿越威胁较小的区域并能在“风险”和“航迹长短”之间进行平衡;三是在飞行器面临多个可行子节点前进方向时依据稀疏度对可行方向二次选择,进一步减少无效搜索和节点计算数量。
附图说明
图1为本发明在直角坐标系下探测概率大于0.1的雷达威胁区示意图。
图2为本发明基于剖分网格的栅格数量关系示意图。
图3为本发明基于剖分网格的第13层级雷达探测威胁三视图。
图4为本发明基于剖分网格的第13层级和第12层级雷达探测威胁层级关联示意图。
图5为本发明剖分组织框架下曼哈顿距离计算要素。
图6为本发明分段变步长中剖分编码位比较示意。
图7为本发明水平面稀疏度适用情景。
图8为本发明水平面稀疏度计算与移动方向之间的对应关系。
图9为本发明空间中同代价的可行方向。
图10为本发明空间中的稀疏度计算及对应方向
图11为本发明仿真实验雷达联合探测威胁区分布示意图。
图12为本发明仿真实验表6中n=30时传统方法及幂指数为1和3的实验结果。
图13为本发明仿真实验n=30时不同雷达威胁区分布情景下算法实现效果。
图14为本发明仿真实验β=1.5时改进算法应用示例。
图15为本发明仿真实验β=1时改进算法应用示例。
具体实施方式
为了使本领域的普通技术人员能更好的理解本发明的技术方案,下面结合附图和实施例对本发明的技术方案做进一步的描述。
基于剖分网格的飞行器航迹规划方法,包括以下步骤,
S1:基于剖分网格,进行多雷达探测区环境建模,对飞行器的飞行物理空间进行表征;
S2:在步骤S1中基于剖分网格组织的物理空间中,利用改进的A*算法进行飞行器航迹规划。
具体的,步骤S1中多雷达探测区环境建模的具体操作包括以下步骤,
S101:计算单雷达探测概率;
S102:根据步骤S101中计算出来的单雷达探测概率,计算多雷达探测区内多雷达联合探测概率;
S103:基于剖分网格组织,依据多雷达联合探测概率,计算任一粒度层级雷达探测概率;
S104:利用剖分网格组织下不同粒度层级之间的关联关系,计算新的粒度层级下剖分体元对应的探测概率大小;
S105:依据相应剖分网格体元内探测概率值大小,对基于剖分网格的雷达探测区进行可视化建模。
在本发明中,单雷达探测概率计算采用的方法为:
S1012:为形成适应于网格计算的概率计算方法,将标准飞行器反射截面积与实际目标反射截面积关联计算,假设目标i的截面积为标准截面积,目标m的截面积为实际截面积,分别表示i和m对应的单雷达平均信号功率,那么在同一物理空间中,其背景噪声和传播条件相同,并且同一雷达的性能一致时,i和m对应的探测概率可分别表示和将两式ln取反后再相除,可得式中,和分别表示i和m对应的探测概率;
S1013:根据雷达基本方程可得,在同一雷达,同一物理空间中(噪声水平N不变)时,式中,为信噪比,为目标飞行器最大平均反射截面积,R表示目标与雷达之间的欧式距离,K0为与雷达发射功率,天线增益,电磁波传播环境等因素相关的常数;在同一物理空间中,可认为这些值都是不变的。
S1014:当一台固定参数雷达的探测概率设为0.1时,将其对应目标反射截面积和最大探测距离分别用σmc和Rmc表示,则理想情况下,在选定的物理空间中,对位于地表、性能已知的探测雷达,其A值不变,为一确知数。
S1015:将步骤S1014中的代入步骤S1013中的可得单雷达探测概率的计算公式为根据该计算公式可知,物理空间内,飞行器平均反射截面积雷达与目标间的欧式距离Ri和相应位置该雷达对目标的探测概率的关系。
在此基础上,多雷达探测区内多雷达联合探测概率的计算方法为:
依据该多雷达联合探测概率计算公式,假设在128×128×32km3的物理空间内,存在3部随机分布、参数相同的雷达,当雷达探测概率为0.1时,设A值为1,最大探测距离为30km。其在直角坐标系下探测概率大于0.1的雷达威胁区如附图1所示(设L_grid代表基础网格体元边长,此处为1km),其中坐标如(1,2,1)指向坐标左下点,右侧颜色区分图代表探测概率值大小及对应的颜色。
在此基础上,基于剖分网格组织,计算任一粒度层级雷达探测概率的具体操作包括以下步骤,
S1031:利用剖分网格对多雷达探测区进行区域网格划分,根据公式可知,假设A值和值已知,在直角坐标系中求取三维坐标差值类似于A*算法中曼哈顿距离的求解,且依据剖分网格划分后的目标区域并不涉及空间曲率的影响,因此视标准体元为飞行器当前参数下所需网格粒度相对应的剖分网格赤道标准体元大小;
S1032:利用剖分编码的二进制三维标识数据和空间体元关系计算方法中的体元位移运算方法,分别求取当前位置坐标编码与雷达所在位置编码的三维剖分体元数量差,再分别乘以基础体元的长、宽、高,得到步骤S102中多雷达联合探测概率计算公式中的坐标差值,如附图2所示;
在此基础上,当需要在不同层级利用探测概率表征雷达威胁区时,可利用公式计算新的粒度层级下剖分体元对应的探测概率大小,式中,code_L代表当前所处粒度层级为L时剖分体元块所对应的编码标识;code_M代表以M层级对空域环境进行划分时每一个该层级剖分体元对应的编码标识;P代表该剖分编码对应的联合雷达探测概率值;Pcode_L代表飞行器在L层级标准体元活动时至少被发现一次的概率,也即飞行器在L层级标准体元活动时的雷达探测概率;8L-M代表三维八叉树组织下,L层级下的基础剖分体元由多少个M层级剖分体元组成;j代表相应的体元序列;T代表当前L层级基础剖分体元结构下包含的M层级剖分体元中含有探测概率信息的体元个数。
该公式表示方式有两个依据:一是依据剖分网格所形成的“剖分网格划分的实际空间—剖分编码的唯一标识—空间包含信息”的表达结构,在编码关系对应基础上,实际空间及相应的信息一一对应;二是为保证飞行器暴露风险可控,在利用剖分网格表征雷达威胁区时,应标明整个剖分体元内可能被雷达发现的概率大小。这就需要将各个细粒度层级的探测概率取反后连乘再用1来减表征新的粒度层级下剖分体元对应的探测概率大小。
在此基础上,步骤S105对基于剖分网格的雷达探测区进行可视化建模采用剖分网格框架,其在可视化中依托于Cesium开源平台进行二次开发,形成数据交互展示平台,可支持航迹规划中的物理空间和雷达威胁区的剖分计算及表征。Cesium平台是基于JavaScript语言编写的使用webGL渲染的地图引擎,用于3D,2D,2.5D形式的地图展示,经过二次开发,其可以交互展示局部区域剖分网格划分结构和雷达威胁区的可视化表征。软硬件实验环境如下表1所示。
表1运行环境与软硬件配置
基于试验展示目的,选取空域范围及雷达数量等参数如下表2所示。
表2参数设定
同时,因计算机性能限制,如附图3所示,其中,(a)为俯视图,(b)为正视图,(c)为侧视图;选用第13层级剖分网格作为最底层并展示了基于剖分网格框架和上述参数信息下,局部空域环境内雷达联合探测威胁的三视图。其中红色深度代表探测概率值得大小,颜色越深代表探测概率值越大。
根据步骤S104中不同层级探测概率关联方式,在附图4中以第12层级网格作为对比例展示了同一环境划分下的网格层级关联,其中,(a)为第12层级网格划分示意图,(b)为第13层级网格划分示意图。
进一步的,在上述剖分网格的基础上,利用改进的A*算法进行飞行器航迹规划的具体操作包括以下步骤,
S201:设置飞行器的约束条件,根据需要粒度确定相应网格层级;
为在基于剖分网格组织的物理空间中规划飞行航迹,对飞行器的约束条件设置如下:
(1)飞行器速度约束:假设飞行器在任务时间段以速度V匀速前进,同时,在剖分网格框架下某次任务中各型飞行器进入下一邻域网格中的时间步长T是一致的。当以length代表某层级标准体元网格的最小边长,下标如L+1代表第L+1层级时,那么这三者的关系为lengthL≥V×T≥lengthL+1。特别地,当飞行器类型变化或本身飞行速度变化时,只需更新V值以寻找下一适应层级网格即可,为只进行一次环境建模也存在V≥Vmin,找出需要规划的各型飞行器或飞行器在整个任务中的最小速度的要求。
(2)在适宜网格粒度下的飞行器位置表达约束:在不同层级对同一飞行器进行位置表达和探测概率计算时,假设飞行器为一质点且总位于该层级网格的中心位置。
(3)飞行器爬升角/俯冲角及最小直飞步长约定:为适应网格移动特点,假设仿真中的网格大小总满足(1)中的飞行器速度约束情况下,飞行器可以在空域网格中以三维空间26邻域搜索前进。
(4)剖分网格组织的物理空间中曼哈顿距离定义:如附图5所示,在附图5中,(a)为二维水平面曼哈顿距离计算三要素,(b)为三维空域下曼哈顿距离计算新增三要素;以每一层级在赤道地表附近的栅格体元为基础体元,曼哈顿距离计算与所在层级的基础体元大小密切联系。在水平二维平面中以其长宽及对角线长的比例(以长为底)乘以10后向下取整作为曼哈顿栅格间距,三维空域中则增加高度维和不同情况下斜对角线长与高度边的比例关系并乘以10后向下取整作为曼哈顿栅格间距。随后依据公式所示,当在不同层级需要进行粒度转换时,采取相应层级基础体元长之比转换为代价比的关系进行关联关系运算。其中,radioL_M代表第L层级剖分网格的M类型边的比例关系,lengthL代表第L层级规则体元的最小边长的长度。
(5)仿真中网格体元比例关系:如下表3所示,7-18层级的最短边范围为217m-44520m,可以满足不同飞行器飞行距离步长约束,同时在这些层级内的基础网格体元边长关系都约等于1:1:1,故在仿真中将基础体元设定为正方体体元进行算法仿真分析。
表3 GeoSOT-3D剖分框架7-18层级体元理论大小与实际大小的对比
(6)飞行器最大飞行高度约定:以歼击机为例,其雷达反射截面积经验值约为1-2m2,常规飞行高度在9-30km。故假设本发明中所做飞行器最高可飞行高度为30km。
从以上约束条件也表明,本发明中的方法弱化了关于飞行器物理限制对算法搜索过程的影响,其核心在于从新的底层表征组织框架原理着手改进算法,从而适应任务级航迹规划的各项需求及飞行器突防背景特点。
S202:确定起点和终点,建立开集和闭集数据表;该步骤操作与传统A*算法相同,本发明中不做赘述。
S203:将飞行器的起点作为父节点,而后依据分段变步长思想通过剖分编码位变化确定子节点,并分别计算每个子节点的f(n)值后放入开集中;
在任务级航迹规划突防背景中,一般有两个特点:(1)涉及物理空间较大,从起飞点到目标点的距离较远,跨度较广;(2)为保护重要目标,雷达威胁区主要分布在目标点附近且较为集中。基于上述两点,本发明将全局航迹规划分为两个航段,第一个航段为距离目标点较远的巡航段,雷达威胁区较少且分布稀疏;另一个航段为距离目标点较近的突防段,威胁数量较多且分布集中。该划分依据以启发式代价因子h(n)的大小进行判断,结合不同任务背景设定判断值D,进行分段变步长设计。
计算父节点的启发式代价因子h(n’),当h(n’)值小于D值时,认为飞行器进入雷达威胁区较为集中的突防段,采用小步长高层级网格搜索所有的子节点;当h(n’)值大于D值时,采用大步长低层级网格搜索所有的子节点。此处需要说明的是,在传统A*算法中,子节点的最终移动代价值f(n)由两部分组成,f(n)=g(n)+h(n),其中g(n)代表从起点到达子节点所花费的实际移动路径代价值,其计算方法是将从起点到当前点所累加的实际移动路径代价值和当前点到子节点所花的实际移动代价值相加而得到的;h(n)与起点无关,是当前位置点距离目标的代价估计值,通常有曼哈顿距离、欧式距离,以及二者权重结合的计算方法,本申请中对h(n’)的计算方法采用现有的曼哈顿距离、欧式距离,或者二者权重结合的计算方法均可,本发明中不做赘述。
虽然本发明中不考虑曲率影响,但结合实际,为避免曲率影响导致前后空间包含关系变化过大,统一采用相邻层级进行分段变步长操作,如附图6所示,其中,以(x,y,z)分别抽象代表纬度维、经度维和高度维的某层级编码数组。
在附图6中,从左至右重复的(x,y,z)分别代表由低层级到高层级的各层级经度、纬度,高度三个维度上的相应编码。在图中所示的位比较操作中,父节点编码和子节点编码分别代表当前所在网格和空域26邻域搜索下的下一邻域剖分编码。假设整个节点编码共有n位,当飞行器在巡航段采用低层级大步长网格编码位进行比较时只需将编码比较位提高到n-3位编码的比较,其利用编码位变化的操作方式简单易行。
相较于传统搜索方式,在八叉树剖分组织结构下,巡航段中改进算法可以由原来的八个小网格中至少搜索两个小网格变为搜索一个大的八叉树网格,减少一倍以上的的搜索资源消耗,且在接近目标点的过程中,因大步长搜索可以减少在分布情况较简单的雷达威胁区的计算节点数量,更快跳出意义不大的局部困境。在飞行器进入突防段时采用小布长高层级网格,剖分编码比较到第n位,在威胁密集分布的场景中找到更合理的可行航迹。结合剖分网格特点所做出的分段变步长设计,可以使飞行器全局航迹规划方法在提高效率的同时也能保证航迹合理性。
在此基础上,为了保证航迹合理性,将子节点对应层级剖分网格体元内的探测概率值进一步统合成为风险控制系数,使飞行器的航迹选择与子节点所包含的探测风险密切关联。在此基础上,通过作用于该系数的幂指数控制因子形成可自行定义的风险可控设计。依据当前层级子节点探测概率大小,本发明中子节点的f(n)的计算方法为
式中,Pcode_L代表飞行器在L层级标准体元活动时的雷达探测概率,P(n)为相应子节点安全通行的概率。
因为Pcode_L、P(n)是雷达探测概率值,都在[0,1]区间内,作为权重来影响g(n)值时并不能取得较好的效果,所以设计1-ln(Pn)的系数结构,其将安全通过子节点的概率取对数再用1来减得到对实际移动代价g(n)的正反馈影响,同时使作用区间代表区间扩展为[1,+∞)。最后在l-ln(Pn)外增加一个β次幂作为风险控制因子来控制系数大小,供指挥官或飞行员根据任务需要进一步设定“暴露风险”和“航迹长短”中探测威胁对飞行器最终航迹的综合影响程度。其中,β值越大,算法搜寻的路径越偏向于更安全(绕路)的路径。
S204:在剖分网格组织的雷达威胁区中,当飞行器以空域26邻域网格搜索时,子节点在遇到威胁时在开集中通常会出现多个f(n)值相同且为最小值的可行子节点,因此,需先在开集中找出最小f(n)对应的所有子节点;
S205:若在开集最小f(n)对应的子节点有多个,继续进行子节点前进方向二次选择,找出更优前进方向及其对应的子节点。
具体的,S2051:在基于空间八叉树组织结构的空间划分中,找到包含目标和父节点的最小八叉树结构,采用公式对八个子块的风险进行评估,式中,q代表包含当前位置和目标位置的最小八叉树网格的八个子块,q=1,2,…,8;Nobstacle、Ntotal代表在最小粒度层级网格体元建模时,包含目标和父节点的最小八叉树结构下,当前层级八个剖分体元块内分别包含的最小粒度体元中威胁体元的总数量和应包含的所有最小粒度体元的总数量;
假设包含当前位置和目标点的八叉树结构层级为L,那么Ntotal=2L-1;而Nobstacle可以通过编码比对快速查找出来,从而较为简单的计算出稀疏度的大小。通过在多个可行方向中找到潜在风险更小的方向的方式可以有效降低计算节点的数量,避免无效搜索。
S2052:从公式中可以看出,Psparsity(q)代表最小八叉树结构下八个子体元块的稀疏度,其值越大代表以最小粒度网格大小划分的子块中存在威胁的网格数量越少,因此,在包含目标和父节点的最小八叉树结构中,选择Psparsity(q)值最大的子块,作为可行方向,也即飞行器航迹的目标点。
本发明中对水平面和立体空域中需要稀疏度判断的情景进行梳理,水平面稀疏度判断情景如附图7所示。附图7中,(a)、(b)、(c)三个场景分别代表了当前位置与目标点分别在一条竖直线上、在一条水平线上以及在一条斜线上时的可行子节点及相对应的可行方向。其中红色网格代表当前位置点;蓝色网格代表目标位置点;白色网格代表无探测威胁的可通行网格;灰色网格代表有探测威胁存在的网格。以目标点与当前位置均位于北半球且目标点位于当前位置的右方或右上方为例,如附图8所示(未标示障碍情况),其中向下移动方向的稀疏度由当前位置所在的剖分体元块的稀疏度进行确定;向右,向上以及向右上方的稀疏度对应关系与附图8的(a)中稀疏度的计算位置保持一致。
在附图8的(b)所示图例中,只考虑其中四个方向的邻域稀疏度计算原因是:结合附图7中所示的三种情景,在曼哈顿距离计算方式下,当目标位置位于当前位置右方或右上方时,只存在所列举的四种情况使得两个可行子节点的移动总代价值f(n)相等的情况,故稀疏度计算只需要这四个方向的计算结果进行判断。并且在不同竖直平面或上层水平面中出现如附图7中所示各可行方向情景时,类似的将各个Psparsity(n)值对应到相应位置即可,在此不赘述。
在空间26邻域搜索中,当当前位置与目标位置不在同一平面时,还需要考虑的方向情况如附图9所示,其中,(a)为垂直方向位置及f(n)代价相同的可行方向,(b)为斜向位置及f(n)代价相同的可行方向;图中,灰色网格代表存在雷达探测威胁的区域,蓝色网格代表目标位置,除图中所示关系外,其余可以简化为类似二维平面结构的关系。随后如附图10中所示关系可知(附图10中(a)为稀疏度对应位置,(b)为竖直面对应方向,(c)为斜面对应方向),只需要将稀疏度在八叉树结构中的位置关系对应到相应可行方向即可,逻辑结构依旧不变。
以北半球编码对应关系为例,八叉树下稀疏度对应位置及方向对应关系如所示,其中水平面的方向对应关系如附图8中的(b)所示。并且在附图10中空间八叉树展示的稀疏度在竖直面和斜面所对应方向也并不涉及背离目标点的方向判断,这是因为从理论上可以推知,在同一父节点的26个子节点中,向背离目标点的子节点移动只会加大相应的最终移动路径代价值f(n)。而在A*算法基本原理介绍中,可以发现这种情况下各个子节点的最终代价不可能相同,并且依据雷达类球或半球状的威胁区分布特点,可知总存在向目标点移动的可行子节点。
仿真实验:
一般情况下,雷达威胁区分布较复杂的情景下飞行器尝试穿越风险区的意义更大。为有效对比算法效果,该仿真实验采用10km为边长的正方体网格并假设雷达数量在10-30个的情景下,将传统A*算法(以下简称传统算法)、基于剖分网格改进的A*算法(以下简称改进算法)进行对比分析。其中实验主机的运行环境及软件配置如下表4所示。
表4仿真环境与软硬件配置
结合任务背景和算法对比分析的需要,区域条件设置如下表5所示,其中R代表符合条件的雷达最大探测半径,P代表雷达联合探测概率的大小。
表5区域条件设置
基于上述条件,假设雷达在探测概率为0.1时的最大飞行器反射截面积为1m2,最大探测范围为30km,本仿真实验中设同一反射截面积飞行器在同样的时间步长T内穿越网格,不再做出进一步的飞行器性能条件设置。同时,为能对比传统算法和改进算法的综合性能,在仿真场景中设置目标点不处于雷达威胁覆盖之下,从而使避障策略下的传统A*算法也能找到一条可行航迹。当以10km为体元边距时,假设n=20个,雷达威胁区分布情况如附图11所示。
改进A*算法效果与对比分析
为体现本发明中改进算法最终航迹暴露风险的程度,在规划时间和最终航迹的长度以及计算的闭集节点数量这三个评价指标之外,引入风险系数描述改进算法的最终航迹风险程度。假定风险系数以最终航迹中各途经网格节点的雷达联合探测威胁概率大小累加的方式进行计算,其计算方法为式中,d(n)代表飞行器移动到当前位置n所累加的探测概率风险值的大小,代表飞行器在L层级标准体元活动时至少被发现一次的概率,j代表其在最终航迹节点集合中的序列位置。
1、β幂大小变化分析
为在不同威胁区分布情景下对比分析幂β大小对算法的影响,本仿真实验选取β=1,2,3这三个值并以随机分布雷达的数量n为10,20,30个进行试验,随后与传统算法进行对比分析。主要结果如下表6所示。
表6不同环境复杂度下幂变化算法实验结果(L_grid=10km)
需要注明的是,上述数据中不同雷达数量对应的雷达威胁区分布情景相同,不同数量的雷达威胁区分布之间并无关系。在不同雷达威胁区分布情景中,各方法计算时间主要与闭集计算节点的数量规模呈正相关性,与雷达数量之间并不呈必然相关,故雷达数量增加并不一定导致方法运行时间增加,因为其还受到雷达威胁区分布情况的影响。
依据改进算法原理,幂值越小,飞行器选择直接穿越风险区域不再绕路的可能性越大,通常能更快搜索到终点从而减少闭集计算节点,所以规划时间相应减少。但由n=10时改进算法不同幂的各项性能表现可知存在幂指数较小的情况下飞行器并未快速找到目标点,反而计算了更多节点的情况,所以幂值大小对飞行器规划效率影响的推论并非绝对。
由表中n=20的数据对比可知,β=3时改进算法与传统算法的最终航迹代价值f(n)相同,风险系数都为0。同时,虽然改进算法的闭集计算节点数量虽然减少了近一半,但是算法规划时间并未降低一半,这说明了改进算法虽然减少了计算节点,但增加的算法逻辑结构使算法规划效率提升效果不明显。通过n=20,β=1和β=3时的结果对比可以发现幂增大,航迹偏向于更远更安全的路径,幂调节作用生效。并且在β=1和β=2中实验结果一致也可以得到针对不同场景幂值要合理设置才能发挥出作用的结论。
由表中n=10,n=30的数据对比可知,雷达威胁区分布越复杂的情况下,改进算法的规划时间好于传统算法,更能适应雷达威胁区分布复杂的场景。综上所述,可以看出本发明中改进算法在一定暴露风险下能找到一条更短航迹,幂变化调节作用生效。同时,因为本发明中改进算法逻辑结构更加复杂,存在在分布较为简单的雷达威胁区情景中规划时间上略高于传统算法的现象,但在更复杂的场景中,本发明中改进算法能取得更好的效果。
附图12中展示了表6中n=30时传统方法及幂指数为1和3的结果,其中幂指数β=1,β=2的最终路径相同,只是闭集节点数量有较小差异,在此不重复画出。从附图12的(a)中代表闭集节点的各色圆圈数量可知传统方法计算的闭集节点数量多且涉及计算范围广,(b)、(c)分别展示了β=1,β=3改进算法得到的最终航迹和其计算的的闭集节点。从图(b)、图(c)中可以发现两者最终航迹有所区别,且多数闭集节点集中在障碍附近,数量相似,较传统方法所计算的闭集节点分布更加合理。
2、雷达威胁区复杂分布的情景下改进算法适用性
为避免附图12和表6中所展示的结果存在特例,进一步分析验证改进算法在雷达威胁区复杂分布情景中的适用性,本仿真实验中中以β=1、n=30为条件,随机试验了四个不同雷达威胁区分布情景中改进算法的作用效果,实验结果如附图13和下表7所示。在附图13中,(a),(b),(c)和(d)分别对应表7中的情景一、情景二、情景三和情景四;其中黄色线段代表改进算法最终航迹,红色线段代表传统算法最终航迹,各色圆圈代表改进算法闭集计算节点,传统算法的闭集计算节点数量过多在此不做展示。
表7 n=30时不同雷达威胁区分布情景下算法实验结果
由表7和附图13中的结果可知,在雷达威胁区分布比较复杂情况下,相较于传统算法,基于剖分网格改进的A*算法可以更有效地规划穿越威胁区的可行航迹,闭集计算节点少且运行效率高,并且在如情景二所示环境中传统算法不能找到可行航迹。可以发现,基于剖分网格的改进A*算法能较好地适应于飞行器突防背景。
同时,依据表7中的闭集节点计算数量的对比也可以看出,与传统算法相比,在环境较为复杂的情景中,改进算法极大降低了找到目标点所需的计算节点数量,从而降低了算法运行时间,分段变步长设计和子节点前进方向的二次选择有效发挥了作用。并且从理论上讲,分段变步长设计仅在远距离巡航段中降低了一倍的计算资源消耗,在附图13中显示为分布较为稀疏的部分圆圈。而在这部分区间,传统算法也并不会过度消耗计算资源,故可以得到在遇到障碍时发挥作用的子节点前进方向二次选择设计发挥了有效作用这一结论。
综上所述,相较于传统A*算法,在雷达威胁区分布越复杂的情景中,本章结合剖分网格框架特点改进的A*算法综合性能越好。其中,以表7中的数据为例,在n=30的条件下,改进算法闭集计算节点减少了50%-90%,相应运行时间减少了18%-59%。并且该方法能通过幂的选取能在“距离”和“风险”之间进一步调节,能有效规划出一条暴露风向可控的可行航迹,更能适应不同需求下的任务级规划飞行器突防背景。
3、基于Cesium平台的应用示例
将应用场景如表8所示进行设置。因为剖分网格划分层级越高,网格数量规模呈指数级上升,在电脑性能限制下无法进行实验,所以本仿真实验选定第12、13层级下,随机分布雷达数量n为7进行实验。并且因为直接选定13层级作为最小粒度层级,假设飞行器依旧严格满足各项约束条件,在此处不再设置飞行器相关参数。第12及13层级网格标准体元的经距和纬距相同,高距以1m误差相差一倍,且第12层级的高距(最短边)约14km,第13层级的高距(最短边)约为7km,假设在0-100km的剖分空间内不考虑网格曲率对航迹规划方法的影响。
表8应用场景参数设置
因为Cesium平台下的算法应用示例是一次模拟计算、缓存并展示,故不能有效展现改进算法的子节点前进方向二次抉择设计,故本仿真实验仅选取β=1,1.5两个幂值在相同环境下得到的最终航迹进行对比演示。β=1.5时改进算法应用示例如附图14所示,其中,(a)为正视图,(b)为俯视图,(c)为侧视图;β=1时改进算法应用示例如附图15所示,其中,(a)为正视图,(b)为俯视图,(c)为侧视图,(d)为有效穿越风险区域图示。
如附图14所示,改进算法在设置β=1.5时,其得到的最终航迹规避了雷达威胁区,寻找到一条绝对安全的航迹。而在附图15所示的最终航迹中,当设置β=1时,改进算在对威胁区域的风险和航迹距离的综合评判后,选择了一条更加贴近威胁区域的航迹,并且其穿越了部分探测风险较小的雷达威胁区。综合两个图示结果可以发现,“暴露风险”和“航迹长短”的调节因子β得到有效应用,基于剖分网格的改进A*算法能为飞行器规划出一条暴露风险可控的可行航迹。
以上显示和描述了本发明的基本原理、主要特征和本发明的优点。本行业的技术人员应该了解,本发明不受上述实施例的限制,上述实施例和说明书中描述的只是说明本发明的原理,在不脱离本发明精神和范围的前提下,本发明还会有各种变化和改进,这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其等效物界定。
Claims (10)
1.基于剖分网格的飞行器航迹规划方法,其特征在于,包括以下步骤,
S1:基于剖分网格,进行多雷达探测区环境建模,对飞行器的飞行物理空间进行表征;
S2:在步骤S1中基于剖分网格组织的物理空间中,利用改进的A*算法进行飞行器航迹规划。
2.根据权利要求1所述的基于剖分网格的飞行器航迹规划方法,其特征在于,步骤S1中多雷达探测区环境建模的具体操作包括以下步骤,
S101:计算单雷达探测概率;
S102:根据步骤S101中计算出来的单雷达探测概率,计算多雷达探测区内多雷达联合探测概率;
S103:基于剖分网格组织,依据多雷达联合探测概率,计算任一粒度层级雷达探测概率;
S104:利用剖分网格组织下不同粒度层级之间的关联关系,计算新的粒度层级下剖分体元对应的探测概率大小;
S105:依据相应剖分网格体元内探测概率值大小,对基于剖分网格的雷达探测区进行可视化建模。
3.根据权利要求2所述的基于剖分网格的飞行器航迹规划方法,其特征在于,步骤S101的具体操作包括以下步骤,
S1012:将标准飞行器反射截面积与实际目标反射截面积关联计算,得到式中,假设目标i的截面积为标准截面积,目标m的截面积为实际截面积,分别表示i和m对应的单雷达平均信号功率,和分别表示i和m对应的探测概率;
5.根据权利要求4所述的基于剖分网格的飞行器航迹规划方法,其特征在于,步骤S103的具体操作包括以下步骤,
S1031:利用剖分网格对多雷达探测区进行区域网格划分,视标准体元为飞行器当前参数下所需网格粒度相对应的剖分网格赤道标准体元大小;
S1032:利用剖分编码的二进制三维标识数据和空间体元关系计算方法中的体元位移运算方法,分别求取当前位置坐标编码与雷达所在位置编码的三维剖分体元数量差,再分别乘以基础体元的长、宽、高,得到步骤S102中多雷达联合探测概率计算公式中的坐标差值;
S1033:利用步骤S1032中求得的坐标差值,结合步骤S102中的多雷达联合探测概率计算公式,计算得到该粒度层级的雷达探测概率。
6.根据权利要求5所述的基于剖分网格的飞行器航迹规划方法,其特征在于:步骤S104中利用公式计算新的粒度层级下剖分体元对应的探测概率大小,式中,code_L代表当前所处粒度层级为L时剖分体元块所对应的编码标识;code_M代表以M层级对空域环境进行划分时每一个该层级剖分体元对应的编码标识;P代表该剖分编码对应的联合雷达探测概率值;Pcode_L代表飞行器在L层级标准体元活动时至少被发现一次的概率,也即飞行器在L层级标准体元活动时的雷达探测概率;8L-M代表三维八叉树组织下,L层级下的基础剖分体元由多少个M层级剖分体元组成;j代表相应的体元序列;T代表当前L层级基础剖分体元结构下包含的M层级剖分体元中含有探测概率信息的体元个数。
7.根据权利要求6所述的基于剖分网格的飞行器航迹规划方法,其特征在于,步骤S2的具体操作包括以下步骤,
S201:设置飞行器的约束条件,根据需要粒度确定相应网格层级;
S202:确定起点和终点,建立开集和闭集数据表;
S203:将飞行器的起点作为父节点,而后依据分段变步长思想通过剖分编码位变化确定子节点,并分别计算每个子节点的f(n)值后放入开集中;
S204:在开集中找出最小f(n)对应的所有子节点;
S205:若在开集最小f(n)对应的子节点有多个,继续进行子节点前进方向二次选择,找出更优前进方向及其对应的子节点。
8.根据权利要求7所述的基于剖分网格的飞行器航迹规划方法,其特征在于,步骤S203分段变步长方法搜索所有的子节点的具体操作包括以下步骤,
S2031:计算父节点的启发式代价因子h(n’),并结合任务背景设定判断值D;
S2032:当h(n’)值小于D值时,采用小步长高层级网格搜索所有的子节点;
S2033:当h(n’)值大于D值时,采用大步长低层级网格搜索所有的子节点。
10.根据权利要求9所述的基于剖分网格的飞行器航迹规划方法,其特征在于,步骤S205中子节点前进方向二次选择的具体操作包括以下步骤,
S2051:在基于空间八叉树组织结构的空间划分中,找到包含目标和父节点的最小八叉树结构,采用公式对八个子块的风险进行评估,式中,q代表包含当前位置和目标位置的最小八叉树网格的八个子块,q=1,2,…,8;
Nobstacle、Ntotal代表在最小粒度层级网格体元建模时,包含目标和父节点的最小八叉树结构下,当前层级八个剖分体元块内分别包含的最小粒度体元中威胁体元的总数量和应包含的所有最小粒度体元的总数量;
S2052:在包含目标和父节点的最小八叉树结构中,选择Psparsity(q)值最大的子块,作为可行方向,也即飞行器航迹的目标点。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202210767445.7A CN115016547A (zh) | 2022-07-01 | 2022-07-01 | 基于剖分网格的飞行器航迹规划方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202210767445.7A CN115016547A (zh) | 2022-07-01 | 2022-07-01 | 基于剖分网格的飞行器航迹规划方法 |
Publications (1)
Publication Number | Publication Date |
---|---|
CN115016547A true CN115016547A (zh) | 2022-09-06 |
Family
ID=83078132
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202210767445.7A Pending CN115016547A (zh) | 2022-07-01 | 2022-07-01 | 基于剖分网格的飞行器航迹规划方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN115016547A (zh) |
Cited By (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN115426035A (zh) * | 2022-11-04 | 2022-12-02 | 中国人民解放军战略支援部队航天工程大学 | 一种基于剖分网格的定位初值搜索方法及系统 |
-
2022
- 2022-07-01 CN CN202210767445.7A patent/CN115016547A/zh active Pending
Cited By (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN115426035A (zh) * | 2022-11-04 | 2022-12-02 | 中国人民解放军战略支援部队航天工程大学 | 一种基于剖分网格的定位初值搜索方法及系统 |
CN115426035B (zh) * | 2022-11-04 | 2023-03-24 | 中国人民解放军战略支援部队航天工程大学 | 一种基于剖分网格的定位初值搜索方法及系统 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN109631900B (zh) | 一种无人机三维航迹多目标粒子群全局规划方法 | |
CN110031004B (zh) | 基于数字地图的无人机静态和动态路径规划方法 | |
CN109814598B (zh) | 无人机低空公共航路网设计方法 | |
CN102880186B (zh) | 基于稀疏a*算法和遗传算法的航迹规划方法 | |
CN102901500A (zh) | 基于概率a星与智能体混合的飞行器最优路径确定方法 | |
Li et al. | AUV 3D path planning based on A* algorithm | |
US11650969B2 (en) | Quadtree terrain data compression using distance-based pruning | |
CN112414405A (zh) | 一种顾及dsm的应急任务的无人机集群航迹规划方法 | |
CN114812564B (zh) | 基于城市动态时空风险分析的多目标无人机路径规划方法 | |
CN112947541A (zh) | 一种基于深度强化学习的无人机意图航迹预测方法 | |
CN115016547A (zh) | 基于剖分网格的飞行器航迹规划方法 | |
Small et al. | A UAV case study with set‐based design | |
Wu et al. | A non-rigid hierarchical discrete grid structure and its application to UAVs conflict detection and path planning | |
Xie et al. | A novel adaptive parameter strategy differential evolution algorithm and its application in midcourse guidance maneuver decision-making | |
Khattak et al. | Estimating turbulence intensity along the glide path using wind tunnel experiments combined with interpretable tree-based machine learning algorithms | |
CN114529164B (zh) | 一种基于立体剖分网格的动态低空空域态势图构建方法 | |
CN113051633B (zh) | 一种航空器运行可视化方法 | |
CN115031738A (zh) | 一种无人机航路规划方法、装置、设备及介质 | |
Sun et al. | Research on Airspace Grid Modeling Based on GeoSOT Global Subdivision Model | |
Fan et al. | A Path-Planning Method for UAV Swarm under Multiple Environmental Threats | |
He et al. | Real-time stealth corridor path planning for fleets of unmanned aerial vehicles in low-altitude penetration | |
CN117930871B (zh) | 一种旋翼物流无人机群实时冲突解脱方法 | |
Gao et al. | Route Planning Based on Grid-Point Optimization Lazy Theat* Algorithm | |
US20230222927A1 (en) | Quantitative approach and departure risk assessment system | |
CN116007629A (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 |