From 237fcba6476f50ad7d5c36ad4074a0eecbb332da Mon Sep 17 00:00:00 2001 From: duxin Date: Wed, 8 Jul 2026 14:00:59 +0800 Subject: [PATCH] =?UTF-8?q?perf:=20=E8=87=AA=E9=80=82=E5=BA=94=E5=88=86?= =?UTF-8?q?=E5=9D=97=20+=2040%=20=E9=87=8D=E5=8F=A0=E7=BC=93=E5=86=B2=20+?= =?UTF-8?q?=20n=5Fclosest=3D50=20=E4=BC=98=E5=8C=96=E5=B1=80=E9=83=A8?= =?UTF-8?q?=E5=85=8B=E9=87=8C=E9=87=91?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 用户需求: 500m 固定分块(50块)容易碎片化, 稀疏区域易受异常值污染 修复 (3项): 1. 自适应块大小: 根据 extent 自动切分为 ~4×4 块(16块), 纵横比极端时自动在长边增加 1 块。 每块面积足够大, 囊括更多采样点特征。 2. 40% 重叠缓冲区: 采样点搜索范围为块长宽的 1.4×, 相邻块交界处使用极度重叠的采样点 → 消除拼缝断层。 网格点(grid_x/y)仍使用严格不重叠块范围 → 无冗余计算。 3. n_closest_points=50: 50 个近邻稀释极端异常值影响, 比之前 15 个更稳健, 且协方差矩阵 50×50 仍远小于全局 8242×8242。 --- src/postprocessing/map.py | 108 ++++++++++++++++++++------------------ 1 file changed, 58 insertions(+), 50 deletions(-) diff --git a/src/postprocessing/map.py b/src/postprocessing/map.py index 375c4cb..21dbe39 100644 --- a/src/postprocessing/map.py +++ b/src/postprocessing/map.py @@ -597,18 +597,15 @@ class ContentMapper: grid_y = grid_yy[:, 0] total_cells = len(grid_x) * len(grid_y) - # 局部窗口：每个块使用 500m 空间窗口内的采样点 - _LOCAL_WINDOW = 500.0 # 米 - _n_workers = min(os.cpu_count() or 4, 8) - print(f"正在使用 局部克里金 (Local Kriging)：" - f"网格={total_cells:,} 点, 窗口={_LOCAL_WINDOW}m, " + print(f"正在使用 局部克里金 (自适应分块 + 40% 重叠缓冲)：" + f"网格={total_cells:,} 点, " f"{_n_workers} workers") grid_content = self._local_kriging( points, values, grid_x, grid_y, - window_size=_LOCAL_WINDOW, n_workers=_n_workers, + n_closest_points=50, ) valid_mask = ~np.isnan(grid_content) @@ -695,84 +692,94 @@ class ContentMapper: return grid_content def _local_kriging(self, points, values, grid_x, grid_y, - window_size=500.0, n_workers=4): - """局部克里金：网格按 window_size 分块，每块只用块内+缓冲区的局部采样点 + n_workers=4, n_closest_points=50): + """局部克里金：自适应分块 + 重叠缓冲区 + 保护性近邻限制 - 核心思想：协方差矩阵锁定在局部 n×n（几十到几百），而非全局 8242×8242。 - 千万级网格秒级完成，精度几乎无损。 + 1. 自适应块大小: 根据 extent 自动切分为 ~4×4 块 (16~25块) + 2. 重叠缓冲区: 采样点范围扩展块长宽的 40%，交界处平滑无拼缝 + 3. 保护性近邻: n_closest_points=50，稀释极端异常值 + 4. 网格点仅使用严格不重叠的块范围（Buffer 仅用于筛选采样点） """ import multiprocessing - # 计算网格空间范围 x_min, x_max = float(grid_x[0]), float(grid_x[-1]) y_min, y_max = float(grid_y[0]), float(grid_y[-1]) grid_dx = float(grid_x[1] - grid_x[0]) if len(grid_x) > 1 else 1.0 grid_dy = float(grid_y[1] - grid_y[0]) if len(grid_y) > 1 else 1.0 + extent_x = x_max - x_min + extent_y = y_max - y_min - # 按 window_size 将空间分割为不重叠的块 - # 每个块向四周扩展一个 window_size 作为缓冲区（获取局部采样点） - n_blocks_x = max(1, int(np.ceil((x_max - x_min) / window_size))) - n_blocks_y = max(1, int(np.ceil((y_max - y_min) / window_size))) - block_dx = (x_max - x_min) / n_blocks_x - block_dy = (y_max - y_min) / n_blocks_y + # ── 自适应块数：目标 4×4 = 16 块（最少 2×2，最多 6×6）── + _TARGET_BLOCKS = 4 + n_blocks_x = max(2, min(6, _TARGET_BLOCKS)) + n_blocks_y = max(2, min(6, _TARGET_BLOCKS)) + # 若 extent 纵横比极端，在较长维度上增加块数 + if extent_y > extent_x * 1.5: + n_blocks_y = min(6, n_blocks_y + 1) + elif extent_x > extent_y * 1.5: + n_blocks_x = min(6, n_blocks_x + 1) + + block_dx = extent_x / n_blocks_x + block_dy = extent_y / n_blocks_y + # 缓冲区 = 块长宽的 40%（兼顾平滑与计算量） + buffer_x = block_dx * 0.4 + buffer_y = block_dy * 0.4 total_blocks = n_blocks_x * n_blocks_y - print(f" 空间分块: {n_blocks_x}×{n_blocks_y} = {total_blocks} 块 " - f"(窗口={window_size}m, 缓冲区={window_size}m)") + print(f" 自适应分块: {n_blocks_x}×{n_blocks_y} = {total_blocks} 块 " + f"(块≈{block_dx:.0f}×{block_dy:.0f}m, 缓冲={buffer_x:.0f}×{buffer_y:.0f}m, " + f"n_closest={n_closest_points})") # 构建任务列表 tasks = [] for iy in range(n_blocks_y): for ix in range(n_blocks_x): - # 块范围（不含缓冲区） + # 严格不重叠的块范围 bx_min = x_min + ix * block_dx bx_max = x_min + (ix + 1) * block_dx by_min = y_min + iy * block_dy by_max = y_min + (iy + 1) * block_dy - # 块自身的网格（不含缓冲，保证每个网格点只被计算一次） - grid_mask_x = (grid_x >= bx_min) & (grid_x <= bx_max) - grid_mask_y = (grid_y >= by_min) & (grid_y <= by_max) - if not np.any(grid_mask_x) or not np.any(grid_mask_y): + # 块自身的纯净网格（不含缓冲） + gmask_x = (grid_x >= bx_min) & (grid_x <= bx_max) + gmask_y = (grid_y >= by_min) & (grid_y <= by_max) + if not np.any(gmask_x) or not np.any(gmask_y): continue - sub_grid_x = grid_x[grid_mask_x] - sub_grid_y = grid_y[grid_mask_y] + sub_grid_x = grid_x[gmask_x] + sub_grid_y = grid_y[gmask_y] - # 带缓冲区的搜索范围（仅用于筛选采样点） - search_xmin = bx_min - window_size - search_xmax = bx_max + window_size - search_ymin = by_min - window_size - search_ymax = by_max + window_size + # 缓冲区采样点范围（40% 重叠） + sx_min = bx_min - buffer_x + sx_max = bx_max + buffer_x + sy_min = by_min - buffer_y + sy_max = by_max + buffer_y - # 缓冲区内的采样点 - point_mask = ( - (points[:, 0] >= search_xmin) & (points[:, 0] <= search_xmax) & - (points[:, 1] >= search_ymin) & (points[:, 1] <= search_ymax) + pmask = ( + (points[:, 0] >= sx_min) & (points[:, 0] <= sx_max) & + (points[:, 1] >= sy_min) & (points[:, 1] <= sy_max) ) - local_pts = points[point_mask] - local_vals = values[point_mask] + local_pts = points[pmask] + local_vals = values[pmask] - sub_n = len(sub_grid_x) * len(sub_grid_y) tasks.append(( local_pts, local_vals, sub_grid_x, sub_grid_y, - bx_min, bx_max, by_min, by_max, # 只保留核心区域结果 - grid_dx, grid_dy, + bx_min, bx_max, by_min, by_max, + n_closest_points, ix, iy, total_blocks, )) - print(f" 有效块: {len(tasks)} (含采样点的空间块)") + print(f" 有效块: {len(tasks)}") if len(tasks) <= 1: return self._local_krige_block(tasks[0]) - # 多进程执行 print(f" 启动 {min(n_workers, len(tasks))} 个 worker 进程...") with multiprocessing.Pool(processes=min(n_workers, len(tasks))) as pool: results = pool.map(_local_krige_block_worker, tasks) - # 拼接结果：初始化全 NaN 数组，逐块填充核心区域 + # 拼接：全 NaN 数组，逐块填回 grid_full = np.full((len(grid_y), len(grid_x)), np.nan, dtype=np.float64) for (block_result, bx_min, bx_max, by_min, by_max) in results: if block_result is None: @@ -794,11 +801,10 @@ class ContentMapper: def _local_krige_block(self, local_pts, local_vals, sub_grid_x, sub_grid_y, bx_min, bx_max, by_min, by_max, - grid_dx, grid_dy, + n_closest=50, block_ix=0, block_iy=0, total=1): """单块局部克里金""" n_pts = len(local_pts) - n_cells = len(sub_grid_x) * len(sub_grid_y) if n_pts < 3: return None, bx_min, bx_max, by_min, by_max @@ -811,7 +817,7 @@ class ContentMapper: z, ss = ok.execute( 'grid', sub_grid_x, sub_grid_y, backend='loop', - n_closest_points=min(15, n_pts), + n_closest_points=min(n_closest, n_pts), ) return np.array(z), bx_min, bx_max, by_min, by_max @@ -3285,7 +3291,7 @@ def _local_krige_block_worker(args): """单个局部克里金块任务（独立进程入口，必须为模块级函数）""" (local_pts, local_vals, sub_grid_x, sub_grid_y, bx_min, bx_max, by_min, by_max, - grid_dx, grid_dy, block_ix, block_iy, total) = args + n_closest, block_ix, block_iy, total) = args import numpy as np from pykrige.ok import OrdinaryKriging @@ -3293,8 +3299,10 @@ def _local_krige_block_worker(args): if len(local_pts) < 3: return (None, bx_min, bx_max, by_min, by_max) - print(f" [LocalKrige] 块 ({block_iy},{block_ix}) [{block_iy * 100 + block_ix}/{total}] " - f"采样点={len(local_pts)}, 网格={len(sub_grid_x)}×{len(sub_grid_y)}") + print(f" [LocalKrige] 块 ({block_iy},{block_ix}) " + f"采样点={len(local_pts)}, " + f"网格={len(sub_grid_x)}×{len(sub_grid_y)}, " + f"n_closest={min(n_closest, len(local_pts))}") ok = OrdinaryKriging( local_pts[:, 0], local_pts[:, 1], local_vals, @@ -3305,6 +3313,6 @@ def _local_krige_block_worker(args): z, ss = ok.execute( 'grid', sub_grid_x, sub_grid_y, backend='loop', - n_closest_points=min(15, len(local_pts)), + n_closest_points=min(n_closest, len(local_pts)), ) return (np.array(z), bx_min, bx_max, by_min, by_max)