一种基于地形通视的无人机任务载荷侦察覆盖率计算方法

1.本发明涉及无人机技术领域，具体涉及一种基于高程地形通视分析的无人机任务载荷侦察覆盖率计算方法。


背景技术：

2.无人机任务规划包括航路规划、任务载荷规划、链路规划以及仿真推演等，其规划的结果将作为无人机系统执行任务过程的主要依据。一般情况下，任务规划的流程是基于输入的任务区和无人机起降机场，规划出无人机的航路；在该基础上进一步进行任务载荷和链路的规划；再通过仿真推演对规划结果进行推演验证，以进一步优化规划结果。其中，航路规划与任务载荷规划相互之间具备依赖关系：一方面，任务载荷规划需结合已规划的任务区域航路信息和需完成的任务进行任务载荷在各航路段的工作模式和参数的设置，另一方面，任务载荷规划的结果将作为任务航线修改的输入并且影响无人机任务执行的效能。无人机任务执行完成度的判断标准主要是对既定侦察区域的覆盖率计算。传统的任务载荷侦察覆盖率计算，主要方式是将任务区和待侦察区等效为二维平面，然后计算被侦察到的区域与全部侦察区域的面积比。使用该方式一方面计算精确度较低；另一方面在高山、丘陵等有地形遮挡情况的区域侦察结果与实际结果可能有较大偏差，对无人机执行该类型地形下的任务带来了挑战。
3.目前，对于任务区阶段无人机任务载荷侦察覆盖率计算典型方法主要包括三个方面：
4.1.基于无人机飞行性能计算无人机转弯半径，确定任务区搜索航线的最小间距，同时确认任务区的进入点和离开点以及航线的类型；
5.2.基于无人机高度、任务载荷姿态以及任务载荷的工作探测范围确定单次扫描覆盖范围；
6.3.计算扫描覆盖区域面积与任务区面积比例，从而得出侦察覆盖率。
7.以上的方式，主要结合了无人机的飞行性能、任务载荷探测性能、威胁区域以及任务区域的具体状态进行综合的分析。该方式虽然也可根据飞行器的任务特点规划螺旋式、扫雪式、8字式和光栅式等航路类型，但该方式忽略了不同地形对任务载荷侦察地面目标的遮挡，因此某些情况下的计算可能会遗漏掉需探测的目标。对于平原地形来说基本可满足使用需求，但对于绝大多数非平原地形，使用该方式将存在侦察区域的探测盲区，无法满足基本使用需求，因此需要拓展其他手段进行提升。利用基于地形高程进行通视计算任务载荷侦察覆盖率的方式可适应于大多数地形，为侦察覆盖率的精确计算提供基础条件。
8.基于该思路进行任务载荷侦察覆盖率计算，基于三个方面：
9.1.基于无人机性能计算方面：基于无人机的性能、燃油和飞行速度等得出任务区航线最大长度，同时无人机飞行转弯半径确定在任务区航线最小间距，以此确定任务区的具体航线；
10.2.基于任务载荷性能计算方面：基于无人机高度、任务载荷姿态以及任务载荷侦
察范围，结合地球曲率计算任务载荷单次扫描的收容宽度；
11.3.基于地形高程计算覆盖率方面：利用地形高程数据的通视分析，计算每一次扫描范围内的通视和非通视区域，计算完整航线的侦察范围，并进一步计算侦察覆盖率。
12.在上述技术手段中，无人机性能计算是基础，任务区航线类型选取确定了航线的基本方向，任务载荷性能计算一方面确定单次扫描的理论覆盖范围，另一方面将方向输入到航路进行航路的合理调整；而侦察覆盖率的计算则是规划结果是否满足要求或是否最优的直观判断。
13.在目前情况下，任务区扫描覆盖率的计算主要考虑无人机飞行高度、任务载荷姿态、任务载荷探测距离和角度范围。进行具体计算时多使用将地面等效为平面的方式加以处理，该方式能一定程度上计算出在任务区的侦察覆盖率。但使用该技术主要存在如下的问题：
14.1.未考虑地形遮挡的情况，对于山区、丘陵等非平原地形的任务区，将存在大量的探测盲区；
15.2.将地面等效为平面，则地面威胁区建模为平面，未考虑地球区域的影响，对无人机的可飞区域预估过于保守，使得在多种情况下由于可飞区域不具备条件，无法得出满足任务要求的规划结果。


技术实现要素：

16.针对现有技术存在的问题，本发明提供一种基于地形通视的无人机任务载荷侦察覆盖率计算方法，具体如下：
17.通过引入地形高程数据而进行的通视分析，进行任务区航线条件下任务载荷侦察覆盖范围和覆盖率的准确计算，本技术适用如下条件：
18.1)无人机执行侦察任务时为平飞状态，即其俯仰、滚转均为0
°
的状态；
19.2)任务区地形为高山地形或者平原地形皆可，但有该区域的高程数据；
20.3)不考虑地面建筑、植被等其他地表覆盖物的情况；
21.4)不考虑光电吊舱eo等光学任务载荷在不同天气情况下被云层、雨雪遮挡的情况；
22.5)不考虑复杂电磁环境对于雷达等任务载荷产生干扰的情况；
23.(1)数据准备
24.该阶段获取基于顶层任务目标所形成的任务区数据、任务载荷参数和状态数据、航路点数据；将上述数据进行二维标定以进行下一步计算的数据准备，具体包括获得无人机航路点、任务区域、任务载荷状态数据、任务载荷性能参数；在任务规划阶段，已知的数据如下：
25.a)无人机位置p在wgs84坐标系下的经纬海拔(b1，l1，h)：
26.b)无人机任务载荷视轴中心的俯仰角θ，该值由操作员认为设定参数，值域为0
°
到90
°
，任务载荷水平为0
°
，竖直为90
°
；
27.c)无人机航向角ψ，由计算时实时获取；
28.d)无人机偏流角γ，由计算时实时获取；
29.e)无人机任务载荷视场角水平大小vh和竖直大小vv，该两值为任务载荷的参数，根
据操作员所选择的工作模式和设定的参数能够方便计算得出；无人机任务载荷当前状态下近界俯仰角
30.f)无人机任务载荷当前状态下远界俯仰角
31.g)无人机任务载荷极限俯仰角θ
max
，该值为任务载荷的固有属性，查询任务载荷手册确认，值域为0
°
到90
°
；
32.h)无人机任务载荷最大作用距离为l
max
，该值为任务载荷的固有属性，查询任务载荷手册确认；
33.i)地球近似为球体，其半径r＝6371000，单位为米；
34.(2)扫描边界计算
35.该阶段基于航路点数据、任务载荷参数和状态数据，计算出扫描边界点；计算过程考虑地球曲率的影响，但先不考虑地形高程的影响；由于侦察覆盖率计算应用于任务预规划阶段，对任务载荷的使用较为理想的为无人机平飞阶段，在该阶段无人机姿态中的俯仰和滚转为0
°
，且任务载荷预置模式为左侧视或右侧视；将任务载荷对地信号作用范围为扇形作为基本假设，获得一次采样的侦察覆盖范围；令
36.p点为无人机当前位置；
37.po为其在地面投影；
38.a点为任务载荷视场左边界与任务载荷探测范围近界在地面投影的交点；
39.b点为任务载荷视场左边界与任务载荷探测范围远界在地面投影的交点；
40.c点为任务载荷视场右边界与任务载荷探测范围远界在地面投影的交点；
41.d点为任务载荷视场右边界与任务载荷探测范围近界在地面投影的交点；
42.ad弧线为无人机任务载荷当前位置扫面范围的近界；
43.bc弧线为无人机任务载荷当前位置扫面范围的远界；
44.m点为任务载荷视轴中心线与任务载荷探测范围近界在地面投影的交点；
45.n点为任务载荷视轴中心线与任务载荷探测范围远界在地面投影的交点；
46.mn连线为任务载荷视场中心线，与无人机飞行方向垂垂直，也就是与飞行方向夹角90度；
47.a、b、c、d、m、n点的经纬度由无人机高度、任务载荷最大作用斜距、任务载荷最大俯仰角这些参数设置后得到；具体计算方式如下：
48.①
计算m点的经纬度
49.a)圆心角计算
50.为计算得更加精确，需将地球近似为球体，并考虑地球曲率，已知无人机位置p、地心o、p点到球面的投影点po、任务载荷近界俯仰角θ1＝∠opm、地球半径r、飞行高度h、任务载荷的最大作用距离l
max
、任务载荷的极限俯仰角θ
max
、最大俯仰角时无人机位置发出的射线与球面的交点m’，无人机位置发出的射线与地球面的切点n’；根据余弦定理，优先计算第一圆心角θ3＝∠mop、第二圆心角具体如下：
51.i.计算任务载荷任务载荷视距范围内最大作用斜距
52.53.ii.计算最小俯仰角∠opn
′
[0054][0055]
iii.计算任务载荷所能达到的最短斜距pm’，根据余弦定理带入参数，如式(3)：
[0056]
pm
′2+(r+h)
2-2*pm
′
*(r+h)*cos(90
°‑
θ
max
)＝r2ꢀꢀ
(3)
[0057]
计算得到pm
′
；
[0058]
iv.校验任务载荷近界俯仰角θ1[0059]
为保证计算的有效性，即所计算的最大俯仰角θ1在合理值范围内，需进行校验；如果θ1小于俯仰角最小值∠opn
′
，则θ1＝∠opn
′
，如果大于任务载荷物理最大俯仰值θ
max
，通过顶层输入或预设，则θ1＝θ
max
＝∠opm
′
，反之则值不变；通过该方式控制任务载荷在有效和可用范围内；
[0060]
v.校验最大作用距离l
max
[0061]
如果任务载荷的最大作用距离l
max
大于最大斜距pn
′
，则l
max
＝pn
′
，如果任务载荷的最大作用距离l
max
小于最小斜距pm
′
，则l
max
＝pm
′
，反之l
max
值不变；通过该方式控制任务载荷最大作用距离有效和在可用范围内；
[0062]
vi.计算第一圆心角θ3，pm为预设的最大俯仰角θ
max
对应的斜距；
[0063]
根据下式余弦定理计算得出pm值：
[0064]
pm2+(r+h)
2-2*pm*(r+h)*cos(90
°‑
θ1)＝r2ꢀꢀ
(4)
[0065]
再次使用余弦定理如下式：
[0066]
r2+(r+h)
2-2*r*(r+h)*cos(θ3)＝pm2ꢀꢀ
(5)
[0067]
计算出θ3：
[0068][0069]
vii.同理，由于l
max
为已知值，使用余弦定理直接计算n与p点所对应的第
[0070]
二圆心角
[0071][0072]
viii.比较θ3与的值，如果θ3大于则交换θ3与的值；反之，不做任何
[0073]
改变；
[0074]
b)m点经纬度计算
[0075]
在计算出θ3与的值后，进一步计算m点的经纬度；此时已知无人机位置在球面正投影点po经度、纬度、航向角以及偏流角(l1，b1，ψ，γ)，po点到m点的圆心角＝θ3，q为北极点，po、m、q构成球面三角形，b1、b2、l1、l2分别是po、m两点的纬度和经度坐标；利用球面三角形计算公式计算出m点的经纬度坐标；
[0076]
为便于计算，需对球面三角形进行简化；
[0077]
其中：q、p0、m所对应的球面三角形弧长分别是a、b、c，a
′
、b
′
、c
′
是球面上的角；∠a
′
是∠p0qm的简写，∠b
′
是∠qp0m的简写，∠c
′
是∠qmp0的简写，a是过球心圆上的弧长p0m，b是球心圆上的弧长qm，c是过球心圆上的弧长p0q；
[0078]
在球面三角形p0mq中，已知p0点的纬度b1，经度l1，航向角ψ，偏流角γ，圆心角∠p0om，需计算目标点m的纬度b2和经度l2；
[0079]
i.计算m点经度l2[0080]
已知∠a
′
＝l
2-l1，c＝∠p0oq＝90-b1，b＝∠moq＝90-b2，a＝∠p0om；a、b、c是球面三角形弧长，单位长度半径的球面，其弧长与角度相等，因此用角度表示；
[0081]
因为pm垂直于航线方向，所以∠b
′
＝90
°
+ψ+γ；
[0082]
根据球面正弦公式：代入已知条件，得：
[0083]
sin a sin b
′
＝sin(90
°‑
b2)sin a
′ꢀꢀ
(8)
[0084]
根据球面三角形余切公式：cot a sin c＝cot a
′
sin b
′
+cos b
′
cos c，代入已知条件，得：
[0085]
cot a
·
sin(90
°‑
b1)＝cot(a
′
)
·
sin b
′
+cos b
′
·
cos(90
°‑
b1)
ꢀꢀ
(9)
[0086]
化简后得：
[0087][0088]
即能够计算出经度l2；
[0089]
ii.计算m点纬度b2[0090]
在计算m点纬度b2之前需先校验b
′
，
[0091]
》如果b
′
值为零，则说明目标点m与p点在同一纬度，b2＝b1+∠p0om；
[0092]
》如果b
′
值不为零，则需要通过判断象限来确定b2的值；计算公式如下：
[0093]
cos b2＝cos a
·
sin b1+sin a
·
cos b1ꢀꢀ
(11)
[0094]
由式(7)转换得：
[0095][0096]
判断l
2-l1的取值范围，如果l
2-l1大于零则sin b2值不变，反之则将sin b2的值取个负号；
[0097]
判断cos b2的取值范围，如果cos b2大于等于零则b在第一、第四象限，反之则b在第二、第三象限；在此基础上，通过判断sin b2取值范围确定b所在象限；
[0098]
》当cos b2的值大于零且sin b2大于等于零，则b在第一象限
[0099]
b2＝90
°‑
csc(sin b2)
ꢀꢀ
(13)
[0100]
》当cos b2的值大于等于零且sin b2小于零，则b在第四象限
[0101]
b2＝90
°
+csc(sin b2)
ꢀꢀ
(14)
[0102]
》当cos b2的值小于零且sin b2大于等于零，则b在第二象限
[0103]
b2＝90
°
+csc(sin b2)
ꢀꢀ
(15)
[0104]
》当cos b2的值小于零且sin b2小于零，则b在第三象限
[0105]
b2＝90
°‑
csc(sin b2)
ꢀꢀ
(16)
[0106]
由此，计算出m点的经纬度；重复上述计算步骤，则能够计算出n、a、b、c、d的经纬度；
[0107]
(3)扫描区域处理；
[0108]
扫描区采用栅格化方法处理；
[0109]
1)地形高程获取
[0110]
在获得考虑地球曲率影响下的各点经纬度后，a、b、c、d、m、n各点海拔高度，由其经
纬度带入高程数据后即能够确定；
[0111]
2)扫描范围栅格化
[0112]
根据a、b、c、d四点所组成的区域进行栅格化，ad方向m等分，ab方向n等分，得到m*n个子四边形区域，存储这些子四边形区域的中心点，为用作通视采样点做准备；
[0113]
栅格化所用参数如下：
[0114]
①
m：方位向离散化等分系数，m为整数，最小值为1，其值越小，细分颗粒度越粗，结果偏差越大，但计算效率越高，由规划人员输入设定；
[0115]
②
n：距离向离散化系数，设定方法同方位向，m和n不一定相等；
[0116]
③
m*n：一次任务载荷侦察覆盖范围的离散点数；
[0117]
④
d：方位向最小细分距离，d＝l
ad
∠m，其中l
ad
是a、d水平投影点的连线；
[0118]
(4)栅格小区域通视计算
[0119]
利用地形高程数据、无人机当前位置数据以及栅格化后的数据进行各栅格小区域通视计算，即能够得出各小区域的通视结果；通视分析时根据当前无人机所在位置和目标点的位置，结合高程地形进行两点之间是否可视的计算和判断；无人机飞行剖面纵向切割地形获得的地形剖面，无人机实时位置与目标点连线，如果与地形剖面有交叉，则表示无人机在当前点与目标不通视，反之则通视；
[0120]
由于无人机任务载荷是否通视地面目标与地面目标是否通视无人机判断依据相同；为便于计算，将地面目标作为视点-v点，无人机任务载荷的探测距离作为作用范围，也就是换算为地距，其边界位置作为目标点-t点，先选择某一个方向；然后从视点到目标点逐点分析交点；其中的f(xi，yi)、g(x
i+1
，y
i+1
)
…
这些点均为视线vt在(x，y)平面上的投影同方形格网单元边的交点，视线vt的斜率α由下式算出：
[0121][0122]
其中：
[0123]zt
为t(x
t
，y
t
)点的高度；
[0124]
zv为v(xv，yv)点的高度；
[0125]
为t点到v点水平距离；
[0126]
以基准水平面向东为x轴，以v点铅锤面向上为y轴构建坐标系，视点v与各交点连线与水平面夹角βi斜率如下：
[0127][0128]
其中：
[0129]
zi为f(xi，yi)点的高度；
[0130]
zv为v(xv，yv)点的高度；
[0131]
为f点到v点水平距离；
[0132]
通视性判断通过比较tan α与tan βi的值完成：
[0133]
如果tan βi大于tan α，则表示不能通视，计算结束；
[0134]
如果tan βi不大于tan α，则f点通视，计算下一个点g(x
i+1
，y
i+1
)；如果可一直推进到目标点t，则视点与目标点之间可以通视；
[0135]
在根据公式(17)和公式(18)以及上述通视判断计算方式完成一次基于航路点数据、航路点等分数据进行栅格小区域通视计算为一次采样计算，在完成一次采样计算并得出结果后，按照栅格的参数按照方位向和距离向以此类推迭代计算，即能够得出一次采样覆盖范围内的通视和非通视栅格；
[0136]
(5)合并计算
[0137]
基于航路点数据、航路点等分数据进行各通视计算点的栅格小区域通视计算，并将通视结果合并，获得通视和非通视栅格的图；此外，该步骤还同步计算扫描区域的合并结果；具体地，在一次采样计算完成后，获得通视和非通视栅格的图，以此类推，按照航路点数据和切分的间隔计算所有航线点的扫描区域，获得通视和非通视栅格的图，并进行合并计算；设定所有航线点扫描区域为u1，不通视区域为u2，扫描区与任务区交集为u3，不通视区域与任务区交集为u4；合并计算的具体算法步骤为：
[0138]
①
令航点1、航点2表示无人机任务航线上的任意2个航点，其他航点的计算以此类推；按照采样间隔，在航点1和航点2中间平滑出所有插值航点，获取所有航点生成的扫描区域数据，扫描区域是多个扇形，每个扇形为航点1和航点2之间某个航点生成的扫描区域数据；
[0139]
②
将多边形区域边界点按照逆时针方向排列存储；
[0140]
③
将这些扇形区域合并，获得由扫描区域边界点形成的多边形，遍历所有扫描区域形成的多边形区域，将当前被合并的区域与多边形区域中各多边形逐个求并集；
[0141]
④
基于遍历和逐个求并集的结果，计算出并集后的多边形区域边界点数据；
[0142]
(6)计算侦察覆盖率
[0143]
设顺着航向方向递推，根据上述步骤
①‑④
计算出的通视区域与任务区交集u3和任务区u4，即能够获得整条航向的侦察覆盖范围；
[0144]
令：
[0145]
p1是第一次计算时无人机所在点，也是任务载荷所在点；
[0146]
a1、b1、c1、d1分别是无人机在p1点时，无人机任务载荷当前位置扫面范围的边界点；a1d1弧线为近界，a1在左边界，d1在右边界，b1c1弧线为远界，b1在左边界，c1在右边界；m1和n1分别为a1d1弧线以及b1c1弧线的中心点；
[0147]
同理，p2点为下一次计算无人机所在点，也是任务载荷所在点；其他点的信息描述和定义与p1点对应的点一致；
[0148]
令：
[0149]
u1为所有航线点扫描区域；
[0150]
u2为所有的不通视区域；
[0151]
u3为扫描区与任务区交集；
[0152]
u4为不通视区域与任务区交集；
[0153]
对该二维图形进行算和处理，计算出侦察覆盖率。
[0154]
在本发明的一个实施例中，在步骤(4)中采用dyntacs通视性算法。
[0155]
在本发明的一个具体实施例中，在步骤(5)中，采样间隔为方位向最小细分距离d的整数倍。
[0156]
在本发明的另一个具体实施例中，在步骤(6)中，递推采样间隔设置为方位向最小
细分距离d的整数倍，d最小值为1，最大值为m。
[0157]
在本发明的另一个实施例中，在步骤(6)中，对该二维图形进行算和处理，具体步骤如下：
[0158]
①
计算所有任务区域集合；
[0159]
②
合并所有航点扫描的通视区域；
[0160]
③
在任务区集合中取出单个任务区；
[0161]
④
将任务区中的边界点数据逆序排列；
[0162]
⑤
计算区域数据与通视区域的交集，得到区域边界点；
[0163]
⑥
通过交集区域边界点，计算出交集面积；
[0164]
⑦
用交集面积除以任务区面积，即得出此任务区的扫描覆盖率；
[0165]
⑧
重复
③‑⑦
即能够计算出所有的任务区的覆盖率。
[0166]
本发明针对于除平原地形外的高山、丘陵等多种复杂地形，使用统一的侦察覆盖率计算方式，在不改变任务规划系统原系统的软硬件架构下，通过引入高程地形数据库并进行通视分析的方式计算侦察覆盖范围和侦察覆盖率，提出一种用于实现符合多种地形地貌下更接近真实任务目标要求的侦察覆盖率计算方法，为任务规划中任务区阶段航路的规划提供一个新的思路。
[0167]
本发明的优点如下：
[0168]
(1)实现航线规划对于非平面地形的支持。本发明通过使用基于地形高程模型的通视计算方式，将原本二维的任务区三维化，充分考虑到不同地形下任务载荷侦察覆盖范围的地形遮挡因素，将航路规划的适用从原来的平原地形提升到对多种地形地貌的支持，提升了航线规划的可用范围。
[0169]
(2)实现既定航线侦察覆盖率的定量描述。在基于地形进行通视分析的基础上，通过标绘和统计，计算出既定航线可完成侦察范围的覆盖率数据计算。通过该方式，可作为规划航路优化的量化依据，能够避免操作员人为和主观判断的误差。
附图说明
[0170]
图1示出侦察覆盖范围示意图；
[0171]
图2示出侦察覆盖范围水平投影示意图；
[0172]
图3示出侦察覆盖范围垂直投影示意图；
[0173]
图4示出侦察覆盖范围中心m、n点计算示意图(考虑地球曲率)；
[0174]
图5示出球面三角形示意图；
[0175]
图6示出球面三角形简化示意图；
[0176]
图7示出栅格化示意图；
[0177]
图8示出通视计算示意图；
[0178]
图9示出通视算法示意图；
[0179]
图10示出通视计算结果示意图；
[0180]
图11示出扫描区域合并示意图；
[0181]
图12示出合并计算示意图；
[0182]
图13示出非通视区域合并计算示意图。
具体实施方式
[0183]
针对于越来越精细化的无人机任务规划需求，在充分考虑常用航路规划约束条件、常用规划算法的基础上。通过针对任务区预规划航路的任务区覆盖率计算，能够有效提升对规划航路的有效性评估，为航路规划任务阶段航线的合理性优化提供依据。在进行任务区覆盖率计算时，使用等效地面模型进行具体覆盖范围的计算，该方式虽然也能计算出大致的覆盖范围结果，但没能考虑到真实地形对于覆盖范围的影响。为此，本发明提出基于高程地形通视分析进行侦察覆盖率的计算。其总体思路为：对于任务区域，除获取其边界经纬度属性外，还通过高程数据获取范围内的地形起伏情况，将二维的地图三维化，再结合无人机位置和任务载荷工作状态进行三维空间内的通视分析，以得出覆盖范围和覆盖盲区，由此提升覆盖范围和覆盖率计算的真实度和准确度。该计算方法过程阐述如下：
[0184]
通过引入地形高程数据而进行的通视分析，进行任务区航线条件下任务载荷侦察覆盖范围和覆盖率的准确计算，本技术可适用如下条件(基本假设)：
[0185]
6)无人机执行侦察任务时为平飞状态，即其俯仰、滚转均为0
°
的状态；
[0186]
7)任务区地形为高山地形或者平原地形皆可，但有该区域的高程数据；
[0187]
8)不考虑地面建筑、植被等其他地表覆盖物的情况；
[0188]
9)不考虑光电吊舱(eo)等光学任务载荷在不同天气情况下被云层、雨雪遮挡的情况；
[0189]
10)不考虑复杂电磁环境对于雷达等任务载荷产生干扰的情况。
[0190]
(1)数据准备
[0191]
该阶段获取基于顶层任务目标所形成的任务区数据、任务载荷参数和状态数据、航路点数据等。将上述数据进行二维标定以进行下一步计算的数据准备，具体包括获得无人机航路点、任务区域、任务载荷状态数据、任务载荷性能参数等。对上述数据进行二维标定的方法为本领域技术人员熟知，不再累述。在任务规划阶段，已知的数据如下：
[0192]
j)无人机位置p在wgs84坐标系下的经纬海拔(b1，l1，h)：
[0193]
k)无人机任务载荷视轴中心的俯仰角θ，该值由操作员认为设定参数，值域为0
°
到90
°
，任务载荷水平为0
°
，竖直为90
°
；
[0194]
l)无人机航向角ψ，由计算时实时获取；
[0195]
m)无人机偏流角γ，由计算时实时获取；
[0196]
n)无人机任务载荷视场角水平大小vh和竖直大小vv，该两值为任务载荷的参数，根据操作员所选择的工作模式和设定的参数可方便计算得出；无人机任务载荷当前状态下近界俯仰角
[0197]
o)无人机任务载荷当前状态下远界俯仰角
[0198]
p)无人机任务载荷极限俯仰角θ
max
，该值为任务载荷的固有属性，可查询任务载荷手册确认，值域为0
°
到90
°
；
[0199]
q)无人机任务载荷最大作用距离为l
max
，该值为任务载荷的固有属性，可查询任务载荷手册确认；
[0200]
r)地球近似为球体，其半径r＝6371000，单位为米。
[0201]
(2)扫描边界计算
[0202]
该阶段基于航路点数据、任务载荷参数和状态数据，计算出扫描边界点。计算过程考虑地球曲率的影响，但先不考虑地形高程的影响。由于本发明侦察覆盖率计算应用于任务预规划阶段，对任务载荷的使用较为理想的为无人机平飞阶段(无人机姿态中的俯仰和滚转为0
°
)，且任务载荷预置模式为左侧视或右侧视。将任务载荷对地信号作用范围为扇形作为基本假设，可获得一次采样的侦察覆盖范围示意如附图1所示，水平投影绘制如附图2所示，垂直视角绘制如附图3所示。其中，
[0203]
p点为无人机当前位置；
[0204]
po为其在地面投影；
[0205]
a点为任务载荷视场左边界与任务载荷探测范围近界在地面投影的交点；
[0206]
b点为任务载荷视场左边界与任务载荷探测范围远界在地面投影的交点；
[0207]
c点为任务载荷视场右边界与任务载荷探测范围远界在地面投影的交点；
[0208]
d点为任务载荷视场右边界与任务载荷探测范围近界在地面投影的交点。
[0209]
ad弧线为无人机任务载荷当前位置扫面范围的近界；
[0210]
bc弧线为无人机任务载荷当前位置扫面范围的远界；
[0211]
m点为任务载荷视轴中心线与任务载荷探测范围近界在地面投影的交点；
[0212]
n点为任务载荷视轴中心线与任务载荷探测范围远界在地面投影的交点；
[0213]
mn连线为任务载荷视场中心线，与无人机飞行方向垂垂直(与飞行方向夹角90度)；
[0214]
a、b、c、d、m、n点的经纬度由无人机高度、任务载荷最大作用斜距、任务载荷最大俯仰角等参数设置后可得到。具体计算方式如下：
[0215]
①
计算m点的经纬度
[0216]
b)圆心角计算
[0217]
为计算的更加精确，需将地球近似为球体，并考虑地球曲率，如附图4所示，已知无人机位置p、地心o、p点到球面的投影点po、任务载荷近界俯仰角θ1(∠opm)、地球半径r、飞行高度h、任务载荷的最大作用距离l
max
、任务载荷的极限俯仰角θ
max
、最大俯仰角时无人机位置发出的射线与球面的交点m’，无人机位置发出的射线与地球面的切点n’。根据余弦定理，可优先计算第一圆心角θ3(∠mop)、第二圆心角(nop)。具体如下：
[0218]
i.计算任务载荷任务载荷视距范围内最大作用斜距
[0219][0220]
ii.计算最小俯仰角∠opn
′
[0221][0222]
iii.计算任务载荷所能达到的最短斜距pm’，根据余弦定理带入参数，如式(3)：
[0223]
pm
′2+(r+h)
2-2*pm
′
*(r+h)*cos(90
°‑
θ
max
)＝r2ꢀꢀ
(3)计算可得pm
′
。式中各参数已前面已经描述，故不赘述。
[0224]
iv.校验任务载荷近界俯仰角θ1[0225]
为保证计算的有效性，即所计算的最大俯仰角θ1在合理值范围内，需进行校验。如果θ1小于俯仰角最小值∠opn
′
，则θ1＝∠opn
′
，如果大于任务载荷物理最大俯仰值θ
max
(顶层输入或预设)，则θ1＝θ
max
＝∠opm
′
，反之则值不变。通过该方式控制任务载荷在有效和可用
范围内。
[0226]
v.校验最大作用距离l
max
[0227]
如果任务载荷的最大作用距离l
max
大于最大斜距pn
′
，则l
max
＝pn
′
，如果任务载荷的最大作用距离l
max
小于最小斜距pm
′
，则l
max
＝pm
′
，反之l
max
值不变。通过该方式控制任务载荷最大作用距离有效和在可用范围内。
[0228]
vi.计算第一圆心角θ3，pm为预设的最大俯仰角θ
max
对应的斜距。
[0229]
根据下式余弦定理可以计算得出pm值：
[0230]
pm2+(r+h)
2-2*pm*(r+h)*cos(90
°‑
θ1)＝r2ꢀꢀ
(4)
[0231]
再次使用余弦定理如下式：
[0232]
r2+(r+h)
2-2*r*(r+h)*cos(θ3)＝pm2ꢀꢀ
(5)
[0233]
可计算出θ3：
[0234][0235]
vii.同理，由于l
max
为已知值，可使用余弦定理直接计算n与p点所对应的第二圆心角
[0236][0237]
viii.比较θ3与的值，如果θ3大于则交换θ3与的值；反之，不做任何改变(此步骤的目的为将坐标系统一到正北模式，便于计算)。
[0238]
b)m点经纬度计算
[0239]
在计算出θ3与的值后，进一步计算m点的经纬度，计算示意如附图5所示。此时已知无人机位置在球面正投影点po经度、纬度、航向角以及偏流角(l1，b1，ψ，γ)，po点到m点的圆心角(即上面求出的θ3)，q为北极点，如附图5所示，po、m、q构成球面三角形，b1、b2、l1、l2分别是po、m两点的纬度和经度坐标。利用球面三角形计算公式可计算出m点的经纬度坐标。
[0240]
为便于计算，需进行简化。简化后的球面三角形如附图6所示。
[0241]
其中：q、p0、m所对应的球面三角形弧长分别是a、b、c，a
′
、b
′
、c
′
是球面上的角。∠a
′
是∠p0qm的简写，∠b
′
是∠qp0m的简写，∠c
′
是∠qmp0的简写，a是过球心圆上的弧长p0m，b是球心圆上的弧长qm，c是过球心圆上的弧长p0q。
[0242]
如附图6所示，在球面三角形p0mq中，已知p0点的纬度b1，经度l1，航向角ψ，偏流角γ，圆心角∠p0om，需计算目标点m的纬度b2和经度l2。
[0243]
i.计算m点经度l2[0244]
已知∠a
′
＝l
2-l1，c＝∠p0oq＝90-b1，b＝∠moq＝90-b2，a＝∠p0om。a、b、c是球面三角形弧长，单位长度半径的球面，其弧长与角度相等，因此可以用角度表示。
[0245]
因为pm垂直于航线方向，所以∠b
′
＝90
°
+ψ+γ。
[0246]
根据球面正弦公式：代入已知条件，可得：
[0247]
sin a sin b
′
＝sin(90
°‑
b2)sin a
′ꢀꢀ
(8)
[0248]
根据球面三角形余切公式：cot a sin c＝cot a
′
sin b
′
+cos b
′
cos c，代入已知条件，可得：
[0249]
cot a
·
sin(90
°‑
b1)＝cot(a
′
)
·
sin b
′
+cos b
′
·
cos(90
°‑
b1)
ꢀꢀ
(9)
[0250]
化简后得：
[0251][0252]
即可计算出经度l2。
[0253]
ii.计算m点纬度b2[0254]
在计算m点纬度b2之前需先校验b
′
，
[0255]
》如果b
′
值为零，则说明目标点m与p点在同一纬度，b2＝b1+∠p0om；
[0256]
》如果b
′
值不为零，则需要通过判断象限来确定b2的值。计算公式如下：
[0257]
cos b2＝cos a
·
sin b1+sin a
·
cos b1ꢀꢀ
(11)
[0258]
由式(7)转换可得：
[0259][0260]
判断l
2-l1的取值范围，如果l
2-l1大于零则sin b2值不变，反之则将sin b2的值取个负号。
[0261]
判断cos b2的取值范围，如果cos b2大于等于零则b在第一、第四象限，反之则b在第二、第三象限。在此基础上，通过判断sin b2取值范围确定b所在象限(cos值只能判断出其所在的两个象限，在此基础上再加上sin值，可确定唯一象限)。
[0262]
》当cos b2的值大于零且sin b2大于等于零，则b在第一象限
[0263]
b2＝90
°‑
csc(sin b2)
ꢀꢀ
(13)
[0264]
》当cos b2的值大于等于零且sin b2小于零，则b在第四象限
[0265]
b2＝90
°
+csc(sin b2)
ꢀꢀ
(14)
[0266]
》当cos b2的值小于零且sin b2大于等于零，则b在第二象限
[0267]
b2＝90
°
+csc(sin b2)
ꢀꢀ
(15)
[0268]
》当cos b2的值小于零且sin b2小于零，则b在第三象限
[0269]
b2＝90
°‑
csc(sin b2)
ꢀꢀ
(16)
[0270]
由此，计算出m点的经纬度。重复上述计算步骤，则可计算出n、a、b、c、d的经纬度。
[0271]
(3)扫描区域处理。扫描区采用栅格化方法处理，栅格化的方法为本领域技术人员熟知，不再累述。
[0272]
1)地形高程获取
[0273]
在获得考虑地球曲率影响下的各点经纬度后，a、b、c、d、m、n各点海拔高度，由其经纬度带入高程数据后即可确定。
[0274]
2)扫描范围栅格化
[0275]
栅格化后的侦察覆盖范围如附图7所示，根据a、b、c、d四点所组成的区域进行栅格化，ad方向m等分，ab方向n等分，就可以得到m*n个子四边形区域，存储这些子四边形区域的中心点，为用作通视采样点做准备。
[0276]
栅格化主要用到的几个参数如下：
[0277]
①
m：方位向离散化等分系数，m为整数，最小值为1，其值越小，细分颗粒度越粗，结果偏差越大(计算效率越高)，由规划人员输入设定；
[0278]
②
n：距离向离散化系数，设定方法同方位向，m和n不一定相等；
[0279]
③
m*n：一次任务载荷侦察覆盖范围的离散点数；
[0280]
④
d：方位向最小细分距离，d＝l
ad
∠m，其中l
ad
是a、d水平投影点的连线。
[0281]
(4)栅格小区域通视计算
[0282]
利用地形高程数据、无人机当前位置数据以及栅格化后的数据进行各栅格小区域通视计算，即可得出各小区域的通视结果。通视分析时根据当前无人机所在位置和目标点的位置，结合高程地形进行两点之间是否可视的计算和判断。结合地形高程的视点和目标点通视分析示意如附图8所示。图中所示为无人机飞行剖面纵向切割地形获得的地形剖面，无人机实时位置与目标点连线，如果与地形剖面有交叉，则表示无人机在当前点与目标不通视(图中的不通视区)，反之则通视。
[0283]
通视算法有多种，在本发明中采用dyntacs通视性算法，该方法为本领域技术人员熟知，不再累述。如附图9所示。
[0284]
由于无人机任务载荷是否通视地面目标与地面目标是否通视无人机判断依据相同。为便于计算，将地面目标作为视点(v点)，无人机任务载荷的探测距离作为作用范围(换算为地距)，其边界位置(先选择某一个方向)作为目标点(t点)。然后从视点到目标点逐点分析交点。其中的f(xi，yi)、g(x
i+1
，y
i+1
)
…
等点均为视线vt在(x，y)平面上的投影同方形格网单元边的交点，视线vt的斜率α由下式算出：
[0285][0286]
其中：
[0287]zt
为t(x
t
，y
t
)点的高度；
[0288]
zv为v(xv，yv)点的高度；
[0289]
为t点到v点水平距离。
[0290]
以基准水平面向东为x轴，以v点铅锤面向上为y轴构建坐标系，视点v与各交点连线与水平面夹角βi斜率如下：
[0291][0292]
其中：
[0293]
zi为f(xi，yi)点的高度；
[0294]
zv为v(xv，yv)点的高度；
[0295]
为f点到v点水平距离。
[0296]
通视性判断可以通过比较tan α与tan βi的值完成(βi值通过公式(18)得到)：
[0297]
如果tan βi大于tan α，则表示不能通视，计算结束；
[0298]
如果tan βi不大于tan α，则f点通视，计算下一个点g(x
i+1
，y
i+1
)。如果可一直推进到目标点t，则视点与目标点之间可以通视。f点不一定是视线vt在(x，y)平面上的投影同方形格网单元边的第一个交点，g(x
i+1
，y
i+1
)也不一定是第二个或者紧挨着f的后面那个点。
[0299]
在根据公式(17)和公式(18)以及上述通视判断计算方式完成一次基于航路点数据、航路点等分数据进行栅格小区域通视计算为一次采样计算，在完成一次采样计算并得出结果后，按照栅格的参数按照方位向和距离向以此类推迭代计算，即可得出一次采样覆盖范围内的通视和非通视栅格。计算出的通视区域显示如附图10所示。图10栅格内的每个
点都与当前无人机所在的点进行通视计算，不通视的点，将其所在栅格填充为黑色，通视的点不做标记(图中为小圆点)。
[0300]
(5)合并计算
[0301]
基于航路点数据、航路点等分数据进行各通视计算点的栅格小区域通视计算，并将通视结果合并，获得如图13那样的通视和非通视栅格的图。此外，该步骤还同步计算扫描区域的合并结果。具体地，在一次采样计算完成后(公式(17)和公式(18))，获得如图10那样的通视和非通视栅格的图，以此类推，按照航路点数据和切分的间隔计算所有航线点的扫描区域，获得如图12那样的通视和非通视栅格的图，并进行合并计算。设定所有航线点扫描区域为u1，不通视区域为u2，扫描区与任务区交集为u3，不通视区域与任务区交集为u4。合并计算的具体算法步骤为：
[0302]
①
如图11所示，航点1、航点2表示无人机任务航线上的任意2个航点，其他航点的计算以此类推。按照采样间隔(建议为方位向最小细分距离d的整数倍)，在航点1和航点2中间平滑出所有插值航点，获取所有航点生成的扫描区域数据，扫描区域是如图11所示的多个扇形(每个扇形为航点1和航点2之间某个航点生成的扫描区域数据，在航点1和航点2中间平滑出所有插值航点，获取所有航点生成的扫描区域数据相当于插值采样，每个插值点都会计算出一个扇形)。
[0303]
②
将多边形区域边界点按照逆时针方向排列存储；
[0304]
③
将这些扇形区域合并，可获得由扫描区域边界点形成的多边形，遍历所有扫描区域形成的多边形区域，将当前被合并的区域与多边形区域中各多边形逐个求并集(区域求并集的方式、遍历方法为本领域技术人员熟知，不再累述)。
[0305]
④
基于遍历和逐个求并集的结果，计算出并集后的多边形区域边界点数据，合并后的区域如图11最外圈粗线表示(该方法为本领域技术人员熟知，不再累述)。
[0306]
(6)计算侦察覆盖率
[0307]
该步骤将合并计算得出的通视区域与任务区求交集，以此得出侦察覆盖范围，并进一步计算得出侦察覆盖率。
[0308]
设顺着航向方向递推(递推采样间隔设置为方位向最小细分距离d的整数倍，d由上文计算方式可知其最小值为1，最大值为m)，根据上述步骤
①‑④
计算出的通视区域与任务区交集u3和任务区u4，即可获得整条航向的侦察覆盖范围。
[0309]
图12中：
[0310]
p1是第一次计算时无人机所在点，也是任务载荷所在点；
[0311]
a1、b1、c1、d1分别是无人机在p1点时，无人机任务载荷当前位置扫面范围的边界点。a1d1弧线为近界，a1在左边界，d1在右边界，b1c1弧线为远界，b1在左边界，c1在右边界。m1和n1分别为a1d1弧线以及b1c1弧线的中心点。
[0312]
同理，p2点为下一次计算无人机所在点，也是任务载荷所在点；其他点的信息描述和定义与p1点对应的点一致。
[0313]
图13中：
[0314]
u1为所有航线点扫描区域；
[0315]
u2为所有的不通视区域；
[0316]
u3为扫描区与任务区交集；
[0317]
u4为不通视区域与任务区交集。图13中最下面的一个卷曲面是u4。
[0318]
该部分主要是二维图形的计算和处理，其数学复杂度较为简单，故不详细介绍。二维图形计算部分主要描述算法((1)魏许青《计算多边形交集、并集面积的算法》.湖南师范大学数学与计算机科学学院；(2)赫恩《计算机图形学》)步骤如下：
[0319]
①
计算所有任务区域集合；
[0320]
②
合并所有航点扫描的通视区域；
[0321]
③
在任务区集合中取出单个任务区；
[0322]
④
将任务区中的边界点数据逆序排列；
[0323]
⑤
计算区域数据与通视区域的交集，得到区域边界点；
[0324]
⑥
通过交集区域边界点，计算出交集面积；
[0325]
⑦
用交集面积除以任务区面积，即可得出此任务区的扫描覆盖率；
[0326]
⑧
重复
③‑⑦
即可计算出所有的任务区的覆盖率。
[0327]
⑨
返回任务区的覆盖率，任务区的交集区域边界点。
[0328]
针对更加广泛的任务载荷探测应用环境，通过本发明主要解决如下问题：
[0329]
(1)在多种地形下，如何实现对单点目标侦察监视的计算；
[0330]
(2)在多种地形下，如何实现对区域的覆盖侦察；
[0331]
(3)基于区域侦察的要求，如何规划合理的航线。
[0332]
为此，通过本发明可实现：
[0333]
(1)在不改变原有无人机规划业务流程的条件下，通过引入基于高程数据进行通视分析计算，从而实现更加准确的侦察覆盖率的计算；
[0334]
(2)在不改变原有无人机使用方式条件下，通过精确计算的覆盖率结果，从而实现更加合理的无人机航线规划；
[0335]
(3)基于该技术采用的方法，可扩展应用到类似的领域(如，无人机链路规划)，使该方法具有更好的适用性。
[0336]
基于上述应用需求，结合地形高程数据通视计算任务区侦察覆盖率的方式可极大提升任务规划结果合理性和可用性。