基于MATLAB：成像系统中的像差，zernike泽尔尼克多项式以及像差校正仿真

公众号：光学小太阳

在几何光学中，从给定物点发出并通过理想成像系统的光线会被聚焦到另一个点。物体的每个点发射（或反射）一个光学场，该光学场在成像系统的入瞳处变成发散的球面波。为了将光线聚焦到像平面上的一点，成像系统必须施加一个球面相位延迟，以将发散的球面波转换为会聚的球面波。

像差是相对于球面相位延迟的偏差，导致从给定物点发出的光线无法正确聚焦，形成有限大小的光斑。当整体观察图像时，像差表现为模糊。来自不同物点的光在像平面上可能会经历不同的像差，这取决于它们与光轴的距离。然而，我们假设这些与视场角相关的像差是恒定的。像差可以用波前𝑊(𝑥,𝑦)来表示，其单位是波长，或者用光学相位 𝜑(𝑥,𝑦)=2𝜋𝑊(𝑥,𝑦)来表示，其单位是弧度。然后，我们可以通过将光瞳孔径调制和像差的影响结合到一个复函数中，来写出一个广义的光瞳函数 𝑃(𝑥,𝑦)：

P(x,y)=P(x,y)exp(iφ(x,y))

通常将任意波前像差表示为多项式展开的形式，如下所示：

式中r和θ是极坐标系中的径向和角坐标，A0,A1,A2,…是多项式的系数，表示不同类型的像差。当像差变得复杂时，使用具有完备性和正交性的表示方法会更好。大多数情况下，我们处理的是圆形孔径，因此我们利用泽尼克圆多项式来表示：

m和n为非负整数，且m ≤ n。(n, m) 到 i 的映射关系比较复杂，部分泽尼克圆多项式的顺序在表1中给出。径向因子 Rnm(r) 和方位角因子 Gm(θ) 由以下公式给出

图1展示了三个不同阶的x方向像散。在图(a)中，n=2 且 m=2；在图(b)中， n=4 且m=2；在图(c)中，n=6 且 m=2。因此，这三个图都具有相同的方位角依赖性，即 cos(2θ)，而径向依赖性则各不相同。每级像散的最高幂次分别为2、4和6，分别对应主散光、次级散光和三级散光。当我们从光瞳中心到边缘跟踪每个模式的径向部分时，高阶模式具有更多的峰值、谷值和零交叉点。

相关matlab代码：

clear;clc;close;
% 不同阶的泽尔尼克多项式像散
%% 设置输入面网格
N = 512; % 网格数
L = 2; % 网格边长 [m]
d1 = L / N; % 网格间距 [m]
D=2;%圆孔径大小
[x1,y1] = meshgrid((-N/2 : N/2-1)*d1);%创建网格
[phi,r1] = cart2pol(x1,y1);
%i=6;12;24
z=circ(x1,y1,D).*zernike(6,r1,phi);
figure
imagesc(z);axis off ;colormap("parula");axis square

真实的像差W(r,θ) 可能包含大量的模式，但我们只对有限模式W′(r,θ)感兴趣，定义有限模式为：

则

这个框架适用于部分校正像差。自适应光学系统中波前传感器用于快速感知像差（有时每秒超过10,000帧），并调整可变形镜面的形状以补偿这些像差。如今许多天文望远镜使用自适应光学来补偿通过地球湍流大气层成像时引起的相位像差。可变形镜只能再现有限数量的泽尼克模式，因此总有一些残余像差是镜子无法校正的。图2给出了生成一个湍流像差的随机样本并产生一个模式受限的W′(r,θ) 的示例，显示了原始屏幕以及限制为3、15、30和100个模式的版本。我们发现随着泽尼克级数表示中包含的模式数量增加，模式受限版本越来越接近原始像差。

相关matlab代码：

% 生成一个湍流像差的随机样本
clear;clc;
N = 50; % 网格数
L = 2; % 网格长度 [m]
delta = L / N; % 网格间距 [m]
%创建网格
[x,y] = meshgrid((-N/2 : N/2-1)*delta);
[theta,r] = cart2pol(x, y);

% 创建圆形孔径索引
ap = circ(x, y, 2);
idxAp = logical(ap);

% 创建原始湍流屏
r0 = L / 20;
screen = ft_phase_screen(r0, N, delta, inf, 0)/ (2 * pi) .* ap;

%%% 利用泽尔尼克模式叠加像差
W = screen(idxAp); 
nModes = 3; % number of Zernike modes
% create matrix of Zernike polynomial values
Z = zeros(numel(W), nModes);
for idx = 1 : nModes
temp = zernike(idx, r, theta);
Z(:,idx) = temp(idxAp);
end
% compute mode coefficients
A = Z \ W;
% 得到有限模式像差
W_prime = Z * A;
% 得到有限模式像差图像
scr = zeros(N);
scr(idxAp) = W_prime;

%% 画图
figure
subplot(1,2,1)
%原始像差
imagesc(screen);axis off ;colormap("parula");axis square
subplot(1,2,2)
%有限模式像差
imagesc(scr);axis off ;colormap("parula");axis square

更多相关完整代码进t店：博士生牛马赚点小钱