CN113703483B - 多uav协同轨迹规划方法及系统、设备、存储介质 - Google Patents

多uav协同轨迹规划方法及系统、设备、存储介质 Download PDF

Info

Publication number
CN113703483B
CN113703483B CN202111011675.2A CN202111011675A CN113703483B CN 113703483 B CN113703483 B CN 113703483B CN 202111011675 A CN202111011675 A CN 202111011675A CN 113703483 B CN113703483 B CN 113703483B
Authority
CN
China
Prior art keywords
uav
collaborative
ith
time
track
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
CN202111011675.2A
Other languages
English (en)
Other versions
CN113703483A (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.)
Hunan Cangshu Aerospace Technology Co ltd
Original Assignee
Hunan Cangshu Aerospace Technology Co ltd
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 Hunan Cangshu Aerospace Technology Co ltd filed Critical Hunan Cangshu Aerospace Technology Co ltd
Priority to CN202111011675.2A priority Critical patent/CN113703483B/zh
Publication of CN113703483A publication Critical patent/CN113703483A/zh
Application granted granted Critical
Publication of CN113703483B publication Critical patent/CN113703483B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • GPHYSICS
    • G05CONTROLLING; REGULATING
    • G05DSYSTEMS FOR CONTROLLING OR REGULATING NON-ELECTRIC VARIABLES
    • G05D1/00Control of position, course or altitude of land, water, air, or space vehicles, e.g. automatic pilot
    • G05D1/10Simultaneous control of position or course in three dimensions
    • G05D1/101Simultaneous control of position or course in three dimensions specially adapted for aircraft
    • G05D1/104Simultaneous control of position or course in three dimensions specially adapted for aircraft involving a plurality of aircrafts, e.g. formation flying
    • YGENERAL 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
    • Y02TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02TCLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T10/00Road transport of goods or passengers
    • Y02T10/10Internal combustion engine [ICE] based vehicles
    • Y02T10/40Engine management systems

Abstract

本发明公开了一种多UAV协同轨迹规划方法及系统、设备、存储介质,所述多UAV协同轨迹规划方法先建立协同轨迹优化的约束模型和目标函数,将协同轨迹规划问题描述为参数最优化问题,然后,结合飞行器性能模型,采用B样条曲线参数化表示UAV轨迹,最后采用进化方法对参数最优化问题进行求解,得到满足空间协同要求和时间协同要求的协同轨迹,实现了多UAV轨迹间的时空协同。

Description

多UAV协同轨迹规划方法及系统、设备、存储介质
技术领域
本发明涉及UAV控制技术领域,特别地,涉及一种多UAV协同轨迹规划方法及系统、设备、计算机可读取的存储介质。
背景技术
多架无人飞机(Unmanned aerial vehicle,UAV)协同执行对地打击任务时,为了提高作战效率,实现预定效果,常常需要对多UAV的作战过程进行规划,生成详细的协同作战计划。其中,协同轨迹规划是多UAV协同规划的关键过程。在对地攻击阶段,为了使成功完成任务的概率最大,要求所有UAV分别在同一时刻从相同或不同的起点出发,同时或按照一定的时间顺序到达特定攻击阵位执行任务,这样即使一部分UAV损毁,其它UAV还可以完成任务。协同轨迹规划的目标是得到一组可行、可飞、连接起点和目标点的近似最优轨迹,在规划过程中不仅要综合考虑敌方威胁、地形、气象等环境因素影响,以及平台机动性能的限制,还要考虑多UAV之间的空间和时间协同,是一个复杂的、非线性、带有较强约束的协同轨迹优化问题,具有较大的挑战性。
针对此问题,目前国内外开展了较为广泛的研究,主要分为多机协同规划和轨迹优化两个层面。在协同层规划中主要研究多机之间的时空协同问题,为了描述多UAV间的信息交换,Brigham Young University的McLain和Beard等学者在文章《"CooperativeControl of UAV Rendezvous,"presented at the Proceedings of the AmericanControl Conference,Arlington,2001》、《"Cooperative Path Planning for Timing-Critical Missions,"in Proceedings of the American Control Conference,Denver,Colorado 2003,pp.296-301.》、《"Coordination Variables,Coordination Functions,and Cooperative Timing Missions,"Journal of Guidance,Control,and Dynamics,vol.28,pp.150-161,2005.》中给出了协调变量(Coordination Variables,CV)和协调函数(Coordination Function,CF)的定义,将要交换的信息定义为协调变量,各UAV在规划时只需要保证协调变量的一致性即可保证多UAV协同并将其应用于协同航线规划问题研究中。其中,文章《"Coordination Variables,Coordination Functions,and CooperativeTiming Missions,"Journal of Guidance,Control,and Dynamics,vol.28,pp.150-161,2005.》将协调变量方法应用于具有时间约束的UAV编队轨迹规划问题研究中,以到达时间作为协同航迹规划的协调变量,针对同时到达、严格时序和松散时序三种时间约束关系,分别提出了相应的协同策略。为了充分发挥UAV的自治性,降低对通信的依赖性,提高系统的鲁棒性和容错性,文献《袁利平,陈宗基,周锐,孔繁峨,"多无人机同时到达的分散化控制方法,"航空学报,vol.31,pp.797-805,2010.》建立了基于“协调变量”和“协调函数”的分布式求解框架,提出了一种多UAV同时到达任务区集结问题的分布式协同控制方法。
但是,以上这些文献主要解决航迹层面的协同问题,例如整体协调多个UAV的到达时间,或者协调多UAV到达同一目标,其均没有考虑轨迹层面的协同问题。而在轨迹协同层面,由于每个轨迹点包括位置和时间信息,在时间协同方面,不能简单地协调到达时间,而需要协调每个轨迹点的时间,在空间协同方面,也不能只涉及协同到达同一目标,而需要协调每个轨迹点的位置,防止多UAV在轨迹飞行过程中发生相撞。
发明内容
本发明提供了一种多UAV协同轨迹规划方法及系统、设备、计算机可读取的存储介质,以解决现有技术的上述缺陷。
根据本发明的一个方面,提供一种多UAV协同轨迹规划方法,包括以下内容:
步骤S1:对每个UAV在协同轨迹规划中的约束条件和代价函数进行建模;
步骤S2:采用B样条曲线表示每个UAV的轨迹;
步骤S3:采用进化算法对多UAV协同轨迹规划问题进行求解,得到满足空间协同要求和时间协同要求的协同轨迹。
进一步地,所述步骤S1中的约束条件包括飞机机动性能约束、敌方威胁约束、飞行禁飞区约束、地形约束、终端位置约束、空间协同约束和时间协同约束;
其中,飞机机动性能约束表示为:h(t)为飞行高度,V(t)为真空速,γ(t)为俯仰角,ψ(t)为航向角,μ(t)滚转角;
敌方威胁约束表示为:||·||2表示两点之间的距离,/>和/>分别为第i个威胁的中心坐标及作用半径,(x(t),y(t),h(t))为UAV在t时刻的位置坐标;
飞行禁飞区约束表示为: 和/>分别为第i个飞行禁飞区的中心坐标及半径,NNFZ为飞行禁飞区的数量;
地形约束表示为:h(t)-hij(t)≥Δh,h(t)为飞机的飞行高度,hij(t)为点(i,j)的地形高度,Δh为最小安全高度阈值;
终端位置约束表示为:(xf,yf,hf)表示武器投放点位置,(Δx,Δy,Δh)表示给定的允许偏差,(xAAR,yAAR,hAAR)为武器可投放区的中心点坐标;
空间协同约束表示为:||vpi(k)-vpj(k)||2≥dsafe,i,j=1,2,…,Nv,i≠j,vpi(k)为第i个UAV的k时刻轨迹点,dsafe为UAV间的最小安全间隔距离,Nv为UAV数量;
时间协同约束表示为:Tsi≤Ti≤Tsii,i=1,…,N,Ts为第一个UAV的到达时间,Ti为第i个UAV的到达时间,Δi表示第i个UAV与第一个UAV之间的时间窗,Δ1=0,τi表示第i个UAV的飞行持续时间。
进一步地,所述步骤S1中的代价函数包括UAV的飞行距离代价函数、飞行高度代价函数、威胁代价函数、协同代价函数和综合代价函数;
其中,飞行距离代价函数表示为:PLRi为第i个UAV的飞行距离代价,/>为第i个UAV的第j个坐标点,N为轨迹点数量,lmin为最小飞行距离,用第i个UAV的起点到终点的直线距离表示;
飞行高度代价函数表示为:MFDi为第i个UAV的飞行高度代价,/>为第i个UAV的第j个坐标点,/>为该点的地形高度,hsafe表示飞行的最小安全高度,N为轨迹点数量;
威胁代价函数用雷达探测概率表示: PRDi为组网雷达系统对第i个UAV的探测概率,Pd(r)为第r个雷达对目标的探测概率,Rpr为第i个UAV的轨迹点p到雷达r的距离,Pf表示雷达系统的虚警概率,/>为雷达r的最大探测距离,K表示与雷达威力相关的归整化因子,σ为RCS值,LoS(p,r)为点p到点r的通视性判断方法,当两点间能够通视时为正,否则为负;
协同代价函数表示为:
ACi为第i个UAV与其它UAV的协同代价,/>为第i个UAV的第k个轨迹点,/>为第j个UAV的第l个轨迹点,dsafe为UAV间的最小安全间隔距离,/>为第i个UAV达到第k个轨迹点的时间,/>为第j个UAV达到第l个轨迹点的时间,tmin为两架UAV到达同一个轨迹点的最小安全间隔时间,Ni和Nj分别为第i个UAV及第j个UAV的轨迹点数量;
综合代价函数表示为:Ji为第i个UAV的综合代价,/>为权系数,/>
进一步地,所述步骤S2中采用3次4阶B样条曲线表示每个UAV的轨迹,3次4阶B样条曲线的表达式为:
Zj(u)为UAV的第j段轨迹,bi,4(u)为对于4次输出zj的第i个B样条基函数,u∈(0,1),为B样条曲线的节点,为控制点序列,每个控制点/>表示为(xi,yi,hi),i=0,1,2,3。
进一步地,所述步骤S3包括以下内容:
步骤S31:初始化协进化算法的种群P0,将每个UAV作为一个子种群P0 i,令Nv表示UAV数量,2m表示每个子种群的规模,则初始化后的种群Pn的规模为N=2Nv*m,采用进化算法对B样条曲线的控制点序列进行编码;
步骤S32:对初始化后的种群Pn执行进化操作,产生下一代种群Qn,种群规模为N,合并初始化后的种群Pn和产生的下一代种群Qn,产生新种群Rn,新种群的规模为2N;
步骤S33:对新种群Rn的每个子种群进行并行搜索,并结合3次4阶B样条曲线生成每个UAV的多条轨迹;
步骤S34:基于每个UAV的约束条件并结合综合代价函数对每个UAV的每条轨迹进行评价,得到每条轨迹的综合评价值;
步骤S35:采用同时到达的时间作为协调变量,并基于集中协调-分布式求解的方式进行时间协同,得到协同到达时间满足要求且综合评价值满足要求的3m组轨迹,后执行步骤S36,否则,从新种群Rn中选择2m组综合评价值满足要求的个体再次组成新的种群Pn+1,再返回执行步骤S32;
步骤S36:对每组轨迹的多UAV协同轨迹规划问题进行综合评价,筛选出综合评价值满足要求的2m组轨迹;
步骤S37:将筛选出的2m组轨迹构成新的种群Pn+1,重复执行步骤S32~步骤S37,不断迭代,满足算法终止条件后停止迭代,从最终筛选出的综合评价值满足要求的2m组轨迹中选择综合评价值最小的一组轨迹作为协同轨迹输出。
进一步地,所述步骤S36中采用以下公式对每组轨迹的多UAV协同轨迹规划问题进行综合评价:
J为多UAV协同轨迹规划的综合评价值,d为UAV的个数。
进一步地,所述步骤S35中采用同时到达的时间作为协调变量,第i个UAV的到达时间的取值范围根据UAV的速度范围[Vmin,Vmax]确定:
其中,表示第i个UAV的第q段轨迹长度,Q为第i个UAV的轨迹段数量,t0为同时出发时刻,为了使各个UAV的飞行时间尽量小,取协同到达时间集合中的最小值作为多UAV的协同到达时间:
另外,本发明还提供一种多UAV协同轨迹规划系统,采用如上所述的规划方法,包括:
建模模块,用于对每个UAV在协同轨迹规划中的约束条件和代价函数进行建模;
轨迹表征模块,用于采用B样条曲线表示每个UAV的轨迹;
分析模块,用于采用同时到达的时间作为协调变量并基于集中协调-分布式求解的架构进行多UAV协同轨迹规划问题求解,得到满足空间协同要求和时间协同要求的协同轨迹。
另外,本发明还提供一种设备,包括处理器和存储器,所述存储器中存储有计算机程序,所述处理器通过调用所述存储器中存储的所述计算机程序,用于执行如上所述的方法的步骤。
另外,本发明还提供一种计算机可读取的存储介质,用于存储进行多UAV协同轨迹规划的计算机程序,所述计算机程序在计算机上运行时执行如上所述的方法的步骤。
本发明具有以下效果:
本发明的多UAV协同轨迹规划方法,先建立协同轨迹优化的约束模型和目标函数,将协同轨迹规划问题描述为参数最优化问题,然后,结合飞行器性能模型,采用B样条曲线参数化表示UAV轨迹,最后采用进化方法对参数最优化问题进行求解,得到满足空间协同要求和时间协同要求的协同轨迹,实现了多UAV轨迹间的时空协同。
另外,本发明的多UAV协同轨迹规划系统、设备、计算机可读取的存储介质同样具有上述优点。
除了上面所描述的目的、特征和优点之外,本发明还有其它的目的、特征和优点。下面将参照图,对本发明作进一步详细的说明。
附图说明
构成本申请的一部分的附图用来提供对本发明的进一步理解,本发明的示意性实施例及其说明用于解释本发明,并不构成对本发明的不当限定。在附图中:
图1是本发明优选实施例的多UAV协同轨迹规划方法的流程示意图。
图2是本发明优选实施例中的武器可投放区的示意图。
图3是本发明优选实施例中的3次B样条曲线的示意图。
图4是本发明优选实施例中的多UAV协同打击多目标的轨迹规划求解框架的示意图。
图5是图1中步骤S3的子流程示意图。
图6是本发明优选实施例中的多UAV协同轨迹规划方法在实际场景仿真中两架UAV协同攻击的对地攻击轨迹示意图。
图7是本发明优选实施例中的多UAV协同轨迹规划方法在实际场景仿真中两架UAV的估计到达时间随时间变化的曲线的示意图。
图8是本发明另一实施例的多UAV协同轨迹规划系统的模块结构示意图。
具体实施方式
以下结合附图对本发明的实施例进行详细说明,但是本发明可以由下述所限定和覆盖的多种不同方式实施。
如图1所示,本发明的优选实施例提供一种多UAV协同轨迹规划方法,包括以下内容:
步骤S1:对每个UAV在协同轨迹规划中的约束条件和代价函数进行建模;
步骤S2:采用B样条曲线表示每个UAV的轨迹;
步骤S3:采用进化算法对多UAV协同轨迹规划问题进行求解,得到满足空间协同要求和时间协同要求的协同轨迹。
可以理解,本实施例的多UAV协同轨迹规划方法,先建立协同轨迹优化的约束模型和目标函数,将协同轨迹规划问题描述为参数最优化问题,然后,结合飞行器性能模型,采用B样条曲线参数化表示UAV轨迹,最后采用进化方法对参数最优化问题进行求解,得到满足空间协同要求和时间协同要求的协同轨迹,实现了多UAV轨迹间的时空协同。
可以理解,所述步骤S1中的约束条件包括飞机机动性能约束、敌方威胁约束、飞行禁飞区约束、地形约束、终端位置约束、空间协同约束和时间协同约束。
其中,飞机机动性能影响每一个阶段的任务执行,必须严格满足,否则生成的计划将不可执行,根据飞机的机动性能限制,建立约束式c1如下:
在公式(1)中,h(t)为飞行高度,V(t)为真空速,γ(t)为俯仰角,ψ(t)为航向角,μ(t)滚转角。
敌方威胁的作用范围可近似为半球,故敌方威胁约束表示为约束式c2
其中,||·||2表示两点之间的距离,和/>分别为第i个威胁的中心坐标及作用半径,(x(t),y(t),h(t))为UAV在t时刻的位置坐标。
另外,飞行禁飞区包括气象禁飞区、危险高度禁飞区、未知区域,本发明采用无限长圆柱体模型表示,所规划的轨迹与该圆柱不能相交,故飞行禁飞区约束的约束式c3表示为:
其中,和/>分别为第i个飞行禁飞区的中心坐标及半径,NNFZ为飞行禁飞区的数量。
另外,地形因素对规划的影响主要体现在飞机的生存概率上,地形约束的约束式c4表示为:
h(t)-hij(t)≥Δh (4)
其中,h(t)为飞机的飞行高度,hij(t)为点(i,j)的地形高度,Δh为最小安全高度阈值。
可以理解,为了完成攻击任务,UAV需要到达指定的攻击阵位,即满足终端位置约束。因此,需要计算攻击某一地面目标的武器可投放区(AAR),如图2所示,在投放点处按设定的投放条件投放制导炸弹,其所有可能的落点组成了可达区,通过在水平面内平移这些落点及其对应的弹道至与目标点重合,这时弹道的起点就组成了可投放区。另外,指定投放点获得可达区与指定目标点获得可投放区本质上是一致的。
设(xAAR,yAAR,hAAR)为武器可投放区的中心点坐标,则终端位置约束表示为:
其中,(xf,yf,hf)表示武器投放点位置,(Δx,Δy,Δh)表示给定的允许偏差。
另外,多UAV间的空间协同约束主要是指UAV间的相撞约束,即在飞行过程中,多UAV间要时刻保持一定的安全距离,则空间协同约束表示为:
||vpi(k)-vpj(k)||2≥dsafe,i,j=1,2,…,Nv,i≠j (6)
其中,vpi(k)为第i个UAV的k时刻轨迹点,dsafe为UAV间的最小安全间隔距离,Nv为UAV数量。
另外,时间协同约束包括同时到达时间约束和时序约束,为了完成有效攻击,多UAV必须按照指定时间或指定时序到达攻击阵位。时间协同约束表示为:
Tsi≤Ti≤Tsii,i=1,…,N (7)
其中,Ts为第一个UAV的到达时间,Ti为第i个UAV的到达时间,Δi表示第i个UAV与第一个UAV之间的时间窗,Δ1=0,τi表示第i个UAV的飞行持续时间。对于同时到达时间约束问题,Δi=τi=0;对于时序约束,Δi和τi为正的恒定值,当τi=0,该时序约束退化为严格时序约束(即各UAV的到达时间间隔相同),上式(7)可进一步分解为:
可以理解,在多UAV协同飞行到达指定攻击阵位过程中,需要期望每架UAV能够以最小的代价到达目标点,同时还要使多UAV的整体代价最小。具体地,所述步骤S1中的代价函数包括UAV的飞行距离代价函数、飞行高度代价函数、威胁代价函数、协同代价函数和综合代价函数。
其中,为了减少燃油消耗,降低飞行风险,减少飞机在敌防空区域内的滞留时间,应最小化UAV的飞行轨迹长度。为使目标函数标准化和可接纳,本发明采用实际飞行距离和最小飞行距离的比值描述UAV的飞行距离代价,则飞行距离代价函数表示为:
其中,PLRi为第i个UAV的飞行距离代价,为第i个UAV的第j个坐标点,N为轨迹点数量,lmin为最小飞行距离,用第i个UAV的起点到终点的直线距离表示。
另外,UAV的飞行高度越低,地形遮蔽的效果越好,越有助于规避未知雷达的探测,但飞行高度较低时,UAV的撞地概率会较大。因此,飞行高度代价函数表示为:
其中,MFDi为第i个UAV的飞行高度代价,为第i个UAV的第j个坐标点,为该点的地形高度,hsafe表示飞行的最小安全高度,N为轨迹点数量。
另外,当UAV不被敌防空雷达探测到时,敌防空威胁就无法对UAV造成毁伤,而且雷达探测概率越小,UAV越安全,因此UAV的威胁代价可用雷达探测概率表示。
其中,Pd为雷达r对目标的探测概率,Rpr为第i个UAV的轨迹点p到雷达r的距离,为雷达r的最大探测距离,Pf表示雷达系统的虚警概率,K表示与雷达威力相关的归整化因子,包含了反映雷达目标探测威力的内在因素,σ为RCS值,即雷达r的散射截面积,LoS(p,r)为点p到点r的通视性判断方法,当两点间能够通视时为正,否则为负。
从而,由n个雷达构成的组网雷达系统对第i个UAV的探测概率表示为:
其中,PRDi为组网雷达系统对第i个UAV的探测概率,Pd(r)为第r个雷达对目标的探测概率。
可以理解,多UAV协同任务计划的评价函数除了包含单个UAV自身的代价外,还应包含各UAV间的协同代价。当多UAV间不满足协同约束时,协同任务将难以执行,因此在进行协同轨迹规划时,可将多UAV对协同约束的满足程度通过UAV的协同代价进行描述。本申请主要考虑空间协同问题,即多UAV避碰代价。假设在给定UAVi和UAVj的轨迹情况下,将每条轨迹的逐个轨迹点进行比较,如果轨迹点对的距离小于最小安全距离,检查到达该点的时间,如果到达时间间隔小于最小安全时间间隔,则该轨迹点对相撞。故结合公式(6)~(8),协同代价函数表示为:
其中,ACi为第i个UAV与其它UAV的协同代价,为第i个UAV的第k个轨迹点,/>为第j个UAV的第l个轨迹点,dsafe为UAV间的最小安全间隔距离,/>为第i个UAV达到第k个轨迹点的时间,/>为第j个UAV达到第l个轨迹点的时间,tmin为两架UAV到达同一个轨迹点的最小安全间隔时间,Ni和Nj分别为第i个UAV及第j个UAV的轨迹点数量。
因此,第i个UAV的综合代价函数表示为:
其中,Ji为第i个UAV的综合代价,为权系数,可根据偏好进行设置。
进一步地,多UAV协同轨迹规划问题的综合代价函数则可表示为:
其中,J为多UAV协同轨迹规划的综合评价值,d为UAV的个数。
至此,所述步骤S1中已经将协同轨迹规划问题转化为参数最优化问题,待后续步骤进行参数最优化求解即可得到满足时空协同要求的协同轨迹。
可以理解,B样条曲线是一种参数曲线,一段B样条曲线是由多段Bezier曲线构成的,并且在其交点处保持一定的连续性。曲线相交点被称为断点,是严格递增的实数序列。k阶B样条曲线的表达式为:
其中,Ci为控制点,作为曲线的基函数的系数,Bi,k为第i个k次B样条基函数,可用考克斯—德布尔递推公式(Cox De Boor recursion)给出如下:
式中,Bi,k(u)的下标i表示序号,k表示阶数。该递推公式表明,若要确定第i个k次B样条基函数Bi,k(u),需要用到ui,ui+1,…,ui+k+1共k+2个节点,并称区间[ui,ui+k+1]为Bi,k(u)的支撑区间,在此区间之外,对应的基函数值为0。
可以理解,B样条曲线采用控制点表示,其形状完全由控制点的位置决定,确定了每个控制点就相当于确定了B样条曲线,移动控制点就可灵活地调整曲线的形状。而采用B样条理论求解轨迹优化问题的思路是将轨迹表示成B样条曲线的控制点序列,从而将轨迹优化问题转换成参数优化问题,并采用数值优化方法进行求解。其中,将轨迹优化问题转换成参数优化问题的关键步骤如下:
首先,提出伪时间概念,将时间要素t引入到B样条理论中,可表达为
其中,t0和tf分别是初始时刻和结束时刻,u∈(0,1)作为B样条曲线的节点。则UAV的轨迹可用时间节点分段的形式表示为:
对应节点序列u1
对应节点序列u2
…;
对应节点序列uq
其中,是对于kj次输出zj的第i个B样条基函数,/>是控制点,也是B样条基函数的系数,pj=lj·(kj-mj)+mj表示控制点数量,其中lj是节点插值的数量,mj是在节点处的平滑条件。
为了减少计算量,期望B样条曲线的次数越少越好,但二次曲线是一条抛物线,不能反映曲线的拐点,所以本发明选用3次B样条曲线进行参数化,曲线的阶k=4。从空间n+1个点pi(i=0,1,…,n)中每次取相邻的四个点作为控制点,即可得到一段3次B样条曲线。B样条曲线基函数可表示为:
故第i段3次B样条曲线可表示为:
Bi,4(u)=b1,4(u)Ci-1+b2,4(u)Ci+b3,4(u)Ci+1+b4,4(u)Ci+2 (20)
相应的矩阵表示形式为:
由于本发明研究的是三维空间中的轨迹规划问题,因此将普通的B样条曲线扩展到三维空间,则每个控制点Ci可表示为(xi,yi,hi)。
因此,第j段的UAV轨迹用3次4阶B样条曲线表示为:
其中,Zj(u)为UAV的第j段轨迹,bi,4(u)为对于4次输出zj的第i个B样条基函数,u∈(0,1),为B样条曲线的节点,为控制点序列,每个控制点/>表示为(xi,yi,hi),i=0,1,2,3。其中,3次B样条曲线的示意图如图3所示。
由公式(22)可知,由于3次B样条基函数bi,4(u)已知,只要确定控制点序列便可由公式(22)生成第j段的UAV轨迹,因此,通过3次B样条曲线将UAV的轨迹映射为控制点序列,从而将轨迹规划问题形式化为参数优化问题,便于采用数值优化方法进行求解。
可以理解,针对多UAV协同轨迹规划的特点,本发明采用集中协调-分布求解的思路,引入协调变量(CV)和协调函数(CF)的思想,构建基于协调变量的协同规划框架,多UAV协同打击多目标的轨迹规划求解框架如图4所示。这种分解方法的优点在于:通过协调变量和协调函数把一个高维的优化问题分解成一个低维、计算量小的问题,集中协调单元只需知道协调变量的取值区间和每架飞机对应的协调函数即可,减少信息和数据传输量,提高了协同的效率。
协同轨迹规划主要包括两个方面,单个UAV轨迹的规划和多个UAV轨迹间的时空协同。本发明采用协调变量和B样条理论进行协同轨迹规划的思路,是在单个UAV轨迹规划中,选用B样条方法表示轨迹,使用群体智能优化算法生成多组满足空间协同要求的轨迹,并利用基于协调变量的方法进行时间协同。而在群体智能优化算法中,协进化算法(Coevolutionary Algorithm)是一种基于群体进化的算法,源于生物界中协同进化的思想:两个或多个物种群体通过相互间的作用,促使双方向前进化。该算法通过种群内个体间的合作与竞争来实现对优化问题的求解,各个物种群体都采用进化算法(EvolutionaryAlgorithm,EA)实现进化过程,其原理简单,易于实现,得到广泛应用,因此本发明选用该算法进行问题求解。
具体地,如图5所示,所述步骤S3包括以下内容:
步骤S31:初始化协进化算法的种群,将每个UAV作为一个子种群,并采用进化算法对B样条曲线的控制点序列进行编码;
步骤S32:对初始化后的种群执行进化操作,产生下一代种群,合并初始化后的种群和产生的下一代种群,以产生新种群;
步骤S33:对新种群的每个子种群进行并行搜索,并结合3次4阶B样条曲线生成每个UAV的多条轨迹;
步骤S34:基于每个UAV的约束条件并结合综合代价函数对每个UAV的每条轨迹进行评价,得到每条轨迹的综合评价值;
步骤S35:采用同时到达的时间作为协调变量,并基于集中协调-分布式求解的方式进行时间协同,得到协同到达时间满足要求且综合评价值满足要求的3m组轨迹,后执行步骤S36,否则,从新种群中选择2m组综合评价值满足要求的个体再次组成新的种群,再返回执行步骤S32;
步骤S36:对每组轨迹的多UAV协同轨迹规划问题进行综合评价,筛选出综合评价值满足要求的2m组轨迹;
步骤S37:将筛选出的2m组轨迹构成新的种群,重复执行步骤S32~步骤S37,不断迭代,满足进化算法终止条件后停止迭代,从最终筛选出的综合评价值满足要求的2m组轨迹中选择综合评价值最小的一组轨迹作为协同轨迹输出。
具体地,首先,初始化协进化算法的种群P0,将每个UAV作为一个子种群P0 i,令Nv表示UAV数量,2m表示每个子种群的规模,则初始化后的种群Pn的规模为N=2Nv*m,采用进化算法对B样条曲线的控制点序列进行编码。其中,进化算法的编码原理属于现有技术,故在此不再赘述。
然后,对初始化后的种群Pn执行交叉、变异等进化操作,产生下一代种群Qn,种群规模为N,合并初始化后的种群Pn和产生的下一代种群Qn,产生新种群Rn,则新种群Rn的规模为2N=4Nv×m。
接着,对新种群Rn的每个子种群进行并行搜索,并结合公式(22)生成每个UAV的多条轨迹。
再基于公式(1)至(5),并结合综合代价函数,即基于公式(9)至(14),对每个UAV的每条轨迹进行评价,得到每条轨迹的综合评价值。至此,可以对每个UAV的每条轨迹进行综合评价,综合评价值越小,空间协同越好。
再采用同时到达的时间作为协调变量,并基于集中协调-分布式求解的方式进行时间协同,得到协同到达时间满足要求且综合评价值满足要求的3m组轨迹,即采用协同到达时间较小和综合评级值较低的3m组轨迹,然后执行后续步骤,否则,从新种群Rn中选择2m组综合评价值较小的个体再次组成新的种群Pn+1,再返回执行上述内容。
然后,采用公式(15)对3m组轨迹中的每组轨迹进行综合评价,筛选出综合评价值较小的2m组轨迹,此2m组轨迹是第一次迭代优化得到的空间协同和时间协同较好的多组轨迹。
最后,继续迭代优化,将筛选出的2m组轨迹构成新的种群Pn+1,重复执行上述内容,不断迭代,当满足算法终止条件后停止迭代,例如达到最大迭代次数后停止迭代,或者模型收敛后停止迭代,从最终筛选出的综合评价值较小的2m组轨迹中选择综合评价值最小的一组轨迹作为协同轨迹输出,该组协同轨迹具有最佳的空间协同和时间协同效果。
可以理解,所述步骤S35中,为了描述环境的主要元素,定义χi为UAVi的状态空间,xi∈χi为UAVi的状态变量,定义Ui(xi)为状态xi的可行决策变量集,ui∈Ui(xi)为UAVi的决策变量。则可定义fii×UicR,其中Rc表示协同空间。
在状态xi下,UAVi的可行协调变量集为:
在时间协同问题中,fi是轨迹点与可能到达时间集合的一种简单映射,Θi(xi)并不一定是连续集。但是,对于特定的轨迹和速度选择,协调变量能够取得唯一值θi=fi(xi,ui)∈Θi(xi)。
同时到达时间约束可表述为:
时序约束可表达为:
设fj是可逆的,其逆函数为fi -1i×Θi→Ui,ui=fi -1(xi,θ),决策量ui是状态变量xi和协调变量θi的函数。对于特定状态xi,一个确定的协调变量θi对应唯一的决策量ui。如果所有UAV得到了一致的协调变量,即θ1=…=θn=θ*,则在ui=fi -1(xi*)的控制下,多UAV必能完成协同任务。
因此,本发明仅考虑同时到达情况,将同时到达的时间作为协调变量,第i个UAV的到达时间的取值范围可根据飞行器的速度范围[Vmin,Vmax]确定:
其中,表示第i个UAV的第q段轨迹长度,Q为第i个UAV的轨迹段数量,t0为同时出发时刻,为了使各个UAV的飞行时间尽量小,取协同到达时间集合中的最小值作为多UAV的协同到达时间:
/>
在公式(28)中,若则需要对各UAV的轨迹重新进行计算。
可以理解,如图6和图7所示,本发明的多UAV协同轨迹规划方法在实际场景仿真中实现了两架飞机协同攻击两个静态的地面目标,同时满足飞机动力学约束、状态量和控制量约束,特别是同时到达时间约束。
另外,如图8所示,本发明的另一实施例还提供一种多UAV协同轨迹规划系统,优选采用如上所述的规划方法,所述规划系统包括:
建模模块,用于对每个UAV在协同轨迹规划中的约束条件和代价函数进行建模;
轨迹表征模块,用于采用B样条曲线表示每个UAV的轨迹;
分析模块,用于采用同时到达的时间作为协调变量并基于集中协调-分布式求解的架构进行多UAV协同轨迹规划问题求解,得到满足空间协同要求和时间协同要求的协同轨迹。
可以理解,本实施例的多UAV协同轨迹规划系统,先建立协同轨迹优化的约束模型和目标函数,将协同轨迹规划问题描述为参数最优化问题,然后,结合飞行器性能模型,采用B样条曲线参数化表示UAV轨迹,最后采用进化方法对参数最优化问题进行求解,得到满足空间协同要求和时间协同要求的协同轨迹,实现了多UAV轨迹间的时空协同。
可以理解,本实施例的系统中的各个模块分别与上述方法实施例中的各个步骤相对应,故每个模块的工作原理在此不再赘述,参考上述方法实施例即可。
另外,本发明的另一实施例还提供一种设备,包括处理器和存储器,所述存储器中存储有计算机程序,所述处理器通过调用所述存储器中存储的所述计算机程序,用于执行如上所述的方法的步骤。
另外,本发明的另一实施例还提供一种计算机可读取的存储介质,用于存储进行多UAV协同轨迹规划的计算机程序,所述计算机程序在计算机上运行时执行如上所述的方法的步骤。
一般计算机可读取存储介质的形式包括:软盘(floppy disk)、可挠性盘片(flexible disk)、硬盘、磁带、任何其与的磁性介质、CD-ROM、任何其余的光学介质、打孔卡片(punch cards)、纸带(paper tape)、任何其余的带有洞的图案的物理介质、随机存取存储器(RAM)、可编程只读存储器(PROM)、可抹除可编程只读存储器(EPROM)、快闪可抹除可编程只读存储器(FLASH-EPROM)、其余任何存储器芯片或卡匣、或任何其余可让计算机读取的介质。指令可进一步被一传输介质所传送或接收。传输介质这一术语可包含任何有形或无形的介质,其可用来存储、编码或承载用来给机器执行的指令,并且包含数字或模拟通信信号或其与促进上述指令的通信的无形介质。传输介质包含同轴电缆、铜线以及光纤,其包含了用来传输一计算机数据信号的总线的导线。
以上所述仅为本发明的优选实施例而已,并不用于限制本发明,对于本领域的技术人员来说,本发明可以有各种更改和变化。凡在本发明的精神和原则之内,所作的任何修改、等同替换、改进等,均应包含在本发明的保护范围之内。

Claims (9)

1.一种多UAV协同轨迹规划方法,其特征在于,包括以下内容:
步骤S1:对每个UAV在协同轨迹规划中的约束条件和代价函数进行建模;
步骤S2:采用B样条曲线表示每个UAV的轨迹;
步骤S3:采用进化算法对多UAV协同轨迹规划问题进行求解,得到满足空间协同要求和时间协同要求的协同轨迹;
所述步骤S1中的约束条件包括飞机机动性能约束、敌方威胁约束、飞行禁飞区约束、地形约束、终端位置约束、空间协同约束和时间协同约束;
其中,飞机机动性能约束表示为:h(t)为飞行高度,V(t)为真空速,γ(t)为俯仰角,ψ(t)为航向角,μ(t)滚转角;
敌方威胁约束表示为:||·||2表示两点之间的距离,和/>分别为第i个威胁的中心坐标及作用半径,(x(t),y(t),h(t))为UAV在t时刻的位置坐标;
飞行禁飞区约束表示为: 和/>分别为第i个飞行禁飞区的中心坐标及半径,NNFZ为飞行禁飞区的数量;
地形约束表示为:h(t)-hij(t)≥△h,h(t)为飞机的飞行高度,hij(t)为点(i,j)的地形高度,△h为最小安全高度阈值;
终端位置约束表示为:(xf,yf,hf)表示武器投放点位置,(△x,△y,△h)表示给定的允许偏差,(xAAR,yAAR,hAAR)为武器可投放区的中心点坐标;
空间协同约束表示为:||vpi(k)-vpj(k)||2≥dsafe,i,j=1,2,L,Nv,i≠j,vpi(k)为第i个UAV的k时刻轨迹点,dsafe为UAV间的最小安全间隔距离,Nv为UAV数量;
时间协同约束表示为:Ts+△i≤Ti≤Ts+△ii,i=1,L,N,Ts为第一个UAV的到达时间,Ti为第i个UAV的到达时间,△i表示第i个UAV与第一个UAV之间的时间窗,△1=0,τi表示第i个UAV的飞行持续时间。
2.如权利要求1所述的多UAV协同轨迹规划方法,其特征在于,所述步骤S1中的代价函数包括UAV的飞行距离代价函数、飞行高度代价函数、威胁代价函数、协同代价函数和综合代价函数;
其中,飞行距离代价函数表示为:PLRi为第i个UAV的飞行距离代价,/>为第i个UAV的第j个坐标点,N为轨迹点数量,lmin为最小飞行距离,用第i个UAV的起点到终点的直线距离表示;
飞行高度代价函数表示为:MFDi为第i个UAV的飞行高度代价,/>为第i个UAV的第j个坐标点,/>为该点的地形高度,hsafe表示飞行的最小安全高度,N为轨迹点数量;
威胁代价函数用雷达探测概率表示: PRDi为组网雷达系统对第i个UAV的探测概率,Pd(r)为第r个雷达对目标的探测概率,Rpr为第i个UAV的轨迹点p到雷达r的距离,Pf表示雷达系统的虚警概率,/>为雷达r的最大探测距离,K表示与雷达威力相关的归整化因子,σ为RCS值,LoS(p,r)为点p到点r的通视性判断方法,当两点间能够通视时为正,否则为负;
协同代价函数表示为:
ACi为第i个UAV与其它UAV的协同代价,/>为第i个UAV的第k个轨迹点,/>为第j个UAV的第l个轨迹点,dsafe为UAV间的最小安全间隔距离,/>为第i个UAV达到第k个轨迹点的时间,/>为第j个UAV达到第l个轨迹点的时间,tmin为两架UAV到达同一个轨迹点的最小安全间隔时间,Ni和Nj分别为第i个UAV及第j个UAV的轨迹点数量;
综合代价函数表示为:Ji为第i个UAV的综合代价,/>为权系数,/>
3.如权利要求2所述的多UAV协同轨迹规划方法,其特征在于,所述步骤S2中采用3次4阶B样条曲线表示每个UAV的轨迹,3次4阶B样条曲线的表达式为:
Zj(u)为UAV的第j段轨迹,bi,4(u)为对于4次输出zj的第i个B样条基函数,u∈(0,1),为B样条曲线的节点,为控制点序列,每个控制点/>表示为(xi,yi,hi),i=0,1,2,3。
4.如权利要求3所述的多UAV协同轨迹规划方法,其特征在于,所述步骤S3包括以下内容:
步骤S31:初始化协进化算法的种群P0,将每个UAV作为一个子种群P0 i,令Nv表示UAV数量,2m表示每个子种群的规模,则初始化后的种群Pn的规模为N=2Nv*m,采用进化算法对B样条曲线的控制点序列进行编码;
步骤S32:对初始化后的种群Pn执行进化操作,产生下一代种群Qn,种群规模为N,合并初始化后的种群Pn和产生的下一代种群Qn,产生新种群Rn,新种群的规模为2N;
步骤S33:对新种群Rn的每个子种群进行并行搜索,并结合3次4阶B样条曲线生成每个UAV的多条轨迹;
步骤S34:基于每个UAV的约束条件并结合综合代价函数对每个UAV的每条轨迹进行评价,得到每条轨迹的综合评价值;
步骤S35:采用同时到达的时间作为协调变量,并基于集中协调-分布式求解的方式进行时间协同,得到协同到达时间满足要求且综合评价值满足要求的3m组轨迹,后执行步骤S36,否则,从新种群Rn中选择2m组综合评价值满足要求的个体再次组成新的种群Pn+1,再返回执行步骤S32;
步骤S36:对每组轨迹的多UAV协同轨迹规划问题进行综合评价,筛选出综合评价值满足要求的2m组轨迹;
步骤S37:将筛选出的2m组轨迹构成新的种群Pn+1,重复执行步骤S32~步骤S37,不断迭代,满足算法终止条件后停止迭代,从最终筛选出的综合评价值满足要求的2m组轨迹中选择综合评价值最小的一组轨迹作为协同轨迹输出。
5.如权利要求4所述的多UAV协同轨迹规划方法,其特征在于,所述步骤S36中采用以下公式对每组轨迹的多UAV协同轨迹规划问题进行综合评价:
J为多UAV协同轨迹规划的综合评价值,d为UAV的个数。
6.如权利要求4所述的多UAV协同轨迹规划方法,其特征在于,所述步骤S35中采用同时到达的时间作为协调变量,第i个UAV的到达时间的取值范围根据UAV的速度范围[Vmin,Vmax]确定:
其中,表示第i个UAV的第q段轨迹长度,Q为第i个UAV的轨迹段数量,t0为同时出发时刻,为了使各个UAV的飞行时间尽量小,取协同到达时间集合中的最小值作为多UAV的协同到达时间:
7.一种多UAV协同轨迹规划系统,采用如权利要求1~6任一项所述的规划方法,其特征在于,包括:
建模模块,用于对每个UAV在协同轨迹规划中的约束条件和代价函数进行建模;
轨迹表征模块,用于采用B样条曲线表示每个UAV的轨迹;
分析模块,用于采用同时到达的时间作为协调变量并基于集中协调-分布式求解的架构进行多UAV协同轨迹规划问题求解,得到满足空间协同要求和时间协同要求的协同轨迹。
8.一种设备,其特征在于,包括处理器和存储器,所述存储器中存储有计算机程序,所述处理器通过调用所述存储器中存储的所述计算机程序,用于执行如权利要求1~6任一项所述的方法的步骤。
9.一种计算机可读取的存储介质,用于存储进行多UAV协同轨迹规划的计算机程序,其特征在于,所述计算机程序在计算机上运行时执行如权利要求1~6任一项所述的方法的步骤。
CN202111011675.2A 2021-08-31 2021-08-31 多uav协同轨迹规划方法及系统、设备、存储介质 Active CN113703483B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111011675.2A CN113703483B (zh) 2021-08-31 2021-08-31 多uav协同轨迹规划方法及系统、设备、存储介质

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111011675.2A CN113703483B (zh) 2021-08-31 2021-08-31 多uav协同轨迹规划方法及系统、设备、存储介质

Publications (2)

Publication Number Publication Date
CN113703483A CN113703483A (zh) 2021-11-26
CN113703483B true CN113703483B (zh) 2024-04-09

Family

ID=78657921

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111011675.2A Active CN113703483B (zh) 2021-08-31 2021-08-31 多uav协同轨迹规划方法及系统、设备、存储介质

Country Status (1)

Country Link
CN (1) CN113703483B (zh)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
FR3131018B1 (fr) * 2021-12-21 2024-03-29 Thales Sa Procede et dispositif de generation de trajectoires d un appareil mobile
CN114415524B (zh) * 2022-03-31 2022-06-21 中国人民解放军96901部队 一种启发式协同多飞行器的轨迹交叉分析及冲突消解方法
CN116126032B (zh) * 2023-04-17 2023-07-14 华南农业大学 一种基于改进多目标进化算法的无人机群路径规划方法

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103913172A (zh) * 2013-12-06 2014-07-09 北京航空航天大学 一种适用于复杂低空下飞行器的路径规划方法
CN106842926A (zh) * 2017-02-08 2017-06-13 北京航空航天大学 一种基于正实b样条的飞行器轨迹优化方法
US10140875B1 (en) * 2017-05-27 2018-11-27 Hefei University Of Technology Method and apparatus for joint optimization of multi-UAV task assignment and path planning
CN112817330A (zh) * 2021-01-05 2021-05-18 北京联合大学 一种多无人机四维航迹协同规划方法及系统
CN112824998A (zh) * 2019-11-20 2021-05-21 南京航空航天大学 马尔可夫决策过程的多无人机协同航路规划方法和装置
CN113268087A (zh) * 2021-06-04 2021-08-17 郑州航空工业管理学院 多约束复杂环境下基于改进蚁群算法的多无人机协同工作的航迹规划方法

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103913172A (zh) * 2013-12-06 2014-07-09 北京航空航天大学 一种适用于复杂低空下飞行器的路径规划方法
CN106842926A (zh) * 2017-02-08 2017-06-13 北京航空航天大学 一种基于正实b样条的飞行器轨迹优化方法
US10140875B1 (en) * 2017-05-27 2018-11-27 Hefei University Of Technology Method and apparatus for joint optimization of multi-UAV task assignment and path planning
CN112824998A (zh) * 2019-11-20 2021-05-21 南京航空航天大学 马尔可夫决策过程的多无人机协同航路规划方法和装置
CN112817330A (zh) * 2021-01-05 2021-05-18 北京联合大学 一种多无人机四维航迹协同规划方法及系统
CN113268087A (zh) * 2021-06-04 2021-08-17 郑州航空工业管理学院 多约束复杂环境下基于改进蚁群算法的多无人机协同工作的航迹规划方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
多机协同航路规划;池彩虹;中国优秀硕士学位论文全文数据库工程科技Ⅱ辑;1-83页 *

Also Published As

Publication number Publication date
CN113703483A (zh) 2021-11-26

Similar Documents

Publication Publication Date Title
CN113703483B (zh) 多uav协同轨迹规划方法及系统、设备、存储介质
CN108388270B (zh) 面向安全域的集群无人机轨迹姿态协同控制方法
Fu et al. Phase angle-encoded and quantum-behaved particle swarm optimization applied to three-dimensional route planning for UAV
Xu et al. Optimized multi-UAV cooperative path planning under the complex confrontation environment
Cao et al. Multi-base multi-UAV cooperative reconnaissance path planning with genetic algorithm
CN108153328A (zh) 一种基于分段贝塞尔曲线的多导弹协同航迹规划方法
Zheng et al. Real-time route planning for unmanned air vehicle with an evolutionary algorithm
Lee et al. Threat evaluation of enemy air fighters via neural network-based Markov chain modeling
CN111121784B (zh) 一种无人侦察机航路规划方法
Geng et al. UAV surveillance mission planning with gimbaled sensors
Tang et al. Collision avoidance for multi-UAV based on geometric optimization model in 3D airspace
CN105144206A (zh) 多目标优化方法和设备
Fu et al. A two-layer task assignment algorithm for UAV swarm based on feature weight clustering
Serrano A bayesian framework for landing site selection during autonomous spacecraft descent
Li et al. Cooperative conflict detection and resolution and safety assessment for 6G enabled unmanned aerial vehicles
Ali et al. Feature selection-based decision model for UAV path planning on rough terrains
Fan et al. Path planning for a reconnaissance UAV in uncertain environment
Chen et al. An improved spherical vector and truncated mean stabilization based bat algorithm for uav path planning
CN109855629A (zh) 一种任务规划方法、装置及电子设备
Geng et al. Cooperative task planning for multiple autonomous UAVs with graph representation and genetic algorithm
CN116700340A (zh) 轨迹规划方法、装置及无人机集群
Xiong et al. Multi-uav 3d path planning in simultaneous attack
Chen et al. A two-stage method for UCAV TF/TA path planning based on approximate dynamic programming
Fu et al. Air combat assignment problem based on bayesian optimization algorithm
CN113255893B (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