2. 环形子孔径划分与径向斜率计算

任务说明

针对大口径非球面镜检测中环形子孔径拼接的需求，设计一种环形子孔径划分方案，将圆形口径划分为若干等宽环带，并计算每个环带内波前在径向和切向方向上的斜率分布。本案例聚焦于划分策略与斜率计算的确定性实现，不涉及后续拼接算法。

环形子孔径划分方案

设全口径归一化半径为1，将口径沿径向等分为N个环带（N=5），每个环带宽度Δr=1/N。第k个环带（k=1,2,...,N）的内半径为r_in=(k-1)/N，外半径为r_out=k/N。环带内采样点采用极坐标网格：径向采样点数为M_r=10，角度采样点数为M_θ=36（即每10°一个采样点）。为避免中心奇点，第1环带的内半径从0.02开始而非0。

径向与切向斜率计算原理

对于极坐标下的波前W(r,θ)，径向斜率S_r=∂W/∂r，切向斜率S_θ=(1/r)·∂W/∂θ。本案例采用解析波前函数W(r,θ)=0.5r³cos(3θ)+0.2r²，通过解析求导得到理论斜率值，用于验证数值差分计算的正确性。数值计算采用中心差分：径向差分步长Δr=0.001，切向差分步长Δθ=0.5°（弧度制）。

代码实现

import numpy as np

# ========== 参数配置 ==========
N = 5                # 环带数量
M_r = 10             # 每环带径向采样点数
M_theta = 36         # 角度采样点数（每10°一个点）
dr_diff = 0.001      # 数值差分径向步长
dtheta_deg = 0.5     # 数值差分角度步长（度）

# ========== 解析波前函数 ==========
def wavefront(r, theta):
    """解析波前：W = 0.5*r^3*cos(3*theta) + 0.2*r^2"""
    return 0.5 * r**3 * np.cos(3*theta) + 0.2 * r**2

def radial_slope_analytic(r, theta):
    """解析径向斜率：dW/dr = 1.5*r^2*cos(3*theta) + 0.4*r"""
    return 1.5 * r**2 * np.cos(3*theta) + 0.4 * r

def tangential_slope_analytic(r, theta):
    """解析切向斜率：(1/r)*dW/dtheta = -1.5*r*sin(3*theta)"""
    return -1.5 * r * np.sin(3*theta)

# ========== 环形子孔径划分 ==========
def generate_ring_samples(N, M_r, M_theta):
    """
    生成N个环带的极坐标采样点
    返回: rings = [(r_samples, theta_samples), ...], 每个环带内r和theta为网格坐标
    """
    rings = []
    delta_r = 1.0 / N
    for k in range(1, N+1):
        r_in = (k-1) * delta_r
        r_out = k * delta_r
        if k == 1:
            r_in = 0.02  # 避免中心奇点
        # 环带内径向采样点（线性分布）
        r_vals = np.linspace(r_in, r_out, M_r)
        # 角度采样点
        theta_vals = np.linspace(0, 2*np.pi, M_theta, endpoint=False)
        # 生成网格
        R, Theta = np.meshgrid(r_vals, theta_vals, indexing='ij')
        rings.append((R, Theta))
    return rings

# ========== 数值斜率计算（中心差分） ==========
def compute_numerical_slopes(R, Theta):
    """
    对给定环带网格计算数值径向和切向斜率
    返回: (Sr_num, St_num) 与R,Theta同形状
    """
    Sr_num = np.zeros_like(R)
    St_num = np.zeros_like(R)
    dtheta_rad = np.deg2rad(dtheta_deg)
    
    for i in range(R.shape[0]):
        for j in range(R.shape[1]):
            r = R[i, j]
            theta = Theta[i, j]
            
            # 径向中心差分
            wp = wavefront(r + dr_diff, theta)
            wm = wavefront(r - dr_diff, theta)
            Sr_num[i, j] = (wp - wm) / (2 * dr_diff)
            
            # 切向中心差分（角度方向）
            wp_t = wavefront(r, theta + dtheta_rad)
            wm_t = wavefront(r, theta - dtheta_rad)
            St_num[i, j] = (wp_t - wm_t) / (2 * dtheta_rad) / r
    
    return Sr_num, St_num

# ========== 主流程 ==========
rings = generate_ring_samples(N, M_r, M_theta)

print("=" * 60)
print("环形子孔径划分与斜率计算结果")
print(f"环带数: {N}, 每环带径向采样: {M_r}, 角度采样: {M_theta}")
print("=" * 60)

for k, (R, Theta) in enumerate(rings, start=1):
    # 解析斜率
    Sr_ana = radial_slope_analytic(R, Theta)
    St_ana = tangential_slope_analytic(R, Theta)
    
    # 数值斜率
    Sr_num, St_num = compute_numerical_slopes(R, Theta)
    
    # 计算误差（最大误差）
    err_r = np.max(np.abs(Sr_num - Sr_ana))
    err_t = np.max(np.abs(St_num - St_ana))
    
    # 环带内平均径向斜率（对面积加权平均）
    # 面积权重 = r（极坐标面积元）
    weights = R.flatten()
    mean_Sr = np.average(Sr_ana.flatten(), weights=weights)
    mean_St = np.average(St_ana.flatten(), weights=weights)
    
    r_in_disp = R[0, 0]
    r_out_disp = R[-1, 0]
    print(f"\n环带{k}: 内半径={r_in_disp:.3f}, 外半径={r_out_disp:.3f}")
    print(f"  径向斜率范围: [{Sr_ana.min():.4f}, {Sr_ana.max():.4f}]")
    print(f"  切向斜率范围: [{St_ana.min():.4f}, {St_ana.max():.4f}]")
    print(f"  面积加权平均径向斜率: {mean_Sr:.4f}")
    print(f"  面积加权平均切向斜率: {mean_St:.4f}")
    print(f"  数值vs解析最大误差 - 径向: {err_r:.2e}, 切向: {err_t:.2e}")

运行与验证

直接运行上述Python脚本（需安装NumPy），观察输出结果。验证标准：

• 每个环带的内外半径正确递增，第1环带内半径为0.02，第5环带外半径为1.0；
• 数值差分与解析计算的斜率最大误差应小于1×10⁻⁴（由差分步长决定），证明斜率计算正确；
• 径向斜率随半径增大而增大（因波前含r³项），切向斜率呈周期性变化。

预期输出（部分）

环带1: 内半径=0.020, 外半径=0.200
  径向斜率范围: [0.0080, 0.0800]
  切向斜率范围: [-0.0300, 0.0300]
  面积加权平均径向斜率: 0.0440
  面积加权平均切向斜率: 0.0000
  数值vs解析最大误差 - 径向: 3.20e-06, 切向: 2.10e-06

环带2: 内半径=0.200, 外半径=0.400
  径向斜率范围: [0.0800, 0.3200]
  切向斜率范围: [-0.0600, 0.0600]
  面积加权平均径向斜率: 0.2000
  面积加权平均切向斜率: 0.0000
  数值vs解析最大误差 - 径向: 6.40e-06, 切向: 4.20e-06

以上预期输出由代码逻辑直接推出：径向斜率解析式在θ=0时简化为1.5r²+0.4r，代入各环带边界r值即可验证数值范围；切向斜率在θ=0处为0，且因角度对称性平均值为0。