基于FPGA的线阵干涉条纹相位解包裹算法实现

2019-10-18郗洪柱王伟魁乐贵高孔德仁

测控技术 2019年9期

郗洪柱，王伟魁，乐贵高，孔德仁

(1.南京理工大学机械工程学院，江苏南京 210094；2.北京遥测技术研究所，北京 100076)

在非接触式光学测量中，一般通过光干涉原理得到干涉图对被测量进行分析。然而，由于正切函数的周期性，从干涉图提取的相位信息被包裹在(-π，+π］范围内［1］。相位解包裹技术应运而生，在合成孔径雷达干涉测量、数字全息测量、条纹投影测量和数字散斑干涉测量等领域有重要应用［2-4］。

在干涉式光声光谱痕量气体检测中采用互补金属氧化物场效应管(CMOS)线阵图像传感器，以此获得光的干涉条纹。然而，由于干涉仪安装和光路调试误差导致所得干涉条纹图具有相当的噪声，并由于对干涉条纹图采用空间相移的Hariharan五步法［5］，所得相位信息存在断点，造成相位图中存在一系列局部坏点，即残差点，阻碍了真实相位的获取。

现有的相位解包裹算法，主要分成空间和时间两类相位解包裹算法［6-7］，如改进的行列去包裹算法［8］及最小二乘算法［9］等一般应用在二维相位解包裹场合中。对于一维相位解包裹算法，目前仍以K.Itoh首先提出的求包裹、求差分和求积分为主［10-11］。在嵌入式系统中应用传统一维相位解包裹算法时，存在包裹系数不易确定及资源占用较大等问题。

为实现线阵干涉条纹相位快速解包裹，优化资源占用，通过设计激光器驱动频率，以及对传统一维相位解包裹算法进行简化，设计了可变包裹系数的逐点差值比较和累加算法，利用FPGA技术实现了专用相位解包裹模块，消除了相位残差点；为了消除解包裹后相位图中的噪声，利用Matlab中的FDATool工具，生成了二阶椭圆滤波器系数，由于系数是浮点数，在FPGA中进行相应的移位截断处理，并采用了防累积误差溢出的措施，实现了稳定、收敛的IIR滤波器。

1 干涉条纹相位解包裹原理及简化

从干涉条纹图中得到相位信息需要进行反正切运算，而在FPGA的反正切IP核计算结果范围为［-π，+π］，为了得到真实相位，需要进一步处理。针对二维相位解包裹，真实相位与缠绕相位可由关系式(1)确立［12］：

式中，φ(i，j)为波面(i，j)处的真实相位；φ(i，j)为反正切计算得到的包裹相位；κ(i，j)表示与波面位置(i，j)有关的整数，满足式(2)确立的关系式：

式中，unwrapping()是解包裹函数，函数内接收两个参数，其一为包裹相位φ(i，j)，另一个参数为φ(N8ij)，即(i，j)点周围有效的8点邻域的相位计算结果。用数学关系式表示(i，j)点周围的8点邻域为

用N8ij表示(i，j)点周围有效的8点邻域，即

一维相位解包裹时，需对上述式(1)～式(4)中的波面(i，j)改为像素点(i)，且8点邻域N8ij修正为两点邻域N2i，即2点邻域为

若用wrapping()表示包裹算子，则φ(i)=wrapping(φ(i))，式(1)可表示为

移项后式(6)也可表示为

式中，n(i)=-κ(i)。利用符号Δ表示差分，则对式(7)两边施加差分运算，并进行包裹运算后有

当-π≤Δφ(i)≤π时，Δn(i)=0，即不需要解包裹。此时有

所以当真实相位差值在［-π，+π］范围内时，其值同时与包裹相位差值相等。实际计算时，采用对相邻包裹相位值计算差值，对差值超过［-π，+π］范围时进行加或减去2π整数倍的方式控制在［-π，+π］内，对初始相位不敏感的场合，将修正后的包裹相位差值累加即得到真实相位的变化［13］。关键之处在于确定整数n(i)的具体数值。

根据相位包裹原理，当前后相位差值超过［-π，+π］范围时，为了将相位差拉回到［-π，+π］内，需要根据具体差值确定整数n(i)。通过判断前后相位差值的绝对值与2π的大小(倍数)关系，可确定整数n(i)。这种处理手段可以起到平滑相位结果的作用，同时，整数n(i)被用作解算结果的修正因子，即当环境中某些对特定波长敏感的气体浓度突然增加时，对应的前后相位差值会显著加大，根据整数n(i)的大小，相应调低激光器光功率，降低待测敏感气体的光热效应，从而可将相位差拉回到［-π，+π］内，根据当前激光器光功率修正待测气体浓度即可。

2 逐点差值比较和累加——滤波法设计

2.1 算法的硬件架构

对干涉式光声光谱痕量气体检测所需的FPGA资源情况作总体评估，选择了Xilinx公司Virtex-II系列300万门芯片XC2V3000，该芯片内部RAM存储器、乘法器及IP核资源丰富，有利于大量数据的实时运算处理。利用AD9240采集干涉条纹，该款数模转换器具有10MSPS转换速率，能够较快地采集CMOS线阵图像传感器输出的512个像素值。

在FPGA中运行解包裹算法，涉及的运算包括求差、求和、比较大小、求余数、数据格式转换、求反正切、滤波(乘加)和四舍五入等，其框架如图1所示，图中a1、a2、a3、b1、b2 和 b3 为二阶 IIR 滤波器系数，pi为圆周率π。为了降低算法运行时的风险，利用FPGA器件成熟IP核实现求余数和求反正切的功能。反正切IP核要求被操作对象以x/y形式输入，同时要求x和y的范围均需满足［-1，1］。基于反正切计算时的上述要求，需要增加额外硬件资源对前端数据进行处理，使之符合反正切的数据输入要求。比较大小和求余数模块正是为了满足这一要求而设计的。在数据参与求余数前，通过大小比较，保证小数除以大数，此时除法IP核输出的整数部分一定为0，余数部分有效。根据前面小数和大数关系确定余数应作为反正切里的x还是y即可，两者确定了一个，另一个一定是1。在具体设计算法时，首先需要规定好数据格式。

2.2 数据范围的确定

干涉条纹模拟电平经AD9240转换后变为14位数字量。扩展为16位后送入除法IP核中，产生的余数仍为16位。按照反正切IP核的要求，将余数扩展为19位，输出为19位，其中，输入数据格式为1N19型，即最高位为符号位，次高位为整数位，余下17位为小数位，最低位权值为2-17。输出数据格式为2N19型，即最高位为符号位，次高位及下一次高位为整数位，余下16位为小数位，最低位权值为2-16。

图1 解包裹框架图

为累加和预留寄存器空间，对相位角整数位进行扩展(扩展至15位)，小数位保持不变(仍为16位)，最高1位仍为符号位。对π的数值按照前述定点数据格式表示法，表示为15N32，即总位数为32位，低16位固定为小数位，最高1位为符号位，中间15位为整数位。

滤波器系数按照16位设计，最高位为符号位，滤波器输出结果为48位，最高位为符号位。

相位解包裹后的结果设计为32位，最高位为符号位。

2.3 算法组成

解包裹算法开始于反正切得到包裹相位值，执行部分包括四块，分别为：当前相位纠正、相位差累加、二阶IIR滤波、结果的截断及四舍五入处理。数据流关系如图2所示，其中，当前相位纠正模块中的初始相位为0，扫描完毕复位指的是按某一固定波长扫描预定时间，得到足够组待平均原始数据后对算法整体进行复位操作，然后切换至下一个波长，重复算法步骤。

图2 逐点差值比较和累加——滤波法数据流图

2.3.1 当前相位纠正

初始相位按0计，当反正切完毕输出相位结果时，将当前相位结果与上次相位结果做差，以差值范围为［+π，+2π］及［-π，-2π］为例，大于 π(按照15N32数据结构表示为：“00000000000000110010010000 111111”) 时，将当前相位加上 - π，即“11111111111111001101101111000001”；小于 -π(按照 15N32 数据结构表示为：“1111111111111001101101111000001”)时，将当前相位加上 π，即“0000000000000110010010000111111”；否则，当前相位不做处理。纠正模块流程示意如图3所示。图中pi为圆周率π。

图3 当前相位纠正模块流程示意图

2.3.2 相位差累加

初始相位累加和按0计，当当前相位纠正模块输出结果后，将更新后的当前相位值与上次值做差后进行累加求和。当接收到扫描完毕复位信号时对寄存器清零。

2.3.3 二阶IIR滤波

IIR滤波器相对于FIR滤波器在相同阶数下具有更好的滤波效果［14］。为了减小资源占用，采用阶数尽量小的IIR滤波器。IIR滤波器通常按照设计原理可分为巴特沃斯、椭圆、切比雪夫1和2型等，按照椭圆滤波器设计出的滤波器阶数最小［15］，且在相同阶数条件下，有最小的通带和阻带波动。利用式xb(n)=b1·x(n)+b2·x(n-1)+b3·x(n-2)简化二阶 IIR 滤波器实现公式，如式(10)所示：

利用Matlab中的FDATool工具，按照椭圆滤波器法生成二阶IIR滤波器系数，采样频率(FS)为6993 Hz，通带截止频率(Fpass)设置为72 Hz。一般通带纹波(Apass)越小越好，阻带衰减(Astop)越大越好，但在实际设置这两项参数时发现，如果通带纹波过小，阻带衰减过大，Matlab生成的滤波器阶数将远超二阶。考虑到实际硬件资源限制，当通带纹波设置为0.01 dB，阻带衰减设置为100 dB时，生成的是二阶滤波器，此时的通带纹波和阻带衰减较为理想。参数设置完毕后单击“Design Filter”，单击菜单栏“Analysis”中的“Filter coefficients”，显示生成的系数及增益等，然后单击菜单栏“Edit”中的“Convert to Single Section”，转换为与式(10)对应的各项系数。最终的系数及在FPGA内的表示方式如表1所示。

表1 滤波器系数

表1中第3列的转换方法为：首先对生成的浮点数值放大16384倍，按照四舍五入原则截取整数部分，并表示为16位有符号数。

在FPGA中利用移位与乘加模块，以及表1的滤波器系数，实现公式(10)所述的二阶IIR滤波器。然而由于滤波器系数经过放大后进行了有限位截取，经公式(10)运算后会产生误差累积。为避免误差累积太大，导致系统不稳定，采取了防累积误差溢出的措施，即在每个波长扫描完毕后，均将相关寄存器清空，保证了每个波长扫描下滤波器的稳定运行。

2.3.4 截断及四舍五入

前一步IIR滤波器系数扩大了16384倍，需对滤波后的数据做相应的缩小，为了尽量减小FPGA中数据截取时造成的误差，设计四舍五入模块进行补偿。算法流程示意如图4所示。得到滤波结果后，将结果的最高位截取作为新数据的符号位，将结果的第44位至第14位作为新数据的数据位。如果滤波结果的第13位为‘1’，表明新数据的小数位绝对值大于或等于0.5，需要进位，将新数据做加一操作，否则，新数据保持不变即可。

图4 截断及四舍五入模块流程图

3 解包裹性能测试

利用高速A/D采集CMOS线阵图像传感器S9227的输出，得到干涉条纹，横坐标为像素个数，按照从S9227采集到的电平顺序排序，纵坐标为对应像素处的A/D量化数值，用黑色实心圆点标注各像素处对应数值，如图5所示。

图5 线阵干涉条纹图

对干涉条纹按照Hariharan五步法，利用FPGA的反正切IP核求得相位值。按照当前相位纠正小节所述，对相位值进行纠正，然后按照相位差累加小节所述，对相位差进行累加操作。连续敲击干涉仪外壳，使悬臂梁受迫振动，经相位解包裹后得到如图6所示振动情况下的相位解包裹结果。横坐标表示相位解算点数，按时间顺序从左往右排列，纵坐标为对应的解算数值用黑色实心圆点标注。图6振动曲线平稳，符合预期认知，表明算法在外界干扰下仍能稳定可靠的运行。进一步分析可知，图中振动周期约为45个点，每个点表示的时间约为143.5 μs，因此，悬臂梁的受迫振动频率约为150 Hz。

当干涉仪中的量子级联激光器固定扫描1038 cm-1时，相位解包裹结果如图7所示。横坐标表示随时间变化得到的解算结果个数，纵坐标表示对应解算值，用黑色实心圆点标注。由图可知，相位解包裹后整体变化过程连续，但存在高频噪声，如同数字信号处理中常见的毛刺一样，对后续处理有一定的影响。按照2.3.3节二阶IIR滤波所述，设计并在FPGA中实现了稳定、收敛的二阶IIR滤波器。对图7的相位解包裹结果进行滤波处理，得到图8所示结果。图中横纵坐标定义与图7相同，且仍用黑色实心圆点标注位置坐标对应的滤波结果。图8所示滤波结果平滑连续，表明滤波效果良好。

图6 敲击时的相位解包裹结果

图7 相位解包裹结果

图8 滤波结果

从Hariharan五步法处理完毕开始算起，到产生解包裹结果为止，算法所耗用资源情况如表2所示。触发器用了4616个，占总资源的16%，四输入查找表(LUT)用了3606个，占总资源的12%，乘法器用了4个，占总资源的4%。此外还用了一个数字时钟管理器(DCM)用于系统时钟的生成和管理，用了一个除法IP核，一个反正切IP核。

表2 算法资源开销情况

本文所述相位解包裹算法从得到相移结果(ad1_frames_data_ok信号拉高)开始算起，经过144个时钟周期计算完毕，得到第一个解包裹结果(phase_add_flag_out信号拉高)。本系统时钟100 MHz，算法全程耗时为1440 ns(开始时间10300 ns，完成时间11740 ns)，如图9所示。由于Goldstein法解1200*1200像素图像约需100 s，TIE 法仍约需3 s［16］，本文应用于嵌入式系统的解包裹算法运行速度上具有明显优势。