周睿 魏凌李新陽王彩霞李梅沈鋒
1)(中國科學院自適應光學重點實驗室,成都 610209)2)(中國科學院光電技術研究所,成都 610209)3)(中國科學院大學,北京 100049)
點光源哈特曼最優閾值估計方法研究
周睿1)2)3)?魏凌1)2)李新陽1)2)王彩霞1)2)李梅1)2)沈鋒1)2)
1)(中國科學院自適應光學重點實驗室,成都 610209)2)(中國科學院光電技術研究所,成都 610209)3)(中國科學院大學,北京 100049)
(2016年11月8日收到;2017年2月4日收到修改稿)
針對夏克-哈特曼波前傳感器探測系統中噪聲隨時間及空間變化頻率較快的特點,為了準確估計系統的最優閾值,根據高斯光斑與噪聲的分布特性,提出一種以滑動窗口內像素均值及圖像信號的局部梯度作為參數,構造關于噪聲權重函數的方法,由此獲得子孔徑閾值的最優估計值,并詳細分析了算法的基本原理和實現過程.以典型處理方法獲取的閾值與理論最優閾值的誤差作為評價標準,仿真和實驗結果表明本文提出的閾值估計方法在不同信噪比、不同光斑大小的條件下,均能取得優于典型閾值處理方法獲得的結果,且與理論最優閾值的誤差小于10%.
夏克-哈特曼波前傳感器,高斯光斑,最優閾值,權重函數
夏克-哈特曼波前傳感器(SHWFS)是一種基于波前斜率測量的光學檢測裝置,因其具有結構簡單、環境適應能力強、實時性強等優點,廣泛應用于自適應光學、激光波前檢驗、生物光學、光學檢測和裝調等領域[1].SHWFS主要由微透鏡陣列和電荷耦合器件(CCD)探測器組成,利用微透鏡陣列對輸入光束進行分割采樣,每個透鏡作為一個子孔徑,將光束聚焦成一個光斑陣列,子孔徑范圍內的波前畸變將造成光斑的位置偏移,SHWFS對輸入波前的測量精度主要取決于各個子孔徑光斑中心位置的提取精度[2?4].光斑中心位置提取的算法有很多種,包括一階矩法、高斯擬合法、權重法、相關法等,其中一階矩法由于算法簡單,魯棒性強,精度較高而被廣泛使用.一階矩法以光斑強度為權重,計算區域內所有像素對光斑質心的貢獻,直接得到強度分布與空間平均相位的關系[5].
然而一階矩法對系統噪聲比較敏感,為了提高一階矩法質心提取精度,通常需要對CCD的讀出信號進行預處理[6,7].在預處理的方法中,閾值處理的方法由于能夠快速去除噪聲的影響而受到廣泛的關注[8].文獻[9]指出,根據SHWFS的特點和噪聲分布特性,存在一個最優閾值,使得質心的探測誤差最小;同時推導了最優閾值等于噪聲的均值與噪聲的3倍標準差之和.然而在實際工作過程中,由于SHWFS系統中有效信號和噪聲很難被區分,因此最優閾值難以實時地準確估計.
目前,獲取SHWFS閾值的方法通常有基于灰度直方圖統計的經驗值法、迭代法、最大類間方差(Ostu)法、基于圖像局部特征的閾值方法等.文獻[10]提出的基于灰度直方圖統計的經驗值法,僅適用于噪聲隨時間、空間頻率變化緩慢的情況,對于背景隨時間、空間頻率變化比較劇烈且需要對SHWFS的采樣數據實時閉環的場景不適合.文獻[11]指出,Ostu法在信號相對面積占圖像總面積達到30%以上時其分割性能接近最優值.然而,通常情況下SHWFS波前探測系統有效信號的面積小于圖像總面積的20%,采用這種方法得到的閾值偏大,從而影響質心計算的精度.文獻[12,13]提出使用最大值的百分比作為閾值,但是由于百分比的確定與探測系統的信噪比(SNR)有關,在實際系統中,難以實現對百分比的優化.文獻[14]提出,根據SHWFS的特點,利用探測靶面的四個邊角區域估計噪聲的均值和方差,這種方法只適用于背景噪聲空間分布比較均勻的情況,對于背景空間分布不均勻的場合,這種為所有子孔徑設置同一閾值的方法不適用.
本文針對SHWFS系統中光斑及噪聲分布的特點,通過求取子孔徑滑動窗口內像素均值及信號的局部梯度方向,構造噪聲權重函數的方法對最優閾值進行估計,對該算法進行了仿真和實驗,并將該算法得到的閾值與典型處理方法得到的閾值進行了對比和分析.結果表明,本文提出的方法與典型方法相比,能夠得到更接近理論最優閾值的結果,從而有效地提高了SHWFS質心探測精度.
在SHWFS探測波前相位畸變時,采用孔徑分割的方法,將相位面分割為多個子孔徑,探測子孔徑的平均斜率來重構相位.通常子孔徑內的光斑強度可以用高斯函數來擬合[15],其表達式為

其中A為高斯光斑的幅值;xc,yc,σx,σy分別為光斑x,y方向的中心位置及x,y方向的光斑寬度.由于SHWFS光斑的強度呈高斯分布且面積較小,其像素灰度值比較集中、分散性小,而噪聲一般可認為是高斯白噪聲,典型SHWFS光斑沿x軸的剖面如圖1所示.根據噪聲和光斑的分布特性,由噪聲區域對應的窗口(i)和有效光斑區域對應的窗口(ii)可以看出,隨著滑動窗口逐漸靠近有效光斑的中心位置,對應窗口內像素擬合曲線與x軸的夾角、窗口內像素均值均呈現單調增長的趨勢,即像素屬于噪聲的概率與滑動窗口內像素擬合曲線與x軸夾角及窗口內像素均值成反比.
同樣,利用上述分析方法,可以對二維SHWFS光斑進行二維滑動窗口的處理,對二維滑動窗口內像素求取局部梯度方向,具體做法是對窗口內的像素信號進行平面擬合,從而得到其擬合平面法線方程與z軸的夾角,其數學過程可以表達為

其中X,Y為光斑在滑動窗口內的坐標位置,Z為光斑強度信息,(ai,j,bi,j,di,j)T為需要求取的擬合平面法向量方程的參數,其矩陣乘法表達式為

以滑動窗口為3×3為例,可以得到

其中M,v,Z分別為滑動窗口內坐標位置組成的矩陣、待擬合平面的法向量及窗口內像素的強度向量;M+為M的偽逆.對于3×3的窗口,可以得到

因此,滑動窗口內各像素擬合平面的法向量可以表達為如下的卷積形式,其中?表示卷積:

圖1 典型高斯光斑x軸剖面圖Fig.1.The x axis pro fi le of tradition SHWFS spot.

由此,可以得到子孔徑內各像素擬合平面法向量與z軸的夾角,即局部梯度方向可以表示為

同時,由前述分析可知,像素屬于有效光斑的概率不僅與滑動窗口內像素擬合曲線與x軸的夾角有關,也與窗口內像素均值有關,滑動窗口內像素均值可以表示為(11)式,其歸一化的均值可以表示為(12)式:

由此,像素屬于噪聲的概率Pi,j可以表示為θi,j,μ′i,j的權重函數,即

在上述分析的基礎上,根據(13)式可得子孔徑內所有像素的概率,從而得到基于滑動窗口統計的噪聲均值如(15)式:

同時,為了有效估計噪聲的標準差,需要將滑動窗口投影到x,y平面以便實現對有效高斯光斑區域的噪聲標準差的估計,如(16)式所示,其中Z′表示去除傾斜后的光強:

由此,可以求得去除傾斜后滑動窗口中心像素的標準差估計如(17)式所示,其中分別為去除傾斜后滑動窗口各個像素的灰度值和均值,

因此,可以得到子孔徑內關于噪聲標準差的估計如下:

根據上述分析,由(15)和(18)式得到了子孔徑內關于噪聲的均值與標準差的估計值,從而可以得到子孔徑內最優閾值的估計為

為了驗證基于滑動窗口投影最優閾值估計方法的有效性,本節針對該方法與迭代閾值法、最大值百分比閾值法及Ostu法求得的閾值與理論最優閾值進行了對比分析.其中迭代閾值法設定圖像最大、最小灰度值的平均數為初始閾值T0,閾值T(i)將圖像分為前景和后景,u,v分別為前景和后景的平均灰度等級,取T(i+1)=(u+v)/2,當迭代到T(i+1)=T(i)時停止,取該值為子孔徑內圖像的閾值;最大值百分比閾值法設定子孔徑內圖像最大值的百分比作為閾值,根據文獻[12,13]通常選擇最大值的30%;Ostu法將圖像分為背景和目標兩類,使得兩類之間的方差最大,從而求得閾值.
為了便于分析和比較在不同光斑大小、不同高斯光斑中心強度及不同噪聲水平情況下幾種閾值處理方法結果與理論最優閾值的接近程度,可以將SNR定義為高斯光斑的能量與噪聲標準差的比值,其表達式如下:

同時為了評價不同方法得到的閾值與理論最優閾值的接近程度,定義如下式所示的性能指標來衡量各種閾值處理方法的效果,

根據(1)式定義的高斯光斑,針對子孔徑尺寸為24 pixels×24 pixels,高斯光斑的半徑從0.5 pixels到3 pixels變化,分別利用高斯光斑能量保持不變、改變噪聲水平的方法和噪聲水平不變、改變高斯光斑能量的方法,在不同SNR條件下,將典型閾值處理方法與本文提出的方法進行了比較和分析.如(20)式所示,隨著SNR的增加,被測光斑的質量逐漸改善.仿真過程中的典型光斑如圖2所示,其中圖2(a)表示在SNR=16.73 dB條件下光斑和噪聲強度的分布情況,圖2(b)表示受噪聲影響光斑強度的三維分布情況.
當滑動窗口大小為3 pixels×3 pixels時,保持噪聲水平不變,通過調整高斯光斑能量的方法改變SNR,將本文提出的最優閾值估計法、迭代閾值法、最大值百分比閾值法及Ostu法的結果與理論最優閾值進行比較,結果如圖3所示.其中圖3(a)—(d)表示光斑的半徑σA=0.5,1.0,2.0,3.0 pixels條件下的仿真結果.從圖3可以看出:在4種不同高斯光斑大小條件下,當SNR較高時,典型閾值處理方法得到的閾值約為最優閾值的2—10倍,遠遠高于理論的最優閾值;并且光斑半徑越小,使用典型閾值處理方法得到的閾值偏離理論最優閾值的幅度越大;隨著SNR的降低,典型閾值處理方法得到的閾值偏離理論最優閾值的程度有所降低,但與理論的最優閾值相比,仍有較大誤差;在SNR較低的情況下,典型閾值處理方法得到的閾值會小于理論的最優閾值,并且光斑半徑越大,使用典型閾值處理方法得到的閾值偏離理論最優閾值的幅度越大,最大偏差約是理論最優閾值的1/2,因此不能有效抑制噪聲信號;然而,在不同的光斑大小及SNR條件下,本文提出的閾值估計方法的結果與理論最優閾值的誤差均能保持在10%以內,且在各種光斑半徑下一致性較好,明顯優于典型閾值處理方法.
圖4顯示了在滑動窗口大小為3 pixels×3 pixels,保持高斯光斑能量不變,通過調整噪聲水平的方法改變SNR時,本文提出的最優閾值估計法與典型閾值處理方法進行比較的結果.從圖4結果可以得到與圖3相同的結論:在SNR較好的情況下,典型閾值處理方法得到的閾值遠遠高于理論的最優閾值,在SNR較低的情況下,典型閾值處理方法得到的閾值低于理論的最優閾值;而本文提出的方法能夠在不同SNR條件下,對4種高斯半徑光斑得到的最優閾值與理論最優閾值的誤差范圍保持在10%以內.

圖2 (a)SNR=16.73 dB時包含背景噪聲的光斑圖像;(b)SNR=16.73 dB時包含背景噪聲的光斑圖像的三維顯示Fig.2.(a)Spot image with background noise(SNR=16.73 dB);(b)the three-dimensional show of the spot image with background noise(SNR=16.73 dB).

圖3 (網刊彩色)噪聲水平不變,不同高斯光斑強度條件下各種閾值估計方法得到的閾值與理論最優閾值的誤差(a)σA=0.5;(b)σA=1.0;(c)σA=2.0;(d)σA=3.0Fig.3.(color online)The error between the theoretic threshold and the threshold obtained by di ff erent methods with the invariable noise level and the di ff erent spot energies:(a) σA=0.5;(b) σA=1.0;(c)σA=2.0;(d)σA=3.0.

圖4 (網刊彩色)高斯光斑能量不變,不同噪聲水平條件下各種閾值估計方法得到的閾值與理論最優閾值的誤差(a)σA=0.5;(b)σA=1.0;(c)σA=2.0;(d)σA=3.0Fig.4.(color online)The error between the theoretic threshold and the threshold obtained by di ff erent methods with the invariable spot energy and the di ff erent noise levels:(a)σA=0.5;(b)σA=1.0;(c)σA=2.0;(d)σA=3.0.
從上述仿真和分析可以看出,在高斯光斑強度或噪聲水平發生變化引起SNR發生改變時,本文提出的最優閾值估計方法與典型閾值處理方法相比,均能夠較為準確地得到噪聲的理論最優閾值.同時,本文也對不同滑動窗口大小情況下,噪聲權重函數中μ′和cosθ取不同的冪指數進行了仿真和分析,均能獲得與圖3和圖4相類似的結論.
為了驗證本文提出的最優閾值估計算法的效果,建立了如圖5所示的實驗系統.該實驗系統主要包括激光器、由微透鏡陣列和CCD探測器組成的SHWFS及計算處理模塊3個部分.其中激光器采用波長λ=1.064μm的電流可調節的激光器,SHWFS采用方形排布8×8的微透鏡陣列,單個微透鏡的口徑為2.667 mm,焦距為77 mm,SHWFS使用的CCD探測器靶面大小為240 pixels×240 pixels,每個子孔徑大小為24 pixels×24 pixels,像元尺寸為24μm,探測器位數為14位.SHWFS獲取的典型的8×8(其中有效子孔徑數為48)的子光斑圖像陣列如圖6(a)所示,利用本文提出的最優閾值估計方法處理后的子光斑圖像陣列如圖6(b),各個子孔徑的SNR如圖6(c)所示.由圖6(b)和圖6(c)可以看出,在各子孔徑內SNR從10—25 dB大范圍波動的情況下,利用本文提出的方法依然能夠有效地去除各個子孔徑內的噪聲.

圖5 實驗裝置示意圖Fig.5.The schematic diagram of experiment Setup.
針對如圖5所示的SHWFS系統,根據文獻[8,16]在光子噪聲可以被忽略的情況下,將背景噪聲、讀出噪聲、雜散光噪聲等信號統一地看作是系統的探測噪聲,并且該噪聲服從高斯分布.首先在關閉激光器的狀態下,獲取SHWFS系統的噪聲圖像信息,通過多幀平均的方法獲得相對準確的噪聲均值和標準差;在實驗過程中探測器的增益、曝光時間以及實驗環境等條件均保持一致.因此,可以認為噪聲的均值和標準差基本保持不變,即理論最優閾值保持不變,根據(19)式可以得到該實驗條件下的理論最優閾值為410 ADU;然后通過調節激光器電流的方法改變有效信號的能量,從而得到受噪聲影響的圖像信號.對SHWFS系統采集的圖像信號,分別利用本文提出的方法、迭代閾值法、最大值百分比閾值法及Ostu法求得閾值,并與理論最優閾值進行了對比分析,4種方法結果與理論最優閾值的誤差曲線如圖7所示;表1顯示了4種方法在不同SNR條件下得到的計算結果.由圖7和表1可以看出,與典型閾值處理方法相比,利用本文提出的最優閾值估計方法得到的閾值能獲得更高精度最優閾值的估計值,在各種SNR條件下,與理論最優閾值的誤差均小于10%,為高精度的質心計算提供了一種較好的閾值處理方法.

圖6 (a)處理前的SHWFS典型圖像;(b)處理后的SHWFS圖像;(c)各個孔徑光斑圖像的SNRFig.6.(a)The image in the SHWFS before processed;(b)the image in the SHWFS after processed using the presented method;(c)the SNR of all the apertures.

圖7 (網刊彩色)實驗條件下各種閾值估計方法與理論最優閾值的誤差Fig.7.(color online)The error between the theoretic threshold and the threshold obtained by di ff erent methods.

表1 不同SNR條件下各種方法對應的估計閾值比較(單位:ADU)Table 1.The estimated threshold using the di ff erent methods with the di ff erent SNR(unit:ADU).
本文針對SHWFS探測系統中,子孔徑光斑呈現高斯分布、噪聲隨時間變化比較快且空間分布不均勻的特點,提出了一種利用滑動窗口均值、滑動窗口內像素局部梯度方,即像素與x-y平面夾角作為參數,構造噪聲權重函數的最優閾值估計方法.比較了這種方法以及迭代閾值法、最大值百分比閾值法和Ostu法等典型的閾值處理方法結果與理論最優閾值的誤差.仿真和實驗結果表明,與典型的閾值處理方法相比,本文提出的方法能夠得到更高精度最優閾值的估計值,并且在各種SNR條件下與理論最優閾值保持小于10%的誤差.
[1]Li J,Gong Y,Hu X R,Li C C 2014Chin.J.Laser41 0316002(in Chinese)[李晶,鞏巖,呼新榮,李春才2014中國激光41 0316002]
[2]Baik S H,Park S K,Kim C J,Cha B 2007Opt.Laser Technol.39 262
[3]Zhu Z Y,Li D Y,Hu L F,Mu Q Q,Yang C L,Cao Z L,Xuan L 2016Chin.Phys.B25 090702
[4]Gao C Q,Gao M W,Weber H 2004Chin.Phys.Lett.21 2191
[5]Wei L,Shi G H,Lu J,Yang J S,Li X Q,Zhang Y D 2013J.Opt.15 055702
[6]Chen L H,Rao C H 2011Acta Phys.Sin.60 090701(in Chinese)[陳林輝,饒長輝 2011物理學報 60 090701]
[7]Li C H,Xian H,Jiang W H,Rao C H 2007Acta Phys.Sin.56 4289(in Chinese)[李超宏,鮮浩,姜文漢,饒長輝2007物理學報56 4289]
[8]Ares J,Arines J 2001Opt.Lett.26 1831
[9]Ma X Y,Rao C H,Zheng H Q 2009Opt.Express17 8525
[10]Liang C,Liao W H,Shen J X,Zhou Y 2009Chin.J.Laser36 430(in Chinese)[梁春,廖文和,沈建新,周宇2009中國激光36 430]
[11]Ren J F,Rao C H,Li M Q 2002Opto-Electron.Eng.29 1(in Chinese)[任劍峰,饒長輝,李明全2002光電工程29 1]
[12]Thatiparthi C,Ommanib A,Burmanc R,Thapa D,Hutchings N,Lakshminarayanan V 2016Proc.SPIE9693 969321
[13]Thomas S 2004Proc.SPIE5490 1238
[14]Nightingale A M,Gordeyev S 2013Opt.Eng.52 071413
[15]Shen F,Jiang W H 1999High Power Laser and Particle Beams11 27(in Chinese)[沈鋒,姜文漢 1999強激光與粒子束11 27]
[16]Li Y K,Zhang J Z,Zhang F Z 2014Proc.SPIE9242 92421V
PACS:07.07.Df,05.40.Ca,42.25.Bs,95.75.QrDOI:10.7498/aps.66.090701
Shack-Hartmann optimum threshold estimation for the point source
Zhou Rui1)2)3)?Wei Ling1)2)Li Xin-Yang1)2)Wang Cai-Xia1)2)Li Mei1)2)Shen Feng1)2)
1)(Key Laboratory on Adaptive Optics,Chinese Academy of Sciences,Chengdu 610209,China)2)(Institute of Optics and Electronics,Chinese Academy of Sciences,Chengdu 610209,China)3)(University of Chinese Academy of Sciences,Beijing 100049,China)
8 November 2016;revised manuscript
4 February 2017)
Shack-Hartmann wavefront sensor,Gaussian spot,optimum thresh,weighted function
10.7498/aps.66.090701
The Shack-Hartmann wavefront sensor(SHWFS)is an optical detection device based on the measurements of wavefront slopes.It is widely used in an adaptive optics system due to its simple structure and strong environment adaptability.The measuring accuracy of the SHWFS depends mainly on the accuracy of the spot image centroid in each sub-aperture.There are many centroid algorithms including the center of gravity algorithm,Gauss fi tting algorithm,and correlation algorithm.As to the simplicity,robustness,high accuracy and stability,the center of gravity algorithm is more widely used.However,the accuracy of gravity algorithm is sensitive to the noise including discretization,aliasing,photon noise,readout noise,stray light,and direct current bias.To improve the accuracy of centroid,the output signals of SHWFS must be pre-processed to suppress the noise e ff ect by using the method of thresholding in general.Many threshold methods have been presented to reduce the error of centroid and there theoretically exists an optimum threshold which causes the minimum error of centroid based on the characteristics of SHWFS and noise.However,it is difficult to separate the signals from the noises,and the optimum threshold cannot be estimated accurately in real time in the SHWFS systems.In this paper aiming at noises in SHWFS,which vary with time and space rapidly,a method based on the noise weighted function of the mean value of pixels and the local gradient direction of image signals in the moving windows is presented according to the characteristics of the Gaussian spot and noise distributions.Moreover,the theory and parameters determination of the method are analyzed.The method utilizes the probability that the pixels in the moving windows belong to the noise,and the probability is inversely proportional to the mean value of pixels and the local gradient direction of image signals,and so the monotonically reducing probability function of pixels is constructed.Finally,the standard deviation and mean value of noise can be obtained,and the estimation value of optimum threshold is equal to the mean value of noise plus three times the standard deviation of noise.To investigate the e ff ects of the optimum threshold estimation with the di ff erent spot sizes,spot strengths and noise levels,the proposed algorithm is compared with traditional methods.The simulation and experimental results show that the proposed method could achieve higher accuracy,and the error between the threshold obtained by the method presented in this paper and theoretical optimum threshold is less than 10%,which is less than those from the traditional methods.
?通信作者.E-mail:zhourui@ioe.ac.cn
?Corresponding author.E-mail:zhourui@ioe.ac.cn