數學形態學區域分割的快速相位解包裹算法

李夢霞，曹 博，盧佳瑋，崔凱華，劉 乾*

（1. 中國工程物理研究院機械制造工藝研究所，四川綿陽621000；2. 清華大學 精密儀器系，北京 100084）

1 引 言

光學干涉以無損、高精度、高分辨率等優勢而廣泛應用于高精密表面測量。從干涉條紋中計算的相位被包裹在(-π，π]區間，需要對包裹相位進行解包裹運算，以獲得連續的表面高度分布。因此，解包裹是干涉測量中必不可少的關鍵環節，直接決定了干涉測量的精度與效率。

解包裹算法分為時域和空域兩種方式。時域方法獲得不同波長對應的相位，根據相位與波長成反比的特點在時域上對相位進行解包裹運算；空域方法根據相鄰像素之間的相位差對相位的級次進行判斷，進而完成解包裹運算。時域解包裹僅與當前點相位隨波長的變化有關［1-2］，可以避免誤差沿空間傳遞，但至少需要兩種條紋頻率的光源或投影設備［3-4］，而空域解包裹只需要一幅相位圖即可求解。一維情況的空域解包裹相對簡單和成熟，比較有代表性的是Itoh 方法［5］。Itoh 方法逐個計算相鄰像素的相位差，對超出 (-π，π]的相位差加減2π，累加相位差后獲得連續相位。依次在行與列方向上應用Itoh 方法，可以簡單地實現二維相位的解包裹。但此方式極易受到噪聲干擾，誤差會隨著展開的路徑而累積，導致拉絲等錯誤的結果［6］。

為了解決拉絲的問題，研究人員發展出了多種改進的解包裹算法。這些算法可以分為路徑相關和路徑無關兩大類。路徑相關解包裹算法的基本思想是沿著一條選定的連續路徑展開相位，其關鍵在于選擇合適的路徑。如果路徑選擇錯誤，測量噪聲會隨著路徑累加造成較大的誤差。路徑相關算法主要有枝切法［6-8］、質量圖導向法［9］等。枝切法利用連接距離最近的正負殘點的枝使解包裹路徑避開殘點區域，計算速度快，但容易在殘點密集區域出現區域性展開誤差。質量圖導向法按照各點相位質量高低指導解包裹路徑，可以避免殘點的影響，但是耗時較長。路徑無關算法以最小二乘法及其衍生算法為主［10-11］。最小二乘法的基本思想是通過最小化的誤差度量濾除局部噪聲，能夠在一定程度上平滑噪聲，但無法限制不可靠數據點引起的相位誤差在空間中傳播，而且處理時間較長。2018 年Zhang Xiaolei 等人［12］提出了一種基于非二次采樣輪廓變換（NSCT）的路徑無關相位解包裹算法。NSCT 算法［12-14］以非下采樣金字塔濾波器和非下采樣方向濾波器對相位圖進行濾波和分解，計算不同區域的相位差異后抬升相位，獲得連續的相位。NSCT 算法使局部誤差無法影響整體的解包裹效果，有效降低噪聲的影響，但由于邊界提取涉及構造濾波器、多次濾波和分解等復雜運算，耗時過長。

針對干涉測量面臨的噪聲干擾、相位分布復雜、數據點缺失等實際問題，以及工程實際中對運算效率的要求，本文提出了一種路徑無關的快速相位解包裹算法。本方法運用數學形態學工具對包裹相位進行區域分割繼而整體展開，具有一定的魯棒性。本文首先介紹該算法的原理，然后通過仿真和實驗驗證了該方法的精度、適用性、計算效率等。

2 原理與步驟

二維包裹相位圖的特點是：不同級次的相位之間有明顯的包裹分界線，在不同級次（區域）內相位是連續變化的。基于此特點，只需確定相鄰區域之間相位的相對抬升值，在各區域相位的基礎上加上對應的抬升值就能實現整張相位圖的解包裹。這是本文相位解包裹算法的基本思想。

數學形態學區域分割（RSMM）解包裹方法的核心在于包裹相位的邊界識別和區域分割。數學形態學是建立在格論和拓撲學基礎上的圖像分析方法，在結構光投影的三維深度分割［15］中有所應用。數學形態學方法的相關算法已經比較成熟，將之與相位解包裹算法結合起來，可以節省運算時間。RSMM 算法的流程如圖1 所示（其中加粗的節點為相位數據，未加粗部分為實際操作步驟），各步驟詳細描述如下：

圖1 RSMM 解包裹算法原理流程圖Fig.1 Flowchart of the RSMM unwrapping algorithm

2.1 無效點填充

在解包裹之前，需要對相位數據缺失的無效點進行填充處理。填充值不能偏離區域極值（±π），且填充值的選擇應與其附近的數據有一定的關聯，毫無關聯的數據會導致錯誤的結果。另外，填充值用于輔助解包裹，并非完美還原實際數據，因此不用過分追求準確性。根據上述原則，本文采用較為簡單、計算量小的移動中位數法［16］進行線性填充，直至所有的無效點都被賦予有效數值，得到填充后的包裹相位。

2.2 降噪處理

包裹相位會存在一些因系統、算法等因素導致的噪聲，可能會引起邊界的誤提取。為了減小噪聲導致的邊界提取誤差，本文采用自適應中值濾波［17］對包裹相位進行去噪處理。

2.3 區域分割

包裹相位在截斷處的相鄰像素會出現大于π的跳變，表現為梯度的陡變。因此，可以利用包裹相位的梯度識別跳變點以確定相位截斷的邊界。首先，用Prewitt 梯度算子［18］計算包裹相位的梯度，同時設置合理的閾值（梯度最大值的一半）得到鄰位相位差大于π 的點，這些點的集合構成包裹相位的邊界。其次，為避免邊界斷裂，使用膨脹函數將同一條邊界上的斷裂點連接起來形成完整的邊界，進而生成全局的二值圖（邊界點為邏輯1，區域點為邏輯0）。然后，將二值圖中的不同邊界分別標記為不同的灰度值。比如，將檢測到的第一個邊界上所有的點都賦值1，第二個邊界上所有的點都賦值2，以此類推，直到所有的邊界點都被賦值。此時灰度圖像中最大值的大小即代表邊界的個數，不同灰度大小的點集即為不同的邊界（記為B）。區域的分割與邊界分割類似。首先，將全局二值圖取反，邊界點變為邏輯0，區域點變為邏輯1。其次，將區域二值圖中的不同區域分別標記為不同的灰度值，不同灰度大小的點集分別代表不同的區域（記為R）。相位圖的邊界類似于等高線，本質上是相同高度點的集合。因此，邊界線不會出現交叉，區域數量最多比邊界數量多1 個。實際操作中可根據這一特點將誤提取的區域刪去，融入鄰近的邊界或區域。然后，對邊界進行小系數膨脹，膨脹后的邊界（記為DB）作為后續操作的判別工具。

圖2 為peaks 函數生成包裹相位圖中劃分的邊界和區域示意圖，圖中B表示邊界，R表示區域。

圖2 相位圖分割示意圖（白色表示選中的區域）Fig.2 Diagram of phase segmentation（white represents the selected area）

2.4 抬升值計算

2.4.1 區域抬升值

為了計算不同區域的相對抬升值，需要明確區域之間的位置關系，因為只有相鄰區域間才有明確的相對相位抬升值。借助膨脹的邊界與原始的區域來判斷不同區域是否相鄰，然后計算相鄰區域的相對相位抬升值。

結合圖3 描述的具體做法為：以圖2 中的邊界B1為例，首先，將其膨脹后的邊界DB1分別和所有區域（Rj，j=1、2……）取交集，若DB1只有和R1、R2的 交 集 不 為 空 集 ，則R1和R2是 被B1分 割開的相鄰區域。以圖3 為例，R1和DB1的交集表示在R1內部且靠近B1的點的集合，R2和DB1的交集表示在R2內部且靠近B1的點的集合。將邊界和區域的交界點的集合記為上標i表示區域標號，下標j表示邊界標號。在包裹相位中計算點集的相位中值，以表示。與同一邊界Bj相鄰的兩個區域Rm和Rn，其各自交界點集的相位中值為區域Rm相對Rn的抬升值可以計算為：

圖3 區域相鄰關系判斷示意圖Fig.3 Diagram of determining the relation between adja?cent areas

通過上述方法可以確定所有相鄰區域的相對抬升值，最后以某個區域為標準（例如設R1的絕對抬升值Er1為0），根據各區域的位置關系和相對抬升關系可確定所有區域的絕對抬升值Erj。

2.4.2 邊界點抬升值

在上一步中只計算了各區域的抬升值，而邊界點比較復雜，需要逐點判斷。邊界Bj上的點bj(x，y)的包裹相位值為φ(x，y)，利用它與相鄰邊界相位中值之間的相對大小判斷其歸屬，進而確定其抬升值，計算如下：

通過上述步驟，可以確定不同區域和邊界上的點的絕對抬升值，記錄這些抬升值并建立全局抬升值矩陣。

2.5 相位抬升

為避免降噪處理等步驟改變原始包裹相位的值，對解包裹結果產生干擾，因此在無效點填充后的原始包裹相位的基礎上進行相位抬升。填充后的包裹相位經（2）（3）（4）步驟獲得抬升值矩陣，與其本身相減即可得到抬升后的相位。

2.6 抬升誤差修正

由于依據式（3）計算的相位抬升值在少量邊界點處存在偏差，進而引起邊界抬升誤差。抬升誤差點與周圍正確點的抬升值相差2π 的整數倍，導致抬升后誤差點與周圍點相位差絕對值大于π。理論上抬升后相位的相鄰像素不應有超過π的躍變。根據這一原則，確定誤差點的位置，做局部中值濾波以修正抬升誤差。如果包裹相位存在無效點，則在解包裹后的相位中將填充的數據還原成無效值，得到最終的解包裹相位。

RSMM 算法優勢在于：（1）使用成熟的數學形態學算子對包裹相位進行邊界分割和抬升值判斷，計算量小，準確度高；（2）包裹相位的無效點不影響區域的判斷，能夠實現含缺失數據的相位解包裹；（3）抬升值在區域上整體計算，而只在邊界逐點計算，省去了大量的相位比對時間，且局部區域噪聲不會影響全局的解包裹效果。

3 計算機仿真

通過計算機仿真對RSMM 算法的計算精度、耗時、適用性進行了驗證。仿真所用計算機CPU 主頻為 3.00 GHz，內存頻率為 2 667 MHz。

3.1 計算精度

以peaks 函數作為基礎相位，生成含有密度為0.1、幅值為π/5 rad 的椒鹽噪聲的真實相位，并將其壓縮在 (-π，π]，得到 500×500 pixels 的包裹相位，如圖4 所示。可以看出，包裹相位中存在明顯的噪聲分布。先對包裹相位圖進行自適應中值濾波，然后利用數學形態學算子提取邊界（如圖5 所示），計算解包裹相位與真實相位的點對點殘差，如圖6 所示。殘差的RMS 值為4.4×10-16rad，表明RSMM 方法具有足夠高的精度，能夠獲得精確的解包裹相位。又經過仿真驗證了不同幅度噪聲對結果的影響，所加相位噪聲幅度在0~π/2 范圍內，解包裹殘差RMS 值均在10-16rad 量級。這是因為區域分割過程中采用了降噪處理，噪聲基本不影響邊界處理結果。為比較解包裹的效果，同時采用了NSCT 算法和最小二乘算法對相位進行解包裹。三種算法的殘差RMS 都在10-16rad 量級，表明三種方法均能正確解包裹，且精度相當。然而，本文方法所需要的計算時間僅為 0.3 s，遠小于 NSCT 算法（13.3 s）和最小二乘算法（2.9 s）。

圖4 含噪聲的peaks 包裹相位Fig.4 Noised wrapped phase of peaks

圖5 提取的邊界Fig.5 Extracted boundaries

圖6 RSMM 算法解包裹相位的殘差Fig.6 Residual error map of unwrapped phase with RSMM algorithm

3.2 計算時間

由于RSMM 算法需要對包裹相位進行分區域計算，因此計算時間會受條紋數量的影響，條紋越多、區域越多，則耗時越長。保持數據圖大小為 500×500 pixels 不變，調整 peaks 函數的整體高度，從而改變包裹相位圖的條紋數。記錄不同包裹條紋數相位解包裹所用時間，并對三種算法進行比較，如圖7 所示。圖中可以看出，RSMM 算法用時最短，僅為最小二乘法的10%、NSCT 法的1.7%。由于NSCT 算法將不同頻率子代下不同方向的邊界融合，又采用區域生長法進行邊界膨脹。提取的邊界比實際大，在條紋密集區域出現相位抬升錯誤。當條紋超過一定密度（6 根/100 像素）以后，NSCT 解包裹算法不能準確地求解，因此未記錄N>12 以后的計算時間。

圖7 三種算法耗時隨條紋數量的變化Fig.7 Calculation times of three algorithms with respect to the fringe number

除條紋數量外，圖像大小也會影響解包裹算法的效率。固定條紋數量為8，對包裹相位的像素數量進行調整，測試從200×200 pixels到 2 000×2 000 pixels 不同像素數量的解包裹運算時間。如圖8 所示，三種算法的計算時間均隨像素數增加而增加，其中RSMM 算法計算時間最短，且隨像素數量增加最慢。這是因為RSMM算法只計算邊界處像素的相位，計算量隨整體像素數量增長較緩慢。

圖8 三種算法耗時隨相位圖大小的變化Fig.8 Calculation times of three algorithms with respect to the phase map size

3.3 復雜相位的適用性

在實際測量時得到的相位可能會因為表面區域分布、污染等情況導致相位數據缺失，存在較復雜的無效點分布。為測試算法的適用性，在生成peaks 相位時加入了復雜分布的無效區域（如圖9）。當包裹相位圖中存在無效點時，RSMM 算法能夠準確解包裹，如圖10（a）所示。而NSCT 解包裹算法出現邊界識別錯誤且后續無法計算，最小二乘算法也出現解包裹錯誤，如圖10（b）所示。

圖9 含無效點的包裹相位Fig.9 The wrapping phase with data dropout

圖10 含無效點的包裹相位解包裹效果Fig.10 Unwrapped phase with data dropout

計算機仿真表明，RSMM 算法具有計算精度高、抗噪能力強、運算耗時短、適應性強等優點。這些優點源于：區域分割等數學形態學工具的應用提升了邊界提取的準確性和效率；相位圖進行分區域整體抬升提高了解包裹精度和效率；解包裹前使用移動中位數法線性填充無效點和自適應中值濾波處理噪聲，增強了算法的魯棒性和適應性。

4 實驗驗證

為測試RSMM 算法解包裹的實際效果，采用干涉顯微鏡分別測量超精密車削工件、標準臺階和表面有瑕疵的平面反射鏡。

車削表面的包裹相位如圖11 所示。采用三種方法對車削表面相位解包裹，結果如圖12 所示。RSMM 法和最小二乘法解包裹得到了比較一致的準確結果，而NSCT 解包裹結果中部分區域出現抬升錯誤。標準臺階標稱高度為80 nm，其包裹相位如圖13 所示。雖然臺階表面自身存在梯度劇變，但在邊界提取時設置了閾值，臺階邊緣并未影響到解包裹效果。RSMM 算法對標準臺階相位解包裹并校平后的結果如圖14所示。在商業設備中計算得到臺階高度為85.639 nm，RSMM 算法解得高度為85.686 nm，最小二乘法解得高度為85.698 nm，NSCT 算法解得高度為85.769 nm，三種算法精度相當。對于這兩種測量對象的包裹相位，三種算法解包裹的運算時間如表1。可以看出，RSMM 算法的計算時間僅為最小二乘方法的1/4、NSCT 算法的1/5。

表1 三種解包裹算法用時Tab.1 Calculation times of three unwrapping algorithms

圖11 車削表面包裹相位Fig.11 Wrapped phase of the turned surface

圖12 車削表面解包裹相位Fig.12 Phase unwrapped of the turned surface

圖13 標準臺階包裹相位Fig.13 Wrapped phase of the standard step

圖14 RSMM 解包裹標準臺階相位Fig.14 Phase unwrapped of the standard step with RSMM algorithm

圖15 為表面有瑕疵的平面反射鏡的包裹相位，其邊界和區域內部均有不規則的數據缺失，邊緣毛糙，噪聲較大。RSMM 法和最小二乘法解包裹的相位如圖16 所示，而NSCT 解包裹算法無法處理含無效點的數據故未給出結果。從圖中可以看出，雖然包裹相位含有較多不規則分布的無效數據，但RSMM 算法仍然能夠正確解包裹且達到較好的效果，而最小二乘法則出現了錯誤結果。鏡面的表面粗糙度經商業儀器測得為9.4 nm，RSMM 算法求得的粗糙度為10.1 nm ，誤差可能源于數據處理過程中的濾波、閾值等處理環節的差異。

圖15 平面反射鏡的包裹相位Fig.15 Wrapped phase of the mirror

圖16 平面反射鏡的解包裹相位Fig.16 Phase unwrapped of the mirror

5 結 論

本文提出的RSMM 算法能夠快速解算包裹相位，精度與傳統算法相當，并且能夠有效處理隨機噪聲干擾等實際問題。該算法適應性強，對于條紋密度高的包裹相位、部分區域數據缺失的復雜包裹相位，都能實現準確的全局解包裹。仿真和實驗結果證明：RSMM 算法殘差RMS 值可以達到 10-16rad 量級，對 1 000×1 000 pixels 的包裹相位解算時間在1 s 以內，僅為經典的最小二乘算法用時的1/4。本文只展示了RSMM 算法應用于干涉測量的效果，這種算法也可以擴展到其他需要相位解包裹的領域，如數字全息、光柵條紋投影輪廓測量等。