变形镜的Matlab仿真实现

1.理论简介

        波前校正器对入射的畸变波前进行校正补偿畸变，是自适应光学系统的核心器件。常见的波前校正器一般有空间光调制器和变形镜两种。

        变形镜是通过改变反射面的形状对波前进行校正（即光程差中的传播路径，空间光调制器改变的是光程差中的折射率），反射的过程中能量损失较少，因此变形镜更适用于光强较低的情况。根据反射面的类型，可以将变形镜分为分段式和连续式两种。

        连续式DM由一个连续的反射镜薄膜和薄玻璃组成，镜面的背面与驱动器相连，以推、拉的方式使与其相连的薄镜面产生形变。连续式DM缺点是相邻驱动器间存在耦合（coupling），从而使致动器之间的AO控制复杂化。所谓“耦合”指的是某一个驱动器的运动作用到该位置处的镜面时，相邻位置处的驱动连接的镜面也会有一定的形变。

        分段式DM每个驱动器的末端由可独立控制的微小反射镜组成，每个镜段都可以单独控制，这类DM被认为更适合校正具有更高空间频率的像差。由于相邻镜子之间没有相互作用，故分段式DM具有零耦合的性质。

2.驱动器对应的响应函数的建模

        使用二维高斯函数描述变形镜第j个驱动器的响应函数fi(x, y)：

 其中，第i个驱动器的坐标是(xi, yi)，驱动器的中心是坐标原点，相邻驱动器间的归一化间距d0=1，高斯指数α=1.73，交连值ω=0.2。

function surf_DM= DM_surface(x,y,numActuators,cir,pointsize)
n           = size(x,1);
surf_DM     = zeros(n,n,numActuators*numActuators);
circle_mark = cir;
a           = 2;      % 高斯指数
d           = 2/(numActuators-1);% 驱动器间距
w           = 0.4;    % 交联值
% 驱动器响应函数的建模，对应上面的exp()公式
temp_num = 0;
for i = 1 : numActuators
    for j = 1 : numActuators
        CX = -1 + (i-1)*d;
        CY =  1 - (j-1)*d;     
        surf_DM(:,:,temp_num+1) = exp(log(w).*(sqrt((x - CX).^2.+(y - CY).^2)./d).^a);
        temp_num = temp_num+1;
    end
end

主函数，得到变形镜各个驱动器对应的面形。

pixofCCD     = 300;
x            = linspace(-1, 1, pixofCCD);
[X, Y]       = meshgrid(x, x);
Y            = -Y;
[theta,rho]  = cart2pol(X, Y);% 由(x,y)换算(r,theta)
Mask         = rho;
Mask(Mask<1) = 1; 
Mask(Mask~=1)= 0;
DM_num      = 11 ;              % 每一行驱动器的数目
DM_val_num  = DM_num * DM_num;  % 驱动器的总数目
pointsize   = pixofCCD / DM_num;
surf_DM     = DM_surface(X, Y, DM_num, Mask, pointsize);
figure(4)
for num_DM = 1 : DM_num * DM_num
    imagesc(surf_DM(:,:,num_DM));colormap(jet);colorbar
    drawnow;
end

某一个驱动器对应的面形如下。

ALpao 69驱动器的变形镜对应的响应函数如下：

 每个图像 就是 该位置处驱动器在单位电压下产生的面形。

或者使用两个贝塞尔函数对响应函数建模。

 响应函数即两者各取一部分、然后合并

clc;clear all;close all
P      = zeros(7,2);
P(1,:) = [0,1];      % 上凸部分的起点坐标
P(2,:) = [0.2,1];
P(3,:) = [0.4,0.7];
P(4,:) = [0.6,0.4];  % 上凸部分的终点坐标，也是下凹部分的起点坐标
P(5,:) = (-1)*P(3,:)+(1+1)*P(4,:);
P(6,:) = [1.0,0];
P(7,:) = [2.0,0];    % 下凹部分的终点坐标
t = linspace(0,1,101)';
bezier1 = ((1-t).^3)*P(1,:) + ...
    3.*((1-t).^2.*t)*P(2,:) + ...
    3.*((1-t).*t.^2)*P(3,:) + ...
    (t.^3)*P(4,:);   % 通过起点和终点坐标绘制贝塞尔曲线
figure(1);plot(bezier1(:,1),bezier1(:,2));title("上凸部分")
t(1) = [];
bezier2 = [((1-t).^3)*P(4,:) + ...
     3.*((1-t).^2.*t)*P(5,:) + ...
     3.*((1-t).*t.^2)*P(6,:) + ...
     (t.^3)*P(7,:) ]; % 绘制下凹部分的贝塞尔曲线
figure(2);plot(bezier2(:,1),bezier2(:,2));title("下凹部分")
bezierr = [bezier1;bezier2];

figure(3);
plot(bezierr(1:101,1),bezierr(1:101,2),'r');title("两条曲线衔接在一起，可导、连续")
hold on
plot(bezierr(102:201,1),bezierr(102:201,2),'b');title("两条曲线衔接在一起，可导、连续")

% plot(bezierr(:,1),bezierr(:,2));title("两条曲线衔接在一起，可导、连续")
p_xScale = spline(bezierr(:,2),bezierr(:,1),0.25); % obj.mechCoupling=25/100=0.25

bezierr(:,1) = bezierr(:,1)/p_xScale; % 横坐标扩大1.4倍数
u = [-flipud(bezierr(:,1)) ; bezierr(2:end,1)]; % flipud(X)实现了矩阵X的上下翻转
v = [ flipud(bezierr(:,2)) ; bezierr(2:end,2)];
figure(4);plot(u,v,'r','LineWidth',2);title("响应曲线")
splineP = spline(u,v);
xlabel('驱动器间距')

合并后的响应函数曲线如下

其三维视图为：

 仿真出来变形镜、哈特曼波前探测器，两者结合起来 就可以进行波前校正了。即哈特曼测波前的斜率，提高变换矩阵将测得的斜率转成变形镜上驱动器的电压信号，驱动器的电压信号在过渡到其对应的面形，最后将原始像差与变形镜拟合的面形作差得到残余波前，多次迭代实现波前校正的过程。