Reed's News
← 返回精选

用几何与CUDA定位一座孤岛

Tech 81 yassa9 2026/8/19 1921 字 原文 ↗

2026年8月16日

说明：本文为人工原创，未使用大语言模型生成。

我撰写此页面，是为了记录Sofia Santos 制作的 gralhix 004 挑战 | Gralhix的解题过程。

你可以查看、克隆并本地运行所有代码文件，以及包含完整操作说明的最终报告， 访问此GitHub仓库即可获取。

任务简介：

主图

这是一张海岛度假村的照片，需完成三个任务： a) 度假村的名称是什么？ b) 该岛屿的坐标是多少？ c) 拍摄照片时，相机朝向哪个基本方位？

在我看来，用谷歌镜头（google lens）解题太浪费乐趣，因此决定用数学和编程方法完成挑战。

a] 元数据分析

当然，首先要检查的是图片元数据。我在Void Linux系统上运行了以下命令：

> exiftool main.png
文件类型：WEBP（无损压缩）
MIME类型：image/webp
图像宽度：736
图像高度：515

不出所料，这里没有任何有用信息——既无EXIF数据、GPS定位，也没有相机品牌或型号。

b] 构建特征指纹

01_00

从图片中可以看到三块陆地：

  • P0：度假村所在的小岛本身
  • P1：右侧的岛屿
  • P2：左前方带有山峰的岛屿

由于照片由无人机拍摄，且元数据中没有海拔信息，我完全无法估算拍摄高度，也就无法建立准确的俯视透视模型。

因此我只能凭直觉估算，只需要获取三座岛屿间的相对距离和三角形夹角即可。

01_01

我编写了一个小型点击式GUI程序01_triangle_gui.py，可以记录每个点的像素坐标并计算三角形的几何参数。

考虑到人工点击无法做到完全精准，我在后续搜索时为所有数值设置了±20%的容差范围。

c] 全球搜索

特征指纹确定后，下一步就是将其与地球上所有真实陆地进行比对！

我选用的数据集是 OpenStreetMap的陆地多边形分割数据集land-polygons-split-4326，这是一套采用WGS84坐标系的全球海岸线矢量数据，大小为882 MB。

我基于直觉设计了一系列启发式过滤规则，花了整整几天时间反复调整参数、不断试错😭，最终确定了一套可行的过滤流程。

01] 热带纬度范围筛选

$$ -30° \le 纬度 \le 30° $$

照片中的小岛具有热带特征，因此我直接排除了热带以外的所有区域，避免后续不必要的复杂几何计算。

经过这一步筛选，剩余的陆地多边形数量为 141,131。

02] 局部密度过滤

$$ N_{5\text{km}}(p) \le 10 $$

$ N_{5\text{km}}(p) $指的是点(p)周边5公里范围内的其他陆地质心数量。我将上限设为10：如果一座小岛周边5公里内有超过10个同类岛屿，说明它位于密集的礁盘、拥挤的海岸线或群岛区域，不符合照片中仅有3-4座孤立岛屿的特征。

这一步将候选数量降至 51,576。

03] 聚类筛选

对每个剩余的点，寻找其20公里范围内的所有其他点（此范围是根据图片视觉判断的启发式值）。如果一个点至少有2个这样的近邻（即形成至少3个点的集群），则保留该集群；无法形成3点及以上集群的点将被排除，因为它们无法构成照片中的三角形布局。

tree = cKDTree(f_coords)
neigh = tree.query_ball_point(
f_coords,
CLUSTER_RADIUS_KM / 111.0)
clusters = set(tuple(sorted(n)) for n in neigh if len(n) >= 3)

$$ \left|{q : \text{两点距离}(p,q) \le 20,\text{km}}\right| \ge 3 $$

筛选后剩余的集群数量为 23,500。

04] 生成三点组合

对每个集群，将其中任意3个点组合成候选三角形，即组合数$ C(n, 3) $。集群越大，组合数增长越快：例如一个包含60个点的集群，仅自身就能生成34,220组三点组合。因此我先将每个集群的点数上限设为60，并按岛屿面积分层抽样，而非随机截取。

$$ \binom{n}{3} = \frac{n(n-1)(n-2)}{6} $$

def stratified_sample(idx_arr, area_arr, cap):
order = np.argsort(area_arr[idx_arr])
n_small = cap // 3
n_large = cap // 3
n_mid = cap - n_small - n_large
mid_start = max(0, (len(idx_arr) - n_large - n_mid) // 2)
keep = np.unique(np.concatenate([
order[:n_small],
order[-n_large:],
order[mid_start:mid_start + n_mid],
]))
return idx_arr[keep]
def gen_cluster_triples(idx_arr):
local = np.array(list(
itertools.combinations(range(len(idx_arr)), 3)),
dtype=np.int64)
return idx_arr[local]

抽样规则为：选取三分之一的小岛、三分之一的大岛，以及三分之一中等面积的岛屿，而非使用完整集群或随机抽样。

最终，23,500个集群共生成80,690,777组三点组合！

05] GPU并行匹配

我为每组三点分配一个CUDA线程。每个线程会先按陆地面积对三个点排序，选出面积最小的P0（即度假村所在小岛），再通过另外两个点的环绕方向确定P1和P2：

long long i = blockIdx.x * (long long)blockDim.x + threadIdx.x;
if (i >= n_triples) return;
int pos[3] = {0, 1, 2};
for (int a1 = 1; a1 < 3; a1++)
{
int key = pos[a1];
double keyval = a[key];
int j = a1 - 1;
while (j >= 0 && a[pos[j]] > keyval)
{
pos[j + 1] = pos[j];
j--;
}
pos[j + 1] = key;
}

P1和P2的区分通过二维叉积实现，无需判断三点所属集群，仅根据叉积符号即可：

$$ \text{叉积} = x_a y_b - x_b y_a $$ $$ P1 = \begin{cases} a & \text{叉积} > 0 \ b & \text{叉积} \le 0 \end{cases} $$

从P0走到a再走到b，若叉积>0，则为左转（逆时针方向）；若叉积<0，则为右转（顺时针方向）。这是判断三点弯曲方向的常用符号法则。

随后，每个线程独立计算P0点的夹角和距离比，公式与构建特征指纹时一致：

$$ \theta_0 = \arccos\left(\frac{\vec{d_1} \cdot \vec{d_2}}{|\vec{d_1}||\vec{d_2}|}\right), \qquad r = \frac{|\vec{d_1}|}{|\vec{d_2}|} $$

只有当夹角、距离比、P0的面积、P0与P1的间距，以及两条边长全部落在特征指纹的容差范围内时，该组三点才算匹配成功。匹配成功的线程会通过原子计数器将结果写入共享输出数组，避免同时完成的线程互相覆盖：

if (hit)
{
unsigned long long slot = atomicAdd(out_count, 1ULL);
out_p0[slot] = p0idx;
out_p1[slot] = p1idx;
out_p2[slot] = p2idx;
}

以下是内核直接输出到命令行的信息：

gpu: NVIDIA GeForce RTX 3050 (sm_86)
显存使用：5169 MB
内核运行时间：204.1 ms

8070万组三点组合并行输入，每组对应一个线程，最终有158,784组通过筛选。

06] 去重

同一组真实岛屿可能属于多个重叠集群，从而被多个GPU线程匹配到，因此需要先对原始匹配结果按岛屿组合去重：

seen = set()
uniq = []
for i in range(len(p0_all)):
key = (p0_all[i], p1_all[i], p2_all[i])
if key not in seen:
seen.add(key)
uniq.append(i)

去重后剩余8,915组唯一的三点组合。

07] 开阔水域验证

02_00

对每一组剩余的三点组合，还要进行最后一项测试：其周边是否存在照片中所示的开阔水域？我们沿P0→P1的边缘构建一个矩形，位置在P2所在侧的对侧，然后检查该矩形范围内是否存在其他陆地。

width = np.hypot(x1, y1)
u = np.array([x1, y1]) / width
v = np.array([-u[1], u[0]])
# 根据构造，p2位于+v侧，
# 因此检查-v侧
length = 2 * width
corners_local = [
(0, 0), (x1, y1),
(x1 - v[0]*length, y1 - v[1]*length),
(-v[0]*length, -v[1]*length),
]

如果矩形范围内存在除三座候选岛屿外的其他陆地，则该组合被排除——因为这不符合照片中开阔无遮挡的水域特征。

经过这一步，8,915组唯一组合被筛选至948组。

下图是这948组候选地点的分布图：

02_01

d] 珊瑚岛形态验证

这一步我们仅关注P0（度假村所在小岛），验证其形态是否符合珊瑚岛的特征。

1] 紧凑度：衡量形状与圆形的接近程度

波尔斯比-波珀得分（Polsby Popper Score）： $$ PP = \frac{4\pi \cdot \text{面积}}{\text{周长}^2} $$

def compactness(row):
return (4 * np.pi * row.area_km2) / (row.perim_km ** 2 + 1e-12)

03_00

得分1.0代表完美圆形，数值越低说明形状越不规则或狭长。珊瑚岛受海浪沉积作用影响，通常呈圆形，因此得分< 0.5的组合被排除。

2] 微型岛礁环绕检查

def micro_cay_count(gdf, sindex, lon, lat):
dists_km = nearby.geometry.distance(pt) * 111.0
mask = (dists_km > 0)
& (dists_km <= HALO_KM)
& (nearby["area_km2"].values < MICRO_KM2)
return int(mask.sum())

我们统计P0周边1.5公里范围内面积小于0.05平方公里的陆地碎片数量（此为启发式参数）。真实的珊瑚礁系统中，主岛周围通常会散布着一些沙洲，而非孤立存在（我曾在这上面栽过跟头😭）。因此，候选点周边至少需要有1个这样的微型岛礁。

最终，948组候选中有213组通过两项检查。

e] 椭圆形验证

这是针对P0多边形的另一项几何过滤。我们先拟合出P0的最小旋转外接矩形，再计算两个比值：

def aspect_and_fill(geom):
mrr = geom.minimum_rotated_rectangle
coords = list(mrr.exterior.coords)
s1 = math.hypot(coords[1][0] - coords[0][0],
coords[1][1] - coords[0][1])
s2 = math.hypot(coords[2][0] - coords[1][0],
coords[2][1] - coords[1][1])
long_side, short_side = max(s1, s2), min(s1, s2)
return long_side / short_side, geom.area / mrr.area

长宽比指外接矩形的长边与短边之比：

$$ \text{长宽比} = \frac{\text{长边}}{\text{短边}} \in [1.05,\ 2.2] $$

比值过于接近1.0说明形状接近完美圆形，不符合照片中略带狭长的形态；比值过高则说明形状过于细长（超过2:1），同样不符合要求。

填充率指岛屿实际面积占外接矩形面积的比例，这一指标有理论依据：任何椭圆的面积恰好是其最小外接矩形面积的$ \pi / 4 $，与椭圆的拉伸程度无关。

$$ \frac{\text{椭圆面积}}{\text{矩形面积}} = \frac{\pi}{4} \approx 0.785 $$

这是光滑椭圆形的理论上限。由于真实珊瑚岛并非完美椭圆，我们将阈值设为该上限的75%作为安全值：

$$ \text{最小填充率} = 0.75 \times \frac{\pi}{4} \approx 0.589 $$

只有填充率不低于完美椭圆75%的形状才能通过筛选。新月形、环形或有缺口的海岸线填充率远低于此，而坚实圆润的珊瑚岛则能达标。

最终，213组候选中有137组通过验证。

f] NDVI植被验证

我们进入最后一个API调用阶段，我将其放在最后，因为它受网络限制而非计算能力限制。

我使用的是 Element84运营的Earth Search，这是一个公开的STAC API，索引了AWS开放数据项目托管的Sentinel-2卫星影像，免费使用且无需API密钥。

你可以访问https://earth-search.aws.element84.com/v1了解详情。

我们需要验证P0是否覆盖植被（如棕榈树），而非裸露的沙地或岩石。程序会从公开STAC目录中获取P0点上空最新的低云量Sentinel-2影像，采样该点的红光和近红外波段数值，计算归一化植被指数（NDVI）：

$$ \text{NDVI} = \frac{\text{近红外波段值} - \text{红光波段值}}{\text{近红外波段值} + \text{红光波段值}} $$

健康植被对近红外光反射强烈，对红光吸收明显，因此茂密的棕榈林会使NDVI值远高于0；而裸露沙地或开阔水域的NDVI值接近0甚至为负。

04_00

这张示意图来自我很喜欢的Geoawesome博客。

我将阈值设为0.6，确保候选点存在真正的树木覆盖，而非零星植被。

最终，137组候选中有66组通过NDVI验证。

g] 海拔与山峰验证

05_00

这是最终揭晓答案前的最后一项检查，包含两个条件：

  • P0本身必须地势低平，符合小型珊瑚礁岛的特征；
  • P2在相机实际朝向的方向上必须有真正的高地。

拍摄的"正前方"是P0指向P1和P0指向P2两个方位角的角平分线：

$$ \theta(P_0, P_i) = $$ $$ \text{atan2}\Big(\sin(\Delta\lambda)\cos\phi_i,\ \cos\phi_0\sin\phi_i - \sin\phi_0\cos\phi_i\cos(\Delta\lambda)\Big) $$

$$ \theta_{\text{正前方}} = $$ $$ \theta(P_0, P_2) + \frac{\big((\theta(P_0,P_1) - \theta(P_0,P_2) + 180) \bmod 360\big) - 180}{2} $$

由此得到相机朝向的方位角。以此为中心，我们在±50°范围内、2公里至20公里半径内选取一系列采样点：

$$ (\text{纬度}, \text{经度}) = \Big(\text{纬度}_0 + \frac{r\cos\theta}{111},\ \ \text{经度}_0 + \frac{r\sin\theta}{111\cos(\text{纬度}_0)}\Big) $$

每个采样点的海拔数据均来自真实的30米分辨率哥白尼DEM瓦片。

Copernicus DEM GLO-30由欧盟哥白尼计划发布，作为免费的云优化GeoTIFF文件托管在AWS开放数据平台，无需账户或密钥即可访问。

更多信息可查看https://registry.opendata.aws/copernicus-dem/

最终，以下两个简单的启发式条件决定候选是否通过（没错，一切都是凭直觉的规则，哈哈）：

$$ \text{P0海拔} \le 50\text{米} $$ $$ 100\text{米} \le \text{弧形区域最高海拔} \le 500\text{米} $$

05_01

从这张抽象示意图可以看到，虚线代表相机正前方的方位角，楔形区域是±50°的搜索范围，延伸至20公里用于海拔检查。

最终，66组候选中有26组通过海拔验证。

这26个候选点均位于南亚、澳大利亚和大洋洲，仅一个候选点位于巴西附近！

05_02

h] 最终人工核验

最后一步是让候选点便于人工核验。程序会通过点面查询，对照国家边界文件为每个候选点匹配所属国家，然后生成P0、P1、P2三点的Google Maps卫星链接。

输出结果为纯HTML表格，每行包含索引、国家名称，以及三个可点击的坐标对。

06_00

我得到了这份最终列表，接下来逐一人工核验。

这里就不逐个赘述了，但前7个候选点明显不符合。

06_01

直到我打开列表中第8个候选点，它位于密克罗尼西亚（我第一次知道有个国家叫密克罗尼西亚）😍：

06_03

再对照P1和P2的位置确认：

06_04

这就是正确答案🥳……

你可以在谷歌地图上查看该地点。

i] 最终答案 ...

a) 度假村的名称是什么？

$$ \text{Oan} $$

b) 该岛屿的坐标是多少？

$$7^\circ,21^\prime,48.4^{\prime\prime},\text{N} \qquad 151^\circ,45^\prime,20.7^{\prime\prime},\text{E}$$

$$ \text{或} $$

$$7.363444^\circ,\ 151.755750^\circ$$

c) 拍摄照片时，相机朝向哪个基本方位？

$$ \because\quad \theta = \text{atan2}\Big(\sin(\Delta\lambda)\cos\phi_1,\ \cos\phi_0\sin\phi_1 - \sin\phi_0\cos\phi_1\cos(\Delta\lambda)\Big) $$

$$ P_0 = (7.3633,\ 151.755983), \quad P_1 = (7.386573,\ 151.739534) $$

$$ \therefore\quad \theta = 324.97^\circ \implies \textbf{西北方向（NW）} $$

j] 数据与许可证

海岸线多边形数据：

land-polygons-split-4326 © OpenStreetMap贡献者，基于 开放数据库许可证（ODbL）1.0提供。仓库中的候选集和最终报告属于衍生数据库，同样遵循该许可证发布。

海拔数据：

Copernicus DEM GLO-30. © DLR e.V. 2010-2014 及 © Airbus Defence and Space GmbH 2014-2018，由欧盟和欧洲航天局通过哥白尼计划提供；保留所有权利。

卫星影像：

包含经过修改的哥白尼哨兵数据（2025-2026），通过Element 84在AWS开放数据平台上运营的Earth Search获取。

国家边界数据：

Natural Earth 10米分辨率 admin-0数据集，属于公有领域。

挑战与原始照片：

OSINT练习#004 由Sofia Santos（gralhix）制作。

(h)部分的卫星截图来自Google Maps / Google Earth