基于无人机GPS坐标的有效区域覆盖面积计算算法研究.pdf

ofJournalAviationUniversity63ofChinaFlightCivilVol.34No.5飞行学院学报中国民航Sept.2023基于无人机GPS坐标的有效区域覆盖面积计算算法研究康驰张鹏孟超陈俊杰（中电科特种飞机系统工程有限公司四川成都611731)摘要：无人机航线区域内覆盖作业需要统计其事后作业覆盖总面积，即根据无人机GPS航线坐标，累积途经区域的覆盖面积并忽略其重复覆盖区域，计算最终有效覆盖面积。难点在于重复覆盖的航线区域面积应该如何进行抠除处理。本文提出了一种全新的面积像素分割算法用以解决该问题。基本思想是将整块航线区域分割为网格小像素，使问题简化为判断像素是否包含在航线覆盖区域中。此面积像素分割算法可以非常精准地计算无人机区域覆盖面积。关键词：无人机面积覆盖像素坐标系统Research on Algorithm for Calculating Area of Coverage Based on GPS Coordinates of UAV.Kang ChiZhang PengMeng ChaoChen Junjie(CETC Special Mission Aircraft System Engineering Co.Ltd,Chengdu,611731,Sichuan,China)Abstract:In the application of coverage-operation by UAV(such as pesticide spraying)in theUAVs route area,the statistics of the total area of coverage after the operation is often required.Namely,the coverage areas through the path is accumulated based on the GPS coordinates of theUAV and the repeated areas shall be ignored.The difficulty lies in how to deal with the repeated cov-erage areas.A new algorithm named the area pixel partition algorithm is proposed to solve this prob-lem.The basic idea is to divide the whole UAVs area into small pixels,which simplifies the problemto determine whether one given pixel is contained in the whole region.The area pixel partition algo-rithm proposed in this paper can accurately solve the problem of calculating the area of coverage byUAV.Keywords:UAVAreaCoveragePixelCoordinate system1引言近年来，无人机技术发展迅速，服务于诸多领域，如无人机代替人工进行大面积农药喷洒作业。在无人机实施区域内覆盖作业的应用场景下，通过GPS追踪器可以取得其飞行航迹，并需要统计分析其事后作业覆盖总面积，即根据无人机GPS航线坐标，累积其途经区域的面积并忽略其重复覆盖区域，计算最终的有效区域覆盖面积。需要考虑无人机的真实覆盖有效面积，即需要抠除掉无人机在作业过程中没有覆盖的区域和重复覆盖的区域。算法设计的关键点在于如何进行未覆盖区域和重复覆盖区域的抠除处理操作。目前关于无人机的面积测量方法中，一般是根据无人机定位区域的若干边界点的GPS坐标组合为一个无规则多边形，然后利用不规则多边形面积计算算法求出区域面积 1-3 。也有很多基于无人机航拍图像进行面积计算的研究，例如文献 4 提出了一种通过图元分割、目视解译来获取目标样本，继而通过神经网络算法来进行图像识别并累加图片像素进而计算出目标面积的方法，文献 5 根据无人机航拍图像，并基于面向对象的分割和最邻近监督分类方法对图像进行目标区域面积快速提取。但是，以上所述方法有如下缺陷和局限性：a)只能根据区域多边形的边界点坐标计算其面积，应用范围存在局限。b)在无人机处于作业状态时，算法不能提供实时计算的能力，只能进行事后统计计算。JournalAviationVol.34No.5UniversityFlight64ofChinaofCivilSept.2023中国民航飞行学院学报c)目标区域极度不规则或者带有孔洞时，将其抽象为平面多边形存在技术难度。d)多边形边界点坐标的选择以及其精度、准确度对其最终面积计算影响大，不同的算法会导致计算结果偏离度较大。e)利用航拍进行面积估算，高度依赖于所选用的图像识别算法，精度和计算力都不能得到保证，不适用于机载微处理器应用环境。为了计算无人机航线区域覆盖的有效面积，本文提出一种新的基于GPS坐标的航线覆盖有效区域面积计算方法，广泛适用于类似无人机农药喷洒或者无人机扫描监视等应用场景，因其算法设计中的关键点在于将不规则区域分割为规则小像素点的集合，故又命名为面积像素分割算法。该算法根据无人机飞行航迹航线、飞行高度、向下覆盖角度计算其飞行区域覆盖面积，不被无人机飞行区域形状和边界局限，可以提供实时计算能力，并且其计算精度可以根据所选用的平面几何形状以及像素点大小来进行定制化调整。2算法分析2.1算法概述本文提出的算法解决了如何计算无人机飞行途经区域覆盖的总有效面积的问题。算法关键点在于：1)无人机坐标系统由地球坐标系转换为地面平面坐标系；2)无人机飞行航线路径覆盖区域如何抽象为平面几何图形表示；3)多个多边形的覆盖有效面积问题转换为像素分割与覆盖问题。2.2算法流程算法整体流程如图1所示。外部接口负责引接无人机飞行航线中所有航点的经度、纬度、离地高度和无人机向下覆盖角度数据，并将所有无人机航点的经纬度(B,L)通过坐标变换方法转换为地球平面坐标系的坐标值(X,Y)，继而结合无人机向下覆盖角度将整条航线区域抽象提取为平面几何图形的组合。因为在所有平面几何体中，以多边形几何体的面积计算最为通用和简便，所以在平面几何抽象过程中，需要尽量选用多边形作为最终抽象对象。接下来需无人机数据接入：1.经纬B,L2.高度H3.向下覆盖角度坐标变换：地球坐标系(B,L)-平面坐标系(X,Y)几何抽象：航线区域-几何图形计算区域是否大于一个阅值(m2)是区域分割否像素覆盖计算输出：航线覆盖区域总面积图1可面积像素分割算法流程图要计算平面中所有几何体的有效覆盖面积：定义一个1mm1mm的像素点，整块图形可以看成由多个像素点构成，算法需要统计出图形中包含的所有像素点即可计算出最终图形有效覆盖面积。需要特别说明，运行主机内存的容量会限制计算区域的大小，当整块图形区域过大导致像素点太多时，需要依靠区域分割方法将图形进行人工分割，从而使像素点计算量减少。最后累加所有被分割的子区域面积数据，即可得到总区域统计面积数据。AviationJournal65ofChinaUniversityFlightofCivilVol.34No.5飞行学院学报中国民航Sept.20232.3详细说明以下章节针对0 2.22算法流程所描述的步骤进行分别详细阐述。2.3.1无人机数据接入算法提供数据接入接口，完成对无人机航线航点的输入引接。必输入数据包括：1)航点经纬度(B,L)；2)航点高度H；3)无人机向下覆盖角度，该数据可以根据真实应用情况自行定义，主要用途是计算无人机航线占用的宽度值。输入数据需要保证数据可靠且按照时间顺序进行排序。2.3.2地图投影坐标变换获取了无人机三维经纬坐标之后，还需要使用地图投影变换技术将之转换为适宜计算的平面坐标。地图投影的含义就是建立地球经纬坐标(B,L)和平面直角坐标(X,Y)的映射关系 6,7 。由于地球本身是一个不可展开为平面的椭球体，所以不论采用何种数学方法进行投影变换，都将对原来的形状产生数学意义上的变形 8 。现在大多数的投影变换方法都会对地球的长度、面积、角度产生变形和破坏，需要根据地图实际作用来合理选择相应的投影方法。现在比较通用的投影方式有方位投影（又称平面投影）、圆锥投影、圆柱投影、伪圆锥投影、伪圆柱投影、多圆锥投影等等 6.8-9 ，例如著名的墨卡托投影(MercatorProjection)就是一种等角正切圆柱投影 8 ，兰勃特投影(LambertProjection)则是一种等角圆锥投影 9 10 。圆锥投影适用于中纬度地区 8 ，本文中将使用基于圆锥投影的兰勃特投影。其基本原理是用一个正圆锥割于地球两个标准纬线，再使用等角条件将地球投影到该正圆锥面上，然后沿着母线剪开，即是兰勃特投影平面 10 。兰勃特投影公式如下所示：X=rsin(1)Y=ro-rcoso(2)r=rearhFxt(3)0=n(L-Lo)(4)COsB(5)mJi-e?xsinB元B-esinBe12(6)42l+esinBIn(mg/mg2)(7)n=In(t/t2)mBlF(8)上式 11,12 中，(B,L)是无人机经度值和纬度值，单位为度；(X,Y)是转换后的平面直角坐标，X表示南北方向，Y表示东西方向，单位是m；(Bo,Lo)是自定义的平面直角坐标原点的经纬度，单位为度；Bi是标准纬线1，选择为北纬2 5 度，B2是标准纬线2，选择为北纬4 7 度；rearth为地球赤道半径，选择为6 3 7 8 13 7 m；r o 为原点纬度处的r值；mBl和mB2分别是标准纬线1和标准纬线2的m值；tB1和tB2分别是标准纬线1和标准纬线2 的t值。2.3.3航线覆盖区域几何图形抽象变换获取了平面坐标系下的航线航点坐标值后，需要根据航点的顺序集合转换为区域面积的集合。我们在后续讨论中认为无人机在作业时是垂直向下覆盖的，其覆盖范围由无人机的离地高度值H和向下覆盖角度决定。如果无人机向下覆盖不是垂直的，则可相应做覆盖几何修正，我们不在此处讨论该问题。算法采用如下方法：相邻两个航点扩展为一个矩形，矩形的长为两航点距离加上无人机向下覆盖直径，矩形的宽为无人机向下覆盖直径，如图2 所示。d/2图2 无人机航线覆盖区域示意图图2 中，H是无人机离地高度，是无人机向下覆盖角度，d是无人机向下覆盖的直径，L是相邻航点之间距离。有下式：JournalVol.34No.5AviationofChinaUniversityFlight66ofCivilSept.2023飞行学院学报中国民航d=H xtan(9)V(x2-x1)+(y2-yl)2(10)上式中，H是无人机离地高度，单位是m；是无人机向下覆盖角度，单位是弧度；（xl,y1)是前航点的平面直角坐标值，单位是m；（x 2,y 2)是后航点的平面直角坐标，单位是m。由此可以得到所有由相邻航点形成的多边形坐标集合，接下来需要计算所有矩阵的覆盖有效总面积。需要说明的是，可以根据实际使用场景对无人机覆盖区域进行不同的几何图形抽象变换，以达到相应的面积计算要求。2.3.4，区域分割在使用像素分割计算方法之前，需要考虑运行主机的内存。如果采用1mm1mm的像素点，且每个像素点占用一个Byte(当然也可以使用一个Bit代表一个像素点，这样可以节省7/8 的内存空间)，则1m的区域需要使用的内存空间为(10001000)/1024/1024=0.95Mb。如果不采用区域分割，当计算一个10 0 m100m的区域时，仅像素点就会占据9.3 Gb的内存大小，显然是不可接受的。算法采用区域分割的方法，将整块计算区域划分为1m的单位方块区域集合，然后分别计算所有单位方块区域中被无人机覆盖的面积，最后累加所有面积值即可得到总覆盖面积。区域分割方法比较简单，不再详述。2.3.5像素覆盖计算将1m单位方块区域齐整地划分为10 0 0 10 0 0个1mm的方块像素点（像素分割示意图见图3），然后轮巡所有代表航线覆盖区域的矩阵几何图形，统计该矩阵图形占据的像素点，并将已经占据的像素点Byte值设置为true，后续不需要再重复判断该像素点是否被占据。轮巡结束后，统计出所有被置为true的像素点的总数量，即可得到该单位方块区域被覆盖的总面积。在判断一个矩阵图形是否包含某像素点时，如果方块像素任一顶点在矩阵内部，或者矩阵任一顶点在像素内部，则说明矩阵和方块像素相交或者包含，可以判定为该矩阵图形包含此像素方块，否则两者不相交也互不包含。图3 像素分割法示意图在判断某点是否在多边形内部时，可采用向量法求解：设有多边形ABCD（逆时针排序）和点P，若P在ABCD内部，则应该同时满足：ABAP0,BCBP0,CDCP0,DADP0(11)否则点P在ABCD外部。考虑准确度的时候，显然最坏情况发生在当航线仅是一条直线且其无人机覆盖宽度恰好导致刚好在航线边缘多了两排像素点的时候。此时面积计算其准确度为：d(12)Ad+0.002上式中，d是无人机向下覆盖的直径，单位是m。例如，d取值为2 m时，计算准确度A=99.9%。2.3.6输出将所有单位方块区域覆盖面积相加，最后输出总覆盖面积。3算法验证使用C+11编程，采用MicrosoftVisualStu-dio2015作为主要开发工具，并将Windows7(64位)作为开发环境。计算主机硬件处理器使用Intel(R)Core(TM)i5-7500 CPU 3.40GHz，内存8GB。算法验证过程如下所述采用平面直角XY坐标系，构造2 2 个虚拟坐标点，每个点都采用(x,y,h)的格式用以记录，其中x和y分别是X方向和Y方向的坐标值，单位是m；h 是无人机离地高度值，单位m。2 2 个航点按照先后顺序分别为：JournalAviationUniversity67ofChinaFlightofCivilVol.34 No.5中国民航飞行学院学报Sept.2023(0,3.1234,1),(10,3.1234,1),(20,3.1234,1),(30,3.1234,1),(40,3.1234,1),(50,3.1234,1),(60,3.1234,1),(70,3.1234,1),(80,3.1234,1),(90,3.1234,1),(100,3.1234,1),(100,2.1234,1),(90,2.1234,1),(80,2.1234,1),(70,2.1234,1),(60,2.1234,1),(50,2.1234,1),(40,2.1234,1),(30,2.1234,1),(20,2.1234,1),(10,2.1234,1)，(0,2.12 3 4,1)，覆盖角度为向下9 0 度。注意，此处没有再验证兰勃特投影算法，故而略过坐标变换这一步骤，直接采用平面直角坐标系来验证面积像素分割算法。该航线理论上覆盖面积为3 0 6.0(m)，该算法计算出来的覆盖面积数据为3 0 6.10 2(m)，本次理论计算准确率为3 0 6/3 0 6.10 2=9 9.9 7%。由计算过程可以看出，如果想进一步提高准确率，可以将像素点设定为更小的方块，代价是需要付出更大的计算量。4结论a)本文提出的面积像素分割算法切实可行，可以解决计算无人机飞行途经区域总有效面积的问题，且具有极高的精度和准确度，若想进一步提高计算准确度，可以将像素点设定为更小的方块。b)通过进一步研究认为，算法还可以继续优化，以节省CPU计算开支，提高计算速度。需要找到一种合适的优化算法，以避免重复计算所有像素点。参考文献1袁文俊，郑武略，贺敏恒，等.基于多旋翼无人机GPS坐标计算的树障面积快速测定法 J.计算技术与自动化,2 0 19(4):4 2-4 62 Xiang Yang,Xiong Ying,Tang Chuzhou.Devel-opment and Experimental Analysis of SCM&GPS-based Farmland Area Measuring Device C/Interna-tional Conference on Intelligent Computation Tech-nology&Automation.IEEE Computer Society,201l1:58-613 Mesas-Carrascosa,Notario-Garcia,MariaDolores,et al.Validation of measurements of landplot area using UAV imagery J.International Jour-nal of Applied Earth Observation and Geoinfor-mation,2014,33:270-2794闵文芳，江朝晖，李想，等。基于无人机平台和图像分析的田间作物检测 J.湖南农业大学学报（自然科学版),2 0 17(2):2 12-2 165吴方明，张淼，吴炳方无人机影像的面向对象水稻种植面积快速提取 J.地球信息科学学报，2019(5):789-7986胡毓钜，龚剑文，地图投影（2 版）M.北京：测绘出版社，19 9 27 Lee Seong Kon.Summarized Reviews on GeodeticCoordinate System and Map Projection for Practition-ers in Exploration Geophysics JJ.2016(4):236-2488梁丽新.浅析地图投影的选择与应用 .西部资源,2 0 12(5)：18 1-18 39李英奎，吕肖庆，李敬多投影间地图投影变换实现的途径与优化 J.地理与地理信息科学，2000(2):79-8410龚清萍，许宇.地图投影在无人机航迹规划中的应用 J.航空电子技术,2 0 0 9(2):4 8-5 21l Savric Bojan,Jenny Bernhard.Automating theselection of standard parallels for conic map projec-tions J.Computers&Geoences,2016,90:202-21212夏兰芳，胡鹏，黄梦龙.地图投影解析变换的数值实现方法 J.测绘科学，2 0 0 7(3):6 9-7 1