CN102749080B - 一种基于流体力学的无人机三维航路生成方法 - Google Patents

Abstract

本发明公开了一种基于流体力学的无人机三维航路生成方法，属于航路规划领域，本发明首先计算得到沿X轴来流方向，位于原点的圆球障碍绕流流线，利用旋转平移矩阵得到平面中任意来流方向与任意障碍位置的绕流流线，根据流线是否满足部分重合条件采用两种方式将所有的流线合并成一条流线，之后根据无人机约束对三维流线进行处理得到无人机三维航路。本发明借鉴了自然界流水能够避开岩石的现象，将流体计算与航路规划相结合，同时考虑无人机飞于约束，在地形较为复杂时，能够规划出光滑且易于飞行的三维飞行避障航路。本发明地形建模简单，计算量小，满足无人机约束，实现方便。

Description

一种基于流体力学的无人机三维航路生成方法

技术领域

本发明属于航路规划领域，具体地说是指一种基于流体力学的无人机三维航路生成方法。

背景技术

航路规划是影响无人机使用和飞行安全的关键技术，一直受到各方面的高度重视，经过几十年的研究和发展，取得了大量研究成果，为目前无人机的大发展奠定了基础。随着无人机的使用空域从中、高空不断向低空、甚至超低空的拓展，复杂地形下的无人机航路规划能力日益成为制约无人机发展的瓶颈。随着飞行高度的不断降低，地形的影响将成为航路规划需要考虑的重要因素，对于复杂地形情况还将成为关键因素。

目前见诸文献的航路规划方法有：基于图形的规划方法、决策型搜索方法、随机搜索方法和人工势场法等。这些航路规划方法通常处理的障碍数目较少且相对简单，但却需要复杂的网格或数据结构进行建模。三维航路规划方法相较于二维更复杂，许多方法在复杂地形条件下存在组合爆炸的问题，并不适合复杂地形下的无人机三维航路规划。

另一方面无人机自身具有一定的飞行性能约束，现有方法规划出的航路对于无人机的实际飞行存在困难。究其原因是航路不够光滑，没有考虑无人机约束。复杂地形条件约束下的无人机航路规划不但要考虑复杂地形对路径的约束，还要考虑地形对飞行安全的潜在影响，以及无人机自身的性能约束，并在垂直面内的地形跟踪和水平面的地形规避之间实现综合，最终规划出三维空间的可供无人机实际飞行的平滑航路。

人工势场法中的电荷法与流函数法分别利用了电磁理论和流体计算中的概念，并将其引入到航路规划中。虽然这两种方法能够快速的规划出相对平滑的航路，但是只能作为二维航路规划方法。如何能够在三维情况下，用一种简单的地形建模方法规划出相对平滑，又适合无人机飞行的三维航路是设计无人机航路规划系统时需要面对的一个问题。

发明内容

针对现有技术中存在的问题，本发明提出一种基于流体力学的无人机三维航路生成方法，通过该方法能够在地形较为复杂的情况下，快速的生成一条满足无人机约束的光滑三维避障飞行航路。

应用本发明提供的航路规划方法，首先计算出理想情况下沿X轴来流方向，位于原点的圆球障碍绕流流线，利用旋转平移矩阵得到XY平面中任意来流方向与任意障碍位置的绕流流线，根据障碍物间距离的不同采用两种方式将所有的流线合并成一条流线，之后根据无人机约束对三维流线进行处理得到无人机三维航路。具体包括如下步骤：

步骤一：单个障碍时，根据定常忽略粘性的不可压缩理想流体，沿空间直角坐标系X轴来流方向作用于位于原点的圆球障碍绕流问题的解析解，在该直角坐标系下求得流速沿X、Y和Z轴的分量；

定常忽略粘性的不可压缩理想流体沿X轴来流方向作用于位于原点的圆球障碍绕流问题，由于流体的运动存在对称性，描述为平面二维极坐标系下的扰流问题，对于空间中任意一点P，极坐标系建立在点P与X轴所构成的平面上，空间直角坐标系的坐标原点与平面二维极坐标系的极点重合，为点O，空间直角坐标系的X轴与平面二维极坐标系的极轴X′轴重合，在该极坐标系下流体的速度势为

其中OP与极轴之间的夹角为θ角，点P与点O的距离为|OP|＝r，V_∞为流体的流速，a为圆球障碍的半径，在空间直角坐标系下，P的坐标为(x,y,z)， 用u、v和w代表三维情况下流速沿空间直角坐标系X、Y和Z轴的分量，则

步骤二：单个障碍时，采用解微分方程组方法，对流速沿X、Y和Z轴的分量所构成的速度势微分方程组进行求解计算，得到沿X轴来流方向的定常忽略粘性的不可压缩理想流体作用于位于原点的圆球障碍绕流流线；

步骤三：单个障碍时，利用旋转平移矩阵，使计算的流线适于空间直角坐标系下XY平面内任何来流方向，并且使多个障碍在空间直角坐标系下XY平面内任意分布，在直角坐标系的XY平面，流线的起点为(x₀,y₀)，终点为(x_c,y_c)，球心坐标为(x_r,y_r)，旋转平移矩阵的计算过程包括以下几个步骤：

(1)计算通过点(x₀,y₀)与点(x_c,y_c)的直线方程y＝k₁x+b₁；k₁为该直线的斜率，b₁为该直线在Y轴上的截距；

(2)计算过点(x_r,y_r)并垂直于直线y＝k₁x+b₁的直线方程；

(3)求得垂足的坐标，并计算L₁，L₂与L₃，其中L₁是垂足到点(x_c,y_c)的距离，L₂是垂足到(x₀,y₀)的距离，L₃是垂足到(x_r,y_r)的距离；

(4)当转换以后点(x_c′,y_c′)在Y轴正半轴，有(x₀′,y₀′)=(-L₂,L₃)，(x_c′,y_c′)=(L₁,L₃)；否则(x₀′,y₀′)=(L₂,-L₃)，(x_c′,y_c′)=(L₁,-L₃)；

定义点(x₀,y₀)到点(x_c,y_c)的向量为 点(x₀′,y₀′)到点(x_c′,y_c′)的向量为向量和间夹角α为如果障碍的高度，则旋转平移矩阵T为

$T=\left[\begin{array}{cccc}\cos \mathrm{\α}& \sin \mathrm{\α}& 0& 0\\ -\sin \mathrm{\α}& \cos \mathrm{\α}& 0& 0\\ 0& 0& 1& 0\\ x_r& y_r& z_r& 1\end{array}\right]$

z_r为障碍的高度，设障碍在原点O时的流线数据为矩阵A，其中矩阵A的每一行代表一个航路点，矩阵A的列数为4列，第一列为航路点的X轴方向坐标，第二列为航路点的Y轴方向坐标，第三列为航路点的Z轴方向坐标，第四列为单位1，任意位置的流线矩阵A′为A′＝A×T；

步骤四：两个障碍时，分别计算两个障碍单独避障的各自流线，判断两障碍单独避障的流线是否满足部分重合条件，计算两个障碍的总流线；

(1)当两个障碍单独避障的流线能够满足部分重合时，此时两个障碍之间的两条流线重合直线线段部分长度大于等于其中O_i和O_i+1分别代表第i个和第i+1个障碍的球心，|O_iO_i+1|表示两个障碍之间的直线距离，|S_i′S_i+1′|表示两个障碍之间的两条流线重合直线线段部分长度，找到两条流线中间相重合的部分，进行流线数据合并作为总流线，合并的过程具体为：

(1)找到流线S_i和S_i+1在L_i区域内的两个流线部分；

(2)流线S_i和S_i+1在L_i区域内的两个流线部分中，每两个相邻的航路点组成一个航路段，分别计算每个航路段相较于来流方向的斜率；

(3)将流线S_i和S_i+1在L_i区域内的两个流线部分中，各航路段的斜率之差的绝对值最小的部分作为合并的位置，并剔除其余的部分；

(2)当两个障碍单独避障的流线不存在部分重合的直线线段部分时，此时两个障碍之间的两条流线重合直线线段部分长度小于计算与两个障碍同时内切的虚拟障碍，对该虚拟障碍利用步骤一～步骤三的方法计算流线，该流线为两个障碍的总流线，该流线具有步骤三中矩阵A的形式，为n行4列的矩阵，n根据流线航路点个数的不同而不同；

Q_i和O_i+1分别代表第i个和第i+1个障碍的球心，各自的避障流线分别为S_i与S_i+1，R_i和R_i+1分别表示Q_i和O_i+1的半径，在线段Q_iQ_i+1上存在点O_new，以该点为圆心，做出与圆Q_i和O_i+1同时相内切的圆O_new，将圆O_new作为虚拟障碍进行躲避，O_new满足：

|O_newO_i|+R_i＝|O_newO_i+1|+R_i+1

步骤五：针对所有障碍，将每两个相邻障碍作为一组，分别计算每个障碍各自的流线，再采用步骤四计算总流线的方法，分别计算每一组中两个障碍的总流线，得到每一组的总流线；然后再采用步骤四中的计算总流线的方法，计算相邻两组总流线再合并后的总流线，顺次采用步骤四中计算总流线的方法，不断合并流线，直至最后产生一条流线最为最终三维流线；

步骤六：步骤一中计算的流线是一种平行来流，在三维直角坐标系下，在YZ平面的坐标范围内改变来流的起始位置，经过步骤一～步骤五的方法，获得多条三维流线；

步骤七：将无人机约束表示为流线的航路点和航路段约束，满足无人机约束的流线作为无人机三维可飞航路；

分别将无人机最大爬升高度(H)约束、最大爬升角(α)约束、最大水平转弯角(Φ)约束转化为对流线的约束如下：

(1)最大爬升高度(H)约束

提取出每条流线中各点的高度，将其与最大爬升高度(H)进行比较，所有的Z轴的坐标z均满足z≤H，则该条流线满足最大爬升高度约束；

(2)最大爬升角(α)约束

计每条流线中相邻两点构成的航路段与空间直角坐标系的水平面的夹角θ_z，将其与最大爬升角(α)进行比较，流线中所有θ_z均满足θ_z≤α，则该条流线满足最大爬升角约束；

(3)最大水平转弯角(Φ)约束

计算每条流线中相邻两航段在空间直角坐标系的水平面内投影的夹角φ_xy，将其与最大转弯角(Φ)进行比较，流线中所有φ_xy均满足φ_xy≤Φ，则该条流线满足最大水平转弯角约束；

对步骤六中所有三维流线按照上述三种约束进行筛选，满足(1)~(3)条所有约束条件的流线作为无人机的三维可飞航路。

本发明的优点在于：

(1)本发明提出一种基于流体力学的无人机三维航路生成方法，借鉴了自然界流水能够避开岩石的现象，将流体计算与航路规划相结合，同时考虑无人机飞行约束，在地形较为复杂时，能够规划出光滑且易于飞行的三维飞行避障航路。

(2)本本发明提出一种基于流体力学的无人机三维航路生成方法，地形建模简单，计算量小，满足无人机约束，实现方便。

附图说明

图1是本发明中旋转平移矩阵求解示意图；

图2是本发明中旋转平移矩阵求解示意图；

图3是本发明中障碍间距较大时流线的计算示意图；

图4是本发明中障碍间距较小时流线的计算示意图。

具体实施方式

下面将结合附图对本发明的具体实施步骤作进一步说明。

本发明提出一种基于流体力学的无人机三维航路生成方法。首先计算得到沿X轴来流方向，位于原点的圆球障碍绕流流线，利用旋转平移矩阵得到XY平面中任意来流方向与任意障碍位置的绕流流线，根据流线是否满足部分重合条件采用两种方式将所有的流线合并成一条流线，之后根据无人机约束对三维流线进行处理得到无人机三维航路。

本发明提出一种基于流体力学的无人机三维航路生成方法，具体包括以下几个步骤：

步骤一：单个障碍时，根据定常忽略粘性的不可压缩理想流体，沿空间直角坐标系X轴来流方向作用于位于原点的圆球障碍绕流问题的解析解，在该直角坐标系下求得流速沿X、Y和Z轴的分量。

定常忽略粘性的不可压缩理想流体沿X轴来流方向作用于位于原点的圆球障碍绕流问题，由于流体的运动存在对称性，可以描述为平面二维极坐标系下的扰流问题，具体如图1所示。对于空间中任意一点P，图1中的极坐标系建立在点P与X轴所构成的平面上。空间直角坐标系的坐标原点与平面二维极坐标系的极点重合，为点O，空间直角坐标系的X轴与平面二维极坐标系的极轴X′轴重合。在该极坐标系下流体的速度势为

OP与极轴之间的夹角为θ角，点P与点O的距离为|OP|＝r，V_∞为流体的流速，a为圆球障碍的半径，如图1。在空间直角坐标系下，P的坐标为(x,y,z)，则有 用u、v和w代表三维情况下流速沿空间直角坐标系X、Y和Z轴的分量，则

步骤二：单个障碍时，采用常用的解微分方程组方法(比如四阶或五阶的定步长或变步长Runge-Kutta方法)，对流速沿X、Y和Z轴的分量所构成的速度势微分方程组进行求解计算，得到沿X轴来流方向的定常忽略粘性的不可压缩理想流体作用于位于原点的圆球障碍绕流流线。微分方程组在MATLAB软件中建立后，可用调用Ode45函数或Ode23函数求解。

步骤三：单个障碍时，利用旋转平移矩阵，使步骤二中计算的流线适于直角坐标系下XY平面内任何来流方向，并且使多个障碍可以在直角坐标系下XY平面内任意分布。因此，直角坐标系下XY平面中任意来流方向与任意障碍物位置的流线计算问题，被转化为沿X轴来流方向，障碍物位于原点的圆球障碍绕流流线计算问题。

如图2所示直角坐标系的XY平面，流线的起点为(x₀,y₀)，终点为(x_c,y_c)，球心坐标为(x_r,y_r)，旋转平移矩阵的计算过程包括以下几个步骤：

(1)计算通过点(x₀,y₀)与点(x_c,y_c)的直线方程y＝k₁x+b₁；k₁为该直线的斜率，b₁为该直线在Y轴上的截距；

(2)计算过点(x_r,y_r)并垂直于直线y＝k₁x+b₁的直线方程；

(3)求得垂足的坐标，并计算L₁，L₂与L₃，其中L₁是垂足到点(x_c,y_c)的距离，L₂是垂足到(x₀,y₀)的距离，L₃是垂足到(x_r,y_r)的距离；

(4)当转换以后点(x_c′,y_c′)在Y轴正半轴，有(x₀′,y₀′)=(-L₂,L₃)，(x_c′,y_c′)=(L₁,L₃)；否则(x₀′,y₀′)=(-L₂,-L₃)，(x_c′,y_c′)=(L₁,-L₃)。

定义点(x₀,y₀)到点(x_c,y_c)的向量为 点(x₀′,y₀′)到点(x_c′,y_c′)的向量为向量和间夹角α为如果障碍的高度z_r≠0，则旋转平移矩阵T为

$T=\left[\begin{array}{cccc}\cos \mathrm{\α}& \sin \mathrm{\α}& 0& 0\\ -\sin \mathrm{\α}& \cos \mathrm{\α}& 0& 0\\ 0& 0& 1& 0\\ x_r& y_r& z_r& 1\end{array}\right]$

如果障碍的高度z_r＝0，则旋转平移矩阵T的形式不变，其中的第四行第三列的项Z_r为零。设障碍在原点O时的流线数据为矩阵A，其中矩阵A的每一行代表一个航路点。由于每条流线的长度可能不同，具有的航路点数目也不同，因此矩阵A的行是不固定的。但是矩阵A的列数是固定的，为4列。第一列为航路点的X轴方向坐标，第二列为航路点的Y轴方向坐标，第三列为航路点的Z轴方向坐标，第四列为单位1，则任意位置的流线矩阵A′为A′＝A×T。

步骤四：两个障碍时，分别计算两个障碍单独避障的各自流线，判断两障碍单独避障的流线是否满足部分重合条件，计算两个障碍的总流线。

(1)当两个障碍单独避障的流线能够满足部分重合时，此时两个障碍之间的两条流线重合直线线段部分长度大于等于即满足优选为如图3所示。其中O_i和O_i+1分别代表第i个和第i+1个障碍的球心，|O_iO_i+1|表示两个障碍之间的直线距离，|S_i′S_i+1′|表示两个障碍之间的两条流线重合直线线段部分长度，找到两条流线中间相重合的部分，进行流线数据合并作为总流线，如图3所示。其中O_i和O_i+1分别代表第i个和第i+1个障碍的球心，他们各自的避障流线分别为S_i与S_i+1，虚线为来流的方向。将中间区域部分L_i作为数据合并的区域，该区域为两个障碍彼此之间的间隔部分，即左侧障碍物的右边缘与右侧障碍物的左边缘之间的平面区域，则合并的过程为：

(1)找到流线S_i和S_i+1在L_i区域内的两个流线部分；

(2)流线S_i和S_i+1在L_i区域内的两个流线部分中，每两个相邻的航路点组成一个航路段，分别计算每个航路段相较于来流方向的斜率；

(3)将流线S_i和S_i+1在L_i区域内的两个流线部分中各航路段的斜率之差的绝对值最小的部分作为合并的位置，并剔除其余的部分。

(2)当两个障碍单独避障的流线不存在部分重合的直线线段部分时，即不符合(1)的条件时，此时两个障碍之间的两条流线不在同一条直线上，如图4。计算与两个障碍同时内切的虚拟障碍，对该虚拟障碍利用步骤一～步骤三的方法计算流线，该流线为两个障碍的总流线，该流线具有步骤三中矩阵A的形式，为n行4列的矩阵，n根据流线航路点个数的不同而不同。

O_i和O_i+1分别代表第i个和第i+1个障碍的球心，各自的避障流线分别为S_i与S_i+1，R_i和R_i+1分别表示O_i和O_i+1的半径。在线段O_iO_i+1上可以找到一点O_new，以该点为圆心，可以做出与圆O_i和O_i+1同时相内切的圆O_new，在航路规划时我们可以将圆O_new作为虚拟障碍进行躲避，在三维情况下这种方法仍然是可行的，图4是二维情况的简化示意图。其中O_new必须满足下式条件

|O_newO_i|+R_i＝|O_newO_i+1|+R_i+1

步骤五：针对所有障碍，将每两个相邻障碍作为一组，分别计算每个障碍各自的流线，再采用步骤四计算总流线的方法，分别计算每一组中两个障碍的总流线，得到每一组的总流线；然后再采用步骤四中的计算总流线的方法，计算相邻两组总流线再合并后的总流线。采用步骤四中计算总流线的方法，不断合并流线，直至最后产生一条流线作为最终的三维流线。

步骤六：步骤一中计算的流线是一种平行来流。因此在三维直角坐标系下，在YZ平面很小的坐标范围内改变来流的起始位置(来流方向不改变)，经过步骤一～步骤五的方法，获得多条三维流线。

步骤七：将无人机约束表示为流线的航路点和航路段约束，满足无人机约束的流线作为无人机三维可飞航路。

分别将无人机最大爬升高度(H)约束、最大爬升角(α)约束、最大水平转弯角(Φ)约束转化为对流线的约束如下：

(1)最大爬升高度(H)约束

提取出每条流线中各点的高度(即空间直角坐标系中Z轴的坐标)，将其与最大爬升高度(H)进行比较，所有的Z轴的坐标z均满足z≤H，则该条流线满足最大爬升高度约束。

(2)最大爬升角(α)约束

计每条流线中相邻两点构成的航路段与空间直角坐标系的水平面(XY平面)的夹角θ_z，将其与最大爬升角(α)进行比较，流线中所有θ_z均满足θ_z≤α，则该条流线满足最大爬升角约束。

(3)最大水平转弯角(Φ)约束

计算每条流线中相邻两航段在空间直角坐标系的水平面(XY平面)内投影的夹角φ_xy，将其与最大转弯角(Φ)进行比较，流线中所有φ_xy均满足φ_xy≤Φ，则该条流线满足最大水平转弯角约束。

对步骤六中所有三维流线按照上述三种约束进行筛选，满足(1)~(3)条所有约束条件的流线作为无人机的三维可飞航路。

Claims (2)

1.一种基于流体力学的无人机三维航路生成方法，其特征在于：具体包括以下几个步骤：

步骤一：单个障碍时，根据定常忽略粘性的不可压缩理想流体，沿空间直角坐标系X轴来流方向作用于位于原点的圆球障碍绕流问题的解析解，在该直角坐标系下求得流速沿X、Y和Z轴的分量；

定常忽略粘性的不可压缩理想流体沿X轴来流方向作用于位于原点的圆球障碍绕流问题，由于流体的运动存在对称性，描述为平面二维极坐标系下的扰流问题，对于空间中任意一点P，极坐标系建立在点P与X轴所构成的平面上，空间直角坐标系的坐标原点与平面二维极坐标系的极点重合，为点O，空间直角坐标系的X轴与平面二维极坐标系的极轴X′轴重合，在该极坐标系下流体的速度势为

其中OP与极轴之间的夹角为θ角，点P与点O的距离为|OP|＝r，V_∞为流体的流速，a为圆球障碍的半径，在空间直角坐标系下，P的坐标为(x,y,z)，用u、v和w代表三维情况下流速沿空间直角坐标系X、Y和Z轴的分量，则

步骤二：单个障碍时，采用解微分方程组方法，对流速沿X、Y和Z轴的分量所构成的速度势微分方程组进行求解计算，得到沿X轴来流方向的定常忽略粘性的不可压缩理想流体作用于位于原点的圆球障碍绕流流线；

步骤三：单个障碍时，利用旋转平移矩阵，使计算的流线适于空间直角坐标系下XY平面内任何来流方向，并且使多个障碍在空间直角坐标系下XY平面内任意分布，在直角坐标系的XY平面，流线的起点为(x₀,y₀)，终点为(x_c,y_c)，球心坐标为(x_r,y_r)，旋转平移矩阵的计算过程包括以下几个步骤：

(1)计算通过点(x₀,y₀)与点(x_c,y_c)的直线方程y＝k₁x+b₁；k₁为该直线的斜率，b₁为该直线在Y轴上的截距；

(2)计算过点(x_r,y_r)并垂直于直线y＝k₁x+b₁的直线方程；

(3)求得垂足的坐标，并计算L₁，L₂与L₃，其中L₁是垂足到点(x_c,y_c)的距离，L₂是垂足到(x₀,y₀)的距离，L₃是垂足到(x_r,y_r)的距离；

(4)当转换以后点(x′_c,y′_c)在Y轴正半轴，有(x′₀,y′₀)＝(-L₂,L₃)，(x′_c,y′_c)＝(L₁,L₃)；否则(x′₀,y′₀)＝(-L₂,-L₃)，(x′_c,y′_c)＝(L₁,-L₃)；

定义点(x₀,y₀)到点(x_c,y_c)的向量为点(x′₀,y′₀)到点(x′_c,y′_c)的向量为向量和间夹角α为如果障碍的高度z_r≠0，则旋转平移矩阵T为

$T=\left[\begin{array}{cccc}\cos \mathrm{\α}& \sin \mathrm{\α}& 0& 0\\ -\sin \mathrm{\α}& \cos \mathrm{\α}& 0& 0\\ 0& 0& 1& 0\\ x_r& y_r& z_r& 1\end{array}\right]$

z_r为障碍的高度，设障碍在原点O时的流线数据为矩阵A，其中矩阵A的每一行代表一个航路点，矩阵A的列数为4列，第一列为航路点的X轴方向坐标，第二列为航路点的Y轴方向坐标，第三列为航路点的Z轴方向坐标，第四列为单位1，任意位置的流线矩阵A′为A′＝A×T；

步骤四：两个障碍时，分别计算两个障碍单独避障的各自流线，判断两障碍单独避障的流线是否满足部分重合条件，计算两个障碍的总流线；

(1)当两个障碍单独避障的流线能够满足部分重合时，此时两个障碍之间的两条流线重合直线线段部分长度大于等于其中O_i和O_i+1分别代表第i个和第i+1个障碍的球心，|O_iO_i+1|表示两个障碍之间的直线距离，|S_i′S_i+1′|表示两个障碍之间的两条流线重合直线线段部分长度，找到两条流线中间相重合的部分，进行流线数据合并作为总流线，合并的过程具体为：

(1)找到流线S_i和S_i+1在L_i区域内的两个流线部分；

(2)流线S_i和S_i+1在L_i区域内的两个流线部分中，每两个相邻的航路点组成一个航路段，分别计算每个航路段相较于来流方向的斜率；

(3)将流线S_i和S_i+1在L_i区域内的两个流线部分中，各航路段的斜率之差的绝对值最小的部分作为合并的位置，并剔除其余的部分；

(2)当两个障碍单独避障的流线不存在部分重合的直线线段部分时，此时两个障碍之间的两条流线重合直线线段部分长度小于计算与两个障碍同时内切的虚拟障碍，对该虚拟障碍利用步骤一～步骤三的方法计算流线，该流线为两个障碍的总流线，该流线具有步骤三中矩阵A的形式，为n行4列的矩阵，n根据流线航路点个数的不同而不同；

O_i和O_i+1分别代表第i个和第i+1个障碍的球心，各自的避障流线分别为S_i与S_i+1，R_i和R_i+1分别表示O_i和O_i+1的半径，在线段O_iO_i+1上存在点O_new，以该点为圆心，做出与圆O_i和O_i+1同时相内切的圆O_new，将圆O_new作为虚拟障碍进行躲避，O_new满足：

|O_newO_i|+R_i＝|O_newO_i+1|+R_i+1

步骤五：针对所有障碍，将每两个相邻障碍作为一组，分别计算每个障碍各自的流线，再采用步骤四计算总流线的方法，分别计算每一组中两个障碍的总流线，得到每一组的总流线；然后再采用步骤四中的计算总流线的方法，计算相邻两组总流线再合并后的总流线，顺次采用步骤四中计算总流线的方法，不断合并流线，直至最后产生一条流线最为最终三维流线；

步骤六：步骤一中计算的流线是一种平行来流，在三维直角坐标系下，在YZ平面的坐标范围内改变来流的起始位置，经过步骤一～步骤五的方法，获得多条三维流线；

步骤七：将无人机约束表示为流线的航路点和航路段约束，满足无人机约束的流线作为无人机三维可飞航路；

分别将无人机最大爬升高度H约束、最大爬升角α约束、最大水平转弯角Φ约束转化为对流线的约束如下：

(1)最大爬升高度H约束

提取出每条流线中各点的高度，将其与最大爬升高度进行比较，所有的Z轴的坐标z均满足z≤H，则该条流线满足最大爬升高度约束；

(2)最大爬升角α约束

计算每条流线中相邻两点构成的航路段与空间直角坐标系的水平面的夹角θ_z，将其与最大爬升角进行比较，流线中所有θ_z均满足θ_z≤α，则该条流线满足最大爬升角约束；

(3)最大水平转弯角Φ约束

计算每条流线中相邻两航段在空间直角坐标系的水平面内投影的夹角φ_xy，将其与最大水平转弯角进行比较，流线中所有φ_xy均满足φ_xy≤Φ，则该条流线满足最大水平转弯角约束；

对步骤六中所有三维流线按照上述三种约束进行筛选，满足(1)～(3)条所有约束条件的流线作为无人机的三维可飞航路。

2.根据权利要求1所述的一种基于流体力学的无人机三维航路生成方法，其特征在于：所述的步骤二中解微分方程组方法包括四阶或五阶的定步长或变步长Runge-Kutta方法。