利用几何和CUDA编程对随机岛屿进行地理定位

# gralhix #004 来源: https://yassa9.github.io/osint/gralhix-004/ ## gralhix004 | 使用几何与CUDA GPU编程定位随机小岛图像 > 16-08-2026 注意：这是一项真实的人工作业，未使用LLM生成。我撰写本页作为Sofia Santos | Gralhix（https://gralhix.com/list-of-osint-exercises/osint-exercise-004/）制作的挑战gralhix 004的解题报告。 > 你可以在此处查看、克隆并本地运行所有代码文件及包含完整说明的最终报告：[GitHub仓库](https://github.com/yassa9/geoint/tree/main/gralhix_004)。 --- ## 任务简介：主文件 这是一张度假村的照片，拍摄于某个岛屿上。 > a) 度假村的名称是什么？ b) 岛屿的坐标是什么？ c) 拍摄照片时相机朝向哪个基本方向？ 在我看来，用`Google Lens`解决这个挑战是浪费一个有趣的机会，因此决定用数学和编程来解决。 --- ## a] 元数据 当然，首先查看的是`元数据`。在`Linux Void`系统上运行： ``` exiftool main.png File Type : WEBP (lossless) MIME Type : image/webp Image Width : 736 Image Height : 515 ``` 如预期，这里没有有用信息。没有EXIF、GPS数据，也没有相机制造商或型号信息。 --- ## b] 构建指纹 01_00 从图片中可以看到3个陆地： - P0：小岛本身 - P1：右侧岛屿 - P2：左前方岛屿（有山峰） 我无法构建正确的鸟瞰视角透视模型，因为图片显然是由无人机拍摄的，且无法估算海拔高度（元数据中也未提供）。因此只能凭直觉估算，只需要三个岛屿的相对距离和构成的三角形角度。 01_01 我构建了一个小型点击GUI `01_triangle_gui.py`，用于按顺序记录每个点的像素坐标并计算三角形几何属性。由于通过视觉精确点击中心点并不完美，我在搜索时为两个值都添加了`±20%容差`范围。 --- ## c] 搜索 指纹锁定后，下一步是将地球上每一个真实陆地与之对比！ > 我使用`OpenStreetMap的分割陆地多边形数据集`作为数据源（[land-polygons-split-4326](https://osmdata.openstreetmap.de/data/land-polygons.html)），这是完整的WGS84全球海岸线矢量数据，大小为`882 MB`。 我创建了启发式过滤器（全部基于直觉和非实证证据），花费数天（确实是整天）调整数值并进行大量试错😭，直到获得可行的过滤器组合。 ### 01] 热带纬度范围 $$ -30° \le \text{纬度} \le 30° $$ 照片中的小岛呈现热带特征，因此我决定在任何昂贵的几何计算之前，立即排除热带以外的区域。 > 正好`141,131`个多边形通过了该范围过滤。 ### 02] 局部密度过滤器 $$ N_{5\text{km}}(p) \le 10 $$ $ N_{5\text{km}}(p) $ 统计在点$(p)$5公里范围内落入的其他质心数量。`上限为10`：如果一个小岛有超过10个如此近的邻居，说明它位于密集的珊瑚礁区、拥挤的海岸线或群岛群中，而非如照片所示的小型孤立3-4岛群。 > 此过滤将候选数量降至`51,576`。 ### 03] 聚类 对每个存活点，查找20公里内的所有其他点（基于图像估计的启发值）。如果至少有2个邻近点（共3个点），则形成一个簇。没有3个以上邻近点的点被丢弃，因为它们无法构成三角形。 ```python 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{dist}(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} $$ ```python 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线程。每个线程按陆地面积对其3个点排序以选出P0（最小，即度假村小岛），然后使用另外两个点的缠绕方向分配P1和P2： ```c 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{cross} = x_a y_b - x_b y_a $$ $$ P1 = \begin{cases} a & \text{cross} > 0 \\ b & \text{cross} \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的间距以及两条边长均落在指纹的容差窗口内，则视为匹配。通过的线程使用原子计数器将结果写入共享输出数组，确保同时完成的线程不会互相覆盖： ```c if (hit) { unsigned long long slot = atomicAdd(out_count, 1ULL); out_p0[slot] = p0idx; out_p1[slot] = p1idx; out_p2[slot] = p2idx; } ``` 现在从内核直接打印在CLI中： ``` gpu: NVIDIA GeForce RTX 3050 (sm_86) vram used: 5169 MB kernel time: 204.1 ms ``` > `80.7百万`个三元组输入，每个线程一个，并行处理。`158,784`个通过了掩码筛选。 ### 06] 去重 由于同一个物理三元组可能属于多个重叠聚类，从而被多个GPU线程命中，因此原始匹配首先按标识进行去重： ```python 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所在侧的对面，然后检查陆地数据集是否有其他物体位于该矩形内。 ```python 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), ] ``` 如果除了3个候选岛屿本身之外有任何物体与该矩形相交，则该候选被丢弃。有陆地存在意味着它不是照片实际显示的开阔无阻碍水域。 > `8,915`个唯一三元组降至`948`个。 下方是948个候选地点的地图。 02_01 --- ## d] 珊瑚环礁形状检查 此阶段仅查看P0（度假村小岛），检查其形状是否确实像珊瑚环礁。 ### 1]`紧凑度`，形状接近圆形的程度： `Polsby-Popper评分：` $$ PP = \frac{4\pi \cdot \text{面积}}{\text{周长}^2} $$ ```python def compactness(row): return (4 * np.pi * row.area_km2) / (row.perim_km ** 2 + 1e-12) ``` 03_00 `1.0`是完美圆形，值越低表示轮廓越参差或细长。珊瑚环礁通常因波浪沉积而呈圆形，因此任何`< 0.5`的都被丢弃。 ### 2] 微型环礁光晕检查： ```python 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个。 > `213/948`个候选通过了两项检查。 --- ## e] 椭圆形状检查 对P0自身多边形进行另一项几何过滤。拟合其最小旋转矩形并从该矩形测量两个比率。 ```python 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{aspect} = \frac{\text{长边}}{\text{短边}} \in [1.05,\ 2.2] $$ 过于接近1.0意味着几乎是完美圆形，而非照片中的略微细长形状。比值过高则形状过于细长（超过2:1）。 `填充比率`是形状实际填充该外接矩形的程度，这背后有一个特性：任何椭圆恰好填充其最小面积外接矩形的$ \pi / 4 $，无论其拉伸程度如何。 $$ \frac{\text{面积}_{\text{椭圆}}}{\text{面积}_{\text{矩形}}} = \frac{\pi}{4} \approx 0.785 $$ 这是完美光滑椭圆的理论上限。真实珊瑚环礁并非完美椭圆，因此阈值设为该上限的一个启发性安全比例： $$ \text{FILL\_RATIO\_MIN} = 0.75 \times \frac{\pi}{4} \approx 0.589 $$ 形状需要保留至少75%的完美椭圆填充度才能存活。月牙形、环形和缺口海岸线远低于此值，而坚固的圆形环礁则高于此值。 > `137/213`个候选存活。 --- ## f] NDVI植被检查 我们到达了最终的API阶段，我将其放在最后，因为它是网络受限而非计算受限的。 > 我们将连接由Element84运行的`Earth Search`，这是一个公共STAC API，索引了托管在AWS开放数据计划上的Sentinel-2影像，免费且无需API密钥。你可以在https://earth-search.aws.element84.com/v1查看它。 我们现在检查P0是否真的有植被覆盖（棕榈树），而非裸露的沙子或岩石。它从公共STAC目录中拉取该点上方最近的低云量Sentinel-2场景，并在该精确像素处采样红光和近红外波段。 $$ \text{NDVI} = \frac{\text{近红外} - \text{红光}}{\text{近红外} + \text{红光}} $$ 活体植被在近红外波段强烈反射并吸收红光，因此健康的棕榈覆盖会使NDVI远高于0，裸沙或开阔水域则接近0或为负值。 04_00 你可以查看我从这篇出色的Geoawesome博客获得的[图像](https://geoawesome.com/eo-hub/understanding-aerial-data-normalized-difference-vegetation-index-ndvi/)。 `阈值设为0.6`，足够高以要求真实的树木覆盖，而不仅仅是零星分布。 > `66/137`个候选通过了NDVI检查。 --- ## g] 海拔与山脉检查 05_00 最终揭晓前的最后检查。有两个条件： - P0本身必须低平，与小环礁一致 - P2必须在相机实际朝向的方向上具有真实的高地地形 “前方”是指向P1方位与指向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{front}} = \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米Copernicus DEM图幅`进行采样。 > `Copernicus DEM GLO-30`，由欧盟哥白尼计划发布，作为免费的云端优化GeoTIFF托管在AWS Open Data上，无需账户或密钥。更多信息可查看[https://registry.opendata.aws/copernicus-dem/](https://registry.opendata.aws/copernicus-dem/)。 最终，这两个简单的启发式条件决定存活（是的，我知道，一切都变得启发式了哈哈）： $$ \text{海拔}(P_0) \le 50\text{米} $$ $$ 100\text{米} \le \max_{\text{弧线}}(\text{海拔}) \le 500\text{米} $$ 05_01 从这张抽象图中可以看到，虚线是相机的正前方航向，楔形区域是向外延伸至20公里的±50°搜索弧，用于海拔检查。 > `26/66`个候选通过了海拔检查。你可以看到这26个幸存者，除一个位于巴西附近外，其余均位于南亚、澳大利亚和大洋洲！ 05_02 --- ## h] 最终报告 最后阶段，只是让最终候选者可供目视检查。每个幸存者通过`点在多边形内查找`（针对国家边界文件）获取其国家名称，然后提供P0、P1、P2的直接Google Maps卫星链接。输出为纯HTML表格，包含索引、国家、每行三个可点击的坐标对。 06_00 我得到了这个最终列表，让我们逐一目视检查。不会在此逐一列举，但对我来说前7个完全不对。 06_01 直到我打开表格中第8个属于密克罗尼西亚的选项😍（第一次知道有密克罗尼西亚这个国家）： 06_03 并通过P1和P2确认： 06_04 那就是答案🎉...你可以在[Google Maps](https://www.google.com/maps/@7.3633,151.755983,50m/data=!3m1!1e3)上查看它。 --- ## 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{西北} $$ --- ## j] 数据与许可 **海岸线多边形**：[land-polygons-split-4326](https://osmdata.openstreetmap.de/data/land-polygons.html) © OpenStreetMap贡献者，根据[开放数据库许可（ODbL）1.0](https://opendatacommons.org/licenses/odbl/)提供。仓库中的候选集和最终报告是衍生数据库，并在同一许可下发布。 **海拔数据**：Copernicus DEM GLO-30