仿大疆司空2面状航线生成——凸多边形区域航线生成算法详解

仿大疆司空2面状航线生成——凸多边形区域航线生成算法详解

一、前言

去年，在针对大疆上云API进行二次开发的过程中，有一个需求是实现大疆司空2中的面状航线功能。在经过上网搜索后，在github上找到了一个开源项目cpRPA（植保无人机凸多边形地块工作路线规划），可以实现面状航线的生成。

参考项目github地址：https://github.com/Char-Ten/cpRPA

经过对其中算法的研究，感觉其实现方式非常值得学习，遂用Java重新实现了相关的逻辑，将其中的关键算法思路整理到本篇博客中。

本篇博客仅作为学习整理使用，由AI总结生成，实际代码请查看原项目。

下面将详细整理这一算法的实现思路，涵盖坐标旋转、扫描线求交、面积计算等核心环节。

二、算法总体流程

输入：凸多边形顶点坐标、飞行间隔(space)、旋转角度(rotate)
                          │
                ┌─────────▼──────────┐
                │ 1. 计算多边形边界框  │
                └─────────┬──────────┘
                          │
              ┌───────────▼────────────┐
              │ 2. 将多边形旋转 -rotate │
              │   （使扫描线平行于纬线）  │
              └───────────┬────────────┘
                          │
              ┌───────────▼────────────┐
              │ 3. 计算旋转后的边界框    │
              └───────────┬────────────┘
                          │
              ┌───────────▼────────────┐
              │ 4. 按间隔生成平行纬线    │
              └───────────┬────────────┘
                          │
              ┌───────────▼────────────┐
              │ 5. 每条纬线与多边形求交  │
              │   奇偶行交替方向(锯齿形) │
              └───────────┬────────────┘
                          │
              ┌───────────▼────────────┐
              │ 6. 将航线旋转回原始方向   │
              └───────────┬────────────┘
                          │
              ┌───────────▼────────────┐
              │ 7. 计算面积等统计信息     │
              └───────────┴────────────┘
                          │
输出：航线航点列表、多边形面积、航线覆盖面积

先旋转再扫描原因： 因为扫描线逻辑固定为水平方向（平行纬线），如果用户希望航线方向是任意角度，就先把多边形反向旋转，在旋转后的坐标系中做水平扫描，最后再旋转回来。这样只需要实现水平扫描线一种逻辑。

三、核心算法详解

3.1 坐标旋转 — 等距柱状投影 + 2D仿射变换

经纬度是球面坐标，不能直接做平面旋转。算法采用等距柱状投影（Equirectangular Projection）将经纬度近似转为平面坐标，再进行旋转。

等距柱状投影的关键： 经度方向的实际距离与纬度有关，纬度越高，同样1度经度对应的实际距离越短。因此需要乘以 cos(纬度) 进行校正。

// 等距柱状投影：经度方向乘以 cos(纬度)
double cosLat = Math.cos(Math.toRadians(center.getLat()));

// 经纬度 -> 平面坐标
double px = point.getLng() * cosLat;
double py = point.getLat();

// 2D旋转变换
double rad = Math.toRadians(deg);
double newX = (px - cx) * cos(rad) - (py - cy) * sin(rad) + cx;
double newY = (px - cx) * sin(rad) + (py - cy) * cos(rad) + cy;

// 平面坐标 -> 经纬度
double newLng = newX / cosLat;
double newLat = newY;

原理： 标准的2D旋转公式，绕点 (cx, cy) 旋转角度 θ：

3.2 扫描线生成 — 按间隔划分平行纬线

根据旋转后多边形的南北边界和飞行间隔，计算需要多少条扫描线：

// 南北两端距离（米）
double distance = haversineDistance(nw, sw);

// 扫描线数量（除以2是因为锯齿形来回算两条航线宽度）
int steps = (int) (distance / space / 2);

// 每条扫描线之间的纬度差
double latStep = (nw.getLat() - sw.getLat()) / steps;

3.3 扫描线与多边形求交 — 核心扫描逻辑

对每条水平扫描线，遍历多边形的每条边，求交点：

/**
 * 计算纬线 y 与线段 (p1, p2) 的交点
 * p1, p2 格式: [lng, lat]
 * 返回交点 [lng, lat]，无交点返回 null
 */
private double[] createInlinePoint(double[] p1, double[] p2, double y) {
    double s = p1[1] - p2[1];
    if (s == 0) return null;  // 水平边，无交点

    // 线性插值求交点的经度
    double x = (y - p1[1]) * (p1[0] - p2[0]) / s + p1[0];

    // 检查交点是否在线段范围内
    if (x > p1[0] && x > p2[0]) return null;
    if (x < p1[0] && x < p2[0]) return null;

    return new double[]{x, y};
}

几何原理： 已知线段两端点 (x1, y1) 和 (x2, y2)，水平线 y = Y，由相似三角形得：

3.4 锯齿形路径生成

对于凸多边形，每条扫描线恰好得到两个交点。将这两个交点按奇偶行交替排列，形成锯齿形路径：

偶数行（第0, 2, 4...行）：从左到右  →
奇数行（第1, 3, 5...行）：从右到左  ←

for (int i = 0; i < latLine.len; i++) {
    // ... 求出当前扫描线与多边形的两个交点 line[0], line[1] ...

    if (i % 2 != 0) {
        // 奇数行：先右后左
        polyline.add(new LatLng(y, Math.max(lng1, lng2)));
        polyline.add(new LatLng(y, Math.min(lng1, lng2)));
    } else {
        // 偶数行：先左后右
        polyline.add(new LatLng(y, Math.min(lng1, lng2)));
        polyline.add(new LatLng(y, Math.max(lng1, lng2)));
    }
}

生成的航线示意：

    ┌──────────────────────┐
    │  ·──────────────→·   │
    │  ·←──────────────·   │
    │   ·──────────→·      │
    │    ·←────────·       │
    └──────────────────────┘

3.5 Haversine公式 — 球面距离计算

计算地球上两点间的距离，使用经典的 Haversine 公式：

private double haversineDistance(LatLng p1, LatLng p2) {
    double lat1 = Math.toRadians(p1.getLat());
    double lat2 = Math.toRadians(p2.getLat());
    double dLat = Math.toRadians(p2.getLat() - p1.getLat());
    double dLng = Math.toRadians(p2.getLng() - p1.getLng());

    double a = sin(dLat/2) * sin(dLat/2)
             + cos(lat1) * cos(lat2) * sin(dLng/2) * sin(dLng/2);
    double c = 2 * atan2(sqrt(a), sqrt(1-a));

    return EARTH_RADIUS * c;  // EARTH_RADIUS = 6371000 米
}

公式推导：

3.6 面积计算 — Shoelace公式

多边形面积采用 Shoelace（鞋带）公式，先将经纬度转换为米制平面坐标：

private double getPolygonArea(List<LatLng> polygon) {
    double s = 0;
    for (int i = 0; i < polygon.size(); i++) {
        LatLng curr = polygon.get(i);
        LatLng next = polygon.get((i + 1) % polygon.size());

        // 经纬度转米制坐标
        double x1 = curr.getLng() * lng2m(curr);  // 1度经度 = 多少米
        double y1 = curr.getLat() * lat2m(curr);   // 1度纬度 = 多少米
        double x2 = next.getLng() * lng2m(next);
        double y2 = next.getLat() * lat2m(next);

        s += x1 * y2 - y1 * x2;  // 叉积累加
    }
    return Math.abs(s) / 2;
}

Shoelace公式：

四、前端凸多边形检测

后端算法仅适用于凸多边形，因此前端在用户绘制时实时检测凸性，使用叉积符号一致性判断：

function isConvex(pts) {
    var n = pts.length;
    if (n < 3) return true;
    var sign = 0;
    for (var i = 0; i < n; i++) {
        var a = pts[i], b = pts[(i + 1) % n], c = pts[(i + 2) % n];
        // 相邻三点的叉积
        var cross = (b.lng - a.lng) * (c.lat - b.lat)
                  - (b.lat - a.lat) * (c.lng - b.lng);
        if (cross !== 0) {
            if (sign === 0) {
                sign = cross > 0 ? 1 : -1;
            } else if ((cross > 0 ? 1 : -1) !== sign) {
                return false;  // 叉积符号不一致，非凸多边形
            }
        }
    }
    return true;
}

原理： 沿多边形顶点依次取相邻三个点，计算向量叉积。如果所有叉积符号一致（全正或全负），说明多边形始终朝同一方向转弯，即为凸多边形；一旦出现符号不同，说明存在"凹陷"。

五、算法复杂度分析

步骤	时间复杂度	说明
计算边界框	O(n)	遍历所有顶点
坐标旋转	O(n)	逐点变换
扫描线求交	O(s * n)	s条扫描线，每条遍历n条边
面积计算	O(n)	Shoelface公式遍历顶点

其中 n 为多边形顶点数，s 为扫描线数量（= 区域高度 / 飞行间隔）。

实际场景中 n 通常 < 20，s 通常 < 1000，计算量很小，可以实时响应。

六、总结

本算法的核心思想是 "旋转-扫描-逆旋转"：

1. 通过坐标旋转，将任意角度的扫描问题统一为水平扫描
2. 利用凸多边形的性质（每条扫描线恰好两个交点），简洁地生成锯齿形路径
3. 逆旋转还原到真实坐标系

这种方法简洁高效，适用于无人机需要区域全覆盖飞行的场景。