4. 区域法波前重构之Fried模型

任务说明

本案例实现基于Fried模型的区域法波前重构。输入为离散采样点上的波前斜率（x方向与y方向梯度），输出为重构得到的波前相位分布。Fried模型假设斜率采样点位于网格单元中心，而重构相位定义在网格节点上，相邻四个斜率点围成一个相位节点，通过最小二乘拟合建立斜率与相位之间的线性关系。

算法原理

设波前相位定义在(M+1)×(N+1)的网格节点上，斜率测量点位于每个网格单元中心，共M×N个。对于节点(i,j)，其周围四个单元中心的x斜率平均值与该节点x方向相位差分相关，y方向同理。构建稀疏矩阵方程Ax = b，其中A为斜率-相位关系矩阵，b为斜率测量值，x为待求相位。由于A秩亏（整体平移不变），需添加约束固定一个参考节点相位为零，再通过稀疏最小二乘求解。

代码实现

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

def fried_reconstructor(sx, sy, dx=1.0, dy=1.0):
    """
    Fried模型区域法波前重构
    参数:
        sx: 二维数组 (M, N)，每个单元中心的x方向斜率
        sy: 二维数组 (M, N)，每个单元中心的y方向斜率
        dx, dy: 网格间距
    返回:
        phi: (M+1, N+1) 数组，重构的波前相位
    """
    M, N = sx.shape
    num_nodes = (M + 1) * (N + 1)
    num_slopes = 2 * M * N

    # 构建稀疏矩阵 A (num_slopes+1, num_nodes)
    A = lil_matrix((num_slopes + 1, num_nodes))
    b = np.zeros(num_slopes + 1)

    row = 0
    # x方向斜率方程: 节点(i,j)与(i,j+1)的差分
    for i in range(M):
        for j in range(N):
            n1 = i * (N + 1) + j          # 节点(i,j)
            n2 = i * (N + 1) + j + 1      # 节点(i,j+1)
            A[row, n1] = -1.0 / dx
            A[row, n2] = 1.0 / dx
            b[row] = sx[i, j]
            row += 1

    # y方向斜率方程: 节点(i,j)与(i+1,j)的差分
    for i in range(M):
        for j in range(N):
            n1 = i * (N + 1) + j          # 节点(i,j)
            n2 = (i + 1) * (N + 1) + j    # 节点(i+1,j)
            A[row, n1] = -1.0 / dy
            A[row, n2] = 1.0 / dy
            b[row] = sy[i, j]
            row += 1

    # 固定参考节点: 左上角相位 = 0
    A[row, 0] = 1.0
    b[row] = 0.0

    # 最小二乘求解
    A = A.tocsr()
    phi_flat, _ = lsqr(A, b, iter_lim=1000, atol=1e-10, btol=1e-10)[:2]
    phi = phi_flat.reshape((M + 1, N + 1))
    return phi

# 示例：构造已知相位并生成斜率，验证重构精度
if __name__ == "__main__":
    M, N = 8, 8
    x = np.linspace(0, 1, N + 1)
    y = np.linspace(0, 1, M + 1)
    X, Y = np.meshgrid(x, y)

    # 真实相位（用于生成斜率）
    phi_true = 0.5 * X**2 + 0.3 * Y**2 + 0.2 * X * Y

    # 计算单元中心坐标
    xc = (x[:-1] + x[1:]) / 2
    yc = (y[:-1] + y[1:]) / 2
    Xc, Yc = np.meshgrid(xc, yc)

    # 解析求导得到斜率
    sx_true = Xc + 0.2 * Yc
    sy_true = Yc + 0.2 * Xc

    # 重构
    phi_rec = fried_reconstructor(sx_true, sy_true, dx=x[1]-x[0], dy=y[1]-y[0])

    # 去除整体平移后比较
    phi_rec -= phi_rec[0, 0]
    phi_true_ref = phi_true - phi_true[0, 0]
    error = np.max(np.abs(phi_rec - phi_true_ref))
    print(f"最大重构误差: {error:.2e}")

运行与验证

运行上述代码需要安装numpy和scipy。执行后程序会构造一个二次型相位分布，解析计算其斜率作为输入，调用fried_reconstructor完成重构，并输出最大误差。由于斜率由解析函数精确生成且无噪声，重构结果应接近机器精度（误差在1e-10量级），这验证了重构矩阵构建和求解过程的正确性。若误差显著大于1e-8，说明矩阵构建或约束添加存在错误。

实际使用时，将波前传感器测得的x、y方向斜率分别填入sx和sy二维数组即可。注意斜率数组维度为(M,N)，对应M行N列个采样单元；重构相位维度自动为(M+1,N+1)。若斜率存在噪声，重构误差会相应增大，这是最小二乘方法的正常特性。