• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    UAV Online Path Planning Algorithm in a Low Altitude Dangerous Environment

    2015-08-09 02:01:57NaifengWenLinglingZhaoXiaohongSuandPeijunMa
    IEEE/CAA Journal of Automatica Sinica 2015年2期

    Naifeng Wen,Lingling Zhao,Xiaohong Su,and Peijun Ma

    UAV Online Path Planning Algorithm in a Low Altitude Dangerous Environment

    Naifeng Wen,Lingling Zhao,Xiaohong Su,and Peijun Ma

    —UAV online path-planning in a low altitude dangerous environment with dense obstacles,static threats(STs) and dynamic threats(DTs),is a complicated,dynamic,uncertain and real-time problem.We propose a novel method to solve the problem to get a feasible and safe path.Firstly STs are modeled based on intuitionistic fuzzy set(IFS)to express the uncertainties in STs.The methods for ST assessment and synthesizing are presented.A reachability set(RS)estimator of DT is developed based on rapidly-exploring random tree(RRT)to predict the threat of DT.Secondly a subgoal selector is proposed and integrated into the planning system to decrease the cost of planning,accelerate the path searching and reduce threats on a path.Receding horizon(RH)is introduced to solve the online path planning problem in a dynamic and partially unknown environment.A local path planner is constructed by improving dynamic domain rapidly-exploring random tree(DDRRT)to deal with complex obstacles.RRT*is embedded into the planner to optimize paths.The results of Monte Carlo simulation comparing the traditional methods prove that our algorithm behaves well on online path planning with high successful penetration probability.

    Index Terms—Online UAV Path planning,threat assessment, IFS,rapidly-exploring random tree(RRT),DDRRT.

    I.INTRODUCTION

    UAV path planning in low altitude spaces is a rather complicated issue[1-4].It is required to overcome several dif fi culties,e.g.,dense and irregular obstacles,multiple dynamic constraints,partially unknown environmental information and limited planning time.The challenge is further exacerbated when planning paths in dangerous environments since a safe penetration path is required to avoid both STs and DTs under limited environmental information[5-7].RH is proved to be suitable for online path planning in dynamic and partially unknown environments[2].Thus,we propose a planning algorithm based on the framework of RH.Because DDRRT is ef fi cient in path planning under complex obstacles and multiple dynamic constraints[8-9],we construct a local path planner based on DDRRT.A threat assessor is integrated into the planning system to provide the information of imminent threats.

    The traditional discretization method and potential fi eld method are dif fi cult to be applied in high-dimensional and complex environments.Alternatively,RRT is chosen to be the basis of our algorithm.It is not nearly encumberedby the complexities of environments[10-11].It guides the growth of path trees towards unexplored continuous areas by the uniformly distributed samples.It employs the piece-wise constraint checking to deal with constraints[11-12].It can be promoted by problem-speci fi c heuristics[13-14].It is effective in creating dynamic feasible paths[10].The RS estimator of DT is also derived from RRT.

    Several background assumptions are made to constrain the focus on path planning.We assume that UAV cannot avoid threats which are ready for intercepting it,simply by increasing its speed or height.DT is de fi ned as an air defense vehicle with a slower speed than UAV,and its intention which means its most willing action is given to intercept UAV.Meanwhile, DT has perfect knowledge of UAV path planning while UAV has no prior knowledge of DT,and no traf fi c rule is available for UAV or DT.Thus,UAV needs to plan a path to avoid a DT based on the estimation of the DT motion and the higher speed of UAV than DT.The planning space is divided into the obstacle-occupied no- fl y zone,the threat area and the free fl ying space.

    The research is organized as follows.Section II brie fl y reviews relevant researches.Section III introduces the online path planning algorithm based on threat assessment.Section IV details the implementation of the algorithm.Section V presents the computational time and approximation accuracy analyses of the algorithm.Section VI presents the simulation results.In Section VII,conclusions are drawn and future works are given.

    II.OVERVIEW OF RELATED APPROACHES

    To solve UAV path planning problem in environments involving various types of STs,Miller et al.proposed a multisteps method[5].They claimed that the probabilistic model was not able to express the uncertainties in planning.The method is still a 2D planning approach,so,it may not be fully applicable to the 3D planning problem.Hanson et al. expressed threats by the fuzzy logic theory,and assessed threats by the belief network[15].But the method needs many rules and reasonings.Liu et al.presented an index based threat modeling method[7].However it is too complicated to be utilized online.Kabamba et al.provided one off-line optimal UAV path planning method based on the probabilistic threat models of missile and radar[16].

    The DT model is different from the ST model because the RS of DT is required to involve both the threatened candidate waypoints of a Host vehicle(HV,represents UAV in this study) and the threat degrees on them.The RS estimator should be able to express intentions of DT.

    Aoude et al.indicated that RRT had advantages in both accuracy and computational ef fi ciency when estimating the RSof an errant DT[17].They designed one threat assessor to help HV avoiding DT at intersections.The method combines the classifying and fi ltering techniques to identify intentions of DT[18].The RS is computed by a variant of RRT named as RRT-Reach of which the samples are created according to the predicted motion model of DT at a bias probability. In[19],the authors modeled the DT-avoiding problem by a zero-sum game.In[20],they modeled the DT-avoiding problem by the ecological recognizer method.In[21],the DT intentions were predicted online through fi ltering by learning the relationship between the states and the intentions of DT beforehand.However,the background of the methods is the intelligent transportation which is different from the problem in this study.The intentions of DT are usually unknown and HV is assumed to obey the traf fi c rules limiting its available strategies to a fi nite set.There could be collisions between DT and HV.

    There are mainly two kinds of methods for avoiding DT given the RS,i.e.,path replanning and occupied space identi fi cation.The fi rst type of method fi rstly plans paths without considering DT,then prunes the DT threatened states or sub paths from the roadmap of HV and reconstructs the roadmap for solutions.Kant et al.investigated a velocity adjustment approach[22],but the planner has no guarantee of completeness.Hsu et al.raised one probabilistic completeness algorithm by setting the extending probabilities of vertexes to be inversely proportional to the number of their adjacent vertexes[23].Jaillet et al.utilized RRT to accelerate the roadmap replanning[24].Ferguson et al.introduced a method to optimize the roadmap progressively[25].This kind of method is too complicated to be utilized online due to the boundary value problems.

    The second type of method is required to identify the occupation situations of the discrete subspaces.Fiorini et al.removed the unavailable velocity set by DT from the state space of HV[26].Frazzoli et al.measured the collision risks of subspaces by the time to collision(TTC)method[27]. Hillenbrand et al.argued that the TTC and its variants were not typically suitable for intersection scenarios where DT approached from various angles.They maintained that the time To react(TTR)was better than TTC[28].However,both TTC and TTR do not consider the intention of DT.The planner embedded with a threat assessor can compensate this shortage[28].Phillips et al.investigated one safe interval path planning method[29].It is able to reduce the computational complexity by dividing the continuous occupied time horizon of each subspace into fi nite discrete safe time intervals.This kind of method supposes that the DT motion model and the environment map are given.The assumption is hard to satisfy. It is dif fi cult to solve the high-dimensional complex problems.

    The studies of path planning online in dynamic environments are as follows.Zhang et al.utilized the model predictive control theory to fully use locally sensed data[2].Petti et al. constructed a path planning iteration model to allocate time among sub planning tasks[30].

    III.PATH PLANNING BASED ON THREAT ASSESSMENT

    A.Algorithm Model

    Fig.1 is the close-loop model of the online path planning algorithm based on threat assessment.It mainly incorporates the threat assessor(TA)and the path planner(PP).TA incorporates the submodules of ST assessment,DT assessment and threats synthesizer.Its approximated quantized results are used to assist PP in identifying safer candidate waypoints online to guarantee the safety of the planned path.The function of PP includes selecting subgoals according to the imminent threats and environmental information,subdividing the global planning task,and planning local paths.Path actuator aims to store and execute the current intermediate paths as well as feeding back the information in the current sensing horizon.

    B.ST Assessment

    1)ST Model Based on the IFS Theory:Two methods are usually applied to model STs.One regards the ST spaces as impassable areas which predigests STs excessively[31].The other models a ST by its damage probability.It is coarse since ST is related to complex factors which cannot be considered thoroughly[5].It cannot represent the uncertainties of STs[5]. Alternatively,an IFS-based method is proposed to model STs.

    A={〈x,μA(x),vA(x)〉|x?X}is called an IFS whereXis a nonempty set,μA(x)andvA(x)denote the membership and non-membership,respectively.μA:X→[0,1],x?X→μA(x)?[0,1],vA:X→[0,1],x?X→vA(x)?[0,1],0≤μA(x)+vA(x)≤1.πA(x)=1-μA(x)-vA(x) denotes the hesitancy ofx.

    We de fi neμA(x)as the active threats of STs,vA(x)as the threat-free degree of UAV,πA(x)as the uncertainties of STs.μA(x)is derived from the classic threat damage probability model to guarantee its validity.

    Fig.1.The overview of the online path planning algorithm based on threat assessment.

    The radar membership is de fi ned as

    whereRis the distance between UAV and the radar center, Rmaxis the maximum detecting radius of the radar.The height of radar detection boundaryhB=KR·L2,KRis the radar characteristic coef fi cient,andLis the horizontal distance from the radar center to the horizontal projections of objects.The area belowhBis the blind spot of radar.

    The membership of surface to air missile(SAM)μM=

    whereRis the distance between UAV and the SAM center.Rmaxis the maximum valid threat radius of SAM whileRminis the minimum one,ΔR=Rmax-Rmin.

    The membership of air defense installation(ADI)is de fi ned as

    whereRis the distance between UAV and the ADI center,Rmaxis the maximum threat radius.μCis one in the area satisfyingR<Rmin.

    Fig.2(a)shows the general appearance and the vertical section of radar membership model.Fig.2(b)and Fig.2(c) show the horizontal and vertical sections of the memberships of SAM and ADI respectively.

    Fig.2.The appearances of membership models.

    The non-membership is calculated by the real-time states of UAV.The threat-free functionNTLT,αT=

    whereHTis the height of UAV andαT(0≤αT≤π) is the absolute value of the angle between the orientation of velocity and the line from UAV to the center of ST. Denote the low altitude fl ying height asHla,the low height limitation asHlmand the maximum height limitation asHhm,LT=Llm=,Lhm=,Lla=The nonmembership contributes to reducing threats by adjustingandαT.Fig.3 shows the nonmembership function curve that varies along withLTandαT.For the requirement of IFS that 0≤μA+vA≤1,we setvA=·NTvT,αT,= 1-μA(x),vA?[0,ˉt].

    Fig.3.The function curve ofNTLT,αTthat varies along withLTandαT.

    2)Synthesizing STs:Fig.4 is an illustration of synthetic threat onWayPointi.Threatimeans a ST.Threat-Bound denotes the valid boundary of a threat.Obs denotes an obstacle.The threat onWayPointiis calculated aswhereδis the step length of RRT,ξis the function of the branch(Pathi-1,i)fromWayPointi-1toWayPointi,ζis a point onξ,dsis the unit arc length ofξ.a(ζ,u(ζ))is the threat atζ,weighted byω(ζ).u(ζ)is the UAV state atζ.The integral utilizes the IFS operator⊕.The threat onWayPointialso stands for the threat onPathi-1,i.

    Threats are regarded as continuous,thus the integral is approximated by a discretization method,as Fig.4 shows.Pathi-1,iis discretized into secondary waypoints denoted asSecWayPointj(1≤j≤s).The threats onSecWayPointjandWayPointiare computed.Then thes+1 values are synthesized by the intuitionistic fuzzy weighted averaging operator(IFWA),IFW∏Aω(a1,a2)=ω1a1⊕ω2a2⊕···⊕ωmam=whereaistands for the threat on a discrete point,weighted byωi.μaiandvaiare the membership and non-membership ofairespectively.ωi=d/δwheredis the distance between a point andWayPointi-1.P=1 are gained after normalization.The process of membership synthesizing by IFWA is consistent with the general damage probability calculation in (5).Thus the utilizing of IFWA is reasonable.

    Fig.4.An illustration of synthetic threats on UAV path tree nodes or branches.

    wherePkis the damage probability after UAV passingkSTs,andPjis the damage probability of thej-th ST.We suppose that STs are not independent.After UAV exposing ink-1 STs,the membership of thek-th STμwkis increased to=1-(1-μwk)(1-μi).μiis the synthetic threat of theith ST on the waypoints of the candidate path from the UAV path tree root to the current node.

    The threats in the intersection of STs are synthesized by IFWA.The weight of ADI or SAM is set to be twice of that of radar.The synthetic threat is larger than any single ST.IFWA is also employed to synthesize the threats of the waypoints on a path.To emphasize big threats,the weight is de fi ned asωi=Threati/Threatmax,Threatiis the threat degree of theith waypoint,Threatmaxis the maximum threat degree of all waypoints.Then the weights are normalized.

    The score function of IFS is employed to compute the threat degreeThreat(x)=μA(x)-vA(x).The threat certainty function isP(x)=1-πA(x).IfμA(x)=1,thenThreat(x)=1. IfvA(x)exceedsμA(x),thenThreat(x)=0.If the degrees of two STs are equal,the threat with biggerP(x)is more dangerous.

    C.DT Assessment

    1)Model of the DT-avoiding Problem:The DT-avoiding problem is formulated as a zero-sum differential game with perfect information of UAV path planning,constraints and free fi nal time[32].The formulation is similar to the one presented in [19]and[33].Lettfbe the earliest time DT to UAV.tf=+∞indicates UAV reaches the goal without any DT whiletf=0 means UAV has already been threatened by DT.The objective of UAV is to maximizetfwhile that of DT is to minimize it. To simplify the problem,we further assume that only one DT is involved.The game terminates with UAV reaching the goal with one feasible and low risk path.

    Let subscripts 1 and 2 denote DT and UAV respectively. The payoff function of the gameGisJ(s(t),u(t),t)=tf. Suppose thatGadmits a saddle point?U1×U2in feedback strategies.If both players execute their feedback saddle-point strategies,the value of the function isJ(s(t),t)?={J(s(t),u(t),t)}.U1andU2represent the set of admissible control functions for the two vehicles.The state of the vehicle iss(t)?S,the input isu(t)?U,with constraints.

    However,the computation of full RS online is dif fi cult.We propose an open-loop sampling-based solution to ef fi ciently approximate feasible?U1for DT.?U2is computed using a close-loop RRT-based method.The payoff function is further simpli fi ed into

    2)DT Assessment Algorithm Model:The DT assessor is subdivided into the DT estimator(DTE)and the threat assessor (TA).DTE is composed by the DT states perceptor,activity parameter calculator(APC)and reachability set estimator(RSE). RSE estimates the RS based on the DT motion parameters from APC.RS is composed of the UAV path tree nodes threatened by DT and their threat states.The states involve the estimated earliest time DT reaching the nodes and the estimated shortest distance between DT and the nodes.

    To incorporate the time-based decisions,time stamps(TSs) are added to the path tree nodes of UAV or DT.TS is the earliest time of DT or UAV arriving at their tree nodes,based on the Dubins distance.The results of RSE using TS is a more ef fi cient approximation of the RS[17-21,34].TS also leads to a low time UAV path.The differential constraints are considered in the processes of path planning of UAV or DT.The TS and other motion parameters are calculated on the path.Thus,paths have the features of trajectory,and the DT-avoiding problem is reasonably resolved.

    3)DT Reachability Set Estimator:Fig.5 shows the interpretation of DT being reachable toNodei.The tree on the horizon plane denotes the 2D DT path tree,and the 3D tree denotes the UAV path tree.Firstly,the horizontal projection point(HP(Nodei))ofNodeiis computed.Then the horizontal threat boundary(TB(Nodei))of DT around HP(Nodei) is calculated by the maximum threat radiusRof the static DT threat model.Finally,the pointBNodekon TB(Nodei) meeting the following conditions is searched.BNodekis the nearest point ofTNodejon TB(Nodei).BNodekand the line betweenBNodekandNodeiare obstacle-free.If the path betweenTNodejandBNodekis also obstacle-free,and the earliest time DT reachingBNodekis less than the TS ofNodei,then DT is reachable toNodei.Due to the continuity of threat,if DT is reachable toNodei,RSE will continue checking other nodes in the neighbors ofNodei.

    DT can adjust its motion fl exibly.Thus the DT path tree is not used to identify the path reaching a goal.Rather the entire tree is analyzed to fi nd the possible threat along each candidate path.RRT is suitable for estimating the long term DT path facing the de fi ciency of the data of its motion model[17].We extend three operators for the RS estimation based on RRT, as Fig.6 shows.The operators do not have a completeness guarantee because of the high burden of computing the full RS[17-21].

    Connect and Greedy-Connect are the basic RRT extending operators.Connect(Tree,qnear,sample)exploresTreeone control step further tosamplefromqnear.qnearis the nearest neighbor ofsample.Connect returns the new tree node and its state as well as the new branch and the control law on it. Greedy-Connect(Tree,qnear,sample)exploresTreemultisteps towardssamplefromqneargreedily,untilsampleis reached or the exploring is blocked by obstacles.

    Fig.5.The de fi nition of reachability of DT toNodei.

    Goal-bias-connect(Nodei,TNode,TTree)is a bias operator whereTTreeis the DT path tree,TNodeis the nearestTTreenode of the UAV path tree(PathTree)nodeNodei.If DT is reachable toNodeithroughTNode,execute Greedy-Connect(TTree,TNode,HP(Nodei)),otherwise utilize Connect(TTree,TNode,HP(Nodei))to explore the space.

    Reach(Nodei,TTree)is a many-to-one operator.The operator fi rstly sorts the k-nearestTTreenodes denoted asTNodej(1≤j≤k)ofNodeiby their TSs.Then it executes Goal-bias-connect(Nodei,TNodej,TTree)in turn,until DT is reachable toNodeithroughTNodejor allTNodejare extended.

    Multi-reach(TNodej,PathTree,TTree)is a one-tomany operator.It increases the opportunities of capturingPathTreenode by DT throughTNode.The k-nearest neighbors denoted asNodei(1≤i≤k)ofTNodeare searched onPathTree,then sorted by their TSs.Goal-bias-connect (Nodei,TNode,TTree)is executed iteratively,until DT is reachable toNodeithroughTNodeor allNodeiare tried. Reach,Multi-reach are utilized randomly according to a certain probability threshold.Goal-bias-connect does not run by itself but run as a sub operator.

    D.Online Path Planning Algorithm

    Because the penetration problem is dynamic with partially unknown threat information,no path can be immediately generated[35].Accordingly,the problem is subdivided into a series of more simple but realistic subproblems according to real-time states of UAV.RH is derived from the method of model predictive control,it is a feasible method to solve the online path planning problem[2].It is utilized as the basic planning frame.

    Fig.6.Three novel operators with target-bias in terms of RRT.

    The methodology of RH is similar to the human decision making activity in complex and dynamic environments[2].RH includes three sub steps,i.e.,scenario prediction(SP),rolling window optimization(RWO)and reinitialization(RI)[2].We construct the planning system according to the framework of RH as follows.In RI,we reinitialized and corrected the environmental information in the current Sensing Horizon (SH)of UAV by the real-time information.In SP,we predict the scenario and dif fi culty of planning if a discrete point on the current SH boundary is chosen as the subgoal.Accordingly, a feasible subgoal is selected.In RWO,a local path reaching the subgoal is planned.The current optimization window is de fi ned by the present SH of UAV.DDRRT is utilized for path searching while RRT*is applied to improve the path until the local planning fi nishes.The subtasks allocator aims to control the local planning for stepping forward until the fi nal goal is reached.

    The periodic RH method is based on optimization and feedback,thus it keeps the advantages of both methods[2].RH plans local paths instead of the global path to fully use the local environmental information.It is also able to accelerate planning processes and make better local controls of paths.

    1)Subgoal Selection:The subgoal of the current window is the mapping of the global goal named asGoal[2].It is required to be feasible and optimized to guide the planner to respond to UAV in limited control time horizon.

    It is assumed that the obstacles distribution is noise-free, but no threat information is available.The threats in SHs are perfectly known.It also assumes that the number of DT is 1. Subgoals are selected for sub planning tasks according to the prior and posteriori information.

    Heuristic subgoals are selected to accelerate path searching and improve path solutions.The boundary of the current SH is discretized into points.The obstacle-occupied,high risk or fl ight-limitation-violated ones are removed.The remaining points are regarded as the candidate subgoals named asSG. The costs ofSGorCostSGare calculated,and the minimum one is chosen as the current subgoal named asSubGoal.The high risk points are the ones on which the threat degrees exceed a certain thresholdTf.The selecting process ofSubGoalis denoted as

    wherePudenotes the current location of UAV,?win(Pu) denotes the boundary of the current SH.SGObsdenotes the set of points occupied by obstacles.SGThreatdenotes the set of high risk points.SGFLVdenotes the set of fl ight-limitationviolated points.“-”is the set minus operator.IfGoalis inside?win(Pu),setSubGoal=Goal.

    The cost-to-go function ofSGis de fi ned as

    whereCγiare weights,=1.Note the line from the end of the last valid sub path(Endi)toSGasLineRToSG, the line fromSGto theGoalasLineSGToG.LRToSGpartially represents the complexity of obstacles onLineRToSG.LSGToGrepresents the parameters onLineSGToG.NSTs,NDTsare the numbers of STs and DTs onLineRToSGrespectively.RSTiis the maximum threat radius of theith ST.RDTjis the maximum in fl uence radius of thejth DT.ωiandσjare the weights related to the types of STs and DTs respectively,ωi=and we setNDTs=1,σj=1.According to the length ofLineRToSGand the maximum velocity of UAV,the estimated shortest time (TSG)that UAV takes in reachingSGis obtained.The current average velocityVDT_of DT can be approximated.De fi neRDTj=VP is the projection function for the vertical direction ofLineRToSGto represent the immediate blocking effect.RDSis the static threat radius of DT.DT in fl uences one direction instead of an area.Thus the square root of the in fl uence distance is extracted.

    The heuristic inSubGoalselection is partially based on the potential fi eld method.Obstacles and threats are regarded as repulsive forces against planning while theSGorGoalprovides positive attraction for planning.Therefore,the space occupancy rate by obstacles and threats is utilized to measure the quality ofSG.LRToSGis de fi ned as

    whereHosiis the height of the obstacles onLineRToSG.Wosiis the width of the projection of the obstacles in the vertical direction ofLineRToSG.Lγjis a weight, P=1.SLRToSGis the length ofLineRToSG.Nos1is the number of the obstacles onLineRToSG. AnalogoustoLRToSG,LSGToG=Lγ1SLSGToG+Lγ2+Lγ3whereSLSGToGis the length ofLineSGToG.GHosiandGWosiare the height and width of the obstacles onLineSGToGcalculated asHosiandWosi.The weights are the same as those of formula 9.LRToSGandLSGToGare approximated parameters according to the minimum cuboid bounding boxes of obstacles.

    Fig.7 is an illustration of the process ofSubGoalselecting. The cross denotes the location of UAV when local planning starts.The cross on the top denotes the initial position while the one below denotes theGoal.The outer rectangle bounds the planning range.The star denotesSubGoal.The circle denotes SH.OneSubGoaland one SH together express an iteration of local planning.SubGoal6and the dashed frame show that the sixth local planning does not end successfully in the given time.The seventh iteration starts atEnd6.The path fromSubGoal5toEND6is reserved.TheGoalis in the ninth SH,thus“SubGoal9=Goal”.

    Fig.7.An illustration of the process ofSubGoalselecting.

    2)Subtasks Allocation:A parallel task allocation model is constructed in a manner of“planning while fl ying”to make good use of available time and sensed data,as Fig.8 shows.

    Fig.8.A parallel task allocation model.

    The time horizon(TH)is equal to the time allocated for local planning,it is set to a fi xedδc.TH must be big enough to ensure that most processes of local path-searching could be fi nished[30].After a local planning times out,the current planning is interrupted and the available intermediate solution is gained.But,TH should not exceed the execution time of any local path.Otherwise,the planning is unstable because one path may be executed before its planning is fi nished.Thus the minimum TH value of stable planning processes represents the convergence situation of the planning algorithm.TH should be less than the available time of the predicted environmental model to ensure the validity of paths[2].TH does not represent an absolutely local partition,since the tree extension in one window may bene fi t the planning in following windows.If the TH declines to 0,the algorithm will degenerate to aconventional heuristic method without fully using the sensed information.The algorithm classi fi es the tree nodes into groups by their TSs,the scale of the planning problem is limited by neglecting the nodes outside the current TH.

    The periodic iterative subtasks allocation in RH accounts for both the time constraints and the validity duration of the local path and environment[30].The rationale of our subtasks allocation scheme is the same as the traditional RH except that our subtasks incorporate local path searching with optimization.That improvement is made based on the requirement of the threat-avoiding problem.

    The task allocation model in Fig.8 is derived from[30]. The path searching aims to fi nd a safe and fl yable initial path. The optimization aims to optimize the acquired path.δscis the initial path planning time before UAV takes off.Iiis the feedback information at the beginning of theith iteration.δpiis theith local path searching time.The algorithm optimizes the obtained path inδc-δpi.δhiis the validity duration of theith local path.φiis theith intermediate solution.tfhis the fi nishing time of planning.At the beginning of each loop,following works should be done.The environmental map is updated according toIi.SubGoaliis selected.The start location for planning is set.Setδc<δhi,if no path reachingSubGoaliis found inδc,the UAV can fl y on the acquired trajectory,and a new iteration will be started,as the second iteration shows.

    3)UAV Local Path Planning Algorithm:DDRRT accelerates the path-searching process by using the environmental information[8-9].It adds dynamic domains(DD)to a node which extended unsuccessfully to reduce the extending probability of the node.It is utilized for rapidly searching paths in obstacle-complex low altitude spaces to use the defects of threats and the masking areas of obstacles.A kind of threat-based heuristic is employed to adjust the exploration probability of a tree node by its threat situation.Thus the safety of a UAV path is improved.

    A goal-bias is also added in the sampling process of DDRRT.Samples are set toSubGoalat a certain probability. The more highly the environment is constrained,the lower the probability should be.k-neighbors of samples are searched by the distance cost.The lower cost neighbors have priorities over the higher ones to extend to samples.The sample is required to be at least inside the DD of one of its neighbors to meet the requirement of DDRRT.Distance costCDistis designed as:

    whereDγiare weights,Distis the Euclidean distance betweenNodeandsample.The orientation ofNodeis de fi ned as the direction of its in-edge.Turnis the angle in radians between the orientation ofNodeand the direction of the line fromNodetosample.TSis the TS ofNode.Threatis the threat degree on the line betweenNodeandsample.DistandTurnare utilized to decrease the control complexity of the path.TSandThreatare used to reduce the duration and risk of a path.

    A threat assessment-based RRT*(TARRT*)is embedded into the planning process to optimize paths.It selects the current topological optimal extending style for path tree by excluding non-optimal solutions[36].The upper bound of cost of theNodei(UB(Nodei))is the actual cost of the path fromNodeito the currentSubGoal.If no such path exists, UB(Nodei)is set to∞.The lower bound(LB(Nodei))is the estimated lowest cost fromNodeitoSubGoal.The cost betweenNodeiandNodejin TARRT*is de fi ned as

    whereCαiare weights,Distis the length ofPathijfromNodeitoNodej.Turnis the angle in radians between the direction ofPathijand the orientation ofNodei.Threatis the threat degree onPathij.

    When a direct path fromNodetoSubGoalis found,set UB(Node)=LB(Node).The following iteration is executed tracing back to the tree root along the path.IfNodeiand its parentNodei-1satisfy UB(Nodei-1)>UB(Nodei)+ C(Nodei-1,Nodei),we set UB(Nodei-1)=UB(Nodei)+ C(Nodei-1,Nodei),and we check all childrenchjofNodei-1exceptNodei.If LB(chj)+C(Nodei-1,chj)>UB(Nodei-1),the subtree rooted atchjis removed because it cannot provide a better solution than the newly found one[36]. If UB(Nodei-1)<UB(Nodei)+C(Nodei-1,Nodei),the iteration is terminated,and a check is needed whether the subtree rooted atNodeican be removed.The computational ef fi ciency is improved,since the tree contains less nodes after being pruned[36].

    4)Constraint Checking and Path Smoothing:The feasibility of a path is ensured if differential and nonholonomic constraints are satis fi ed[37].The parameterization of the allowable state transitions via controls is dif fi cult[11].Thus RRT is applied to check constraints during the planning process. Before extending a node,whether a curvature-continuous path between the node and the possible new tree node exists or not is checked by the existence conditions in[38].If the path exists,the current control strategy(u)is determined. Then constraints are checked onu.If a constraint is violated, the extension will be canceled.Otherwise the tree node is extended.All sub paths including the connecting points are guaranteed to be feasible by the piece-wise constraint checking.After replanning or path smoothing,a fi nal fl yable path can be obtained.The state including motion parameters of every waypoint on the path can be calculated.Thus a trajectory can be planned by combining the state of each waypoint with time.

    After a local path reaching its subgoal is planned,it will be smoothed.The path is required to have smooth secondorder derivative to guarantee the continuities of the velocity and lateral acceleration of UAV.The constant-curvature arcs of the traditional Dubins path are replaced with the clothoid arcs,providing smooth curvature transitions between linear path segments.The curvature of the clothoid arc varies linearly over the path length from zero at the boundaries to a maximum at the arc mid-point.We utilize the method proposed in[38]to calculate the curvature-continuous Dubins path based on the differential geometry approach using the Frenet frame.

    To construct the 3D path,we de fi ne the orientation of the tangent of a tree node by the direction of the edge ending at the tree node.A tree node is coplanar with its father node.We use the approach proposed in[38]to construct the 3D path directlyby composite 2D paths.Firstly,the initiative and terminative maneuver planes are de fi ned by the initiative and terminative poses respectively.The straight-line maneuver of the 3D path is de fi ned by the intersection of the two planes.Two clothoid maneuvers are combined in the initiative maneuver plane to transit the tangent vector of the initiative pose to the straightline.The transited initiative pose is rotated about its tangent vector to make the binormal vector of the initiative pose to be parallel to that of the terminative pose.The new initiative pose is coplanar with the terminative pose and the fi nal curve can be constructed.The curve parameters can be calculated based on the differential geometry using the Frenet frame.

    IV.ALGORITHM IMPLEMENTATION

    Fig.9 shows the DT path treeTTreeupdating subject to the new observation of DT.The crosses denote the observed DT positions.Firstly,the Root ofTTreeis relocated to the new position of DT.A nodeqnewa step further from the current DT position is made along the predicted moving direction of DT.The direction is de fi ned from the last position to the new one.Connect(TTree,Root,qnew)is executed.TheTTreenodes with less TSs than the TS ofqneware pruned as well as the corresponding edges,as the dotted lines in Fig.9 shows.TTreeis converted into a forest.The isolating subtrees are tried to be reconnected byqnew.The subtrees which cannot connect toTTreeare kept.Finally,the TSs ofTTreenodes are updated.

    Fig.9.An illustration for updating process of DT path tree after the location of DT is changed.

    Algorithm 1.RS estimation algorithm

    1)A uniformly distributed rand numberRand?[0,1] is made.A probability thresholdptwhich decides executing reach or multi-connect is de fi ned.

    2)IfRand≤pt,reach(qnew,TTree)will be executed whereqnewis the newly generatedPathTreenode;otherwise,a random sampleqrandwill be generated,its nearestTTreenodeTNodewill be found,then multi-connect (TNode,PathTree,TTree)will be called.

    3)If DT is reachable to aPathTreenodeqr,then the algorithm will search for other DT reachable nodes in the neighborhood ofqr.qris called a DT directly reachable node. The other DT reachable nodes in the neighborhood ofqrare called indirectly reachable nodes.The estimated nearest distances between DT and the reachable nodes are calculated by the predicted time DT takes in reaching the nodes.

    Algorithm 2.Path planning algorithm

    1)If the planning time is available,Step 2 will be executed iteratively untilGoalis reached.Then Step 3 will be executed. IfGoalis not reached when the planning times out,some strategies will be utilized for the UAV navigation.

    2)If the current local planning timeLPTime≤TH and the current subgoal is not reached,then the following works will be done.Firstly,the environment andPathTreeare updated. SecondlyPathTreeis extended by the improved DDRRT and a new tree nodeqnewis generated.Finally the static threat onqnewis assessed.The RS is estimated and the dynamic threats on all the reachable nodes are computed.The DT static threat model is the same as the radar model excepthB=0. The dynamic and static threats on the reachable nodes are synthesized by IFWA.

    3)If TH-LPTime>0 and the current subgoal is reached, then the algorithm will search for a waypointRootron the acquired path closest to the current UAV location while UAV cannot reach in“TH-LPTime”.ARootris also the root of a subtree to be improved with no need for considering the in fl uence of the path actuator.TARRT*is applied to optimize the subtree in“TH-LPTime”.A local path is further obtained.

    V.COMPUTATIONAL TIME AND APPROXIMATION ACCURACY ANALYSES

    The randomized techniques do not have any running-time upper-bound:a complete motion is not guaranteed to be computed within an acceptable time[30].In the processing phase,the asymptotic time complexities of RRT and RRT* are O(nlogn)in terms of the number of simple operations, such as comparisons,additions,and multiplications,wherenis the times of sampling[36].Their time complexities are O(n)in the query phase[36].DDRRT reduces the number of collision detection calls by a sampling space reduction to make its path-searching time lower than RRT in obstaclecomplex spaces[8].Thus the upper bound of the asymptotic processing time complexity of DDRRT is O(nlogn).In the path searching process,we add the k-neighbors method into DDRRT.The k-neighbors searching can be realized without increasing the time complexity of DDRRT[36].We use RRT* in the path-optimization process.Thus the processing time complexity of our path planner is not higher than RRT.The time spent on subgoal selection is less than path planning,thus the asymptotic time complexity of our algorithm is not higher than RRT.The resulting path tree of our method has the same style as RRT,thus the query time complexity of our method is also O(n).

    DDRRT are sure to converge to any solution,this property is called the probabilistic completeness[8-9].Our method also does not destroy the probabilistic completeness.The planning path cannot be guaranteed to be optimal,because DDRRT cannot improve the result and RRT*is executed in a limited time.However some problem-speci fi c heuristics are introduced into our local planner to generate a low-cost path under complex obstacles and threats.After an initial path is obtained based on DDRRT,the path is further improved by RRT*in the remaining local planning time.

    VI.SIMULATION RESULTS

    We performed experiments to verify the effectiveness of our algorithm.We made 100 Monte Carlo experiments for different methods in the 2D and 3D scenarios respectively. The experiments are done in the computer with Intel Pentium 4,2.5GHz CPU,1GB RAM,Windows XP OS and Matlab R2010a.The simulation results verify that our planner is effective and robust under complicated dangerous environments. Our heuristicSubGoalselection is proved to be superior to the traditional Euclidean distance-based method(EDBM).The effectiveness of our IFS-based models and ST-avoiding methodare veri fi ed.The RS estimator is certi fi ed to work well in improving the ef fi ciency and accuracy of the proposed model.

    A.Simulation Environment Description

    The meanings of the simulation elements are marked in the simulation fi gures.The shadowed objects denote obstacles.The shadowed areas denote ST areas.The masking of obstacles by STs can be observed.

    The durationσof one extending step of path trees is 0.01s.The related UAV path tree branch is dispersed into 10 secondary waypoints when synthesizing STs.DT starts at[10,10]in the 2D scenario,and at[950,950,0]in the 3D scenario.The maximum static threat radii of DT are 10 and 100 in the 2D and 3D environments respectively.KRis set to 0.003 in the radar model.In the light of[8]and[9],theDDof the planning algorithm is set to the optimal value of 10δ.δis one extending step length,δ=vσwherevis the current velocity of UAV or DT.

    In the 2D scenario,the maximum velocity of UAV is 5 unit/s while the minimum one is 1 unit/s.The maximum velocities of DT is set to 1,2,3 unit/s respectively,and the minimum one is 0.2 unit/s.The results are gained when the maximum velocity of DT is 3 unit/s.The radars are located at[55,95],[12,55],[56,15]and[90,85],and the corresponding maximum detectable radii are 30,30,30 and 22.The maximum threat radius of the SAM located at[50,90] is 10,and the maximum radiuses of the other two are 8 located at[22,50]and[88,78].The minimum threat radiuses of the SAMs are 1.The maximum threat radius of ADI located at [56,10]is 8.In the path planning process,we setδsc=5s,and the radius of SH to 20.

    In the 3D scenario,the maximum velocity of the UAV is 15 unit/s while the minimum one is 5 unit/s.The maximum velocity of DT is 5 unit/s and the minimum one is 2 unit/s. TheHlain(4)is set to 100,Hlmis set to 50,Hhmis set to 300.The radars are located at[800,800,0],[300,200,0] and[400,900,0],and the maximum threat radiuses are 230, 300,300 respectively.One SAM and two ADIs locate at [400,900,0]and[720,800,0],[300,200,0]respectively.The maximum threat radius of the SAM is 150 while the minimum radius is 20.The threat radius of ADIs is 100.We setδsc=8s, and the radius of SH to 200.

    Experiments are performed for settingTH.Table I shows that whenTHis bigger or equal to 3s in the 2D scenario or bigger or equal to 5s in the 3D scenario,the Path Searching Time(PST)tends to be stable.Therefore theTHs are set to 3s and 5s in the 2D and 3D scenarios respectively.

    TABLE I COMPARISON ON EFFICIENCY INDICATORS OF THE ALGORITHM WITH DIFFERENTTHIN 2D AND 3D _________________ENVIRONMENTS

    The weights in equations are given by experts by the analytic hierarchy process.In(8),Cγ1=Cγ3=0.3,Cγ2=Cγ4=0.2.In(10),Dγ1=Dγ2=0.3,Dγ3=Dγ4= 0.2.In(11),Cα1=Cα3=0.4,Cα2=0.2.In(9),Lγ1=0.6,Lγ2=0,Lγ3=0.4 in the 2D scenario,Lγ1=0.5,Lγ2=0.2,Lγ3=0.3 in the 3D scenario.

    B.UAV Motion Model

    The UAV dynamics model extended from Dubins is denoted as

    The UAV stateXis denoted as[x,y,z,θ,φ,v]Twherex,y,zindicate the location of UAV.θandφare the yaw angle and pitch angle respectively.v?[vmin,vmax]denotes the speed.σdenotes one extending step duration.The control input is denoted asu=(a,ωθ,ωφ)wherea?[amin,amax]represents the feasible acceleration,ωθandωφare angular velocities bounded by the minimum turning radiusRmin,the turning radiusr≥Rmin.Rmin=(vmax2)/(g×(nymax2-1)1/2) wherenymaxis the maximum normal overload.Whenφis set to zero,the model is utilized by DT.

    C.Results Analyzing

    Fig.10 shows the simulations of our Heuristic DDRRT (HDDRRT)in the complicated 2D dangerous environment. The effectiveness and fl exibility of the planner are veri fi ed by the following observations.UAV avoids STs by going around them or going through the low threat areas far away from the centers of STs,as the path in Fig.10(a)shows. UAV also avoids DT by fl ying around it.The avoidances become more obvious as the maximum DT velocity increases, as Fig.10(b),(c),(d)show.The ef fi ciency and robustness of the RS estimator are illustrated by the following observations. DT path trees grow towards UAV path trees with strong bias, satisfying the requirement of(6).More candidate waypoints are pursued as the maximum velocity of DT increasing.The pursued nodes with low maximum DT velocity are probable to be found again when DT velocity is higher,meaning the DT degree is higher with larger velocity.Meanwhile,the comprehensiveness of the planner is demonstrated.Fig.10 shows the planner not only takes the threat-avoiding problem into account,but also tries to satisfy the dynamic constraints. Thus the planner not only tries to decrease the threats on a path,but also tries to reduce the length and the number of turns of the path.

    Fig.11 is the simulations of HDDRRT in the complex 3D dangerous environment.The effectiveness and fl exibility of the planner are illustrated by the following observations.The resulting paths avoid threats in the similar way as in the case of 2D results.The algorithm is inclined to search paths in low altitude to avoid STs by utilizing the defect of radar and the masking of obstacles by threats,as Fig.11(a)shows.UAV avoids DT by the masking of obstacles or going around it,as the path in Fig.11(b)shows that UAV turns obviously and fl ies low at the middle part of its path to avoid DT.Low cost paths are found according to(10)and(11).The solutions have multistyles due to the random characteristic of RRT-based methods and the real-time situations.The ef fi ciency of the RS estimator is demonstrated by the following observations.The DT path tree grows towards the UAV path tree with strong bias.The cyan squares in Fig.11(c)show the DT path.The DT path follows the UAV path reasonably.

    Fig.10.Simulations on path planning in the 2D environment.

    Fig.11.Simulations on path planning in the 3D environment.

    HDDRRT is compared with the methods of No-Fly Zone DDRRT(NFZ-DDRRT),Threat Assessment based DDRRT (TADDRRT)and TARRT*for ef fi ciency.The results are shown in Table II.TADDRRT has the same planning process as our algorithm except that it does not add any heuristic and optimize paths.The NFZ-DDRRT deletes the samples in affected areas by DT as traditional methods do.Since DDRRT cannot deal with threats itself,we add NFZ-DDRRT with the same ST-avoiding module as our method.All approaches share the same subgoal selection.

    In Table II,the indicators include:number of the UAV path tree nodes(NTN);number of collision detection callbacks (NCDC);extending success rate of the path tree nodes(ESR), ESR=NTN/NCDC;path length(PL);total threat amount (TTA),TTA is calculated by adding the threat degrees on all waypoints and secondary waypoints of a path;path cost (PC),PC is computed subject to(11),PC=Cα1·PL+Cα2· Turn+Cα3·TTA;failure rate of planning(FR),an execution fails if no global path is returned after the planning times out. The accepted planning time is set to 50s and 150s in the 2D and 3D environments respectively.NTN and NCDC repress the time complexities of algorithms.PST,NTN,NCDC and ESR are the results in the path searching stage,representing the ef fi ciencies of our method.PL,PC and FR are the results when the planning is fi nished.

    Our HDDRRT is veri fi ed to be the most effective and robust among all the aforementioned methods,and TADDRRT is better than NFZ-DDRRT and TARRT*.Table II shows the results of the 2D simulation are consistent with that of the 3D simulation.Our method has the lowest PST,NTN and NCDC,verifying its best path-searching capability.The PL and PC of our method are the lowest,certifying the ef fi ciency of our heuristic.The lowest PC also certi fi es that the paths created by our algorithm are the most fl yable and safe.The lowest PR shows our method is the most robust.The ESR of our approach is a bit lower because it tries to fi nd shorter and safer paths in the obstacle-dense areas.The high PR demonstrates NFZDDRRT is not robust,because it avoids DT only by fl ying around the DT threatening areas,but it is hard to be realized. Since long optimization time requirement is hardly satis fi ed online and the process of path- fi nding is the same as RRT,the ef fi ciency of TARRT*is the lowest.

    TABLE II COMPARISONS ON EFFICIENCY OF VARIOUS METHODS

    Table III gives the comparisons of our Heuristic Subgoal Selection method(HSS)with the traditional EDBM.The results show that our method is better than EDBM.Our method is proved to be able to speed up path-searching,reduce the path cost and avoid falling into the local optimum.The conclusion is certi fi ed by the lower PST,NTN,NCDC,PC and FR.Our method is able to assist the planner to avoid obstacles,which is demonstrated by the higher ESR.Our method not only considers PL but also takes the complexity of path-searching into account.EDBM just considers the PL,thus its PL is a bit lower than our method.EDBM is inef fi cient in these environments with multi-traps composed of a series of narrow hallways,irregular obstacles(either convex or nonconvex)and complex threats.

    TABLE III COMPARISONS ON EFFICIENCY OF VARIOUS METHODS __WITH DIFFERENT SUBGOAL SELECTING METHOD__S___

    Table IV shows the results related to the ST-avoiding abilities of algorithms.Both our method and TADDRRT can avoid STs effectively while our method is superior to TADDRRT. The ST-avoiding ability of TARRT*is the lowest.The related indicators include:synthetic threat of the solution Path(STP); path duration(PD);time of UAV fi rstly detected by STs on the solution path(TFD),since the PD of different methods varies a lot,TFD/PD is used;number of high risk waypoints (NHRW),if the threat degree of a waypoint exceeds 0.5,the waypoint will be regarded as a high risk waypoint.

    TABLE IV COMPARISONS ON THE THE ST-AVOIDING CAPABILITIES ______________OF VARIOUS ALGORITHMS

    HDDRRT and TADDRRT have better STP,TTA,TFD/PD and NHRW,because they consider STs explicitly in their heuristics.The better results of our method than those of TADDRRT also prove that the improvement of our method plays a positive role.TARRT*does not work well,because its performance relies on the optimization time which is limited in the online planning problem.

    If the RS is given,the DT-avoiding is the same as the ST-avoiding problem.Accordingly,the ef fi ciencies of different RS estimators are compared rather than the DT-avoiding abilities. If the intention of DT is given,the newly proposed RRT-reach in[19]will turn into the goal-bias-connect(GBC)method. Hence our multi-exploring(ME)method is compared to GBC as well as RRT.

    Table V shows the results of different RS estimators.The related indicators include:Synthetic DT on the solution path (SDTP),it is calculated in the same way as STP;total actual DT amount on a path(TDTA),it is calculated as TTA; estimated Total DT amount(ETDTA);accuracy of estimation (AE),AE=ETDTA/TDTA represents the accuracy of estimators;Time of one UAV candidate waypoint is fi rstly pursued by DT(TFP),because the PD of different methods varies a lot, TFP/PD is applied to repress the ef fi ciency of the estimator; RS size(RSS),it is related to the effectiveness of estimators.

    Our method has the least SDTP certifying that it is the most effective.The highest ETDTA/TDTA of our method veri fi es its highest accuracy.The earliest TFP/PD and the most RSS prove that our method is the most ef fi cient.GBC is much better than the totally random RRT estimator.

    TABLE V COMPARISONS ON THE EFFICIENCY OF VARIOUS DT _________REACHABILITY SET ESTIMATORS

    Table VI presents the results of the algorithms based on the IFS models of STs and the probabilistic models.It is proved that the IFS model is more ef fi cient than the probabilistic one in helping UAV to avoid threats and rapidly fi nd paths.The conclusion is veri fi ed by the less STP,PST,NTN,PC and the higher Penetration Success Rate(PSR)of our method based on the IFS model.PSR is computed subject to(5). The reason of the observations are that the nonmembership of the IFS model assists the path tree quickly and bypass the threat areas.The nonmembership considers the threat-freeeffect of the UAV fl ying strategy.The probabilistic model cannot consider this effect.The FR of our method is a bit lower in the 2D environment.The PL of our method is lower in the 2D environment but a little higher in the 3D environment while the PC of our method is lower.That is because the 3D environment is less complex than the 2D environment.Our method not only considers the PL as the probabilistic model but also takes the threat-avoiding strategies into account.

    TABLE VICOMPARISONS ON EFFICIENCY OF THE ALGORITHMS WITH THE IFS THREAT MODEL AND THE PROBABILISTIC ___________________THREAT MODEL

    VII.CONCLUSIONS AND FUTURE WORK

    We propose a novel method for UAV online path planning in low altitude dangerous environments.Firstly we model STs based on IFS to express the uncertainties in STs.A RRT-based RS estimator is proposed to solve the DT-avoiding problem. Secondly a heuristic subgoal selector is proposed to decrease the cost of planning,accelerate the path searching and reduce threats on a path.Finally a threat assessment-based online path planning method is designed according to RH,to plan path in a dynamic,partially unknown and dangerous environment. The simulation results certify our ST model is more reasonable than the probabilistic model.The RS estimator is veri fi ed to be more effective.The subgoal selector is certi fi ed to be superior to EDBM.Our integrated planning algorithm is veri fi ed to be more ef fi cient and robust than the traditional methods.Its good threat-avoiding ability is also certi fi ed.

    The main contributions of this study include:

    1)The IFS-based ST models are proposed to deal with the dif fi culties in probabilistic model and assess threats according to the real-time states of UAV.A RS estimator is constructed to assess the threat of DT.

    2)A heuristic subgoal selector is proposed and integrated into the planning system to guide the planner to reduce the cost of planning and rapidly obtain a low-threat path.

    3)RH is introduced to solve the online planning problem in a dynamic and partially unknown environment.DDRRT is improved for path searching in an obstacle-complex dangerous environment.

    The future work deals with the uncertainties of UAV motion and observation in the planning process.

    REFERENCES

    [1]Ye Wen,Fan Hong-Da,Zhu Ai-Hong.Mission Planning for Unmanned Aerial Vehicles.Beijing:National Defense Industry Press,2011. 1-201(in Chinese)

    [2]Zhang Chun-Gang,Xi Yu-Geng.Rolling path planning of mobile robot in dynamic unknown environment.Robot,2002,24(1):71-75(in Chinese)

    [3]Liu Gang,Lao Song-Yang,Yuan Can,Hou Lv-Lin,Tan Dong-Feng. OACRR-PSO algorithm for anti-ship missile path planning.Acta Automatica Sinica,2012,38(9):1528-1537(in Chinese)

    [4]Liu Gang,Lao Song-Yang,Tan Dong-Feng,Zhou Zhi-Chao.Research status and progress on anti-ship missile path planning.Acta Automatica Sinica,2013,39(4):347-359(in Chinese)

    [5]Miller B,Stepanyan K,Miller A,Andreev M.3D path planning in a threat environment.In:Proceedings of the 50th IEEE Conference on Decision and Control and European Control Conference.Orlando,USA: IEEE,2011.6864-6869

    [6]Aoude G S.Threat Assessment for Safe Navigation in Environments with Uncertainty in Predictability[Ph.D.dissertation],Massachusetts Institute of Technology,USA,2011.

    [7]Liu J,Li X B,Lei Y L,Wang W P.An index based threat modeling method for path planning.In:Proceedings of the 1st International Conference on Advances in System Simulation.Washington,DC:IEEE, 2009.1-5

    [8]Yershova A,Jaillet L,Simeon T,LaValle S M.Dynamic-Domain RRTs: ef fi cient exploration by controlling the sampling domain.In:Proceedings of the 2005 IEEE International Conference on Robotics and Automation. Barcelona,Spain:IEEE,2005.3856-3861

    [9]Jaillet L,Yershova A,LaValle S M,Simeon T.Adaptive tuning of the sampling domain for Dynamic-Domain RRTs.In:Proceedings of the 2005 IEEE/RSJ International Conference on Intelligent Robots and Systems.Edmonton,Canada:IEEE,2005.2851-2856

    [10]LaValle S M.Planning Algorithms.Cambridge:Cambridge University Press,2006.482-580

    [11]LaValle S M,Kuffner J J.Randomized kinodynamic planning.In: Proceedings of the 1999 IEEE International Conference on Robotics and Automation.Detroit,USA:IEEE,1999.473-479

    [12]Yang G,Kapila V.Optimal path planning for unmanned air vehicles with kinematic and tactical constraints.In:Proceedings of the 41st IEEE Conference on Decision and Control.Las Vegas,USA:IEEE,2002. 1301-1306

    [13]Urmson C,Simmons R.Approaches for heuristically biasing RRT growth.In:Proceedings of the 2003 IEEE/RSJ International Conference on Intelligent Robots and Systems.Las Vegas,USA:IEEE,2003. 1178-1183

    [14]Lee J H,Pippin C,Balch T.Cost based planning with RRT in outdoor environments.In:Proceedings of the 2008 IEEE/RSJ International Conference on Intelligent Robots and Systems.Nice,France:IEEE, 2008.684-689

    [15]Hanson M L,Sullivan O,Harper K A.On-line situation assessment for unmanned air vehicles.In:Proceedings of the 40th International Florida Arti fi cial Intelligence Research Society Conference.Florida, USA:AAAI,2001.44-48

    [16]Kabamba P T,Meerkov S M,Zeitz F H.Optimal UCAV path planning under missile threats.World Congress,2005,16(1):2002-2008

    [17]Aoude G,Joseph J,Roy N,How J P.Mobile agent trajectory prediction using bayesian nonparametric reachability trees.In:Proceedings of the 2011 AIAA Infotech Aerospace Conference.St.Louis,USA:AIAA, 2011.1587-1593

    [18]Aoude G S,Luders B D,Lee K K,Levine D S,How J P.Threat assessment design for driver assistance system at intersections.In: Proceedings of the 13th International IEEE Conference on Intelligent Transportation Systems.Funchal,Portugal:IEEE,2010.1855-1862

    [19]Aoude G S,Luders B D,How J P,Pilutti T E.Sampling-based threat assessment algorithms for intersection collisions involving errant drivers. In:Proceedings of the 2010 Symposium on Intelligent Autonomous Vehicles.Lecce,Italy:IEEE,2010.

    [20]Aoude G S,Luders B D,Levine D S,How J P.Threat-aware path planning in uncertain urban environments.In:Proceedings of the 2010IEEE/RSJ International Conference on Intelligent Robots and Systems. Taipei,China:IEEE,2010.6058-6063

    [21]Aoude G S,Luders B D,Joseph J M,Roy N,How J P.Probabilistically safe motion planning to avoid dynamic obstacles with uncertain motion patterns.Autonomous Robots,2013,35(1):51-76

    [22]Kant K,Zucker S W.Planning collision-free trajectories in time-varying environments:a two-level hierarchy.The Visual Computer,1988,3(5): 304-313

    [23]Hsu D,Kindel R,Latombe J C,Rock S.Randomized kinodynamic motion planning with moving obstacles.The International Journal of Robotics Research,2002,21(3):233-255

    [24]Jaillet L,Simeon T.A PRM-based motion planner for dynamically changing environments.In:Proceedings of the 2004 IEEE/RSJ International Conference on Intelligent Robots and Systems.Sendai,Japan: IEEE,2004.1606-1611

    [25]Ferguson D,Stentz A.Anytime,dynamic planning in high-dimensional search spaces.In:Proceedings of the 2007 IEEE International Conference on Robotics and Automation.Roma,Italy:IEEE,2007.1310-1315

    [26]Fiorini P,Shiller Z.Motion planning in dynamic environments using velocity obstacles.The International Journal of Robotics Research,1998, 17(7):760-772

    [27]Frazzoli E,Dahleh M A,Feron E.Real-time motion planning for agile autonomous vehicles.JournalofGuidance,Control,andDynamics,2002, 25(1):116-129

    [28]Hillenbrand J,Kroschel K,Schmid V.Situation assessment algorithm for a collision prevention assistant.In:Proceedings of the 2005 Intelligent Vehicles Symposium.Las Vegas,USA:IEEE,2005.459-465

    [29]Phillips M,Likhachev M.Sipp:safe interval path planning for dynamic environments.In:Proceedings of the 2011 IEEE International Conference on Robotics and Automation.Shanghai,China:IEEE,2011. 5628-5635

    [30]Petti S,Fraichard T.Safe motion planning in dynamic environments.In: Proceedings of the 2005 IEEE/RSJ International Conference on Intelligent Robots and Systems.Edmonton,Canada:IEEE,2005.2210-2215

    [31]Anderson S J,Peters S C,Pilutti T E,Iagnemma K.An optimalcontrol-based framework for trajectory planning,threat assessment,and semi-autonomous control of passenger vehicles in hazard avoidance scenarios.International Journal of Vehicle Autonomous Systems,2010, 8(2):190-216

    [32]Isaacs R.Differential Games:A Mathematical Theory with Applications to Warfare and Pursuit,Control and Optimization.New York:Dover Publications,1999.15-126

    [33]Ehtamo H,Raivio T.On applied nonlinear and bilevel programming for pursuit-evasion games.Journal of Optimization Theory and Applications, 2001,108(1):65-96

    [34]Ardiyanto I,Miura J.Real-time navigation using randomized kinodynamic planning with arrival time fi eld.Robotics and Autonomous Systems,2012,60(12):1579-1591

    [35]Kim Y,Gu D W,Postlethwaite I.Real-time path planning with limited information for autonomous unmanned air vehicles.Automatica,2008, 44(3):696-712

    [36]Karaman S,Frazzoli E.Sampling-based algorithms for optimal motion planning.The International Journal of Robotics Research,2011,30(7): 846-894

    [37]Jiang Yan,Gong Jian-Wei,Xiong Guang-Ming,ChenHui-Yan. Research on differential constraints-based planning algorithm for autonomous-driving vehicles.Acta Automatica Sinica,2013,39(12): 2012-2020(in Chinese)

    [38]Tsourdos A,White B A,Shanmugavel M.Cooperative Path Planning of Unmanned Aerial Vehicles.West Sussex:Wiley&Sons,2011.1-185

    Naifeng Wen Ph.D.candidate at the School of Computer Science and Technology,Harbin Institute of Technology.His research interests include information fusion,UAV navigation and online path planning,as well as UAVs cooperation.Corresponding author of this paper.

    Lingling Zhao Ph.D.,lecturer at the School of Computer Science and Technology,Harbin Institute of Technology.Her research interests include information fusion and target tracking.

    Xiaohong Su Ph.D.,professor at the School of Computer Science and Technology,Harbin Institute of Technology.Her research interests include information fusion,software engineering and neural networks.

    Peijun Ma Ph.D.,professor at the School of Computer Science and Technology,Harbin Institute of Technology.His research interests include information fusion,software engineering and pattern recognition.

    t

    September 26,2013;accepted May 15,2014.This work was supported by National Natural Science Foundation of China (61175027).Recommended by Associate Editor Bin Xian.

    :Naifeng Wen,Lingling Zhao,Xiaohong Su,Peijun Ma.UAV online path planning algorithm in a low altitude dangerous environment.IEEE/CAA Journal of Automatica Sinica,2015,2(2):173-185

    Naifeng Wen,Lingling Zhao,Xiaohong Su and Peijun Ma are with the School of Computer Science and Technology,Harbin Institute of Technology,Harbin 150001,China(e-mail:wennaifeng@126.com;Zhaoninglinghit@126.com;sxh@hit.edu.cn;ma@hit.edu.cn).

    中文字幕精品亚洲无线码一区| 精品久久久久久久久久久久久| 成熟少妇高潮喷水视频| 亚洲av中文字字幕乱码综合| 夜夜看夜夜爽夜夜摸| 日本色播在线视频| 男女边吃奶边做爰视频| 狠狠狠狠99中文字幕| 免费在线观看成人毛片| 天堂网av新在线| 成人亚洲欧美一区二区av| 亚洲精品成人久久久久久| 亚洲最大成人av| 高清毛片免费观看视频网站| 欧美高清性xxxxhd video| 欧美人与善性xxx| 91午夜精品亚洲一区二区三区| 亚洲精品色激情综合| 18禁在线无遮挡免费观看视频| 免费看日本二区| av黄色大香蕉| 亚洲无线观看免费| 两个人的视频大全免费| 亚洲av男天堂| 好男人视频免费观看在线| 女人被狂操c到高潮| 亚洲三级黄色毛片| 国产伦在线观看视频一区| 国产高潮美女av| 日韩精品有码人妻一区| 亚洲av第一区精品v没综合| 99热6这里只有精品| 欧美最黄视频在线播放免费| h日本视频在线播放| 亚洲成人精品中文字幕电影| 在线观看美女被高潮喷水网站| 热99re8久久精品国产| 啦啦啦观看免费观看视频高清| 熟女电影av网| 色吧在线观看| 一区二区三区免费毛片| 亚洲第一电影网av| 三级男女做爰猛烈吃奶摸视频| 看非洲黑人一级黄片| 在线观看美女被高潮喷水网站| 国产视频内射| 久久亚洲国产成人精品v| 亚洲激情五月婷婷啪啪| 久久精品综合一区二区三区| 国产成人aa在线观看| 啦啦啦韩国在线观看视频| 你懂的网址亚洲精品在线观看 | 边亲边吃奶的免费视频| 欧美日韩一区二区视频在线观看视频在线 | 成人午夜精彩视频在线观看| 欧美最黄视频在线播放免费| av免费在线看不卡| 免费av毛片视频| 亚洲欧美精品专区久久| 你懂的网址亚洲精品在线观看 | 看片在线看免费视频| 精品人妻一区二区三区麻豆| 99热精品在线国产| 色视频www国产| 亚洲欧美精品综合久久99| 国产久久久一区二区三区| 亚洲天堂国产精品一区在线| 国内精品久久久久精免费| 成人综合一区亚洲| 看片在线看免费视频| 1000部很黄的大片| 一区福利在线观看| 深夜精品福利| 欧美xxxx黑人xx丫x性爽| av黄色大香蕉| 我的老师免费观看完整版| 色综合亚洲欧美另类图片| 中国美女看黄片| 亚洲国产精品国产精品| 国产亚洲精品久久久久久毛片| 久久久久久伊人网av| 国产成人a∨麻豆精品| 亚洲av电影不卡..在线观看| 波野结衣二区三区在线| 亚洲激情五月婷婷啪啪| 熟女电影av网| 日韩精品青青久久久久久| 国产精品日韩av在线免费观看| 久久久午夜欧美精品| 亚洲成人av在线免费| 欧美最黄视频在线播放免费| 午夜福利成人在线免费观看| 国产色爽女视频免费观看| 国产综合懂色| 亚洲欧美日韩高清在线视频| 啦啦啦韩国在线观看视频| 亚洲欧美成人精品一区二区| 日韩,欧美,国产一区二区三区 | 波多野结衣巨乳人妻| 欧美激情久久久久久爽电影| 欧美区成人在线视频| 男人舔奶头视频| 亚洲18禁久久av| 久久6这里有精品| 免费观看的影片在线观看| 欧美zozozo另类| 日本一本二区三区精品| 看十八女毛片水多多多| 成人亚洲欧美一区二区av| 精品国内亚洲2022精品成人| 欧美区成人在线视频| 悠悠久久av| 深夜精品福利| 国产成人福利小说| 亚洲18禁久久av| 久久人人精品亚洲av| 久久久久久久久久久免费av| 美女脱内裤让男人舔精品视频 | 嫩草影院新地址| 亚洲不卡免费看| 国产精品人妻久久久久久| 乱码一卡2卡4卡精品| 高清在线视频一区二区三区 | 毛片一级片免费看久久久久| 精品久久久久久久末码| а√天堂www在线а√下载| 51国产日韩欧美| av黄色大香蕉| 国产精品乱码一区二三区的特点| 小蜜桃在线观看免费完整版高清| 黄色一级大片看看| 欧美最黄视频在线播放免费| 人人妻人人澡人人爽人人夜夜 | 日韩欧美三级三区| 18禁在线播放成人免费| 午夜老司机福利剧场| 插阴视频在线观看视频| 国产又黄又爽又无遮挡在线| 黄色日韩在线| 成人亚洲欧美一区二区av| 日本色播在线视频| 国产精品三级大全| 联通29元200g的流量卡| 久久精品久久久久久久性| 六月丁香七月| 国产国拍精品亚洲av在线观看| 亚洲av一区综合| 超碰av人人做人人爽久久| 久久久久网色| 又爽又黄无遮挡网站| 99热全是精品| 国产精品美女特级片免费视频播放器| 国产视频内射| 在线播放国产精品三级| 99久久精品一区二区三区| a级一级毛片免费在线观看| 亚洲图色成人| 亚洲色图av天堂| 看免费成人av毛片| 亚洲av男天堂| av女优亚洲男人天堂| 欧美变态另类bdsm刘玥| 成人国产麻豆网| 久久久久九九精品影院| 国产精品野战在线观看| 精品久久久久久久久亚洲| 精品久久久久久久末码| 久久久久久久久中文| 亚洲精品成人久久久久久| 亚洲精品国产av成人精品| av在线老鸭窝| 久久久午夜欧美精品| 久久韩国三级中文字幕| 欧美最黄视频在线播放免费| 日本色播在线视频| 精品欧美国产一区二区三| 午夜免费男女啪啪视频观看| av免费观看日本| 一本久久精品| 精品久久久久久久人妻蜜臀av| 中文欧美无线码| 六月丁香七月| 97超碰精品成人国产| eeuss影院久久| 美女内射精品一级片tv| 少妇的逼水好多| 久久精品国产自在天天线| 99久久精品国产国产毛片| 国产午夜精品一二区理论片| 亚洲美女视频黄频| 高清毛片免费观看视频网站| 国产v大片淫在线免费观看| 午夜激情欧美在线| 1000部很黄的大片| 日韩中字成人| 最近视频中文字幕2019在线8| 国产片特级美女逼逼视频| 美女 人体艺术 gogo| a级毛色黄片| 99久久人妻综合| 久久这里有精品视频免费| 国产女主播在线喷水免费视频网站 | 亚洲av成人av| 精品久久久久久久末码| 久久精品国产99精品国产亚洲性色| 国产成人精品一,二区 | 最近中文字幕高清免费大全6| 欧美激情在线99| 国产一级毛片在线| 日韩欧美在线乱码| 欧美bdsm另类| 插阴视频在线观看视频| 亚洲在线观看片| 日日摸夜夜添夜夜添av毛片| 高清日韩中文字幕在线| 韩国av在线不卡| 欧美区成人在线视频| 听说在线观看完整版免费高清| 特大巨黑吊av在线直播| 国产成人精品久久久久久| 久久综合国产亚洲精品| 国产又黄又爽又无遮挡在线| av免费观看日本| 国产老妇伦熟女老妇高清| 伦精品一区二区三区| 色吧在线观看| 一级av片app| av.在线天堂| 一级二级三级毛片免费看| 国内久久婷婷六月综合欲色啪| 麻豆久久精品国产亚洲av| 深夜a级毛片| 久久综合国产亚洲精品| 国产久久久一区二区三区| 此物有八面人人有两片| 午夜亚洲福利在线播放| 啦啦啦韩国在线观看视频| 国内揄拍国产精品人妻在线| 综合色丁香网| 亚洲在线观看片| 少妇被粗大猛烈的视频| 久久久久性生活片| 欧美日韩一区二区视频在线观看视频在线 | 能在线免费看毛片的网站| 免费一级毛片在线播放高清视频| 九九久久精品国产亚洲av麻豆| 狂野欧美白嫩少妇大欣赏| 亚洲经典国产精华液单| 成人漫画全彩无遮挡| 国产综合懂色| 久久精品久久久久久噜噜老黄 | 国产激情偷乱视频一区二区| 亚洲一级一片aⅴ在线观看| 国产亚洲欧美98| av在线蜜桃| 91在线精品国自产拍蜜月| 精品国内亚洲2022精品成人| 日日撸夜夜添| 如何舔出高潮| 人人妻人人澡欧美一区二区| 亚洲欧美日韩高清专用| 白带黄色成豆腐渣| 国产黄色小视频在线观看| 99热这里只有是精品50| 午夜a级毛片| 老司机影院成人| 黄色配什么色好看| 18禁在线播放成人免费| 欧美潮喷喷水| 深夜精品福利| 精品久久久久久成人av| 免费av观看视频| 99视频精品全部免费 在线| 亚洲人与动物交配视频| 日韩一区二区三区影片| 久久久久免费精品人妻一区二区| 国产爱豆传媒在线观看| 综合色丁香网| 久久精品国产亚洲网站| 欧美最黄视频在线播放免费| 日韩国内少妇激情av| 全区人妻精品视频| 九九在线视频观看精品| 99热只有精品国产| 一边摸一边抽搐一进一小说| 欧美成人一区二区免费高清观看| 日本黄色片子视频| 精品久久久久久久人妻蜜臀av| 91午夜精品亚洲一区二区三区| 我要看日韩黄色一级片| 亚洲高清免费不卡视频| 熟女人妻精品中文字幕| 51国产日韩欧美| 久久久久久久久久成人| 国内少妇人妻偷人精品xxx网站| 国产视频首页在线观看| 卡戴珊不雅视频在线播放| 又爽又黄无遮挡网站| 网址你懂的国产日韩在线| 亚洲国产欧美人成| 中文字幕免费在线视频6| 免费人成在线观看视频色| 久久久久久久久大av| 亚洲一区高清亚洲精品| 久久精品人妻少妇| 亚洲成人久久爱视频| 国产精品国产高清国产av| 深爱激情五月婷婷| 国产色婷婷99| 国产男人的电影天堂91| 精品人妻一区二区三区麻豆| 菩萨蛮人人尽说江南好唐韦庄 | 一进一出抽搐动态| ponron亚洲| 一本久久中文字幕| 国产毛片a区久久久久| 久久综合国产亚洲精品| av在线亚洲专区| a级毛片免费高清观看在线播放| 黄色欧美视频在线观看| 91麻豆精品激情在线观看国产| 人妻久久中文字幕网| 日产精品乱码卡一卡2卡三| 国产亚洲精品久久久com| 国产日韩欧美在线精品| 美女国产视频在线观看| 亚洲国产精品国产精品| 色综合亚洲欧美另类图片| 午夜福利成人在线免费观看| 黄片wwwwww| 桃色一区二区三区在线观看| 91午夜精品亚洲一区二区三区| 精品一区二区三区视频在线| 日韩欧美 国产精品| 一本久久精品| 在线a可以看的网站| 两个人视频免费观看高清| 欧美变态另类bdsm刘玥| 欧美性猛交黑人性爽| av专区在线播放| АⅤ资源中文在线天堂| 精品久久久久久久久久免费视频| 一本—道久久a久久精品蜜桃钙片 精品乱码久久久久久99久播 | 麻豆久久精品国产亚洲av| 春色校园在线视频观看| 永久网站在线| 国产淫片久久久久久久久| 亚洲人成网站在线播放欧美日韩| 男人的好看免费观看在线视频| 国产白丝娇喘喷水9色精品| 亚洲av第一区精品v没综合| 麻豆成人av视频| 青春草亚洲视频在线观看| 国产精品野战在线观看| 九九热线精品视视频播放| 成人一区二区视频在线观看| 哪里可以看免费的av片| 联通29元200g的流量卡| 国产精品三级大全| 亚洲人成网站在线观看播放| 亚洲国产欧美人成| 精品久久久久久久人妻蜜臀av| 色哟哟·www| 亚洲丝袜综合中文字幕| 综合色av麻豆| 免费av观看视频| 噜噜噜噜噜久久久久久91| 精品国内亚洲2022精品成人| 丰满乱子伦码专区| 国产午夜精品论理片| 欧美一区二区亚洲| 欧美+日韩+精品| 日本免费a在线| 日韩成人伦理影院| 免费看美女性在线毛片视频| 麻豆成人av视频| 日本免费a在线| 一个人看视频在线观看www免费| or卡值多少钱| 国产伦理片在线播放av一区 | 老师上课跳d突然被开到最大视频| 精品国产三级普通话版| 成年av动漫网址| 久久精品国产自在天天线| 国产精品久久久久久久久免| 久久这里只有精品中国| 真实男女啪啪啪动态图| 国产成人福利小说| 国产三级中文精品| 麻豆乱淫一区二区| 国产精品福利在线免费观看| 亚洲av不卡在线观看| 免费观看a级毛片全部| 成年免费大片在线观看| www日本黄色视频网| 国产精品蜜桃在线观看 | 欧美色欧美亚洲另类二区| 国产精品一区二区性色av| 最近手机中文字幕大全| 成人毛片a级毛片在线播放| 99在线视频只有这里精品首页| 少妇高潮的动态图| 亚洲va在线va天堂va国产| 人体艺术视频欧美日本| 精品一区二区三区视频在线| 99国产极品粉嫩在线观看| 成年免费大片在线观看| 久久精品久久久久久噜噜老黄 | 国产精品一区二区在线观看99 | 美女国产视频在线观看| 精品久久久噜噜| 亚洲18禁久久av| 成人一区二区视频在线观看| 国产精品精品国产色婷婷| av在线亚洲专区| 亚洲av第一区精品v没综合| 免费无遮挡裸体视频| 毛片一级片免费看久久久久| 精品熟女少妇av免费看| 亚洲欧美中文字幕日韩二区| 国产三级中文精品| 极品教师在线视频| 久久精品国产亚洲av香蕉五月| 乱码一卡2卡4卡精品| 国产精品一区二区三区四区久久| 国产精品野战在线观看| 免费观看人在逋| 国产精品.久久久| 欧美成人a在线观看| 在线a可以看的网站| 久久九九热精品免费| 我的老师免费观看完整版| 热99re8久久精品国产| 国产极品天堂在线| 久久草成人影院| 亚洲美女视频黄频| 十八禁国产超污无遮挡网站| 日韩欧美一区二区三区在线观看| 色视频www国产| 两个人的视频大全免费| 国产黄a三级三级三级人| 少妇人妻精品综合一区二区 | 亚洲乱码一区二区免费版| 日韩av在线大香蕉| 成人永久免费在线观看视频| 国国产精品蜜臀av免费| 久久久久久久亚洲中文字幕| 久久亚洲精品不卡| 国产精品一二三区在线看| 国产高潮美女av| 1024手机看黄色片| 亚洲第一区二区三区不卡| 亚洲精品影视一区二区三区av| 亚洲,欧美,日韩| 欧美性猛交黑人性爽| 欧美潮喷喷水| 亚洲国产欧美在线一区| 国产色爽女视频免费观看| 丰满的人妻完整版| 又爽又黄无遮挡网站| 国产极品天堂在线| 国产麻豆成人av免费视频| av.在线天堂| 美女cb高潮喷水在线观看| 色综合色国产| 成人漫画全彩无遮挡| 国产精品人妻久久久久久| 国产精品日韩av在线免费观看| 黄色欧美视频在线观看| 99久久中文字幕三级久久日本| 久久这里有精品视频免费| 午夜福利在线观看免费完整高清在 | 久久久久国产网址| 精品久久久久久久末码| 国产爱豆传媒在线观看| 国产精品蜜桃在线观看 | 午夜福利成人在线免费观看| 日韩三级伦理在线观看| 欧美变态另类bdsm刘玥| 日日摸夜夜添夜夜爱| 日本黄色视频三级网站网址| 可以在线观看的亚洲视频| 久久99蜜桃精品久久| 精品99又大又爽又粗少妇毛片| 国产私拍福利视频在线观看| 精品久久久久久久人妻蜜臀av| 美女xxoo啪啪120秒动态图| 黄色欧美视频在线观看| 精品午夜福利在线看| 亚洲aⅴ乱码一区二区在线播放| a级毛色黄片| 男女啪啪激烈高潮av片| 最近的中文字幕免费完整| 日韩欧美 国产精品| 久久九九热精品免费| 国产精品精品国产色婷婷| 午夜福利成人在线免费观看| 欧美人与善性xxx| 免费大片18禁| 中文字幕av在线有码专区| 爱豆传媒免费全集在线观看| 色尼玛亚洲综合影院| 亚洲欧美精品专区久久| 一级黄片播放器| 精品人妻偷拍中文字幕| 最近2019中文字幕mv第一页| 久久热精品热| 麻豆成人午夜福利视频| av.在线天堂| 亚洲最大成人手机在线| 永久网站在线| 国产老妇伦熟女老妇高清| 精品熟女少妇av免费看| 国产老妇伦熟女老妇高清| 国产高潮美女av| 国产精品美女特级片免费视频播放器| 22中文网久久字幕| 中文字幕精品亚洲无线码一区| 在线观看av片永久免费下载| 国产人妻一区二区三区在| 国产三级中文精品| 日本黄色片子视频| 亚洲成人中文字幕在线播放| 精品久久久久久成人av| 亚洲乱码一区二区免费版| 久久人人爽人人片av| 国产精品久久久久久亚洲av鲁大| 韩国av在线不卡| 国产 一区精品| 可以在线观看的亚洲视频| 国产精品久久久久久亚洲av鲁大| 啦啦啦啦在线视频资源| 欧美日韩在线观看h| 精品久久久久久成人av| 日韩欧美一区二区三区在线观看| 又黄又爽又刺激的免费视频.| 亚洲欧美精品专区久久| 晚上一个人看的免费电影| 亚洲熟妇中文字幕五十中出| 成人高潮视频无遮挡免费网站| 春色校园在线视频观看| 精品人妻视频免费看| 成年av动漫网址| 欧美一区二区国产精品久久精品| .国产精品久久| 久久人人爽人人片av| 亚洲av男天堂| 亚洲,欧美,日韩| 久久午夜亚洲精品久久| 亚洲欧美日韩无卡精品| 啦啦啦观看免费观看视频高清| av在线播放精品| 日本成人三级电影网站| 一级黄片播放器| 久久久久久久久大av| 欧美zozozo另类| 22中文网久久字幕| 国产精品永久免费网站| 久久久久久国产a免费观看| 国产探花极品一区二区| 美女xxoo啪啪120秒动态图| 久久精品91蜜桃| 国产v大片淫在线免费观看| 哪里可以看免费的av片| 三级经典国产精品| 国产精品久久久久久精品电影小说 | 三级经典国产精品| 非洲黑人性xxxx精品又粗又长| 亚洲综合色惰| 久久亚洲国产成人精品v| 午夜老司机福利剧场| 蜜臀久久99精品久久宅男| 99久久成人亚洲精品观看| 欧美一区二区亚洲| av在线蜜桃| 内地一区二区视频在线| 久久久久久久午夜电影| 人体艺术视频欧美日本| 一级av片app| 丰满人妻一区二区三区视频av| 亚洲无线观看免费| 国产亚洲精品久久久com| 黄片wwwwww| 久久久久久久久中文| 亚洲va在线va天堂va国产| 成人三级黄色视频| 国产激情偷乱视频一区二区| 午夜老司机福利剧场| 亚洲国产精品sss在线观看| 可以在线观看毛片的网站| 亚洲av熟女| 日本免费一区二区三区高清不卡| 成年av动漫网址| 乱人视频在线观看| 国内精品美女久久久久久| 女的被弄到高潮叫床怎么办| 国产黄a三级三级三级人| 三级男女做爰猛烈吃奶摸视频| 热99re8久久精品国产| 国产精品国产三级国产av玫瑰| 特级一级黄色大片| 国产老妇伦熟女老妇高清| 男女做爰动态图高潮gif福利片| 午夜免费男女啪啪视频观看| 性色avwww在线观看| 婷婷色综合大香蕉| 18+在线观看网站| 插阴视频在线观看视频| 日韩中字成人| 97超碰精品成人国产| 好男人在线观看高清免费视频| 男人舔女人下体高潮全视频| 久久亚洲国产成人精品v| 精品一区二区三区视频在线| 国产成人91sexporn|