中国剩余定理相位解包裹方法的改进∗

赵甜甜 李 锋

（江苏科技大学电子信息学院 镇江 212003）

1 引言

光栅投影测量是三维测量方法的一个重要分支，因其具有非接触、低成本和高精度等优点，被广泛应用在光学三维测量领域［1～3］。由此方法得到的相位是立体匹配的一个特征［4］，结合系统标定参数，再根据三角法进行计算就可以得到被测物体三维信息。该方法的关键问题之一就是如何正确地获取待测物体的相位。目前最常用的测量相位的方法就是相移法，该方法具有较高的测量精度，但是相位被包裹在一个周期内，在整个测量空间呈锯齿状分布，必需将包裹相位展开。目前，相位解包裹技术主要有两种［5～6］，一是空间解包裹技术，二是多频解包裹技术。

空间相位解包裹只需要一幅包裹相位图，但是只有在理想情况下（被测物体轮廓简单，信噪比足够高）才能正确地获得绝对相位图。多频相位解包裹是假设在时间一致的条件下，通过向被测物体表面投影一系列不同的光栅条纹，然后对相机拍摄到的编码图案解包裹确定绝对相位值。其典型的方法有小数重合法、时间解包裹法、外差法和基于中国剩余定理（Chinese Remainder Theorem，CRT）的方法。因解包裹相位的实质为求解编码点的像素坐标对不同频率的包裹相位和频率构成的同余方程组的问题。而中国剩余定理根据系列模数和对应的剩余，采用简单的运算获得被除数，是求解同余方程组的有效方法，且可以实现最大范围解包裹，计算效率高。因而，基于中国剩余定理解包裹相位的方法吸引了国内与国外专家的普遍关注与研究。

Cushov［7］等初次将中国剩余定理应用于光学干涉仪相位轮廓测量，计算绝对相位存在偏差。Takeda［8］等 对Gerohbery-Saxton 算 法 做 出 了 改 进，提出了自适应高度偏移联合误差校正查找表修正较大偏差的方法。Zhong［9］等采用比值为无理数的两个频率在一定范畴内获得所有折叠整数的同余结果，并保留误差值，根据阈值进行选择符合条件的折叠整数作为最优结果。Towers［10］等假设在相位误差存在上、下限的条件下，遍历搜索绝对相位误差，根据测量范围约束确定最优结果。Xia［11］等提出了健壮中国剩余定理，但是该方法需要进行二维搜索，增大了计算成本。Wang［12］等在此基础上进行改进，将该定理从整数域扩充到实数域，可以进行封闭式求解。张旭［13］等通过分析亮度噪声对包裹相位的影响，提出了选频准则。文献［14］提出了基于CRT 选频准则的微相位测量轮廓术。于晓洋［15］等提出了一种CRT 工程化求解方法。但是存在抗干扰能力差，测量误差比较大的缺点。而且测量范围只适用于模互质和模存在最大公约数不为1 的情况，并没有对模非互质又不存在公约数的情况进行求解和分析。

因此，为了解决上述现有关于中国剩余定理解包裹存在的问题，本文提出了一种改进的中国剩余定理工程化求解方法。首先，本文采用标准四步相移的方法获取包裹相位，降低了包裹相位中存在的误差，提高了测量精度。其次，通过采用合并方程组的思想求解同余方程组，此方法既能计算模非互质又不存在公约数的情况，又能计算模互质和模存在最大公约数非互质的情况，扩大了中国剩余定理解包裹的应用范围，同时提高了测量速度，增强了抗干扰能力。

2 基于改进的CRT相位求解

本文采用标准四步相移作为获取包裹相位的方法，其波动范围为［0，2 π］，需要进行相位解包裹。从标准四步相移计算包裹相位的计算公式中，可以看出其具有消除背景光照和常数项影响的优点，对投影仪的非线性起到抑制作用，减小了系统误差，具有测量精度高、抗干扰能力强的优点。

2.1 中国剩余定理原理

设np、ri、pi为整数，i=1，2，…，k ；令ri为np对模数pi取余后的余数。当各pi之间不一定互质时，采用合并方程的思想进行求解同余方程组。设k1、k2、…、ki为整数，满足：

首先，合并公式np=r1+k1p1，np=r2+k2p2整理化简得：

显然，要想有解，必须gcd(p1，p2)|(r2-r1)。令gcd(p1，p2)=d ，r2-r1=C ，则有

依次类推迭代下去，最终求得k 个方程的最小正数解np为r%p。

实际求解绝对相位时，包裹相位和绝对坐标值都不是整数，而是实数。因此需要将CRT应用范围从整数域扩充到实数域。首先，对输入的模数pi即编码周期进行分析，设最大公约数为d ，当pi两两互质时，d=1；当pi存在最大公约数不为1 时，d= gcd (pi，pj)≠1(1 ≤i≠j≤k) ；排除上述两种情况，d=1。其次，设ri、np为实数，ri为余数即包裹相位，np为绝对坐标值。ri可以分解为riod+ric，其中rio为整数、0 ≤ric＜d。在理想情况下，即ri没有误差，所有的ri都是相等的，则riod也是相等的，ric相等，记为rac。有：

式中：rounddown()代表向下取整。最后，将rio代入式（1）替换余数ri，通过合并方程的方式求解出整数解np，与riod对应的整数解为npd。因为ri比riod大了rac，令rc=rac/d，0 ≤rc＜1。所以中国剩余定理的实数解为

通过式（1）～（9），将CRT 应用范围推广到对模数没有任何要求的情况，并且在实数域范围内对同余方程组进行求解。

2.2 小数差判据

实际工程中，设ri为包裹相位，其存在误差不相同，需要考虑CRT求解的准确性问题。

定义包裹相位ri的误差为|，如果，则rio会偏离原位置，只要不全部相等，必然会导致np求解失败。因此，，是求解正确的前提。设变换后的包裹相位误差和包裹相位分别为，即：r̂i的整数部分为bio，小数部分为bi，则

根据ei和rc的不同情况，r̂i=bio+bi有三种情况：1）当0 ≤ei+rc＜1 时，bio=rio，bc=ei+rc；2）当1 ≤ei+rc＜2 时，bio=rio+1 ，bi=ei+rc-1 ；3）当-1 ＜ei+rc＜0 时，bio=rio-1，bi=ei+rc+1。分析每种情况，bi都是(ei+rc)和整数的和，则包裹相位小数部分之差∆bij=bi-bj(1 ≤i≠j≤k)为包裹相位误差之差(ei-ej)与整数的和，所以∆bij可以描绘出包裹相位误差之间的差异。bi、bj与分别对应，且分别有三种情况，所以∆bij=bi-bj有九种组合：（1）和都属于1），∆bij=ei-ej；（2）和都属于2），∆bij=ei-ej；（3）和都属于3），∆bij=ei-ej；（4）属于1），属于（2），∆bij=ei-ej-1 ；（5）属 于1），属 于3），∆bij=ei-ej-1 ；（6）属 于2），属 于1），∆bij=ei-ej-1 ；（7）属 于 2），属 于 3），∆bij=ei-ej-2 ；（8）属 于3），属 于1），∆bij=ei-ej+1 ；（9）属 于3），属 于2），∆bij=ei-ej+2。

第一类情况I，包含组合（1）、（2）和（3），为

第二类情况II，包含组合（4）、（5）、（6）和（8），为

第三类情况III，包含组合（7）和（9），为

可以看出-2γ＜1-2γ＜2-2γ、2γ＜1+2γ＜2+2γ，通过限定γ使式（17）成立。

这样可以通过| ∆bij|来区 分I、II 和III，求得γ＜0.25 ，即 |ei|＜0.25 。又 因 为0 ≤bi＜1 ，所 以0 ≤| ∆bij|＜1，而2-2γ＞1，所以Z 类情况不存在。即：

根据式（10）可得，当编码周期存在最大公约数d时，包裹相位误差限定为。

2.3 取整方式

根据式（18），可以看出实数解的误差为(Ei+Ej)/2，小于误差 |Ei|和 |Ej|之间的最大值。同理，对其余两种情况进行分析也是满足上述结论的。

i 带入求得第i个模数pi的实数解==(np+1)d+ (rc+ei-1)d=npd+(rc+ei)d，同理，+=npd+(rc+ej)d。在根据式（18）获得最终实数解=()/2 =npd+rcd+(Ei+Ej)/2，误差仍然为(Ei+Ej)/2。同理，求解其余三种情况，得到如下结论：当0.5 ＜1，对包裹相位进行四舍五入取整，可以正确求取实数解，不存在求解失败的现象，其误差为(Ei+Ej)/2。

当存在k个模数时：首先判断所有小数部分的差值，如果所有的最大值小于 0.5 ，则 对 所 有r̂i进 行 向 下 取 整，即rounddown(r̂i)，将其作为余数测量值整数按照前面所述的过程进行计算，得到的实数解误差为如果所有的最大值大于0.5，则对所有进行四舍五入取整，即round(r̂i)，获得的实数解误差仍然为

总结上述分析，使用改进的CRT进行准确求解的条件为：当包裹相位测量误差＜d/4 时，即＜0.25，如果所有的最大值小于0.5，则所有的rio=rounddown(̂)，如果所有中最大值大于0.5，则所有的rio=round(r̂i)。将取整结果带入式（1），合并方程组进行求取整数解np，根据式（9）、（18）得到实数解，其误差为

2.4 改进CRT算法步骤

本文改进的CRT 工程化算法可以分为以下几步。

步骤1：计算所有变换后包裹相位̂ 之间的小数部分之差。步骤2：当所有中的最大值小于0.5 时，将所有包裹相位̂向下取整；当所有中的最大值大于0.5 时，对所有的包裹相位̂ 采用四舍五入的方式进行取整。

步骤3：将包裹相位的取整结果带入式（1），采用合并方程组的思想，求取绝对坐标整数解。

步骤4：令rc为包裹相位与其取整结果之差，根据式（10）求解获得第i个模数pi的实数解。

步骤5：利用式（18）解得同余方程组的实数解，即编码图案绝对坐标值。

本文CRT 求解算法既可以求解模数非互质的情况，也可以求解互质的情况，具有解析求解、运算简单、求解快的优点。

3 仿真与实验分析

3.1 仿真分析

为了验证本文方法的可行性，选取三个周期T1=13，T2=15，T3=18，最大展开范围为1170pix⁃el。在Matlab2017a 中采用标准四步余弦相移的方法产生条纹图形，大小为1024pixel×768pixel，并在上面分别添加高斯分布随机噪声，计算得到三个不同周期下的包裹相位a1，a2，a3。根据包裹相位ai和周期Ti，采用本文方法获取编码图案绝对坐标。

表1 不同方差高斯随机噪声下的相位解包裹结果

表1 为对条纹图形归一化后，分别加入不同方差的高斯随机噪声误差，与未加入噪声误差的条纹图形求取编码图案绝对坐标进行比较。为进一步证明实验数据的有效性，本文对无噪声的条纹图形进行30 次求解取其平均值带入计算，同样对加入不同方差的随机噪声得到的条纹图形进行30 次求解取平均值。从表1 中可以看出，在随机噪声方差不大于0.25 时，最大相位误差不大于12.75%，平均相对误差最大为0.08%，从而看出使用本文方法求解出绝对坐标值的精确度较高，抗干扰能力强，验证了本文方法在求解周期即不互质也不存在最大公约数不为1 情况下的可行性。但当随机噪声方差大于0.25时，求解平均相对误差和最大相对误差急剧增大，求解得到的绝对坐标值精确度大幅降低，出现绝对坐标值的跳变。当随机高斯噪声方差为0.25 时，对模拟包裹相位进行求解，其结果如图1（a）为求解得到的绝对坐标面型图，（b）为选取第600 行的绝对坐标图。从图1 中可以看出，存在较大随机噪声时，使用本文方法依然能精确地对包裹相位进行展开。从理论上验证了本文方法在求解非互质周期且最大公约数不为1的可行性。

图1 相位展开结果

3.2 实验分析

为了进一步验证本文方法的有效性，构建了一套由EPSONEB-C760X 投影仪（分辨率为1024pixel×76 8pixel）、单个AVT相机Manta G-125（分辨率为1292pixel×964pixel）和一台个人电脑组成的结构光测量系统（Windows7，Matlab2017a），如图2所示。

图2 结构光测量系统

采用搭建的实验平台，向标准白板投射周期为13、15、18 的光栅图像，采用本文算法得到绝对相位，从而对该平面进行三维重建，重构出来的平面表面光滑，点云分布均匀。对其进行平面拟合，得到结果的最大绝对偏差为1.96mm，平均绝对偏差为0.546mm。

图3 鼠标及其相位展开结果

对图3（a）中的鼠标进行测量。本文算法与多频外差法、基于CRT的三步梯形相移法［16］进行对比实验，在相同环境下，直接对获取的产生畸变的条纹图案进行求解包裹相位。图3（b）、（c）、（d）分别为采用三种方法获得的绝对相位图。图3 表明，多频外差法和基于CRT 的三步梯形相移法抗噪声能力差，解析得到的绝对相位图上出现明显的噪点，轮廓不清晰，且对系统的非线性抑制能力弱。而本文方法求解出来的绝对相位，表面光滑，物体轮廓清晰，除去由物体阴影部分产生的误差，表面没有噪点。本文方法明显优于图3 中其他方法，具有抗干扰能力强，求解精度高的优点。图4 为使用本文方法得到的三维测量结果。

图4 鼠标测量结果

表2 给出了本文方法与多频外差法、基于CRT的三步梯形相移法相位解包裹所用的平均时间，可以看出本文方法在测量时间上也是优于其他方法的。

表2 不同方法的相位解包裹时间

4 结语

现有的基于CRT 相位解包裹方法对输入的包裹相位误差非常敏感，且不能对非互质频率的包裹相位进行求解。由于对包裹相位误差敏感，本文采用标准四步相移法获取包裹相位，减小了系统误差，增加了抗干扰能力。通过包裹相位测量值的小数差选择取整方式，并在计算整数解时采用合并方程组的思想。本方法无需限定输入频率，从而扩展了中国剩余定理解包裹相位应用范围。仿真和实验以及对比分析结果表明，本文方法运算简单，求解速度快，同时抗干扰能力强，求解精确度高。但是该方法限定了包裹相位测量误差，对包裹相位误差大时要采用其他措施