使用 SciPy 实现 NumPy 数组的重叠拼接与加权融合

本文介绍一种基于 `scipy.linalg.block_diag` 的通用方法，将两个二维数组按指定重叠宽度进行对齐拼接，并对重叠区域元素取平均值；支持稀疏结构延展，兼顾内存效率与数值精度。

在科学计算与信号处理中，常需将多个局部数据块（如分块图像、时序片段或有限元子域）沿某一维度拼接，且要求在交界处平滑过渡——而非简单留零或硬截断。scipy.linalg.block_diag 虽能高效构建块对角矩阵，但默认不支持重叠融合。本文提供一种简洁、可扩展、数值稳健的实现方案，核心思想是：通过零填充构造同维数的“占位矩阵”，利用线性叠加与掩模归一化，自然导出重叠平均结果。

基本原理与实现步骤

给定两个形状均为 (n, n) 的方阵 A 和 B，设定重叠宽度为 overlap（整数，需满足 1 ≤ overlap < n），目标输出矩阵尺寸为 (n + n – overlap, n + n – overlap)。关键在于：

将 A “右下对齐”放置于大矩阵左上角区域；将 B “左上对齐”放置于大矩阵右下角区域；重叠区域（大小为 overlap × overlap）内，每个位置的值为 A 与 B 对应位置的算术平均。

该逻辑可通过两次 block_diag 构造、一次逐元素加法与一次条件归一化完成，无需显式循环或全量初始化大数组。

以下为完整可运行示例（适配任意方阵及重叠宽度）：

import numpy as npfrom scipy.linalg import block_diagdef overlap_stitch(A, B, overlap): """ 拼接两个同形方阵 A 和 B，使其在右下/左上方向重叠 ‘overlap’ 行列，并对重叠区取平均。 Parameters ———- A, B : ndarray, shape (n, n) 输入方阵，必须维度相同 overlap : int 重叠宽度（行列数），1 <= overlap < A.shape[0] Returns ——- out : ndarray, shape (2*n – overlap, 2*n – overlap) 重叠拼接后的矩阵 """ n = A.shape[0] if A.shape != B.shape or overlap <= 0 or overlap >= n: raise ValueError("A and B must be square and identical; overlap must satisfy 0 < overlap < n") # 计算非重叠填充尺寸 pad_A = n – overlap # A 在右侧/下方需补零的行/列数 pad_B = n – overlap # B 在左侧/上方需补零的行/列数 # 构造零填充块：使 A 和 B 可对齐到同一目标尺寸 ZerosA = np.zeros((pad_A, pad_A)) ZerosB = np.zeros((pad_B, pad_B)) # 分别生成带零填充的扩展矩阵（A 左上，B 右下） # block_diag(A, ZerosA) → shape (n+pad_A, n+pad_A) = (2n-overlap, 2n-overlap) # block_diag(ZerosB, B) → 同样尺寸，但 B 位于右下 Sum = block_diag(A, ZerosA) + block_diag(ZerosB, B) # 构建计数掩模：每个位置记录有多少个原始矩阵在此处贡献了非零值 Denom = block_diag(np.ones_like(A), ZerosA) + block_diag(ZerosB, np.ones_like(B)) # 安全除法：仅在 Denom > 0 处执行 Sum/Denom，其余置 0 return np.where(Denom > 0, Sum / Denom, 0)# 示例验证A = np.linspace(1, 9, 9).reshape(3, 3)B = np.linspace(10, 18, 9).reshape(3, 3)result = overlap_stitch(A, B, overlap=2)print(result)

输出：

[[ 1. 2. 3. 0. ] [ 4. 7.5 8.5 12. ] [ 7. 10.5 11.5 15. ] [ 0. 16. 17. 18. ]]

注意事项与进阶建议

✅ 内存友好性：本方法仅依赖 block_diag 与逐元素运算，天然兼容 scipy.sparse 矩阵。只需将 np.zeros 替换为 scipy.sparse.csr_matrix((pad_A, pad_A))，并确保 A, B 为稀疏格式，即可获得稀疏输出（需注意 np.where 不直接支持稀疏索引，此时推荐改用 sparse.diags 或手动构造 coo_matrix）。⚠️ 维度限制：当前实现假设 A 和 B 为同形方阵。若需处理矩形阵列（如 (m,n) 与 (p,q)），需分别指定水平与垂直重叠量，并调整 block_diag 的填充策略。? 加权融合扩展：若需非等权平均（如指数衰减权重），可将 Denom 替换为自定义权重矩阵 W，并将 Sum 改为加权和 wA * A_extended + wB * B_extended。? 重叠定位灵活性：当前为“右下–左上”标准重叠；如需其他对齐方式（如中心对齐、顶部对齐），可通过 np.pad 配合切片替代 block_diag 实现更细粒度控制。

综上，该方法以清晰的线性代数视角重构了重叠拼接问题，在保持代码简洁性的同时，提供了良好的可读性、可维护性与工程扩展潜力。