莫爾條紋傾角快速求解算法及應(yīng)用

陳 好，賈 靜，王 詢，何易德，蔣澤偉

（西南技術(shù)物理研究所，成都610041）

引 言

等效熱焦距參量是評(píng)估高能固體激光器熱效應(yīng)的核心參量［1-2］。隨著固體激光器輸出功率的提高，激光材料熱效應(yīng)愈加明顯，等效熱焦距參量的測(cè)量也逐漸受到重視［3-5］?；谔┎?莫爾條紋長(zhǎng)焦距測(cè)量技術(shù)是廣泛應(yīng)用于長(zhǎng)焦距光學(xué)元件焦距參量精確測(cè)量的方法，也是實(shí)現(xiàn)等效熱焦距參量精確測(cè)量的重要途徑之一［6］。它的基本原理是一束平行光束連續(xù)通過(guò)兩個(gè)相同周期的Rochi光柵，利用一個(gè)光柵的衍射自成像與第2個(gè)光柵的疊加形成莫爾條紋，通過(guò)莫爾條紋傾角測(cè)得被測(cè)元件焦距值。此方法具有測(cè)量范圍大、測(cè)量精度高、光路簡(jiǎn)單、操作方便等優(yōu)點(diǎn)［7-9］。莫爾條紋傾角求解是該技術(shù)的核心之一，是測(cè)量系統(tǒng)測(cè)量精度和測(cè)量速度的主要影響因素。其中莫爾條紋傾角測(cè)量誤差引起的焦距誤差占整個(gè)焦距測(cè)量誤差的90%［10］。為提高莫爾條紋傾角的求解精度，浙江大學(xué)HOU等人［11］提出了精確計(jì)算莫爾條紋傾角的傅里葉迭代算法；YU在他的基礎(chǔ)上提出了基于信息熵的迭代算法［12］。這兩種方法均通過(guò)傅里葉變換求取莫爾條紋頻譜，并通過(guò)迭代運(yùn)算逐步逼近1級(jí)頻譜的精確坐標(biāo)，進(jìn)而求解莫爾條紋傾角。算法突破數(shù)字圖像像素的限制，顯著提高莫爾條紋的測(cè)量精度，傾角求解精度可達(dá)10-4rad量級(jí)。

迭代運(yùn)算的引入在提高測(cè)量精度的同時(shí)也降低了測(cè)量速度，當(dāng)?shù)螖?shù)為10次時(shí)，采用常規(guī)配置計(jì)算機(jī)完成一次測(cè)量用時(shí)15s左右；基于信息熵的迭代算法在迭代20次時(shí)完成一次測(cè)量用時(shí)47s［12］。對(duì)于光學(xué)元件焦距的測(cè)量，測(cè)量時(shí)間不是主要問(wèn)題，但激光材料熱效應(yīng)等效熱焦距參量隨抽運(yùn)時(shí)間快速變化，單次測(cè)量時(shí)間過(guò)長(zhǎng)將無(wú)法反應(yīng)熱效應(yīng)的變化過(guò)程，需要以大于1Hz的頻率進(jìn)行測(cè)量［13-15］。為提高激光材料熱效應(yīng)等效熱焦距參量的測(cè)量速度，本文中在對(duì)頻域迭代法測(cè)量速度影響因素分析的基礎(chǔ)上，基于迭代運(yùn)算坐標(biāo)位置的相似性，將離散傅里葉變換公式進(jìn)行變形，將頻譜坐標(biāo)中包含小數(shù)值坐標(biāo)點(diǎn)進(jìn)行公共項(xiàng)分離，公共部分采用傅里葉變換系數(shù)進(jìn)行濾波，不同部分在提取相似項(xiàng)的基礎(chǔ)上進(jìn)行組合運(yùn)算，將2維條紋圖像降為1維數(shù)據(jù)進(jìn)行處理，整體上減少運(yùn)算復(fù)雜度。同時(shí)坐標(biāo)點(diǎn)的分離將不同部分限定在極小的可知范圍，極大地方便了查表運(yùn)算和位運(yùn)算的采用，提高了算法執(zhí)行效率。編寫(xiě)軟件，通過(guò)實(shí)驗(yàn)驗(yàn)證了快速求解算法的可行性和運(yùn)算速度，結(jié)果表明，在保證測(cè)量精度的前提下完成一次測(cè)量耗時(shí)為0.4s，很大地提高了莫爾條紋傾角的求解速度，滿足了等效熱焦距參量測(cè)量對(duì)測(cè)量速度的需求。

1 莫爾條紋傾角快速求解算法分析

莫爾條紋是由兩塊等間隔排列的直線族或曲線族疊加產(chǎn)生的周期性條紋圖像。在基于泰伯-莫爾條紋長(zhǎng)焦距測(cè)量技術(shù)中由一個(gè)光柵的衍射自成像與第2個(gè)光柵的疊加形成，其傾角是莫爾條紋與圖像坐標(biāo)軸的夾角，直接反映被測(cè)元件的焦距值。在莫爾條紋頻譜域中，條紋傾角與其頻譜±1級(jí)的連線存在π／2的夾角，因此，莫爾條紋傾角的求解可通過(guò)求解頻譜的±1級(jí)坐標(biāo)求解。對(duì)莫爾條紋圖像進(jìn)行傅里葉變換即可得到莫爾條紋頻譜。但CCD采集到的莫爾條紋圖像為數(shù)字圖像，受CCD分辨率的限制，僅進(jìn)行快速離散傅里葉變換只能求解整數(shù)像素點(diǎn)坐標(biāo)，若條紋數(shù)量為N，則算法測(cè)量精度僅為（1／N）rad，測(cè)量精度無(wú)法達(dá)到要求。

1.1 頻域迭代原理

為方便莫爾條紋傾角求解，一般CCD采集的莫爾條紋行列方向像素相同，設(shè)為M，則可采用離散傅里葉變換公式在整數(shù)點(diǎn)坐標(biāo)的基礎(chǔ)上求解單個(gè)像素內(nèi)含分?jǐn)?shù)點(diǎn)坐標(biāo)（u，v）的頻譜值。將坐標(biāo)點(diǎn)的求解精度提高到亞像素級(jí)別，進(jìn)而可將莫爾條紋傾角的求解精度提高100倍左右。

式中，f（x，y）為條紋圖像灰度值，F(xiàn)（u，v）為條紋圖像頻譜。

由于坐標(biāo)值（u，v）存在小數(shù)部分，無(wú)法采用快速傅里葉變換算法，運(yùn)用（1）式計(jì)算像素內(nèi)各位置的頻譜就需要龐大的計(jì)算量，計(jì)算速度必然很慢。莫爾條紋圖像是按一定規(guī)律分布的標(biāo)準(zhǔn)余弦圖像，由于系統(tǒng)誤差的存在，其頻譜圖的+1級(jí)在整像素內(nèi)是按以最大值為中心的2維高斯分布。這個(gè)特點(diǎn)說(shuō)明可以通過(guò)迭代逐漸逼近的方法找到其最大值，從而減少運(yùn)算量?；诖祟l域迭代算法采用像素分割，比較定位的方法逼近最大值。首先對(duì)莫爾條紋進(jìn)行快速傅里葉變換得到其頻譜圖，確定頻譜圖的+1級(jí)坐標(biāo)點(diǎn)，也即對(duì)頻域的+1最大值進(jìn)行整數(shù)像素點(diǎn)初定位，得到其序列值，設(shè)其序列值為（m，n）；然后以頻域的+1級(jí)點(diǎn)為中心，把最大值點(diǎn)所在像素按九宮格進(jìn)行分割，九宮格像素分割法如圖1所示。

Fig.1 Pixel segmentation based on nine palace lattice

九宮格內(nèi)各宮格的坐標(biāo)值為：（m+1／3，n+1／3），（m+1／3，n），（m+1／3，n-1／3），（m，n+1／3），（m，n），（m，n-1／3），（m-1／3，n+1／3），（m-1／3，n），（m-1／3，n-1／3）。采用（1）式求解各宮格中心處的頻譜值，比較各個(gè)頻譜值，確定最大值的坐標(biāo)值（u0，v0）；重復(fù)上述方法，以每次得到的最大值的坐標(biāo)值為中心，按步長(zhǎng)1／3k進(jìn)行九宮格分割完成第k次迭代運(yùn)算，再求解最大值，這樣經(jīng)過(guò)多次分割后逐漸逼近+1級(jí)頻譜最大值，達(dá)到提高測(cè)量精度的目的。算法覆蓋了整像素內(nèi)的整個(gè)區(qū)域，也相當(dāng)于求解了整個(gè)區(qū)域頻譜值。

1.2 快速求解算法原理

采用頻域迭代法完成莫爾條紋傾角測(cè)量，當(dāng)?shù)?0次需要計(jì)算90個(gè)頻譜值，減少了很多次計(jì)算，但由于最大值坐標(biāo)點(diǎn)包含小數(shù)部分，頻譜功率值計(jì)算只能采用指數(shù)運(yùn)算，算法執(zhí)行效率很低，完成一次測(cè)量仍然需要花費(fèi)十幾秒的時(shí)間。

為進(jìn)一步提高莫爾條紋傾角求解速度，本文中從頻譜值的求解上尋找方法。分析頻域迭法頻譜值的計(jì)算過(guò)程知，9個(gè)坐標(biāo)點(diǎn)均是由上次迭代頻譜值最大值點(diǎn)坐標(biāo)（u0，v0）與步長(zhǎng) 1／3k的和差，因此對(duì)于本次迭代運(yùn)算，上次迭代最大值點(diǎn)坐標(biāo)是公共部分，步長(zhǎng)為不同部分。將每一個(gè)頻譜點(diǎn)坐標(biāo)（u，v）寫(xiě)為u＝u0+uk，v＝v0+vk，其中（u0，v0）為上次迭代頻譜值最大值點(diǎn)坐標(biāo)，（uk，vk）為本次運(yùn)算的迭代步長(zhǎng)，因此可將（1）式改寫(xiě)為：

若設(shè)：

基于（4）式對(duì)迭代運(yùn)算中各點(diǎn)頻譜值進(jìn)行求解，f′（x，y）為一次迭代公共部分，采用傅里葉變換指數(shù)項(xiàng)對(duì)莫爾條紋圖像進(jìn)行濾波運(yùn)算完成。迭代運(yùn)算中，w（x，y）的坐標(biāo)點(diǎn)（uk，vk）分別為（1／3k，1／3k），（1／3k，0），（1／3k，-1／3k），（0，1／3k），（0，0），（0，-1／3k），（-1／3k，1／3k），（-1／3k，0），（-1／3k，-1／3k），其中 k為迭代次數(shù)。對(duì)于坐標(biāo)點(diǎn)（1／3k，1／3k），（-1／3k，-1／3k），若 x+y相同，則 w（x，y）相等，因此可先計(jì)算相同點(diǎn)f′（x，y）的和，將M×M 2維數(shù)據(jù)降為長(zhǎng)度為2M的1維數(shù)據(jù)參與計(jì)算。同理，對(duì)于（1／3k，0），（-1／3k，0）可先計(jì)算 x相同點(diǎn) f′（x，y）的和；對(duì)于（0，1／3k），（0，-1／3k）可先計(jì)算 y相同點(diǎn) f′（x，y）的和；對(duì)于（1／3k，-1／3k），（-1／3k，1／3k）可先計(jì)算 x-y相同點(diǎn) f′（x，y）的和，將M×M 2維數(shù)據(jù)降為4組1維數(shù)據(jù)參與計(jì)算。因此（4）式可變?yōu)椋?/p>

式中，H（n）為 f′（x，y）滿足坐標(biāo)關(guān)系的條紋灰度值和。E（n）為長(zhǎng)度為2M的指數(shù)項(xiàng)，對(duì)于每次迭代運(yùn)算E（n）的值固定，可建立數(shù)據(jù)表，在計(jì)算 F（u，v）過(guò)程中采用查表運(yùn)算和位運(yùn)算，極大地提高了算法運(yùn)算效率。

Fig.2 Flow chart of fast algorithm

快速求解算法基于迭代運(yùn)算坐標(biāo)位置的相似性，將離散傅里葉變換公式進(jìn)行變形，將頻譜坐標(biāo)分為公共項(xiàng)和不同部分，公共部分采用傅里葉變換系數(shù)進(jìn)行濾波，不同部分在提取相似項(xiàng)的基礎(chǔ)上進(jìn)行組合運(yùn)算，將2維條紋圖像降為1維數(shù)據(jù)，同時(shí)坐標(biāo)點(diǎn)的分離將不同部分限定在極小的可知范圍，極大地方便了查表運(yùn)算和位運(yùn)算的采用。算法整體上提高了莫爾條紋傾角的求解速度?？焖偾蠼馑惴鞒倘鐖D2所示。步驟如下：（1）在測(cè)量開(kāi)始前建立20次迭代運(yùn)算中所有指數(shù)項(xiàng)E（n）的數(shù)據(jù)表；（2）對(duì)莫爾條紋圖像進(jìn)行快速傅里葉變換（fast Fourier transform，F(xiàn)FT），求解零頻坐標(biāo)（m0，n0）和 +1級(jí)頻譜整像素點(diǎn)坐標(biāo)值（m，n）；（3）代入整像素點(diǎn)坐標(biāo)（m，n），求解一次迭代運(yùn)算公共部分f′（x，y）的值，同時(shí)求解滿足 x+y，x，y，x-y相等坐標(biāo)點(diǎn) f′（x，y）的值和 H（n）；（4）分別求解九宮格中各坐標(biāo)點(diǎn)的頻譜值，并求解頻譜最大值的坐標(biāo)（mk，nk）；（5）代入坐標(biāo)（mk，nk），求解 k+1次迭代運(yùn)算公共部分 f′（x，y）的值，重復(fù)步驟（3）和步驟（4），求解第 n次迭代運(yùn)算坐標(biāo)點(diǎn)計(jì)算莫爾條紋傾角。

2 實(shí)驗(yàn)結(jié)果與分析

構(gòu)建了仿真實(shí)驗(yàn)和等效熱焦距測(cè)量實(shí)驗(yàn)，分別對(duì)快速求解算法的測(cè)量精度和測(cè)量速度進(jìn)行驗(yàn)證。其中，通過(guò)仿真實(shí)驗(yàn)?zāi)M不同角度莫爾條紋，采用傳統(tǒng)頻域迭代算法和快速求解算法對(duì)模擬莫爾條紋的傾角進(jìn)行測(cè)量，三者進(jìn)行對(duì)比驗(yàn)證快速求解算法的測(cè)量精度；通過(guò)等效熱焦距測(cè)量實(shí)驗(yàn)采集莫爾條紋，對(duì)兩種測(cè)量算法的測(cè)量速度進(jìn)行對(duì)比，驗(yàn)證快速求解算法的測(cè)量速度，同時(shí)驗(yàn)證快速求解算法對(duì)激光材料熱效應(yīng)等效熱焦距參量測(cè)量的適用性。

2.1 仿真實(shí)驗(yàn)

基于MATLAB軟件構(gòu)建了理想莫爾條紋模擬程序和算法仿真程序。模擬了傾角為 5°，15°，30°，45°，60°，75°，85°，90°的莫爾條紋，其中傾角為 45°的模擬莫爾條紋如圖3所示。

Fig.3 Moiréfringe with 45°angle

分別采用傳統(tǒng)頻域迭代算法和快速求解算法對(duì)莫爾條紋傾角進(jìn)行求解，迭代10次，測(cè)量結(jié)果如表1所示。

由表1可知，傳統(tǒng)算法和快速求解算法的測(cè)量結(jié)果完全一致，而且和模擬角度均小于0.002°。排除圖像離散化過(guò)程引入的測(cè)量誤差，兩種算法均可精確求解莫爾條紋傾角，測(cè)量誤差小于0.002°，證明了快速求解算法的測(cè)量精度。

Table 1 Results with different algorithms

2.2 物理實(shí)驗(yàn)

采用實(shí)驗(yàn)室基于泰伯-莫爾條紋技術(shù)研制的長(zhǎng)焦距測(cè)量?jī)x對(duì)500J釹玻璃激光器激光材料熱效應(yīng)等效熱焦距參量進(jìn)行測(cè)量。測(cè)量算法以長(zhǎng)焦距測(cè)量?jī)x測(cè)量軟件為基礎(chǔ)，采用Visual Studio 2008編碼實(shí)現(xiàn)，計(jì)算機(jī)為CPU Intel i3處理器（主頻3.5GHz）及內(nèi)存4GB的工控機(jī)。通過(guò)對(duì)比兩種測(cè)量算法的測(cè)量時(shí)間驗(yàn)證快速求解算法的測(cè)量速度，并驗(yàn)證測(cè)量算法對(duì)等效熱焦距參量的測(cè)量的有效性。CCD采集的莫爾條紋見(jiàn)圖4。

Fig.4 Moiréfringe

采集任意狀態(tài)下12幅莫爾條紋圖，采用傳統(tǒng)頻域迭代算法和快速求解算法對(duì)條紋傾角進(jìn)行求解，測(cè)量時(shí)間如表2所示。

由表2可知，采用傳統(tǒng)算法完成一次測(cè)量平均測(cè)量時(shí)間為15.0s，采用快速求解算法完成一次測(cè)量平均測(cè)量時(shí)間為0.4s，前者測(cè)量時(shí)間是后者的38倍，直接說(shuō)明快速求解算法測(cè)量速度相對(duì)于傳統(tǒng)算法得到很大的提高，驗(yàn)證了快速求解算法的測(cè)量速度。

對(duì)實(shí)驗(yàn)室500J釹玻璃激光器在完成單次脈沖發(fā)射過(guò)程中激光材料熱效應(yīng)等效熱焦距參量進(jìn)行測(cè)量，測(cè)量間隔時(shí)間設(shè)置為0.5s，測(cè)量結(jié)果分布曲線見(jiàn)圖5。

由圖5可知，等效熱焦距參量隨時(shí)間快速變化，在脈沖發(fā)射后10s內(nèi)，等效熱焦距參量由-158.4m快速減小到-3400m然后變?yōu)?400m后再次快速減小到293.3m；10s以后等效熱焦距變換相對(duì)緩慢，逐漸減小至130m后逐漸增大，到50s時(shí)增大到165m。從測(cè)量結(jié)果知，采用快速求解算法將測(cè)量時(shí)間間隔設(shè)置為0.5s得到的測(cè)量結(jié)果很好地反映了脈沖發(fā)射后等效熱焦距參量的變化過(guò)程，從而證明了快速求解算法對(duì)激光材料熱效應(yīng)等效熱焦距參量測(cè)量的適用性。

Table 2 Measurement time results with different algorithms

Fig.5 Result curve of thermal focal length

3 結(jié) 論

針對(duì)基于泰伯-莫爾條紋長(zhǎng)焦距測(cè)量技術(shù)在測(cè)量激光材料熱效應(yīng)等效熱焦距時(shí)對(duì)測(cè)量速度的需求，提出了一種莫爾條紋傾角快速求解算法，并對(duì)該算法的原理和流程進(jìn)行了分析，整體上減少了算法復(fù)雜度，提高了莫爾條紋傾角的測(cè)量速度。建立了仿真實(shí)驗(yàn)，通過(guò)與模擬傾角和傳統(tǒng)測(cè)量算法測(cè)量結(jié)果進(jìn)行對(duì)比驗(yàn)證了快速求解算法的測(cè)量精度。建立物理實(shí)驗(yàn)，采用相同配置計(jì)算機(jī)，在保證測(cè)量精度的前提下傳統(tǒng)測(cè)量算法平均測(cè)量時(shí)間為15s，快速求解算法的平均測(cè)量時(shí)間為0.4s，測(cè)量速度提高了38倍，證明了該算法的測(cè)量速度。對(duì)實(shí)驗(yàn)室500J釹玻璃激光器在完成單次脈沖發(fā)射過(guò)程中激光材料熱效應(yīng)等效熱焦距參量進(jìn)行了測(cè)量，測(cè)量結(jié)果證明快速求解算法滿足等效熱焦距參量的測(cè)量需求。