基于B 樣條函數(shù)的快速波前復原

陳 浩，魏 凌，李恩德，何 益，楊金生，李喜琪，樊新龍，楊澤平，張雨東

1 中國科學院光電技術研究所，自適應光學重點實驗室，四川 成都 610209；

2 中國科學院大學材料科學與光電子技術學院，北京 100049；

3 江蘇材料光學重點實驗室，江蘇 蘇州 215163；

4 中國科學院蘇州生物醫(yī)學工程技術研究所，江蘇 蘇州 215163

1 引 言

夏克-哈特曼波前傳感器(Shack-Hartmann wavefront sensor，簡稱哈特曼)是一種常用的斜率型波前探測器，被廣泛應用于自適應光學的相關領域[1-3]。它并不能直接探測波前的相位信息，只能通過子孔徑內的斜率重構波前相位。由于傳感器的噪聲、采樣誤差、子孔徑內存在高階像差等原因，斜率的計算易受噪聲影響。盡管可以提高斜率計算的精度和抗噪能力[4-5]，為復原計算提供更精確的輸入，然而，斜率探測誤差是不能完全消除的，這需要波前復原算法具有抗噪能力。因此，波前復原算法需要綜合考慮抗噪性、準確度和計算復雜性。常見的復原算法可以分為兩類：區(qū)域法和模式法。

區(qū)域法在待測波前點與斜率測量點之間建立差分模型，然后通過求解線性的差分方程組復原波前信息，傳統(tǒng)的區(qū)域法模型有Fried[6]和Southwell[7]等。這種方法具有良好的局部性，可以較精確地復原局部區(qū)域的像差，并且可以用于任意形狀的光瞳。然而由于差分模型的精度及差分過程導致的誤差，算法抗噪能力較弱。另外，由于差分模型的限制，復原的波前點數(shù)量與子孔徑數(shù)量相當，需要通過插值來獲得高空間分辨率的波前信息，而這會帶來額外的誤差。

模式法一般采用空間正交多項式復原波前。最常用的是Zernike 多項式，因為它與Seidel 像差有非常直接的對應關系[8]。模式法復原的波前具有解析表達式，可以得到任意空間分辨率的波前信息。由于Zernike多項式只在圓域內正交，因此只能復原圓形光瞳的波前，不少學者采用變換等方式使Zernike 可以適用于方形或其它形狀的光瞳[9]，但仍存在變換復雜，不能完全適用于所有形狀光瞳的缺點。另外，盡管Zernike多項式在圓域內是正交的，但其導數(shù)并不正交，即各個模式的斜率不正交，因此會發(fā)生模式泄露，并且在不同階次的Zernike 空間中，展開系數(shù)是不一樣的。Nam 等[10]提出正交化斜率矩陣的方法來解決這個問題，但此時的解空間已不再是原本的Zernike 空間了。另外，由于正交多項式的非局部性，導致任一子孔徑的斜率可能會對所有的模式產生影響，對邊界上的波前復原精度較差，因此模式法不能很好地復原局部區(qū)域的像差。由于多項式的正交性，模式法一般適用于圓形光瞳或方形光瞳，當光瞳形狀與預設不一致時，復原精度會受到比較嚴重的影響。

在高能激光[3]、光學加工檢測[11]以及基于像差的PSF 計算[12]等應用中，對局部像差的復原能力要求相對較高，而基于Zernike 多項式的模式法由于局部性較差，并且復原的像差容易在邊緣出現(xiàn)較大的值，在這些場合應用有限。

為了保證復原的局部性同時提高抗噪能力，樣條函數(shù)的方法得到了廣泛的應用。Seifert 等提出了一種采用三次B 樣條進行波前復原的方法[13]，該方法用二維三次樣條曲面對波前進行最小二乘擬合。由于B 樣條曲面的平滑性和緊支撐特性，可以在保留局部性特征的同時，具有比傳統(tǒng)區(qū)域法更強的抗噪聲能力。文章仿真結果表明，隨著樣條數(shù)量的增多，局部復原精度越高，但當樣條數(shù)量接近斜率采樣點數(shù)量時，復原矩陣的條件數(shù)會增大，復原精度會有所下降。然而文章并未明確B 樣條數(shù)量的選擇方法，只是建議針對應用場合的像差仿真來確定最優(yōu)樣條數(shù)量。Ares 定量比較了三次B 樣條波前復原算法與Zernike 模式法的性能區(qū)別[14]，證實了對于復雜的波前，三次B 樣條復原方法更具有優(yōu)勢。隨后，Cornelis 等提出了基于多變量的非線性樣條波前重構方法(spline based aberration reconstruction,SABRE)[15]，該方法本質上屬于區(qū)域法，通過特定的三角形分割方式，采用樣條函數(shù)空間模型代替常規(guī)的差分模型，該方法可以適用于任意形狀的光瞳和任意的子孔徑布局。SABRE 方法在模型劃分時，有多種分割模型，不同分割模型的復原精度也不盡相同，而文章中并未提出分割模型的選擇方法。Huang 等提出了一種使用局部化的三次樣條函數(shù)復原波前的方法[16]，盡管能提高局部復原精度，但也更容易受到噪聲的影響。Pant 使用加權的三次樣條函數(shù)復原波前[17]，在提高算法抗噪能力的同時也提高了算法計算復雜度。

從以上文章可以看出，相比普通樣條函數(shù)，B 樣條函數(shù)在控制點處均具有相同的樣條基函數(shù)，在抗噪能力和局部復原能力上均衡，計算復雜度也較小。本文以2 階B 樣條為例，提出了一種基于B 樣條函數(shù)的快速波前復原算法，并將原本的斜率最小二乘化問題轉化為泊松方程，然后使用超松馳迭代法進行求解，大大減少算法的存儲需求和計算量，若使用分塊加速超松馳迭代法，還可以提高算法的并行程度，方便應用到多核PC 系統(tǒng)、FPGA 或DSP 等硬件平臺上。

本文結構如下：第二節(jié)介紹基于B 樣條的波前重構算法；第三節(jié)介紹基于泊松方程的2 階B 樣條快速復原算法；第四節(jié)通過對驅動器響應函數(shù)測試的實驗，將本算法與傳統(tǒng)波前重建算法進行比較；第五節(jié)對全文進行總結和展望。

2 基于B樣條的哈特曼波前快速復原算法

2.1 哈特曼基本原理

哈特曼波前傳感器通過微透鏡陣列、小孔或光柵等光學元件，將待測波前分割為N個小區(qū)域后，每個區(qū)域會在CCD 上成像。當待測量波前與參考波前存在偏差時，這個偏差會轉化為子孔徑內光斑的偏移(Δx,Δy)i,i=1,2,...,N，令波前為φ，微透鏡焦距為f，則有：

波前重構算法即是通過采樣到的{?φi}復原出波前φ。

2.2 B 樣條曲線

B 樣條(Basis spline)最初由Isaac Jacob Schoenberg提出，其具有幾何不變性、凸包性、保凸性、局部支撐性等優(yōu)良特性，是應用最廣泛的樣條函數(shù)之一[18]。

k階B 樣條曲面方程可表示為

其中：i和j為節(jié)點序號，為控制頂點，Ni,k及Nj,k為k次規(guī)范B 樣條基函數(shù)，基函數(shù)為次分段多項式，可由遞推公式得到：

定義u i=i，即均勻采樣的情況，則可以得到1 階和2 階B 樣條基函數(shù)分別為

圖1 繪制了控制節(jié)點間隔為1 的1 階和2 階B 樣條基曲面。

2.3 基于2 階B 樣條基曲面的波前復原

為簡化表述，令：

圖1 B 樣條基曲面。(a) 1 階B 樣條基曲面；(b) 2 階B 樣條基曲面Fig.1 Surfaces of B-spline basis.(a) First-order B-spline surface;(b) Second-order B-spline surface

易知B(x,y) 為關于原點對稱的函數(shù)。由2.2 節(jié)知，B(x,y)具有緊支撐特性，通過將像差波面分解為一系列B(x,y) 波面的疊加，可以在有效抑制噪聲的基礎上更好地復原局部波前。

為使B 樣條復原的空間頻率最高，需要B 樣條基的數(shù)量不少于波前采樣點的數(shù)量，對于方形布局的微透鏡陣列，在傳統(tǒng)的區(qū)域法中，有兩種常用的模型可以借鑒——Fried 模型和Southwell 模型。

圖2 描述了這兩種模型的B 樣條基位置與子孔徑之間的關系。其中灰色圓圈表示樣條基的位置，正方形為子孔徑位置，正方形內的正交箭頭為子孔徑內斜率。

假設子孔徑數(shù)量為(M×N)，為達到最優(yōu)的空間分辨率，在 Fried 模型中，B 樣條基數(shù)量為(M+1）× (N+1)，Southwell 模型中，由于B(0,0)=0，需要在上下左右至少分別補充1 排B 樣條基，即B 樣條基的數(shù)量最少為(M+2）× (N+2)。

以子孔徑邊長為單位，令左上角坐標為(0,0)，則Fried 模型和Southwell 模型排布的B 樣條基分別如式(7)和式(8)所示：

代入式(1)，則對于每個子孔徑有：

圖2 方形排布的B 樣條基位置與子孔徑之間的關系。(a) Fried；(b) SouthwellFig.2 Positional relation between B-spline basis and subaperture with a square layout.(a) Fried model;(b) Southwell model

由于B(x,y) 的緊支撐特性，R為一稀疏矩陣，因此式(12)為稀疏方程。

易知，若B(x,y) 使用1 次B 樣條函數(shù)，式(12)與傳統(tǒng)的Fried 區(qū)域法一致，而Southwell 布局也有類似的關系，因此傳統(tǒng)的區(qū)域法可以視為B 樣條函數(shù)復原算法在1 階條件下的特例。為了描述簡潔，本文后面的討論和仿真均采用Fried 模型，而Southwell 模型推導方式與之類似。

3 基于泊松方程的快速2 階B 樣條函數(shù)波前復原算法

式(12)的直接解法為求解R的廣義逆矩陣，即

盡管R為一個稀疏矩陣，但R的廣義逆矩陣R+并不是稀疏矩陣，在復原的時候即使提前計算好R+，也需要(M+1）2× (N+1)2個數(shù)的存儲空間，并且一次重構需要(M+1）2× (N+1)2次乘法。當子孔徑數(shù)目比較多的時候，使用FPGA 或DSP 等計算時，會帶來巨大的空間開銷，可以預估，當子孔徑數(shù)提升1 倍時，存儲空間和計算量會提升16 倍，這種情況下，這種計算方法將制約B 樣條波前復原算法的應用。

由式(10)知，波前重構是求解如下最優(yōu)化問題：

代入式(10)，可得：

圖3 B 樣條散度積分示意圖Fig.3 B-spline divergence integral diagram

式(19)左邊等價于計算圖3 中編號為5 的B 樣條基處的散度(其中圓圈代表B 樣條基位置，實線方形表示子孔徑區(qū)域，虛線方形表示B 樣條基的積分區(qū)域)，支撐域包含編號5 區(qū)域的B 樣條基為1～9 號，對ΔBi(x,y)(i=1,2,…,9)的解析表達式在編號5 區(qū)域中積分，式(19)可改寫為

改變成式(20)后，即可使用逐次超松馳迭代法(successive over relaxation method,SOR)進行求解[19]：

接下來計算??g。為方便描述，以圖4 所示為1,2,3,4 表示4 個子孔徑，令為子孔徑x和y方向的斜率，為計算黑色圓點處的散度??g，在黑色圓點處對斜率進行1 階Taylor 展開，可得：

進一步可得出：

圖4 散度近似計算示意圖Fig.4 Diagram of divergence approximation

這樣就可以通過子孔徑的斜率估計出B 樣條基處(黑色圓點)的散度。至此，式(19)中的必要元素已具備。通過求解式(19)，則可得到所有B 樣條基的系數(shù)，最后通過式(9)則可得到復原波前的解析式。

對于邊界位置B 樣條基，由于積分區(qū)域改變，式(20)會有所改變，以圖3 中的1 號B 樣條基為例：

而邊界位置處的散度估計，可以通過對斜率的外拓來簡化計算，可采用重復延拓或補0 延拓。

可見，基于泊松方程的2 階B 樣條快速重建算法，存儲空間為(M+1）× (N+1)，乘法需求量為(M+1）× (N+1)×L，其中L為迭代步數(shù)。一般而言，SOC 算法在迭代幾十次左右后將收斂。因此相比廣義逆求解，存儲空間和計算量的需求減少了很多，特別當子孔徑數(shù)目增大時，優(yōu)勢更明顯。

4 實 驗

本文利用自主研制的121 單元六邊形排布的壓電陶瓷變形鏡，測量其中心位置處7 個驅動器(命名為1～7 號，正中心位置為4 號)的響應函數(shù)，分別通過干涉儀和哈特曼波前傳感器采集數(shù)據(jù)，再通過不同的復原算法從哈特曼的數(shù)據(jù)中復原波前，并與干涉儀的數(shù)據(jù)進行對比。

為減少變形鏡初始面型以及光學系統(tǒng)對測試結果的影響，每個驅動器分別加±80 V 后，對測量得到的波前進行對減。基準數(shù)據(jù)采用干涉儀(Zygo GPI-600)測量的結果，4 號驅動器的響應函數(shù)如圖5 所示。

使用哈特曼測量時，將變形鏡的光瞳口徑縮束后與哈特曼波前傳感器成3：4 匹配(即哈特曼波前傳感器的有效探測區(qū)域為變形鏡的75%)，哈特曼波前傳感器的微透鏡陣列為方形排布，陣列數(shù)為22×22，子孔徑大小為0.25 μm，焦距為9.52 mm(使用Yang 等提出方法標定[20])，測量波長為808 nm。

通過區(qū)域法、模式法和本文提出的復原算法(迭代次數(shù)取20 次)復原的波前和殘差波前如圖6 所示，為實現(xiàn)相同數(shù)字采樣率，區(qū)域法的結果采用雙線性插值方法進行處理。為實現(xiàn)對減，首先通過相關找到原始兩個波面的位置偏移，區(qū)域法通過雙線性插值重新生成與干涉儀位置一樣的波前；模式法和本文提出的算法通過偏移后的坐標計算波前，考慮到模式法的邊緣效應，僅計算單位圓域內的波前數(shù)據(jù)，而本文提出的算法僅計算樣條內部的區(qū)域。4 號驅動器的復原波前和殘差如圖6 所示，相關系數(shù)、殘差如表1 所示，所有驅動器的波前殘差PV 值與RMS 值如圖7 所示，對應的數(shù)值結果如表2 所示。

從與干涉儀對比的結果可知，由于Zernike 模式法的局限性，其對局部像差的復原能力較差，因此對這種驅動器產生的響應函數(shù)無法精確復原，復原殘差很大。區(qū)域法和本文提出的方法均可以較好地復原這種局部像差。從復原的數(shù)值上看，本文的方法優(yōu)于區(qū)域法，其原因在于：區(qū)域法采用雙線性插值時會帶來插值誤差；區(qū)域法的抗噪能力有限，因此殘差波前的PV 值顯著大于本文提出的方法。另外，觀察圖6 中驅動器頂端復原的波前細節(jié)，可以看出，區(qū)域法在這個位置由于受采樣誤差和噪聲的影響，未能復原出平滑的“峰”。

圖5 4 號驅動器的干涉儀測量數(shù)據(jù)Fig.5 Measurement data of No.4 actuator obtained with ZYGO interferometer

圖6 不同復原方法復原結果。(a) 本文算法重建波前；(b) 本文算法的殘余波前；(c) 基于Zernike 多項式的模式法重建波前；(d) 模式法的殘余波前；(e) Fried 區(qū)域法重建波前；(f) 區(qū)域法的殘余波前Fig.6 Wavefronts restored by different methods.(a) Wavefront restored by our method;(b) Residual wavefront error of (a);(c) Wavefront restored by the modal method;(d) Residual wavefront error of (c);(e) Wavefront restored by the zonal method;(f) Residual wavefront error of (e)

表1 不同復原方法結果比較(4#驅動器)Table 1 Comparison results of different wavefront reconstruction methods

圖7 不同驅動器復原殘差Fig.7 PV and RMS results of residual wavefront reconstructed by different methods

表2 不同驅動器復原殘差及相關值Table 2 Residual reconstruction error and correlation values of different actuators

5 結 論

本文提出了一種基于2 階B 樣條快速波前復原算法，采用泊松方程改寫斜率最小二乘方程，使其可以通過超松馳迭代法進行快速求解，大大降低了復原算法的存儲空間需求和計算復雜度，容易在DSP 或FPGA 等硬件中實現(xiàn)，在子孔徑數(shù)目增多時其優(yōu)勢更加明顯，并且可以采用分塊的超松馳迭代法等方法[21]提高算法的并行計算能力。該方法可以推廣到其它斜率型波前傳感器，如金字塔波前傳感器[22]，四波剪切干涉儀[23]等。

通過對變形鏡驅動器響應波前的測量實驗可以看出，模式法存在明顯的局限性，無法精確復原局部像差，而本文提出的算法具有與區(qū)域法相當?shù)木植肯癫顝驮芰Γ瑫r有優(yōu)于后者的抗噪能力。由于復原的波前具有解析解，無需插值即可得到非常高的數(shù)值采樣精度的波前信息，這使得復原波前更加平滑。

通過泊松方程改寫的式(19)，左側為對樣條函數(shù)的散度積分，右側為利用斜率對散度的估計。因此，左側可以擴展到不同階次和不同擴展度的B 樣條基，實現(xiàn)數(shù)值采樣精度的靈活選取。而右側根據(jù)式(23)中的處理方法，可以靈活選取計算散度的區(qū)域。一般來說，計算散度利用的斜率點越多，抗噪能力越強，對局部像差的復原能力越弱。因此，通過靈活選取計算散度的區(qū)域，可以有效平衡對局部區(qū)域像差的復原能力和抗噪能力需求。