作者：Steven
著作權宣告：著作權歸作者所有，商業轉載請聯系作者獲得授權，非商業轉載請注明出處

注：本文所講內容為本人碩士畢業論文：《基于干涉影像質量分析的激光干涉儀抗振技術研究》，如有參考需要標注來源哈，如有疑問可以評論回復也可以郵箱1257147469@qq.com聯系我，

前言

自二十世紀末以來，移相干涉技術(Phase Shifting Interferometry, PSI)因其高解析度、高計算精度的優勢，逐漸取代了傳統的波面測量技術并被作為一種標準波面測量技術沿用至今，PSI是以光波波長作單位的納米級超高精度測量手段，通過移相式激光干涉儀實作，被廣泛應用在光學領域中，特別是在光學系統成像質量評價和光學器件面形檢測方面，

本文介紹了移相干涉技術中最基礎卻也非常重要的一步——相位提取，主要闡述了移相干涉測量原理、四步移相法提取相位、多步平均法推導程序、多步解包裹后平均法這四個部分，希望能給同樣從事該領域研究的你帶來一點幫助，

一、移相干涉測量原理

過去的干涉測量技術是通過人的肉眼或者相機拍攝，來直觀判斷干涉圖中條紋特征進而完成測量，該方法的不穩定因素（比如人的主觀意志）很多，其精度誤差在λ/10左右；現代干涉測量技術通過將電子技術、計算機技術、光電影像處理技術、演算法模型和高精度光學元件結合，大大推動了干涉測量領域的發展，使其進入高精度時代，移相干涉技術（PSI）作為現代干涉測量技術的代表，已經成為干涉測量領域學者的研究重點，

移相干涉技術測量的原理：移相式激光干涉儀通過控制壓電陶瓷驅動器（PZT）移動參考鏡位置或者改變激光器波長的方式，令參考光和測驗光的光程差改變，從干涉影像中可看出干涉條紋進行了相應的位移，在移相程序中，通過光電探測器（CCD或者CMOS）采集不同移相量下的干涉圖，圖1所示為一組移相間隔為90°的移相干涉圖組，計算機根據特定的數學演算法和模型對干涉圖組的像素資料進行系列運算，進而可以求得被測元件的面形資訊，并根據計算好的相關引數進行質量評價，

基于PZT實作移相的程序為：當PZT控制模塊發出指定命令時，PZT伸長量進行相應改變，即參考鏡和被測鏡的相對距離改變，進而完成相位調制以達到移相目的，

在上述移相干涉測量程序中，干涉影像中的條紋分布是由參考光束和測驗光束疊加形成的干涉場導致，干涉場的光強資訊與干涉圖中像素點的灰度資訊相關，其分布函式公式(1)可表示為：

$\mathrm{I}_{i}(\mathrm{x}, \mathrm{y})=A(\mathrm{x}, \mathrm{y})+B(\mathrm{x}, \mathrm{y}) \cos \left(\varphi(\mathrm{x}, \mathrm{y})+\delta_{i}\right)$

(1)

其中：A(x,y)是干涉圖的背景光強，B(x,y)是調制度， $\varphi(\mathrm{x}, \mathrm{y})$ 是待測物的波面初始相位資訊， $\delta_{i}$ 是移相量，(x,y)是干涉圖中像素點的坐標資訊，

參考鏡和被測鏡間的面形差P(x,y)與位相分布的關系可表示為公式(2)：

$P(\mathrm{x}, \mathrm{y})=\frac{\varphi(\mathrm{x}, \mathrm{y})}{4 \pi} \lambda$

(2)

λ為激光器波長（本文為He-Ne激光器——632.8nm），在移相干涉測量技術中，測量單位為λ，當通過解相演算法處理干涉圖組的像素資料，可計算出 $\mathrm{I}_{i}(\mathrm{x}, \mathrm{y})$ ，并進一步得到面形資訊 $\varphi(\mathrm{x}, \mathrm{y})$ 以完成檢測作業，

下面開始介紹3種比較常見和基礎的相位提取演算法，

二、四步移相法

作為經典的解相演算法，四步移相法可以有效計算出理想移相圖組的相位分布，先獲取無移相的干涉影像，再進行三次移相間隔為π/2的移相，每次移相完成采集1幀影像，共得到4幀相位相差π/2的干涉影像，其影像光強資訊分別為：

$\begin{array}{l} \mathrm{I}_{1}(\mathrm{x}, \mathrm{y})=A(\mathrm{x}, \mathrm{y})+B(\mathrm{x}, \mathrm{y}) \cos (\varphi(\mathrm{x}, \mathrm{y})) \\ \mathrm{I}_{2}(\mathrm{x}, \mathrm{y})=A(\mathrm{x}, \mathrm{y})+B(\mathrm{x}, \mathrm{y}) \cos (\varphi(\mathrm{x}, \mathrm{y})+\pi / 2) \\ \mathrm{I}_{3}(\mathrm{x}, \mathrm{y})=A(\mathrm{x}, \mathrm{y})+B(\mathrm{x}, \mathrm{y}) \cos (\varphi(\mathrm{x}, \mathrm{y})+\pi) \\ \mathrm{I}_{4}(\mathrm{x}, \mathrm{y})=A(\mathrm{x}, \mathrm{y})+B(\mathrm{x}, \mathrm{y}) \cos (\varphi(\mathrm{x}, \mathrm{y})+3 \pi / 2) \end{array}$

(3)

公式(3)中4個式子進行公式變換消去A和B可得到和四幀影像像素資訊的關系式：

$\varphi(\mathrm{x}, \mathrm{y})=\arctan \left(\frac{\mathrm{I}_{4}(\mathrm{x}, \mathrm{y})-\cdot \mathrm{I}_{2}(\mathrm{x}, \mathrm{y})}{\mathrm{I}_{1}(\mathrm{x}, \mathrm{y})-{ }^{\cdot} \mathrm{I}_{3}(\mathrm{x}, \mathrm{y})}\right)$

(4）

該演算法求解速度快，適合對理想干涉圖組進行求解，當有外部干擾時，相位變化呈無規則態，此時四步移相法的解相效果非常糟糕，無任何抗振能力，與四步移相法類似的單周期定步長演算法還有三步演算法等，

三、多步平均法（5-9步）

為了抑制移相器線性移相誤差或外在環境干擾帶來的影響，相關學者在四步移相法的基礎上引入了平均技術（本文介紹的平均法基于四步移相法，移相間隔均為π/2）：通過獲取T+K（T=4）幀移相間隔為π/2的干涉影像，使干涉圖組的相位變化不再局限于一個2π周期內，避免因單周期內干擾過大而導致求解偏差較大的情況發生；鑒于余弦分布的特點，不同周期下干涉圖的移相量均可化為單周期移相值，即0、π/2、π、3π/2，每個移相值對應多幀不同周期下的干涉圖，對這些干涉圖而言，其光強資訊在理論上應完全一致，考慮到現實環境干擾的復雜性和無規則性，每幀圖的資訊都會出現或多或少的偏差，應用平均技術的思想能一定程度上起到減少誤差的作用，

基于平均技術思想，多步平均法雖能有效抑制振動干擾，但當幀數過多時，也可能會出現因取圖時間過長而引入過量振動的問題，最佳幀數的設定是多步平均法非常關鍵的一步，綜合考慮干涉儀內置相機的性能（如果CCD或CMOS的幀率很高，可以適當增加獲取影像的幀數），一般獲取9幀影像最佳，既可以保證圖組的相位資訊跨過兩個周期以使平均法達到抑制振動影響的效果，又可以避免引入過量的擾動，

接下來推導九步平均演算法的解相公式：

a.將四步移相法公式(4)變形得公式(5)：

$\varphi(\mathrm{x}, \mathrm{y})=\arctan \left(\frac{\mathrm{I}_{4}(\mathrm{x}, \mathrm{y})- \mathrm{I}_{2}(\mathrm{x}, \mathrm{y})}{\mathrm{I}_{1}(\mathrm{x}, \mathrm{y})- \mathrm{I}_{3}(\mathrm{x}, \mathrm{y})}\right)=\arctan \left(\frac{N_{0}}{D_{0}}\right)$

(5）

上式中， $N_{0}$ 和 $D_{0}$ 分別代表從第一幀圖（即移相量從0開始）開始的4幀干涉影像間解相運算的分子分母部分，四步移相法作為后續平均演算法的基本演算法式，基礎步數為T（=4），

如果繼續獲取第五幀圖，該干涉圖的移相量為2π，化為單周期移相值為0，此時將以第二幀圖（即移相量從）開始的4幀干涉圖作四步解相運算，可得公式(6)：

$\varphi(\mathrm{x}, \mathrm{y})=\arctan \left(\frac{\mathrm{I}_{4}(\mathrm{x}, \mathrm{y})- \mathrm{I}_{2}(\mathrm{x}, \mathrm{y})}{\mathrm{I}_{5}(\mathrm{x}, \mathrm{y})- \mathrm{I}_{3}(\mathrm{x}, \mathrm{y})}\right)=\arctan \left(\frac{N_{\pi/2 }}{D_{\pi/2}}\right)$

(6）

將公式(5)和(6)結合可得五步（T+1）平均演算法，公式(7)：

$\begin{aligned} \varphi(\mathrm{x}, \mathrm{y}) &=\arctan \left(\frac{N_{0}+N \frac{\pi}{2}}{D_{0}+D_{\frac{\pi}{2}}}\right) \\ &=\arctan \left(\frac{2 \cdot\left( \mathrm{I}_{4}(\mathrm{x}, \mathrm{y})-\mathrm{I}_{2}(\mathrm{x}, \mathrm{y})\right)}{\mathrm{I}_{1}(\mathrm{x}, \mathrm{y})+\mathrm{I}_{5}(\mathrm{x}, \mathrm{y})-2 \mathrm{I}_{3}(\mathrm{x}, \mathrm{y})}\right) \end{aligned}$

(7）

b.繼續獲取第六幀圖，該干涉圖的移相量為2π+π/2，化為單周期移相值為π/2，同理可獲得公式(8)：

$\varphi(\mathrm{x}, \mathrm{y})=\arctan \left(\frac{\mathrm{I}_{4}(\mathrm{x}, \mathrm{y})-\mathrm{I}_{6}(\mathrm{x}, \mathrm{y})}{\mathrm{I}_{5}(\mathrm{x}, \mathrm{y})-\mathrm{I}_{3}(\mathrm{x}, \mathrm{y})}\right)=\arctan \left(\frac{N_{\pi}}{D_{\pi}}\right)$

(8）

此時運用“擴展平均”思想以獲得進一步的補償，即運用N+2的處理方式：分別對前T+1（5）幀和后T+1（5）幀圖采用五步平均演算法，累加平均得到六步（T+2）平均演算法（也可稱作擴展平均演算法，后文將其統一作平均演算法），公式(9)：

$\begin{array}{c} \varphi(\mathrm{x}, \mathrm{y})=\arctan \left(\frac{N_{0}+N_{\frac{\pi}{2}}+\left(N_{0}+N_{\frac{\pi}{2}}\right) \frac{\pi}{2}}{D_{0}+D_{\frac{\pi}{2}}+\left(D_{0}+D_{\frac{\pi}{2}}\right)_{\frac{\pi}{2}}}\right)=\arctan \left(\frac{N_{0}+2 N_{\frac{\pi}{2}}+N_{\pi}}{D_{0}+2 D_{\frac{\pi}{2}}+D_{\pi}}\right) \\ =\arctan \left(\frac{4 \cdot \mathrm{I}_{4}(\mathrm{x}, \mathrm{y})-3 \cdot \mathrm{I}_{2}(\mathrm{x}, \mathrm{y})-\mathrm{I}_{6}(\mathrm{x}, \mathrm{y})}{\mathrm{I}_{1}(\mathrm{x}, \mathrm{y})-4 \cdot \mathrm{I}_{3}(\mathrm{x}, \mathrm{y})+3 \cdot \mathrm{I}_{5}(\mathrm{x}, \mathrm{y})}\right) \end{array}$

(9）

c.以此類推，繼續擴展可得到T+K步平均演算法，

當K=3時，獲取到第七幀圖，四到七幀圖作四步運算得公式(10)：

$\varphi(\mathrm{x}, \mathrm{y})=\arctan \left(\frac{\mathrm{I}_{4}(\mathrm{x}, \mathrm{y})- \mathrm{I}_{6}(\mathrm{x}, \mathrm{y})}{\mathrm{I}_{5}(\mathrm{x}, \mathrm{y})- \mathrm{I}_{7}(\mathrm{x}, \mathrm{y})}\right)=\arctan \left(\frac{N_{3 \pi / 2}}{D_{3 \pi / 2}}\right)$

(10）

分別對前T+2（6）幀和后T+2（6）幀作六步平均運算并累加可得七步（T+3）平均演算法，公式(11)：

$\begin{aligned} \varphi(\mathrm{x}, \mathrm{y}) &=\arctan \left(\frac{N_{0}+2 N_{\frac{\pi}{2}}+N_{\pi}+\left(N_{0}+2 N_{\frac{\pi}{2}}+N_{\pi}\right) \frac{\pi}{2}}{D_{0}+2 D_{\frac{\pi}{2}}+D_{\pi}+\left(D_{0}+2 D_{\frac{\pi}{2}}+D_{\pi}\right) \frac{\pi}{2}}\right) \\ =& \arctan \left(\frac{N_{0}+3 N_{\frac{\pi}{2}}+3 N_{\pi}+N_{3 \pi / 2}}{D_{0}+3 D_{\frac{\pi}{2}}+3 D_{\pi}+D_{3 \pi / 2}}\right) \cdots \cdots . \\ =& \arctan \left(\frac{4 \cdot\left(2 \cdot \mathrm{I}_{4}(\mathrm{x}, \mathrm{y})-\mathrm{I}_{2}(\mathrm{x}, \mathrm{y})-\mathrm{I}_{6}(\mathrm{x}, \mathrm{y})\right)}{\mathrm{I}_{1}(\mathrm{x}, \mathrm{y})-7 \cdot \mathrm{I}_{3}(\mathrm{x}, \mathrm{y})+\cdot 7 \cdot \mathrm{I}_{5}(\mathrm{x}, \mathrm{y})-\mathrm{I}_{7}(\mathrm{x}, \mathrm{y})}\right) \end{aligned}$

(11）

當K=4時，獲取到第八幀圖，五到八幀圖作四步運算得公式(12)：

$\varphi(\mathrm{x}, \mathrm{y})=\arctan \left(\frac{\mathrm{I}_{8}(\mathrm{x}, \mathrm{y})-\cdot \mathrm{I}_{6}(\mathrm{x}, \mathrm{y})}{\mathrm{I}_{5}(\mathrm{x}, \mathrm{y})-\cdot \mathrm{I}_{7}(\mathrm{x}, \mathrm{y})}\right)=\arctan \left(\frac{N_{2 \pi}}{D_{2 \pi}}\right)$

(12）

分別對前T+3（7）幀和后T+3（7）幀作七步平均運算并累加可得八步（T+4）平均演算法，公式(13)：

$\begin{array}{l} \varphi(\mathrm{x}, \mathrm{y}) \\ =\arctan \left(\frac{N_{0}+3 N \frac{\pi}{2}+3 N_{\pi}+N_{3 \pi / 2}+\left(N_{0}+3 N \frac{\pi}{2}+3 N_{\pi}+N_{3 \pi / 2}\right) \frac{\pi}{2}}{\left.D_{0}+3 D_{\frac{\pi}{2}}+3 D_{\pi}+D_{3 \pi / 2}+D_{0}+3 D_{\frac{\pi}{2}}+3 D_{\pi}+D_{3 \pi / 2}\right) \frac{\pi}{2}}\right) \\ =\arctan \left(\frac{N_{0}+4 N_{\frac{\pi}{2}}+6 N_{\pi}+4 N_{3 \pi / 2}+N_{2 \pi}}{D_{0}+4 D_{\frac{\pi}{2}}+6 D_{\pi}+4 D_{3 \pi / 2}+D_{2 \pi}}\right) \\ =\arctan \left(\frac{15 \cdot \mathrm{I}_{4}(\mathrm{x}, \mathrm{y})-\cdot 5 \cdot \mathrm{I}_{2}(\mathrm{x}, \mathrm{y})-11 \cdot \mathrm{I}_{6}(\mathrm{x}, \mathrm{y})+\mathrm{I}_{8}(\mathrm{x}, \mathrm{y})}{\mathrm{I}_{1}(\mathrm{x}, \mathrm{y})-11 \cdot \mathrm{I}_{3}(\mathrm{x}, \mathrm{y})+\cdot 15 \cdot \mathrm{I}_{5}(\mathrm{x}, \mathrm{y})-5 \cdot \mathrm{I}_{7}(\mathrm{x}, \mathrm{y})}\right) \end{array}$

(13）

當K=5時，獲取到第九幀圖，六到九幀圖作四步運算得公式(14)：

$\varphi(\mathrm{x}, \mathrm{y})=\arctan \left(\frac{\mathrm{I}_{8}(\mathrm{x}, \mathrm{y})-\cdot \mathrm{I}_{6}(\mathrm{x}, \mathrm{y})}{\mathrm{I}_{9}(\mathrm{x}, \mathrm{y})-\mathrm{I}_{7}(\mathrm{x}, \mathrm{y})}\right)=\arctan \left(\frac{N_{5 \pi / 2}}{D_{5 \pi / 2}}\right)$

(14）

分別對前T+4（8）幀和后T+4（8）幀作八步平均運算并累加可得九步（T+5）平均演算法，公式(15)：

$\begin{array}{l} \varphi(\mathrm{x}, \mathrm{y}) \\ =\arctan \left(\frac{N_{0}+4 N_{\frac{\pi}{2}}+6 N_{\pi}+4 N_{3 \pi / 2}+N_{2 \pi}+\left(N_{0}+4 N_{\frac{\pi}{2}}+6 N_{\pi}+4 N_{3 \pi / 2}+N_{2 \pi}\right) \frac{\pi}{2}}{D_{0}+4 D_{\frac{\pi}{2}}+6 D_{\pi}+4 D_{3 \pi / 2}+D_{2 \pi}+\left(D_{0}+4 D_{\frac{\pi}{2}}+6 D_{\pi}+4 D_{3 \pi / 2}+D_{2 \pi}\right) \frac{\pi}{2}}\right) \\ =\arctan \left(\frac{N_{0}+5 N_{\frac{\pi}{2}}+10 N_{\pi}+10 N_{3 \pi / 2}+5 N_{2 \pi}+N_{5 \pi / 2}}{D_{0}+5 D_{\frac{\pi}{2}}+10 D_{\pi}+10 D_{3 \pi / 2}+5 D_{2 \pi}+D_{5 \pi / 2}}\right) \\ \end{array}$

$\varphi(\mathrm{x}, \mathrm{y})=\arctan \left(\frac{26 \cdot \mathrm{I}_{4}(\mathrm{x}, \mathrm{y})-\cdot 6 \cdot \mathrm{I}_{2}(\mathrm{x}, \mathrm{y})-26 \cdot \mathrm{I}_{6}(\mathrm{x}, \mathrm{y})+6 \cdot \mathrm{I}_{8}(\mathrm{x}, \mathrm{y})}{\mathrm{I}_{1}(\mathrm{x}, \mathrm{y})-16 \cdot \mathrm{I}_{3}(\mathrm{x}, \mathrm{y})+\cdot 30 \cdot \mathrm{I}_{5}(\mathrm{x}, \mathrm{y})-16 \cdot \mathrm{I}_{7}(\mathrm{x}, \mathrm{y})+\mathrm{I}_{9}(\mathrm{x}, \mathrm{y})}\right)$

(15）

至此，九步平均演算法的運算式推導完成，

四、多步解包裹后平均法

從上文推導的多步平均演算法公式不難看出，多幀圖的光強資訊權重并不一樣，當每幀干涉圖均有強度接近的振動干擾時，上述演算法能起到一定的誤差平衡效果，一旦某幀影像受擾動影響過大，比如(15)式中第四幀圖，其高權重系數會導致解相程序受到成倍的振動干擾，

為了解決因權重系數差異大而帶來的解相不穩定問題，本文還提出了一種多步解包裹后平均法，設解包裹運算子號為unwrap，對一組幀數為9的干涉圖組，從第一幀圖開始每4幀圖進行一次四步運算得到，K的意義同多步平均法中一致，可得到六組包含不同移相量的相位分布，公式(16)：

$\begin{array}{c} \varphi(\mathrm{x}, \mathrm{y})^{(0)}=\varphi(\mathrm{x}, \mathrm{y})=\arctan \left(\frac{\mathrm{I}_{4}(\mathrm{x}, \mathrm{y})-\cdot \mathrm{I}_{2}(\mathrm{x}, \mathrm{y})}{\mathrm{I}_{1}(\mathrm{x}, \mathrm{y})-\cdot \mathrm{I}_{3}(\mathrm{x}, \mathrm{y})}\right) \\ \varphi(\mathrm{x}, \mathrm{y})^{(1)}=\varphi(\mathrm{x}, \mathrm{y})+\frac{\pi}{2}=\arctan \left(\frac{\mathrm{I}_{5}(\mathrm{x}, \mathrm{y})-\cdot \mathrm{I}_{3}(\mathrm{x}, \mathrm{y})}{\mathrm{I}_{2}(\mathrm{x}, \mathrm{y})-\cdot \mathrm{I}_{4}(\mathrm{x}, \mathrm{y})}\right) \\ \varphi(\mathrm{x}, \mathrm{y})^{(2)}=\varphi(\mathrm{x}, \mathrm{y})+\pi=\arctan \left(\frac{\mathrm{I}_{6}(\mathrm{x}, \mathrm{y})-\cdot \mathrm{I}_{4}(\mathrm{x}, \mathrm{y})}{\mathrm{I}_{3}(\mathrm{x}, \mathrm{y})-\cdot \mathrm{I}_{5}(\mathrm{x}, \mathrm{y})}\right) \\ \end{array}$

$\begin{array}{c} \varphi(\mathrm{x}, \mathrm{y})^{(3)}=\varphi(\mathrm{x}, \mathrm{y})+\frac{3 \pi}{2}=\arctan \left(\frac{\mathrm{I}_{7}(\mathrm{x}, \mathrm{y})-\cdot \mathrm{I}_{5}(\mathrm{x}, \mathrm{y})}{\mathrm{I}_{4}(\mathrm{x}, \mathrm{y})-\cdot \mathrm{I}_{6}(\mathrm{x}, \mathrm{y})}\right) \\ \varphi(\mathrm{x}, \mathrm{y})^{(4)}=\varphi(\mathrm{x}, \mathrm{y})+2 \pi=\arctan \left(\frac{\mathrm{I}_{8}(\mathrm{x}, \mathrm{y})-\cdot \mathrm{I}_{6}(\mathrm{x}, \mathrm{y})}{\mathrm{I}_{5}(\mathrm{x}, \mathrm{y})-\cdot \mathrm{I}_{7}(\mathrm{x}, \mathrm{y})}\right) \end{array}$

$\begin{array}{l} \varphi(\mathrm{x}, \mathrm{y})^{(5)}=\varphi(\mathrm{x}, \mathrm{y})+\frac{5 \pi}{2}=\arctan \left(\frac{\mathrm{I}_{9}(\mathrm{x}, \mathrm{y})-\mathrm{I}_{7}(\mathrm{x}, \mathrm{y})}{\mathrm{I}_{6}(\mathrm{x}, \mathrm{y})- \mathrm{I}_{8}(\mathrm{x}, \mathrm{y})}\right) \end{array}$

(16）

反三角函式arctan通過計算機程式計算得到的區間為單周期（-π，π），因此干涉圖的相位分布會出現斷層式變化，如圖所示是由四幀實際獲取的干涉圖解相后的相位分布情況，其中圖2為二維展示，圖3為三維展示，此時的相位分布雖然攜帶了面形資訊，但并非真實的周期性相位資訊，因此公式(16)的相位分布無法直接疊加，否則疊加后的相位分布將受到跨周期交錯的影響，而呈現例外的波浪式分布，則無法通過解包裹演算法還原真實的相位分布，

為此，將公式(16)的相位分布分別進行解包裹運算，使其還原成真實的相位資訊，即相位變化呈連續性而不是斷崖式，雖然這六組相位分布的初相資訊不一致，但這僅會導致相位的整體平移量變化，而平移量對面形測量而言是無關緊要的，對單組相位分布來說，其變化趨勢是相對的且攜帶面形資訊，將六組資料資訊直接疊加并取平均，即可獲得進行了一定程度誤差平衡的真實相位分布 $\bar{\varphi}(\mathrm{x}, \mathrm{y})^{\prime}$ ：