CN103714516A - 一种自适应光学成像中的点扩散函数估计方法 - Google Patents

Abstract

本发明公开了一种自适应光学成像中的点扩散函数估计方法。该发明基于研究多元样条的B网方法，首先将经过大气降质的畸变波前近似表示为单纯形剖分上的二元样条，进而得到波前传感器斜率测量数据与样条函数系数之间的数学关系；其次，结合相邻三角剖分上样条函数的光滑拼接条件，利用斜率测量数据通过约束最小二乘估计估计样条函数系数进而重构畸变波前；最后将得到的畸变波前结合透镜光瞳函数，即可估计出经大气湍流成像的瞬时点扩散函数。本发明与传统的Zernike模式法估计点扩散函数相比克服了变换矩阵秩不完备的缺陷，提高了估计精度；并且适用于中心被遮挡的波前传感器及微透镜阵列的多种几何排布的情况。

Description

一种自适应光学成像中的点扩散函数估计方法

技术领域

本发明属于空间目标探测识别技术领域，涉及到空间目标自适应光学图像复原理论与方法，具体涉及到一种自适应光学成像中的点扩散函数估计方法。

背景技术

根据空间目标光学成像探测的特性，自适应光学图像复原方法主要可分为两类：盲解卷积和波前解卷积。对于盲解卷积，一个相对准确的点扩散函数估计值作为迭代算法的初始估计将会加快收敛、减少计算量；对于波前解卷积，点扩散函数估计值与真实值之间的偏差更是会直接影响复原的效果。因此快速准确的点扩散函数估计方法对于自适应光学图像的复原来讲意义显著。

点扩散函数与波前存在着理论上的对应关系，在自适应光学***中用波前传感器探测波前信息。目前使用H-S型波前传感器的自适应光学***居多，所获得数据是波前相位的斜率采样值。通过波前重构可以在有噪声及干扰的情况下，利用有限离散观测信息恢复出完整的波前，进而重构短曝光图像的瞬时点扩散函数。波前重构的基本方法大致上可以分为两大类：区域法和模式法。由于Zernike多项式各模式的像差意义明显，故Zernike模式法在面向图像复原的波前重构方法中一直得到广泛的应用。Zernike模式法最终归结为求解线性方程组，但由于Zernike函数偏导数的不完全正交性以及在有限采样点上函数的非正交性，在求最小二乘解时会出现变换矩阵D的秩不完备、D^ΤD奇异或接近奇异的情况，这将使梯度矢量的测量误差被放大，导致波前重构结果精度降低，进而影响点扩散函数的估计精度；其次，许多大型望远镜的中心常常被遮挡，若仍用Zernike多项式展开波前，因为展开域是环形的，所以Zernike展开不再具有正交性，仍应用此法则需要十分繁琐的数学推导与计算，这些因素都限制了Zernike模式法的应用效力。另外，虽然近年来在大型自适应光学设备快速发展的推动下不断涌现出新的面向校正控制的波前重构算法，但把这些方法进一步的应用到点扩散函数估计直至图像复原中的效力还并未得到有效的研究和验证。

发明内容

本发明的目的在于克服现有成像点扩散函数估计方法的不足，提出一种自适应光学成像中的点扩散函数估计方法，该发明基于研究多元样条的B网方法，将经过大气降质的畸变波前近似表示为单纯形剖分上的二元样条，避免了求最小二乘解时不完备或奇异矩阵的出现，提高了估计精度，并可以应对多种波前传感器微透镜阵列的几何排布方式。

为了达到上述目的，本发明所采用的技术方案包括以下步骤：

1）指定样条函数的次数d及刻画光滑程度的连续性阶数r，确定波前区域内三角剖分的形式；

2）建立斜率测量模型的偏微分方程并确定方程中各参数矩阵的值；

3）根据具体的参数值确定变换矩阵D；

4）结合相邻样条函数间的光滑拼接条件确定重构矩阵R；

5）利用斜率测量数据计算各二元样条函数的系数向量

；

6）根据系数

重构出畸变波前φ(x,y)；

7）由重构的降质波前去估计大气瞬时点扩散函数h(x,y)。

所述步骤1）中，根据微透镜阵列排布方式确定三角剖分的形式，确定三角剖分的形式时采用基于几何布局的网格单元直接剖分方法。

所述步骤2）中，斜率测量模型的偏微分方程的建立方法具体为：

首先将波前表示为：

$\mathrm{\φ}\left(x\right)\≈s_d^r\left(b\left(x\right)\right)=\underset{\left|\mathrm{\λ}\right|=d}{\mathrm{\Σ}}{c}_{\mathrm{\λ}}{B}_{\mathrm{\λ}}^d\left(b\left(x\right)\right)=B^d\left(b\left(x\right)\right)\·c$

其m阶方向导数

可表示为 $D_u^m{s}_d^r\left(b\left(x\right)\right)=\frac{d!}{(d-m)!}{B}^{d-m}{P}_u^{d,d-m}\left(a\right)\·c,$ 其中，

为二元单纯形样条，d为样条函数次数，r为连续性阶数，b(x)为一个三角剖分上笛卡尔坐标到面积坐标的变换，

为d次Bernstein多项式基函数，λ=(λ₀,λ₁,λ₂)为多重符号且其长度为|λ|=λ₀+λ₁+λ₂=d，c_λ为多项式系数，B^d(b(x))和c分别为

和c_λ的矩阵形式，u为笛卡尔坐标下方向导数的方向向量，a为与u相对应的面积坐标下的方向坐标，

为分块对角矩阵，每对角块为对应于a和u的de Casteljau矩阵；其次根据相邻三角剖分间的r阶光滑拼接条件，每条边都有

个约束条件，从而构成约束方程Hc=0，H为由E条边构成的E*Z行约束矩阵，且H的每行都包含一个单独的连续性约束，这样，波前传感器的斜率测量偏微分方程即可表示为：

${\mathrm{\σ}}_u=d{B}^{d-1}{P}_u^{d,d-1}c+n$

Hc=0

其中，σ_u表示波前传感器斜率测量矩阵，n表示测量噪声及模型误差。

相邻三角剖分间的r阶光滑拼接条件描述如下：

设

分别是定义在相邻三角形t=<ν₁,ν₂,ν₃>和

上的d次多项式，ν_i(i=1,2,3)为三角形的顶点，c_λ和

分别是S(b(x))和

关于t和

的贝塞尔坐标，则S(b(x))和

之间的r阶光滑拼接条件为：

${\stackrel{\&OverBar;}{c}}_{{\mathrm{\&lambda;}}^s}=c_{{\mathrm{\&lambda;}}^0}^{\left(s\right)}\left(\stackrel{\&OverBar;}{\mathrm{\&tau;}}\right),s=\mathrm{0,1},\&CenterDot;\&CenterDot;\&CenterDot;,r$

其中，

是关于t的面积坐标，λ^s=(s,λ₁,λ₂)，λ⁰=(0,λ₁,λ₂)，λ₁+λ₂=d-s。

所述步骤3）中，确定变换矩阵D的方法具体为：

步骤2）建立了一个约束最小二乘问题σ_u=Dc+n,Hc=0；与Zernike模式法相对应，通过对比此约束最小二乘问题与斜率测量偏微分方程，即可将变换矩阵D表示为：

$D=d{B}^{d-1}{P}_u^{d,d-1}.$

所述步骤4）中，确定重构矩阵R的方法具体为：

将由斜率测量数据去重构样条函数系数向量的矩阵定义为重构矩阵，结合步骤1）和步骤2），通过求解约束最小二乘问题，多项式系数向量c的显式解表达式为：

$\hat{c}=[{\left(D^{T}D\right)}^{-1}{D}^T-{\left(D^{T}D\right)}^{-1}{H}^T{\left[H{\left(D^{T}D\right)}^{-1}{H}^T\right]}^{-1}{\left(D^{T}D\right)}^{-1}{D}^T]\mathrm{\&sigma;}$

重构矩阵R，则表示为：

R=(D^ΤD)^-1D^Τ-(D^ΤD)^-1H^T[H(D^ΤD)^-1H^T]^-1(D^ΤD)^-1D^Τ

其中，Τ表示矩阵转置，σ为波前斜率测量向量，此时由重构矩阵所确定的样条函数系数是在约束Hc=0下的最优线性无偏估计。

所述步骤5）中，根据步骤4），重构矩阵是由斜率测量数据去重构样条函数系数向量的矩阵，二元样条函数的系数向量直接由重构矩阵与波前斜率测量向量相乘得到，则其计算公式为：

$\hat{c}=\mathrm{R\&sigma;}.$

所述步骤6）中，畸变波前φ(x,y)的重构方法为：

根据步骤2），已将畸变波前表示为单纯形剖分上的二元样条，步骤5）又得到了样条函数的系数向量而d次Bernstein多项式基函数B^d(b(x))是确定的，将二者相乘即可得到畸变波前的估计量

即：

$\widehat{\mathrm{\&phi;}}\left(x\right)=B^d\left(b\left(x\right)\right)\hat{c}.$

所述步骤7）中，估计大气瞬时点扩散函数h(x,y)的方法具体为：

已知透镜光瞳函数为A，根据步骤6），由傅里叶光学理论，结合透镜光瞳函数A，得到其波前对应的点扩散函数为：

$h(x,y)={\left|{\mathrm{FT}}^{-1}\right[\mathrm{Aexp}\left(j\widehat{\mathrm{\&phi;}}\left(x\right)\right)\left]\right|}^2$

其中，FT^-1表示傅里叶逆变换，

与现有技术相比，本发明具有以下有益效果：

本发明在重构波前时满足变换矩阵D满秩、D^ΤD非奇异的条件，不会放大斜率测量误差，并且基于二元样条的方法可以选择使用非线性基函数，这都促使波前重构精度的提高，进而可以提高点扩散函数的估计精度；另外，本发明不受波前传感器微透镜阵列几何排布方式的限制，对圆形或方形布局都有效，并且还适用于中心被遮挡的情况，方法的通用性较强。

附图说明

图1为本发明的原理示意图；

图2为本发明的具体流程图；

图3为本发明一种典型54子孔径波前传感器微透镜阵列排布图；

图4为本发明针对图3波前传感器做的一种三角剖分示例图。

具体实施方式

下面结合附图和实施例对本发明作进一步详细的描述：

参见图1和图2，本发明包括以下步骤：

1）指定样条函数的次数d及刻画光滑程度的连续性阶数r，确定波前区域内三角剖分的形式；

根据微透镜阵列排布方式确定三角剖分的形式，确定三角剖分的形式时采用基于几何布局的网格单元直接剖分方法。

2）建立斜率测量模型的偏微分方程并确定方程中各参数矩阵的值；

斜率测量模型的偏微分方程的建立方法具体为：

首先将波前表示为：

$\mathrm{\&phi;}\left(x\right)\&ap;s_d^r\left(b\left(x\right)\right)=\underset{\left|\mathrm{\&lambda;}\right|=d}{\mathrm{\&Sigma;}}{c}_{\mathrm{\&lambda;}}{B}_{\mathrm{\&lambda;}}^d\left(b\left(x\right)\right)=B^d\left(b\left(x\right)\right)\&CenterDot;c$

其m阶方向导数可表示为 $D_u^m{s}_d^r\left(b\left(x\right)\right)=\frac{d!}{(d-m)!}{B}^{d-m}{P}_u^{d,d-m}\left(a\right)\&CenterDot;c,$ 其中，

为二元单纯形样条，d为样条函数次数，r为连续性阶数，b(x)为一个三角剖分上笛卡尔坐标到面积坐标的变换，为d次Bernstein多项式基函数，λ=(λ₀,λ₁,λ₂)为多重符号且其长度为|λ|=λ₀+λ₁+λ₂=d，c_λ为多项式系数，B^d(b(x))和c分别为

和c_λ的矩阵形式，u为笛卡尔坐标下方向导数的方向向量，a为与u相对应的面积坐标下的方向坐标，

为分块对角矩阵，每对角块为对应于a和u的de Casteljau矩阵；其次根据相邻三角剖分间的r阶光滑拼接条件，每条边都有

个约束条件，从而构成约束方程Hc=0，H为由E条边构成的E?Z行约束矩阵，且H的每行都包含一个单独的连续性约束，这样，波前传感器的斜率测量偏微分方程即可表示为：

${\mathrm{\&sigma;}}_u=d{B}^{d-1}{P}_u^{d,d-1}c+n$

Hc=0

其中，σ_u表示波前传感器斜率测量矩阵，n表示测量噪声及模型误差。

相邻三角剖分间的r阶光滑拼接条件描述如下：

设S(b(x))和

分别是定义在相邻三角形t=<ν₁,ν₂,ν₃>和

上的d次多项式，ν_i(i=1,2,3)为三角形的顶点，c_λ和分别是S(b(x))和

关于t和

的贝塞尔坐标，则S(b(x))和

之间的r阶光滑拼接条件为：

${\stackrel{\&OverBar;}{c}}_{{\mathrm{\&lambda;}}^s}=c_{{\mathrm{\&lambda;}}^0}^{\left(s\right)}\left(\stackrel{\&OverBar;}{\mathrm{\&tau;}}\right),s=\mathrm{0,1},\&CenterDot;\&CenterDot;\&CenterDot;,r$

其中，

是关于t的面积坐标，λ^s=(s,λ₁,λ₂)，λ⁰=(0,λ₁,λ₂)，λ₁+λ₂=d-s。

3）根据具体的参数值确定变换矩阵D；

确定变换矩阵D的方法具体为：

步骤2）建立了一个约束最小二乘问题σ_u=Dc+n,Hc=0；与Zernike模式法相对应，通过对比此约束最小二乘问题与斜率测量偏微分方程，即可将变换矩阵D表示为：

$D=d{B}^{d-1}{P}_u^{d,d-1}.$

4）结合相邻样条函数间的光滑拼接条件确定重构矩阵R；

确定重构矩阵R的方法具体为：

将由斜率测量数据去重构样条函数系数向量的矩阵定义为重构矩阵，结合步骤1）和步骤2），通过求解约束最小二乘问题，多项式系数向量c的显式解表达式为：

$\hat{c}=[{\left(D^{T}D\right)}^{-1}{D}^T-{\left(D^{T}D\right)}^{-1}{H}^T{\left[H{\left(D^{T}D\right)}^{-1}{H}^T\right]}^{-1}{\left(D^{T}D\right)}^{-1}{D}^T]\mathrm{\&sigma;}$

这样重构矩阵R即可表示为：

R=(D^ΤD)^-1D^Τ-(D^ΤD)^-1H^T[H(D^ΤD)^-1H^T]^-1(D^ΤD)^-1D^Τ

其中，Τ表示矩阵转置，σ为波前斜率测量向量，此时由重构矩阵所确定的样条函数系数是在约束Hc=0下的最优线性无偏估计。

5）利用斜率测量数据计算各二元样条函数的系数向量?c；

根据步骤4），重构矩阵是由斜率测量数据去重构样条函数系数向量的矩阵，二元样条函数的系数向量

可直接由重构矩阵与波前斜率测量向量相乘得到，则其计算公式为：

$\hat{c}=\mathrm{R\&sigma;}.$

6）根据系数

重构出畸变波前φ(x,y)；

畸变波前φ(x,y)的重构方法为：

根据步骤2），已将畸变波前表示为单纯形剖分上的二元样条，步骤5）又得到了样条函数的系数向量

而d次Bernstein多项式基函数B^d(b(x))是确定的，这样将二者相乘即可得到畸变波前的估计量

即：

$\widehat{\mathrm{\&phi;}}\left(x\right)=B^d\left(b\left(x\right)\right)\hat{c}.$

7）由重构的降质波前去估计大气瞬时点扩散函数h(x,y)。

估计大气瞬时点扩散函数h(x,y)的方法具体为：

已知透镜光瞳函数为A，根据步骤6），由傅里叶光学理论，结合透镜光瞳函数A，则其波前对应的点扩散函数为：

$h(x,y)={\left|{\mathrm{FT}}^{-1}\right[\mathrm{Aexp}\left(j\widehat{\mathrm{\&phi;}}\left(x\right)\right)\left]\right|}^2$

其中，FT^-1表示傅里叶逆变换，

本发明的原理在于：

本发明基于研究多元样条的B网方法，首先将经过大气降质的畸变波前近似表示为单纯形剖分上的二元样条，进而得到波前传感器斜率测量数据与样条函数系数之间的数学关系；其次，结合相邻三角剖分上样条函数的光滑拼接条件，利用斜率测量数据由约束最小二乘估计样条函数系数进而重构畸变波前；最后，畸变波前结合透镜光瞳函数，即可估计出大气湍流的成像点扩散函数。

实施例：

为使本发明的目的、技术方案、优点等更加清楚明白，以下结合具体实例，并参照附图，对本发明做进一步详细说明。本发明是一种基于多元样条的自适应光学成像点扩散函数估计方法，图2为本发明估计点扩散函数的具体流程，以下分别做详细介绍：

1）初始化阶段。依次为，首先依据波前传感器微透镜阵列的布局方式确定波前单纯形区域上三角剖分的形式，这里以54子孔径波前传感器为例，如图3所示，此正六边形微透镜阵列的布局方式本质上是三角的，直观上即可以依据不同的样条函数次数确定出几种三角剖分，不需要复杂的剖分算法，但对于此种波前传感器来讲，传统的基于有限差分的波前重构方法将失效，并且由于其中心被遮挡，Zernike模式法更是失去了正交性，使用时需要十分繁琐的数学修正，此时本发明的优点就显现出来；

其次确定二元样条函数的次数d及相邻三角剖分间光滑拼接的连续性阶数r，如图4所示，针对图3的布局，54个小圆圈表示微透镜阵列的位置，粗实线区域被剖分为3个三角形，整个环形波前区域共分为6个粗实线区域18个三角形，确定样条函数的次数d=2，连续性阶数r=1；

再次，依据以上的描述建立波前斜率测量偏微分方程并确定方程参数矩阵，波前近似表示为：

$\mathrm{\&phi;}\left(x\right)\&ap;s_d^r\left(b\left(x\right)\right)=\underset{\left|\mathrm{\&lambda;}\right|=d}{\mathrm{\&Sigma;}}{c}_{\mathrm{\&lambda;}}{B}_{\mathrm{\&lambda;}}^d\left(b\left(x\right)\right)=B^d\left(b\left(x\right)\right)\&CenterDot;c---\left(1\right)$

其m阶方向导数为 $D_u^m{s}_d^r\left(b\left(x\right)\right)=\frac{d!}{(d-m)!}{B}^{d-m}{P}_u^{d,d-m}\left(a\right)\&CenterDot;c,$ 其中

为d次Bernstein多项式基函数，c_λ为多项式系数，相应的粗体为其矩阵形式，u为笛卡尔坐标下方向导数的方向向量，a为与u相对应的面积坐标下的方向坐标，为分块对角矩阵，每对角块为对应于a和u的de Casteljau矩阵；其次根据相邻三角剖分间的r阶光滑拼接条件

每条边都有个约束条件，从而构成约束方程Hc=0，H为由E条边构成的E*Z行的约束矩阵。这样，波前传感器的斜率测量偏微分方程即可表示为：

${\mathrm{\&sigma;}}_u=d{B}^{d-1}{P}_u^{d,d-1}c+n---\left(2\right)$

Hc=0

其中σ_u表示波前传感器斜率测量值，n表示测量噪声及模型误差。模型建立好其相应的B^d-1、与H阵也已随即确定；

2）确定变换矩阵D。与Zernike模式法类似，在此也将样条函数系数与斜率测量数据之间的变换关系用变换矩阵D来描述。约束最小二乘问题的模型为σ_u=Dc+n；Hc=0，结合步骤1）则有：

$D=d{B}^{d-1}{P}_u^{d,d-1}---\left(3\right)$

与Zernike模式法不同的是这里的D^ΤD是非奇异的，这是由于

的构成方式以及加入了剖分之间光滑拼接条件约束所直接决定的，从而不会放大斜率数据的测量噪声。

3）计算重构矩阵R。将由斜率测量数据去重构样条函数系数向量的矩阵定义为重构矩阵。结合步骤1）和步骤2），由约束最小二乘问题的显式解表达式：

$\hat{c}=[{\left(D^{T}D\right)}^{-1}{D}^T-{\left(D^{T}D\right)}^{-1}{H}^T{\left[H{\left(D^{T}D\right)}^{-1}{H}^T\right]}^{-1}{\left(D^{T}D\right)}^{-1}{D}^T]\mathrm{\&sigma;},$ 从而有：

R=(D^ΤD)^-1D^Τ-(D^ΤD)^-1H^T[H(D^ΤD)^-1H^T]^-1(D^ΤD)^-1D^Τ   （4）

此时由重构矩阵所确定的样条函数系数是在约束Hc=0下的最优线性无偏估计。

4）计算样条函数系数向量

由步骤3）可得：

$\hat{c}=\mathrm{R\&sigma;}---\left(5\right)$

5）重构畸变波前φ(x,y)。由步骤1）：

$\widehat{\mathrm{\&phi;}}\left(x\right)=B^d\left(b\left(x\right)\right)\hat{c}---\left(6\right)$

6）估计点扩散函数h(x,y)。在步骤5）的基础上，由傅里叶光学理论，结合透镜光瞳函数A，其波前对应的点扩散函数为：

$h(x,y)={\left|{\mathrm{FT}}^{-1}\right[\mathrm{Aexp}\left(j\widehat{\mathrm{\&phi;}}\left(x\right)\right)\left]\right|}^2---\left(7\right)$

其中FT^-1表示傅里叶逆变换。至此，则完成了成像点扩散函数的估计。

以上实施例仅用以说明本发明而并非限制本发明所描述的技术方案；因此尽管本说明书参照上述的实施例对本发明已进行了详细的说明，但是本领域的普通技术人员应当理解，仍然可以对本发明进行修改或等同替换；而一切不脱离发明的精神和范围的技术方案及其改进，均应涵盖在本发明的权利要求范围之内。

Claims (9)

1.一种自适应光学成像中的点扩散函数估计方法，其特征在于，包括以下步骤：

1）指定样条函数的次数d及刻画光滑程度的连续性阶数r，确定波前区域内三角剖分的形式；

2）建立斜率测量模型的偏微分方程并确定方程中各参数矩阵的值；

3）根据具体的参数值确定变换矩阵D；

4）结合相邻样条函数间的光滑拼接条件确定重构矩阵R；

5）利用斜率测量数据计算各二元样条函数的系数向量

；

6）根据系数重构出畸变波前φ(x,y)；

7）由重构的降质波前去估计大气瞬时点扩散函数h(x,y)。

2.根据权利要求1所述的自适应光学成像中的点扩散函数估计方法，其特征在于：所述步骤1）中，根据微透镜阵列排布方式确定三角剖分的形式，确定三角剖分的形式时采用基于几何布局的网格单元直接剖分方法。

3.根据权利要求1所述的自适应光学成像中的点扩散函数估计方法，其特征在于：所述步骤2）中，斜率测量模型的偏微分方程的建立方法具体为：

首先将波前表示为：

$\mathrm{\&phi;}\left(x\right)\&ap;s_d^r\left(b\left(x\right)\right)=\underset{\left|\mathrm{\&lambda;}\right|=d}{\mathrm{\&Sigma;}}{c}_{\mathrm{\&lambda;}}{B}_{\mathrm{\&lambda;}}^d\left(b\left(x\right)\right)=B^d\left(b\left(x\right)\right)\&CenterDot;c$

其m阶方向导数可表示为 $D_u^m{s}_d^r\left(b\left(x\right)\right)=\frac{d!}{(d-m)!}{B}^{d-m}{P}_u^{d,d-m}\left(a\right)\&CenterDot;c,$ 其中，

为二元单纯形样条，d为样条函数次数，r为连续性阶数，b(x)为一个三角剖分上笛卡尔坐标到面积坐标的变换，

为d次Bernstein多项式基函数，λ=(λ₀,λ₁,λ₂)为多重符号且其长度为|λ|=λ₀+λ₁+λ₂=d，c_λ为多项式系数，B^d(b(x))和c分别为

和c_λ的矩阵形式，u为笛卡尔坐标下方向导数的方向向量，a为与u相对应的面积坐标下的方向坐标，

为分块对角矩阵，每对角块为对应于a和u的de Casteljau矩阵；其次根据相邻三角剖分间的r阶光滑拼接条件，每条边都有

个约束条件，从而构成约束方程Hc=0，H为由E条边构成的E*Z行约束矩阵，且H的每行都包含一个单独的连续性约束，这样，波前传感器的斜率测量偏微分方程即可表示为：

${\mathrm{\&sigma;}}_u=d{B}^{d-1}{P}_u^{d,d-1}c+n$

Hc=0

其中，σ_u表示波前传感器斜率测量矩阵，n表示测量噪声及模型误差。

4.根据权利要求3所述的自适应光学成像中的点扩散函数估计方法，其特征在于：相邻三角剖分间的r阶光滑拼接条件描述如下：

设

分别是定义在相邻三角形t=<ν₁,ν₂,ν₃>和上的d次多项式，ν_i(i=1,2,3)为三角形的顶点，c_λ和

分别是S(b(x))和关于t和

的贝塞尔坐标，则S(b(x))和

之间的r阶光滑拼接条件为：

${\stackrel{\&OverBar;}{c}}_{{\mathrm{\&lambda;}}^s}=c_{{\mathrm{\&lambda;}}^0}^{\left(s\right)}\left(\stackrel{\&OverBar;}{\mathrm{\&tau;}}\right),s=\mathrm{0,1},\&CenterDot;\&CenterDot;\&CenterDot;,r$

其中，

是关于t的面积坐标，λ^s=(s,λ₁,λ₂)，λ⁰=(0,λ₁,λ₂)，λ₁+λ₂=d-s。

5.根据权利要求3所述的自适应光学成像中的点扩散函数估计方法，其特征在于：所述步骤3）中，确定变换矩阵D的方法具体为：

步骤2）建立了一个约束最小二乘问题σ_u=Dc+n,Hc=0；与Zernike模式法相对应，通过对比此约束最小二乘问题与斜率测量偏微分方程，即可将变换矩阵D表示为：

$D=d{B}^{d-1}{P}_u^{d,d-1}.$

6.根据权利要求1所述的自适应光学成像中的点扩散函数估计方法，其特征在于：所述步骤4）中，确定重构矩阵R的方法具体为：

将由斜率测量数据去重构样条函数系数向量的矩阵定义为重构矩阵，结合步骤1）和步骤2），通过求解约束最小二乘问题，多项式系数向量c的显式解表达式为：

$\hat{c}=[{\left(D^{T}D\right)}^{-1}{D}^T-{\left(D^{T}D\right)}^{-1}{H}^T{\left[H{\left(D^{T}D\right)}^{-1}{H}^T\right]}^{-1}{\left(D^{T}D\right)}^{-1}{D}^T]\mathrm{\&sigma;}$

重构矩阵R，则表示为：

R=(D^ΤD)^-1D^Τ-(D^ΤD)^-1H^T[H(D^ΤD)^-1H^T]^-1(D^ΤD)^-1D^Τ

其中，Τ表示矩阵转置，σ为波前斜率测量向量，此时由重构矩阵所确定的样条函数系数是在约束Hc=0下的最优线性无偏估计。

7.根据权利要求1或6所述的自适应光学成像中的点扩散函数估计方法，其特征在于：所述步骤5）中，根据步骤4），重构矩阵是由斜率测量数据去重构样条函数系数向量的矩阵，二元样条函数的系数向量

直接由重构矩阵与波前斜率测量向量相乘得到，则其计算公式为：

$\hat{c}=\mathrm{R\&sigma;}.$

8.根据权利要求7所述的自适应光学成像中的点扩散函数估计方法，其特征在于：所述步骤6）中，畸变波前φ(x,y)的重构方法为：

根据步骤2），已将畸变波前表示为单纯形剖分上的二元样条，步骤5）又得到了样条函数的系数向量

而d次Bernstein多项式基函数B^d(b(x))是确定的，将二者相乘即可得到畸变波前的估计量

$\widehat{\mathrm{\&phi;}}\left(x\right)=B^d\left(b\left(x\right)\right)\hat{c}.$

9.根据权利要求8所述的自适应光学成像中的点扩散函数估计方法，其特征在于：所述步骤7）中，估计大气瞬时点扩散函数h(x,y)的方法具体为：

已知透镜光瞳函数为A，根据步骤6），由傅里叶光学理论，结合透镜光瞳函数A，得到其波前对应的点扩散函数为：

$h(x,y)={\left|{\mathrm{FT}}^{-1}\right[\mathrm{Aexp}\left(j\widehat{\mathrm{\&phi;}}\left(x\right)\right)\left]\right|}^2$

其中，FT^-1表示傅里叶逆变换，

Images