3. 直接斜率法波前重构

任务定义

本案例实现基于 Shack-Hartmann 波前传感器测量数据的直接斜率法波前重构。输入为子孔径内波前斜率（x、y 方向），输出为离散采样点上的波前相位分布。采用 zonal 重构策略：将波前相位建模为子孔径网格节点上的未知量，利用相邻节点间的相位差与测量斜率之间的线性关系建立方程组，通过最小二乘求解得到波前相位。

算法原理与实现

设子孔径按 M×N 网格排列，节点数为 (M+1)×(N+1)。x 方向斜率测量值对应水平相邻节点间的相位差除以子孔径间距 d，y 方向同理。构建稀疏矩阵 A（尺寸 2MN × (M+1)(N+1)），每行对应一个斜率测量方程，未知量为各节点相位。由于整体平移不影响斜率，需固定一个参考节点相位为 0 以保证解唯一。采用稀疏最小二乘求解。

import numpy as np
from scipy.sparse import lil_matrix
from scipy.sparse.linalg import lsqr

def direct_slope_reconstruction(slope_x, slope_y, d=1.0):
    """
    直接斜率法波前重构 (zonal)
    参数:
        slope_x: 形状 (M, N) 的 x方向斜率 (每子孔径)
        slope_y: 形状 (M, N) 的 y方向斜率
        d: 子孔径间距
    返回:
        phase: 形状 (M+1, N+1) 的节点相位
    """
    M, N = slope_x.shape
    n_nodes = (M+1) * (N+1)
    n_eq = 2 * M * N + 1  # +1 为参考节点约束

    A = lil_matrix((n_eq, n_nodes))
    b = np.zeros(n_eq)

    row = 0
    # x方向斜率: 水平相邻节点差 / d = slope_x
    for i in range(M):
        for j in range(N):
            node_left = i * (N+1) + j
            node_right = i * (N+1) + (j+1)
            A[row, node_left] = -1.0 / d
            A[row, node_right] = 1.0 / d
            b[row] = slope_x[i, j]
            row += 1

    # y方向斜率: 垂直相邻节点差 / d = slope_y
    for i in range(M):
        for j in range(N):
            node_bottom = i * (N+1) + j
            node_top = (i+1) * (N+1) + j
            A[row, node_bottom] = -1.0 / d
            A[row, node_top] = 1.0 / d
            b[row] = slope_y[i, j]
            row += 1

    # 固定参考节点 (0,0) 相位为 0
    A[row, 0] = 1.0
    b[row] = 0.0

    A = A.tocsr()
    phase_flat, _ = lsqr(A, b, atol=1e-10, btol=1e-10)[:2]
    return phase_flat.reshape((M+1, N+1))

测试与验证

构造已知相位分布 phi(x,y) = 0.5x² + 0.3y²，在 4×4 子孔径网格上解析计算各子孔径中心处的斜率（中心差分），作为输入进行重构，并与原始相位比较。

# 生成测试数据
M, N = 4, 4
d = 1.0
x = np.linspace(-2, 2, N+1)
y = np.linspace(-2, 2, M+1)
X, Y = np.meshgrid(x, y)
true_phase = 0.5*X**2 + 0.3*Y**2

# 子孔径中心坐标
xc = (x[:-1] + x[1:]) / 2
yc = (y[:-1] + y[1:]) / 2
Xc, Yc = np.meshgrid(xc, yc)

# 解析斜率
slope_x = Xc          # d/dx (0.5x^2) = x
slope_y = 0.6 * Yc    # d/dy (0.3y^2) = 0.6y

# 重构
recon_phase = direct_slope_reconstruction(slope_x, slope_y, d)

# 误差评估（去除整体平移）
recon_phase = recon_phase - recon_phase[0,0]
true_phase = true_phase - true_phase[0,0]
rmse = np.sqrt(np.mean((recon_phase - true_phase)**2))
print(f"重构 RMSE: {rmse:.2e}")

预期输出：RMSE 在 1e-10 量级（因解析斜率与离散差分模型完全匹配，重构精度仅受数值求解器容差限制）。

验证方法

运行上述代码，观察终端输出的 RMSE 值。若 RMSE < 1e-8，说明重构算法正确实现了斜率到相位的映射；若误差较大，检查子孔径索引映射与符号约定是否正确。

关键说明

• 本实现适用于子孔径斜率数据完整、无缺失的情况；若存在坏点需先插值或引入加权处理。
• 参考节点选择 (0,0) 处节点，消除重构波前的整体平移不确定性，不影响波前形状。
• 矩阵规模为 (2MN+1) × (M+1)(N+1)，对典型 20×20 子孔径约 800×441，稀疏求解效率高。