相移算法

相移算法（精选7篇）

相移算法 篇1

波束形成技术是指将一定几何形状排列的多元基阵的各阵元输出经过加权、延时和求和等处理后形成空间指向性的信号处理方法。它的作用, 一方面是进行空间处理以获取抗噪声和混响干扰的空间增益, 另一方面是为了得到高的目标分辨力, 用以测定目标方位。现役某型航空吊放声纳的信号处理机中, 所使用的处理器处理速度较慢, 且已经停止生产, 所以需要设计基于新型DSP处理器的波束形成系统。

波束形成的基本原理是将基阵的各阵元接收到的信号经过延时、加权、求和以在某个方向上的输出同向叠加成最大, 其他方向相应变小, 这便形成了某一方向上的集中响应。运算表达式如 (1) 式。

$y(t)={\displaystyle \sum _{i=1}^{{\rm N}}w}{}_{i}s{}_i(t-\tau {}_i)$ (1)

(1) 式中, N为阵元总数, si为第i个阵元接收到的信号, τi是为了形成某方向的波束而引入第i个阵元的时间延迟, wi为第i个阵元的幅度加权因子。波束形成基本原理如图1所示 。

1 宽带数字相移波束形成算法

基于上述原理的波束形成, 对任意频率的入射声波信号, 只要施以固定的阵元时延, 则主波束方向不会改变, 而各阵元所需时延只与阵形几何参数以及要求的主波束方向有关。对于单载波信号, 一定的时延等效一定的相移, 对于有一定带宽的信号, 时延可以用相移近似。因此, 时延问题就变成了如何在一个宽带内实现任意恒定的相移。下面介绍的正交调制法可以通过调节两个实数权系数实现恒定相移。

设具有振幅A、载频f0、初相Φ0和相移θ的单频信号为s (t, θ) , 则

s (t, θ) =Acos (2πf0t+Φ0+θ) =Acos (2πf0t+Φ0) cosθ+Asin (2πf0t+Φ0) (-sinθ) =Acos (2πf0t+Φ0) v1+Asin (2πf0t+Φ0) v2 (2)

其中v1=cosθ, v2=-sinθ。

(2) 式表明, 通过对信号和它的正交信号分别乘以v1, v2则可实现对Acos (2πf0t+Φ0) 的任意相移θ, 这就是正交调制相移法的原理, 如图2所示。

图2表明, 该方法需要对原信号实现90°相移, 根据数字信号处理理论, Hilbert变换器具有这种特性[2]。设连续时间信号x (t) , 其Hilbert变换$\widehat{x}$ (t) 定义为

$\begin{array}{l}\widehat{x}(t)=\frac{1}{\text{π}}{\displaystyle {\int }_{-\infty }^{\infty }\frac{x(\tau )}{t-\tau }}\text{d}\tau =\frac{1}{\text{π}}{\displaystyle {\int }_{-\infty }^{\infty }\frac{x(t-\tau )}{\tau }}\text{d}\tau \\ =x(t)*\frac{1}{\text{π}t}(3)\end{array}$

由傅立叶变换的理论可知, jh (t) =j/πt的傅立叶变换是符号函数sgn (Ω) , 因此Hilbert变换器的频率响应

${\rm H}(\text{j}\Omega )=-\text{j}\text{s}\text{g}\text{n}(\Omega )=\{\begin{array}{cc}-\text{j}‚\begin{array}{ccc}& \Omega >0& \end{array}& \\ \text{j}‚\begin{array}{cc}& \Omega <0\end{array}& \end{array}(4)$

若记

H (jΩ) =1 (5)

$\phi (\Omega )=\{\begin{array}{cc}-\text{π}/2‚\begin{array}{cc}& \Omega >0\end{array}& \\ \text{π}/2‚\begin{array}{cc}& \Omega <0\end{array}& \end{array}(6)$

(6) 式说明, Hilbert变换器是全通滤波器, 信号通过Hilbert变换器后, 其负频率成分作+90°相移, 正频率成分作-90°相移, 因此, 可以获得对输入信号的90°相移。根据系统设计要求, 采用切比雪夫最佳一致逼近法设计了一个11阶Hilbert滤波器。Hilbert变换器的权系数和相移系数, 均为事先求好存入程序中待调用。考虑到Hilbert变换器在数字频率0.05-0.45范围内有满意的频率响应, 且解调后回波信号中心频率为840 Hz, 带宽为1000 Hz, 所以取采样频率为3 424 Hz, 则要设计的Hilbert变换器的下限频率fl=340/3 424=0.1, 上限频率fh=1 340/3 424=0.4, 得到Hilbert变换器的权系数:

h1=[-3.796 531×10-2, -1.251 21×10-7, -1.426 754×10-1, -6.908 435×10-8, -6.102 951×10-1, 0, 6.102 951×10-1, 6.908 435×10-8, 1.426 754×10-1, 1.251 21×10-7, 3.796 531×10-2]。

频率响应如图3所示。

下面以圆弧阵为例, 介绍并分析本文的波束形成算法的原理和性能。圆弧阵设有12个基元, 均匀分布在圆阵上。相邻5个基元形成一个波束, 将在水平方向上形成张角为30°的12个波束。图4是圆弧阵及其12个基元分布示意图。

把基元按顺时针方向编号为H1, H2, …H12圆心通过基元H1的方向选作0度方向。为了讨论方便, 把时间的参考点选在圆心O上。设回波信号来自θ方向, 到达O点的信号假定为Acos (2πf0t) 。其中f0为信号中心频率, 第i个基元Hi所接收到的信号为 (为简单起见, 设Φ0=0)

si (t) =Acos[2πf0t+τi (θ) ] (7)

(7) 式中τi (θ) 为Hi相对于圆心O的延时, 即

τi (θ) =r0cos[θ- (i-1) α0]/c (i=1…n) (8)

(8) 式中c表示水中声速, r0表示圆弧阵半径, α0表示相邻两基元夹角。

以1号波束为例, 要形成1号波束, 需要5个基元H11, H12, H1, H2, H3参加定向, 1号波束最大响应轴在0°方向, 我们应将Si (t) 信号加以时延τi (θ0) 。根据 (8) 式分别算出H11, H12, H2, H3在H1方向上的时延, 此时各基元上的信号对应圆心O应延时的时间分别为:

τ11 (θ0) =0.033 3 ms, τ12 (θ0) =0.057 7 ms,

τ1 (θ0) =0.066 6 ms, τ2 (θ0) =0.057 7 ms,

τ3 (θ0) =0.033 3 ms。

进而由

Φi=2πfτi (i=1, 2…12) (9)

得到各基元上的信号对应圆心O应加的相移分别为

Φ11=Φ3=2πf0τ11 (θ0) =2.887=165.4°,

Φ12=Φ2=2πf0τ12 (θ0) =5.003=286.7°,

Φ1=2πf0τ1 (θ0) =5.775=330.9°。

这样, 经过相位加权后的信号为

Si[t-τi (θ) ]=Acos[2πf0t+τi (θ) -τi (θ0) ] (10)

1号波束形成的输出信号为

$y{}_1(t)={\displaystyle \sum _{11}^{3}w}{}_{i}s{}_i[t-\tau {}_i(\theta {}_0)]$ (11)

(11) 式中wi为权系数, 为满足指向性要求和后续处理方便, 取w1=w2=w12, w3=w11=1。

同理, 可以推广到其他11个波束形成的处理中, 形成全部12个主波束。最后, 通过实时检波算法实现对波束的检波。检波算法如下:

$\begin{array}{l}y(n+1)=\sqrt{\frac{1}{n}{\displaystyle \sum _{i=1}^{n+1}2}x{}^2(i)}=\\ \sqrt{\frac{1}{n+1}\left[{\displaystyle \sum _{i=1}^{n}2}x{}^2(i)+2x{}^2(n+1)\right]}=\\ \sqrt{\frac{n}{n+1}y{}^2(n)+\frac{2}{n+1}x{}^2(n+1)}(12)\end{array}$

2 算法的DSP软件实现

DSP的软件设计主要包括以下内容: (1) 从接口分机读取12路数据, 存储到数据缓冲区 (本系统采用双缓冲区设计, 即一个缓冲收满开始处理, 另一个缓冲继续收数, 两缓冲连续交替) ; (2) 用FIR 滤波器实现信号的Hilbert变换; (3) 每路信号与相应相移系数加权求和, 得到12路主波束和24路半波束; (4) 波束输出。

下面以1号波束为例, 进行波束形成的分析计算。

将上文得到的相移分别作正弦、余弦, 得到相应的相移系数分别为:

v1=[-0.967 709 17, 0.287 360 519,

0.873 772 223, 0.287 360 519, -0.967 709 17]

v2=[-0.252 069 358, 0.957 822 494,

0.486 335 38, 0.957 822 494, -0.252 069 358]

波束形成的核心算法的程序流程如图5所示。

程序流程说明: (1) 读入接口分机发送的12路阵元信号; (2) 先以H1作为基准, 将阵元H11、H12、H1接收到的信号分为两路; (3) 其中一路经过Hilbert滤波形成90°相移, 另一路不变; (4) 两路信号分别与相移系数相乘后再求和, 形成H1的左半波束; (5) 用相同的方法, 通过H1、H2、H3形成H1的右半波束; (6) 将H1左、右半波束相加, 得到基于阵元H1的主波束; (7) 循环阵元H1到H12, 依次得到12路主波束、24路半波束。

在这一过程中, 关键的步骤是实现两路信号与相移系数的加权并求和。这里用汇编语言编制了函数mul16 () 专门实现这一运算。

此函数调用如下mul16 (vx, vy, vr, n) 。其中vx为信号, vy为相应的相移系数, vr为运算结果, n为运算长度。汇编程序的主程序如下:

STLM a, ar_x ; //pointer to x

MVDK *sp (arg_y) , * (ar_y) ; //pointer to y

MVDK *sp (arg_z) , * (ar_z) ; //pointer to z

LD *sp (arg_n) , a ; //a = n

SUB #1, a

STLM a, brc ; //brc = n - 1

_START:

LD #0, a

RPTB ELOOP-1

MAC *ar_x+, *ar_y+, a ; //a = X*Y+a

ELOOP

STH a, *ar_z+;

STL a, *ar_z;

_END:

首先, 分别将x, y, z 3个变量的寄存器设置好并读入x, y的数据;其次, 读入运算长度n并存入brc中以用于RPTB循环;最后, 将x, y相乘求和并存入z中。

通过这一函数, 可以快速实现16位的信号加权求和, 使整个波束形成的过程速度提高。

对用本文提出的方法完成的波束形成信号进行检波积分, 并使入射角在0-360°范围内扫描, 得出波束的指向性图, 并把它同理论公式给出的波束指向性图进行比较。 (以180°作为基准方向)

通过分析对比理论指向性图和实际指向性图, 我们看到二者非常接近, 这说明用本文的方法实现的波束形成在工程上可行性, 结果可靠。另外, 通过实测波束形成的时间, 在TMS320VC5416处理板上波束形成程序处理时间约为10 ms, 相比原有的声纳信号处理机的20 ms处理时间, 提高了近一倍, 可见本文的实现方法实现简单, 运算速度快, 提高了现役声纳信号处理机的处理速度。

3 结论

本文主要介绍了宽带数字相移波束形成算法的数学原理和DSP软件实现方法。通过与理论指向性图的对比分析, 该方法形成的波束指向准确, 实现简单, 有较高的工程应用价值。该方法已经在某新型航空吊放声纳的信号处理机上得到应用。新型信号处理机采用TMS320VC5416芯片作为处理器, 处理速度得到了较大幅度提升。

参考文献

[1]赵航芳, 祝献, 宫先仪.混响背景下的信号检测.哈尔滨工程大学学报, 2004;25 (1) :34—37

[2]Ma Y L., Zhang B S., Investigation to the supergain Hydrophone Arrays for wideband application in isotropic ambient noise field.Proceedings of ECUA2000, Lyon, France, 10—13July2000:1043—1048

[3]Zhang Y W., Ma Y L.An efficient architecture for real-time narrow-band beamforming.IEEE Journal of Oceanic Engineering, 2001;10 (4) :635—638

[4]王英哲, 杨益新, 张忠兵, 等.多波束宽带恒定束宽波束形成器的DSP实现.电子技术, 2000; (7) :50—52

[5]胡广书.数字信号处理理论、算法与实现.北京:清华大学出版社, 2003

[6]乔瑞萍, 崔涛, 张芳娟.TMS320C54XDSP原理及应用.西安:西安电子科技大学出版社, 2005

[7]李启虎.声纳信号处理引论.北京:海洋出版社, 2000

同步相移干涉表面精密测量 篇2

表面加工质量的保证, 通常由精密制造设备和合适的制造工艺来实现, 其前提是制造系统必须具有良好的稳定性和可靠性。然而, 高精密加工时, 表面加工质量极易受到制造系统及其单元在静力学、动力学和热力学方面运行状态的影响[1,2]。为保证表面质量, 现有的高精密表面加工 (如表面精密磨削、金刚石单点加工等) 不得不采取经常性的加工表面离线检测, 以便于可能需要的工艺调整或补偿。如此质量控制方式将因为加工表面离线检测和重新置位的过程, 工艺连续性、加工环境和状态的一致性受到破坏而降低工艺调整和补偿的可靠性, 同时也降低效率。表面在线测量能保持工艺连续性, 保持测量前后加工环境和状态的一致性, 从而保证由表面检测所获得的工艺调整和质量控制的依据可靠, 实现更为可靠的质量保证, 并提高加工效率, 因而受到高度重视。然而, 一般精密表面测量方法 (如激光波面干涉测量方法) 由于对外界环境极度敏感, 易受到在线环境下振动干扰的限制, 难以用于在线测量。因此, 表面在线精密测量是高精密制造领域的一个难题。

微电子、微机电系统 (MEMS) 、光电子信息技术和航空航天技术等的快速发展, 对表面质量要求正稳步提高[3], 表面在线精密测量及质量控制问题因而变得更为突出, 需开发适用于在线条件下表面精密测量和有效质量控制的方法[3,4]。

为此, 本文基于激光波面干涉表面测量原理[5,6,7], 提出以偏振相移[5]为基础的同步相移技术, 研究在线在机表面精密测量的解决方案——同步相移干涉精密表面在线测量系统。该系统在传统泰兰-格林激光波面干涉仪[7]基础上, 利用偏振干涉[8], 采用偏振的同步相移技术, 有效地抑制在线环境对波面干涉测量精度和可靠性的影响, 从而实现精密加工平面和球面的表面形状在线测量。

1 测量原理

1.1系统原理

提出的同步相移激光干涉表面精密在线测量系统原理如图1所示。图中, L1为准直透镜;L2、L3为扩束透镜;BS为消偏振分光棱镜;P1、P2、P3为检偏器;M1、M2为反射镜;Q1、Q2、Q3为1/4波片。He-Ne激光器发出的光经空间滤波器、偏振片和准直镜L1后, 成为准直平面波进入偏振分光棱镜, 偏振方向垂直纸面的分量经分光面反射, 经过1/4波片Q1到达参考表面, 并被反射回来穿过分光棱镜成为参考光波, 该光波由于两次通过Q1使偏振方向变为平行于纸面, 从而穿过偏振分光棱镜分光面进入干涉光路;偏振方向平行于纸面的偏振分量穿过分光面, 经过1/4波片Q2后入射到被测表面, 反射回来的光波为带有被测表面偏离信息的准平面波。该准平面波再经过Q2, 偏振方向变为垂直纸面并由偏振分光棱镜分光面反射进入干涉光路。这两束偏振方向正交的光波在干涉光路中共路经过1/4波片Q3后变成旋转方向相反的圆偏振光, 它们一同经过L2和L3组成的扩束镜扩束或收束至合适于后面CCD的孔径后, 进入同步相移干涉单元。在同步相移单元, 多步相移同时发生并被同步采集送入计算机进行分析和干涉相位恢复[9,10], 得到从被测表面返回的波面波前形状, 从而得到被测表面面形。用于球面测量时, 则在偏振分光棱镜后置放一适当聚焦镜, 产生标准球面波入射到被测表面, 返回的波带有被测表面偏离信息, 随后进入测量系统, 得到测量结果。

1.2同步相移原理

相移干涉技术是激光波面干涉实现表面高分辨率高精度测量的重要技术, 其测量不确定度达到λ/100[5,6,9,10] (λ为激光波长) 。

常见的相移干涉技术大多采用参考镜的步进驱动来实现相移, 是分时相移, 以获得多步相移干涉图, 进而基于相移算法进行干涉相位的分析和恢复。分时相移中, 诸如振动等时域变化的环境因素影响将在相移过程中引入额外的相移, 成为相移误差, 从而引入干涉相位的恢复误差, 进而造成测量误差。

为此, 本系统基于偏振干涉原理[5], 提出同步相移技术, 建立同步相移单元, 同时获得产生不同相移的多幅相移干涉特征, 以避免环境振动等时域因素引入相移误差进而引起测量结果误差。

采用的同步相移原理如图2所示。进入的旋转方向相反的两共路相干光波被消偏振分光棱镜BS1、BS2和BS3分为光强相等的三组相干光束 (其中光路中使用的消偏振分光棱镜可以将入射光分为能量相等的两束出射光, 而且对偏振态不影响) , 分别通过检偏器P1、P2、P3后发生三路偏振干涉[8], 偏振干涉相位对应被测表面形状误差及检偏器的检偏角度。

以穿过BS1的一组相干光波为例。在图1中, 假设由参考表面和被测表面反射回来的出射到1/4波片Q3、偏振态分别沿平行于纸面的x向和垂直于纸面的y向正交的两光波电场分量表示为琼斯向量:

式中, E1、E2分别为两光波的电场分量幅值;wck (x, y) 、wbc (x, y) 分别为两光波的波前相位。

1/4波片的琼斯矩阵为

检偏器的琼斯矩阵为

Jp=[cos θ sin θ] (3)

式中, θ为检偏器的检偏角。

这两相干光波穿过Q3、BS2、BS1和检偏器P1后, 幅度减少为原来的1/4, 出射的光波总的电场分量可表示为

$\begin{array}{l}\mathit{E}=\mathit{J}{}_{\text{p}}\mathit{J}{}_{\lambda }(\frac{\mathit{E}{}_1}{4}+\frac{\mathit{E}{}_2}{4})=\\ [\cos \theta \sin \theta ]\left[\begin{array}{cc}\begin{array}{l}\\ 1\end{array}& \begin{array}{l}\\ -\text{j}\end{array}\\ -\text{j}& 1\end{array}\right]\left(\begin{array}{c}\begin{array}{l}\\ \frac{1}{4}E{}_1\text{e}{}^{-\text{j}\text{w}{}_{\text{c}\text{k}}(x,y)}\end{array}\\ \frac{1}{4}E{}_2\text{e}{}^{-\text{j}\text{w}{}_{\text{b}\text{c}}(x,y)}\end{array}\right)=\\ \frac{E{}_1}{4}\text{e}{}^{-\text{j}(w{}_{\text{c}\text{k}}(x,y)+\theta )}+\frac{E{}_2}{4}\text{e}{}^{-\text{j}(\frac{\text{π}}{2}+w{}_{\text{b}\text{c}}(x,y)-\theta )}(4)\end{array}$

那么, 在CCD上接受到的干涉光强分布为

可以看出, 两束相干光的相位差不仅与被测相位有关, 而且还与检偏器的检偏角成线性关系。因此, 三路偏振干涉相位分别受到2倍于P1、P2、P3检偏角度的相位调制。于是, 通过设置检偏器P1、P2、P3合适的检偏角, 使得三路相干光波偏振干涉分别产生相位值为0、2π/3、4π/3的相移, 三步相移干涉特征得以同时获得。

三步相移干涉特征由同步控制采集卡控制的同步CCD1、CCD2、CCD3同时获取, 实现同步相移干涉测量。

1.3抗振特性分析

采用偏振相移的方法获得的三步相移干涉特征分布为

Ii (x, y) =Ia+Ibsin[Δw (x, y) +αi]} i=0, 1, 2 (6)

式中, αi为第i步相移量;wbc (x, y) 、wck (x, y) 为对应被测表面和参考表面的相干波面相位分布;Δw (x, y) 为被测表面和参考表面对应点的相位差。

采用同步偏振相移, 获取α0=0、α1=2π/3、α2=4π/3三步相移干涉特征I0、I1、I2:

由于环境振动干扰, 被测表面在测量过程中存在振动, 每时刻除偏振相移产生的相移外, 还会由于振动引入相位变化, 从而产生一定的相移误差。设被测表面存在振动, V (t) =Vmsin ω t (Vm为最大振幅, ω为振动频率) , 图3为振动引起相移误差示意图。

若采用一般的分时相移, 设相移在时间上间隔为td, 开始相移时间ts是随机的, 每步相移采样时振动引起的相移误差φi等于该时刻振动量Vmsin (ω t) , 即

以三步相移为例, 三步获得的相移分别为

由于时间ts的随机性, 采样间隔td的不一致, 每步的相移误差φi都是不确定的, 可导致大的测量误差。

而采用同步相移技术, α′0=Vmsin ω ts, α′1=2π/3+Vmsin ω ts, α′2=4π/3+Vmsin ω ts。同步采集获得干涉特征为

即由于相移时间上的同步, 使振动引入的相移误差都是Vmsin ω ts。采用相移算法得到干涉相位为[10]

$\begin{array}{l}\text{Δ}{w}^{\prime }(x,y)=\text{Δ}w(x,y)+V{}_{\text{m}}\sin \omega t{}_{\text{s}}=\\ \arctan \frac{{\rm I}{}_3-{\rm I}{}_2}{{\rm I}{}_1-{\rm I}{}_2}(11)\end{array}$

得到表面被测点高度为

$\begin{array}{l}h(x,y)=\frac{\text{Δ}{w}^{\prime }(x,y)}{4\text{π}}=\frac{\lambda [\text{Δ}w(x,y)+V{}_{\text{m}}\sin \omega t{}_{\text{s}}]}{4\text{π}}=\\ \frac{\lambda [\text{Δ}w(x,y)]}{4\text{π}}+\frac{\lambda (V{}_{\text{m}}\sin \omega t{}_{\text{s}})}{4\text{π}}(12)\end{array}$

由于每个被测点高度都包含相同的误差项λ (Vmsin ω ts) / (4π) , 即由于振动引入的相移误差只是将被测表面的每个点提升了同样的高度, 而对最后被测表面面形误差测量结果没有影响, 故同步相移技术消除了振动引入的测量误差。

2 实验装置与实验结果

为了验证理论分析结果, 模拟了在现场振动情况下的同步相移测量和普通分时相移测量实验, 实验装置如图4所示。其中, 图4a所示为模拟被测样品振动的装置。样品固定于压电陶瓷的一端, 压电陶瓷的另一端固定。压电陶瓷电压输入端接入PZT驱动器, 由PZT驱动器产生一定频率、幅值和波形的电压信号, 驱动压电陶瓷产生动态位移, 带动被测样品做振动, 以模拟现场振动。图4b为同步相移测量系统装置图, 其光源采用5mW He-Ne激光源, CCD为敏通公司的MTV-1881EX黑白CCD, 能实现同步控制, 同步采集卡为西安微视公司生产。

将一个边长为20.0mm的被测Wafer样品固定于压电陶瓷上, 压电陶瓷电压输入端由PZT驱动器输入频率为50Hz、电压值30～110V变化的锯齿波信号, 驱动压电陶瓷产生频率为50Hz、峰值为4.1μm、呈锯齿波变化的动态位移, 带动样品振动。在振动的情况下, 由测量系统分别进行同步和分时相移测量, 同步相移按上述同步相移原理实现, 分时相移由另一PZT在被测表面振动的同时步进驱动参考表面实现。得到的相移干涉图见图5。

将图5中的相移干涉图分别取截面灰度显示, 如图6所示, 图中, 3条曲线之间的错位量对应实际相对相移量。可以看出, 图6a中错位量不均等, 实际相移量远离名义相移量;而在图6b中, 错位量均等, 实际相移量近似等于名义相移量。因此, 在确定的名义相移量情况下, 由于分时采样时刻相对被测表面振动引起相移周期的随机性, 三幅相移干涉图间的实际相移量存在不确定的误差;而同步相移中, 三幅相移干涉图在被测表面振动周期中的同一时刻获取, 实际相移量只与偏振相移产生的准确性有关, 不受振动影响。

将分时和同步采集的干涉图分别送入计算机进行处理, 完成相位信息提取, 计算表面形貌, 得到表面形貌结果分别如图7a和图7b所示。

结果显示, 在振动干扰情况下, 相对于分时相移, 同步相移测量方法较好抑制了振动影响, 取得了理想的测量结果, 而分时相移则由于振动引起相移混乱, 得到不可信的测量结果。

3 结论

针对精密表面加工中的在线在机质量监测问题, 基于波面干涉表面形貌精密测量原理, 结合提出的以偏振相移为基础的同步相移技术, 研究了具有抗振特性的精密表面激光波面干涉测量方案。介绍了偏振干涉同步相移表面测量原理, 分析了其抗振特性, 建构了测量系统实验装置, 并完成了测试实验。分析和测试结果验证了研究的同步相移激光波面干涉测量系统用于在线表面精密测量的适应性。

参考文献

[1]Chen Mingjun, Zhao Qingliang, Dong Shen, et al.The Critical Cconditions of Brittle-ductile Transi-tion and the Factors Influencing the Surface Qualityof Brittle Materials in Ultra-precision Grinding[J].Journal of Materials Processing Technology, 2005, 168 (1) :75-82.

[2]Jager V E, Hinzmann G, Just A.The Influence ofThermal Expansion of the Cutting Tool during Ul-traprecision Diamond Turning of Electroless Nickelon the Accuracy of the Workpiece and a SimpleMethod of Correction[C]//Europe Society for Pre-cision Engineering.1st International Conference andGeneral Meeting of the European Society for Preci-sion Engineering and Nanotechnology.Bremen:EU-SPEN, 1999:143-146.

[3]Richter U, Hentsch W, Liepack H, et al.In SituMeasurements of Optical Film Parameters andPlasma Monitoring during Reactive Sputtering forAdvanced In-line Process Control[J].Thin SolidFilms, 2009, 516:4541-4545.

[4]Rosati G, Boschetti G, Biondi A.et al.On-line Di-mensional Measurement of Small Components onthe Eyeglasses Assembly Line[J].Optics and La-sers in Engineering, 2009, 47:320-328.

[5]苏大图.光学测试技术[M].北京.北京理工大学出版社, 1996.

[6]吴震.光干涉测量技术[M].北京.中国计量出版社, 1993.

[7]Malacara D.Optical Shop Testing[M].2ed.NewYork:John Wily&Sons, Inc, 1992.

[8]Shamier J.Optical Systems and Processes[M].Washington D C:SPIE Optical Engineering Press, 1999.

[9]Meng X F, Cai L Z, Wang Y R.Wavefront Recon-struction by Two-step Generalized Phase-shift-ing Interferometry[J].Optics Communications, 2008, 281:5701-5705.

一种全新的波片相移方法 篇3

0 引言

人眼只能探测到的波长范围在400nm到700nm之间。为了能够获得更多的肉眼不可测得信息, 人们采用开发了大量探测设备, 如红外探测[1]、光谱探测[2,3]、超分辨率探测[4]、太赫兹检测[5,6,7]等。然而, 即使是在人眼可测波段, 眼睛所能获得的也只是电磁波的振幅信息, 而不能获得相位信息。为了获得相位信息, 可以采用相移法。相移法是光学干涉计量中普遍应用的直接相位测量技术。相移技术的基本思想是在干涉光路中引入相移器, 从而使得参考光与物光波的相位差增加一个可控参考位相, 通过几步相移, 测出物光波前位相, 提高测量精度。对于同一试件状态, 变化参考相位得到不同的干涉图, 通过对各点灰度进行简单计算, 可以得到全场位相分布。相移技术在干涉计量[8,9]、显微成像[10,11,12]、粒度分析[13]、温度场测量[14]、三维形貌测量[12,15]以及图像加密[16,17]等方面得到了广泛的应用。

目前传统的相移方法主要为时间相移法和空间相移法。时间相移方法包括压电陶瓷 (PZT) 法[18,19]、移动光栅法[20]、拉伸光纤法[21]、液晶相移法[22]、偏振相移法[23]、空气相移法[24]等。但是传统的时间相移方法记录的图像是在同一空间位置不同时刻采集的, 因此时间相移局限于对静态或准静态相位的测量。空间相移与时间相移恰好相反, 干涉图为同一时刻不同空间位置获得的, 因此可以进行动态测量, 主要有“普通分光棱镜分光+偏振相移”[25]、“光栅分光+偏振相移”[26,27]、“光栅分光+光栅相移”[28]三类。空间相移虽然能够进行动态测量, 但是由于干涉图记录的空间位置不同, 要求不同探测器 (若采用多个探测器接收不同相移干涉图) 或探测器不同部分 (若采用单个探测器的不同部分接收不同相移干涉图) 的光电性能一致, 且不同空间位置的干涉图之间要求进行像素级的位置匹配和灰度校正, 因此实际操作过程难度高精度低, 难以被实际应用。

本文提出了一种新的波片相移方法, 将物光波和参考光调制为互相垂直的线偏振光, 在干涉光路中用波片引入相移量, 实现时间相移。基于这一原理制备波片阵列, 还可以实现实时相移。文中第二部分给出了波片相移光路图, 并给出了该相移原理的理论推导;第三部分采用数值模拟的方法对波片相移方法进行验证, 波片相移法重构的物面的复振幅与传统四步相移法重构的物面的复振幅相一致;第四部分采用100×物镜观察叶片气孔细胞, 当图像传感器 (CCD) 未放置在物镜的像平面时, 通过波片相移法重构记录平面的复振幅, 并对其进行菲涅尔变换以获得像平面的复振幅。实验成功获得像平面处物光波的光强分布, 证实了波片相移方法的正确性。第五部分基于波片相移原理提出了波片阵列结构, 通过将其与偏振片集成在图像传感器前可以实现实时相移数字全息, 能够测量动态物体。

1 原理

波片相移光路见图1, 由激光器发出的光经过扩束准直后照射偏振分光棱镜, 经偏振分光棱镜透射和反射的光为两束偏振方向互相垂直的线偏振光, 其中反射光照射到待测物体上, 为物光, 透射光为参考光。物光和参考光分别经过全反镜反射后到达第二块偏振分光棱镜, 被反射的物光和被透射的参考光经过波片和偏振片被CCD接收。经过第二块偏振分光棱镜后物光波振幅Es, 频率ω, 初始相位φ1;参考光振幅Er, 频率ω, 初始相位φ2。位于CCD前的偏振片的透偏振方向与物光波经第二块偏振分光棱镜反射时的偏振方向一致。位于偏振片前的四分之一波片和二分之一波片的快慢轴分别为f轴和s轴。偏振片的透偏振方向与波片快轴夹角为θ。

将物光波与参考光矢量沿波片快慢轴方向分解, 其矢量形式为

通过波片后慢轴相对于快轴相位延迟α, 则

偏振片透偏振方向与f轴夹角θ, 因此透过偏振片后光的振幅为

当不放置四分之一波片和二分之一波片时, 慢轴相对于快轴的相位延迟量α=0, 此时记录平面的振幅E0和光强I0分别为:

当只放置二分之一波片时, 慢轴相对于快轴的相位延迟量α=π, 此时记录平面的振幅Eπ和光强Iπ分别为:

当只放置四分之一波片是, 慢轴相对于快轴的相位延迟量α=π/2, 此时记录平面的振幅Eπ/2和光强Iπ/2分别为:

当同时放置四分之一波片和二分之一波片时, 慢轴相对于快轴的相位延迟量α=3π/2, 此时记录平面的振幅E3π/2和光强I3π/2分别为:

由上式可得,

因此只需要选择三个相位延迟量即可。

联立公式 (7) 、 (9) 、 (11) 可得, 当θ=0, π/4, π/2, 3π/4时, 方程组有多解, 当取其他值时, 有唯一解。因此θ取值不能为0, π/4, π/2, 3π/4。为了计算方便, 取θ=π/8。则

由上式可得

式中只有Ir一个变量, 因此可解出Ir值, 也即解出了Er值。将Er值代入公式 (17) 和 (18) , 可分别得cos (Δφ) 值和sin (Δφ) 值, 进而可得Δφ值。

因此, 采用此种方法可解出Es, Er, Δφ, 也即得到了记录平面上的物光的复振幅u (x, y) 。物面的复振幅U (X, Y) 可以通过菲涅尔变换得到:

2 数值模拟

为了验证理论分析的正确性, 进行了如下数值模拟。给定物光波的振幅和相位, 记录CCD的像素数为1024X1024, 每个像素的尺寸为10um, 物面与CCD距离为10cm, 激光波长为532nm, 物面的光强和相位见图2.a和图2.b。采用四步相移法, 获得记录平面的四幅光强图, 并重建物光波的图像, 见图2.c和2.d。采用波片相移法, 获得记录平面的四幅光强图, 见图2.i~图2.l, 并重构物光波的图像, 见图2.e和2.f。由图可见, 采用波片相移法与四步相移法获得物光波光强和相位信息基本相同。

3 实验验证

为了进行实验来验证上述技术的正确性, 我们采用100×长焦物镜观察叶片表面气孔细胞。光路图如图1所示, 待测物体位置放置叶片表面气孔细胞标本, 由物光波照射标本。图中a处放置100×长焦物镜用来对标本进行显微成像。CCD不在显微物镜的像平面上, 因此CCD上采集到的物光图像是由显微物镜像平面上的物光波经过菲涅尔变换后到达CCD上的物光波的光强图像, 为散焦的图像。通过波片相移法计算CCD上物光波的复振幅, 再进行菲涅尔变换即可获得物镜像平面处物光波的光强分布。实验时, 首先不加四分之一波片和二分之一波片, 采集图像, 获得I0;放入四分之一波片, 采集图像, 获得I1;撤掉四分之一波片, 放入二分之一波片, 采集图像, 获得I2。之后采用非偏振光源对四分之一波片和二分之一波片的透光率进行标定, 分别获得相应的的透光率系数矩阵。用I1和I2分别除以相应的透光率系数矩阵, 即可获得Iπ/2和Iπ。图3给出了叶片气孔细胞显微镜下图像、加入二分之一波片时记录平面干涉光强、记录平面参考光光强以及通过菲涅尔变换获得与记录平面不同距离z处的光强图像, 其中重构距离z选择为140mm到220mm, 间隔为10mm, 可见在距离记录平面180mm处, 物光波的光强图像最清晰, 与显微镜下拍摄的图像最一致, 也即叶片气孔细胞在100×物镜观察下的像平面位于距离记录平面约180mm处。

图2图中的图像都为256X256像素。物面:a.光强分布, b.相位分布;应用四步相移法重构的图像:c.光强分布, d.相位分布;应用波片相移法重构的图像:e.光强分布, f.相位分布;应用集成波片阵列方法重构的图像:g.光强分布, h.相位分布。il分别为应用波片相移法获得的记录平面的四幅光强图像I0, Iπ/2, Iπ, I3π/2。Fig.2Imagesare256X256pixels.Objectplane:a.intensitydistribution, b.phasedistribution;Imagereconstructedbyfourstepphaseshiftmethod:c.intensitydistribution, d.phasedistribution;Imagereconstructedbywaveplatephaseshiftmethod:e.intensitydistribution, f.phasedistribution;Imagereconstructedbywaveplatearraymethod:g.intensitydistribution, h.phasedistribution.i~haretheimagesofI0, Iπ/2, Iπ, I3π/2acquiredbywaveplatephaseshiftmethodintherecordplane.

4 实时相移数字全息

2004年, YasuhiroAwatsuji等人提出了实时准相移数字全息[29], 在参考光路中加入相位延迟器阵列, 该阵列每个单元与图像传感器的像素单元一一对应, 相位延迟器阵列每2×2单元的玻璃厚度不同, 使得物光波与参考光之间的相位延迟量分别为0, π/2, π, 3π/2。该方法实现了实时相移方法, 但相位延迟器位于参考光路中, 延迟器单元与图像传感器像素单元的对准较为繁琐, 且抗振要求高, 实际操作中难以实现。4Dtechnology公司的JamesMillerd等人提出了基于金属纳米光栅的像素偏振片阵列应用于实时全息的方法[30], 将偏振片阵列集成到CCD上, 解决了前述每次实验都需要对准的问题, 并且提高了抗震强度。但是基于金属纳米光栅的偏振片阵列的制作难度有以下几点[31]:a.由于金属纳米光栅制作受紫外光衍射极限的限制 (光栅栅线宽度约为100nm) , 标准的紫外光刻技术达不到要求, 需要对激光器光源倍频, 进行全息曝光或干涉曝光;b.由于纳米光栅的线宽仅有约100nm, 通过刻蚀方法将纳米光栅结构由光刻胶模板层转移到金属层上的成功率很低;c.要想获得具有四个偏振方向的偏振片阵列, 需要四次沉积铝层和二氧化硅层;d.由于偏振片阵列是基于像素级别, 要想获得具有四个偏振方向的偏振片阵列, 需要进行四次全息曝光或干涉曝光 (即a与b步骤要重复四次) , 并进行四次标准的像素尺寸掩膜版紫外光刻。鉴于国内尚且无法制备单一偏振方向的基于金属纳米光栅的偏振片, 因此制备具有四个偏振方向的基于金属纳米光栅的偏振片阵列更是难以做到的。目前掌握偏振片阵列制作技术的国家为美国 (4D technology公司) 和日本 (Photron公司) , 且不对中国出售偏振片阵列芯片, 因此找到一种可以替代偏振片阵列的结构是非常重要的。

基于本文提到的波片相移方法, 设计了一种基于像素的波片阵列, 可以应用于实时相移数字全息中。首先制作如图4.a所示的波片阵列, 其中每2×2个波片单元的慢轴相对于快轴的相位延迟量分别为0, π/2, π, 3π/2, 每个波片单元的尺寸与CCD像素的尺寸相同。偏振片可选用含碘化物的聚乙烯醇薄膜, 成本较低, 厚度约10um, 或选择刻蚀金属纳米光栅, 成本较高, 厚度仅100nm左右。将波片阵列、偏振片集成到CCD上, 波片阵列的每个单元与CCD的像素单元对准, 偏振片的透偏振方向与波片阵列的快轴夹角为π/8。

波片阵列采用双折射材料进行像素尺寸刻蚀, 刻蚀出不同的深度即可获得不同的相位延迟量, 对于双折射率较高的材料 (如方解石) , 波片阵列的厚度可达到仅有几个微米, 但由于刻蚀技术中对深度难以精确控制, 因此本文给出了一种波片阵列制作的可行性方案。如图4.b, 在玻璃基底上用紫外敏感胶粘合第一片真零级四分之一波片, 然后对该四分之一波片进行刻蚀, 得到的第一块四分之一波片中每2×2个单元的相位延迟量分别为0, 0, 0, π/2;在第一片四分之一波片上用紫外敏感胶粘合第二片真零级四分之一波片, 然后对该四分之一波片进行刻蚀, 得到的第二块四分之一波片中每2×2个单元中的相位延迟量分别为0, 0, π/2, π/2;在第二片四分之一波片上用紫外敏感胶粘合第三片真零级四分之一波片, 然后对该四分之一波片进行刻蚀, 得到的第三块四分之一波片中每2×2个单元的相位延迟量分别为0, π/2, π/2, π/2;因此其组合后的波片阵列中每2×2个单元的相位延迟量分别为0, π/2, π, 3π/2。

实验时, 采集单帧图像即可计算物光波的复振幅。图像中每2×2个单元的光强值分别为相位延迟量0, π/2, π, 3π/2时的光强值。计算记录平面物光波的复振幅方法如下:

计算机生成四幅空值图像, 图像的分辨率与CCD采集到的图像分辨率相同;

将CCD采集到的图像中所有相位延迟量为0的点的灰度值复制到计算机生成的第一幅空值图像中对应位置;将CCD采集到的图像中所有相位延迟量为π/2的点的灰度值复制到计算机生成的第二幅空值图像中对应位置;将CCD采集到的图像中所有相位延迟量为π的点的灰度值复制到计算机生成的第三幅空值图像中对应位置;将CCD采集到的图像中所有相位延迟量为3π/2的点的灰度值复制到计算机生成的第四幅空值图像中对应位置;

计算机生成的四幅图像中的空值像素的灰度值采用周围非空值像素灰度插值平均的方法补齐, 得到四幅完整的图像, 即I0, Iπ/2, Iπ和I3π/2;应用这四幅图像即可算出记录平面上的物光复振幅。

该方法的模拟效果图见图2.g和2.h, 分别为物面处物光波的光强和相位。可见, 采用该方法能够很好的还原物光的光强和相位信息。同时由于本方法只需要一幅记录图像即可计算物光复振幅, 因此可以应用于实时相移全息中, 可测量动态的物体的光强和相位。

与基于金属纳米光栅的偏振片阵列相比, 波片阵列的制作工艺:1.不需要步骤a和步骤b, 2.粘贴四层四分之一波片的工艺难度与成本远低于步骤c的四次沉积铝层和二氧化硅层的工艺难度和成本, 3.波片阵列需要四次像素尺寸的掩膜版紫外光刻, 这与步骤d里的难度与单元的对准精度要求是一样的, 4.波片阵列的制作不需要除以上3步额外的复杂工序, 因此波片阵列的制作难度和成本远低于基于像素的偏振片阵列, 成品率高于基于像素的波片阵列。

5 结论

本文提出了一种新的时间相移数字全息技术, 这项技术操作简单, 容易实现, 具有传统时间相移技术干涉图不需要位置匹配等优点。为了验证波片相移方法的正确性, 进行了数值模拟, 成功重构了物光波的复振幅, 与传统四步相移结果一致。另外还进行了实验验证, 通过波片相移方法显微观察叶片气孔细胞, 获得记录平面物光波的复振幅, 并通过菲涅尔变换获得显微物镜像平面处气孔细胞的光强图像。基于波片相移技术原理, 提出基于像素的波片阵列结构, 并将其与偏振片集成到图像传感器前, 既继承了传统时间相移方法不需要位置匹配等优点, 又具有传统空间相移能够进行实时相位测量的优点。本文对波片相移方法作了详细的理论推导, 并且采用数值模拟和实验验证了该方法的正确性。

6 附录

公式10到公式11推导:

基于数字光栅相移法的三维重构 篇4

文中基于数字光栅相移的三维测量技术,结合最新的相位高度映射关系和棋盘标定方法,通过空间变换减少对数字光栅和CCD几何约束要求,可操作性强,标定的过程简单。实现了对被测物高度等三维信息的准确重构。该方法快捷高效,适用性强。

1 数字相移测量

将仿真的正弦光栅条纹投射到被测物体表面,光栅条纹经被测物体表面调制后会发生变形,被测物体的高度信息就包含在这些发生变形的条纹相位信息中,只要找到相位和高度的映射关系就可以重构出被测物体的高度[5,6,7]。

数字相移法一般有三步相移或四步相移,相移N次需要采集N+1帧图像,文中采用三步相移法,三步相移法[6]计算被测物的包裹相位原理如下。

一幅理想的数字条纹图像的灰度;表达式如下

其中,I'(x,y)是平均灰度,或者是背景灰度,I'(x,y)是灰度调制值;φ(x,y)是要求解的相位场。

理想的三幅相移条纹图像的灰度表达式

联立式(2)、式(3)和式(4)可得到式(5)。把由式(5)得到的相位解包裹即可得到经被测物调制后的相位[8,9]和高度的映射关系

系统测量示意图如图1,使用的测量系统减少了CCD和投影仪之间约束关系,CCD和投影仪可以不在同一个平面,CCD和投影仪的光心可以不相交在参考面上,通过计算可得到相位和高度的关系式[10]

h(x,y)为(x,y)点的高度;Δϕ(x,y)为被测物相对参考平面的相位差;a(x,y)和b(x,y)为关系因子,可通过系统标定求得。

2 三维重构的棋盘法标定

当测量系统确定后,式(6)的a(x,y)和b(x,y)按理论是可以确定的[11,12],只需利用一个各角点空间坐标已知的棋盘格标定板,拍摄两幅光栅条纹图。

由式(7)和式(8)可得到

已知棋盘格大小,假设位置1时角点空间坐标Z坐标为0,根据CCD标定得到的2个位置的旋转和平移矩阵R1、T1、R2、T2,由式(11)可以得到位置2时角点在O1X1Y1Z1中的高度;利用三步相移法得到的相位图就可以计算出相应于2个不同位置角点高度差对应的相位差,得到高度和相位之间映射系数a(x,y)和b(x,y)。

具体步骤如下:

(1)CCD和投影仪位置固定后,使棋盘格标定板平行于CCD成像平面放置在位置1如图2所示;

(2)以标定板左上方的一个角点为原点建立如图3所示的世界坐标系0w1Xw1Yw1Zw1;

(3)用CCD标定工具栏进行标定,可以得到摄像机内部参数,以及位置1世界空间坐标系和摄像机坐标系XcYcZc之间的旋转矩阵R和平移矩阵T;

(4)旋转标定板一个角度到位置2,重复步骤(2)和(3),建立此位置的世界坐标系Ow2Xw2Yw2Zw2,计算位置2世界空间坐标系和摄像机坐标系之间的旋转和平移矩阵R1、T1以及光栅相位分布图。

理论上根据上面得到的光栅相位图和标定板各角点的高度信息就可以计算当前系统的被测物高度和光栅相位分布的映射关系。

3 实验结果

实验以截面为梯形高度为20 mm实物如图4所示,作为测量对象,用文中方法进行测量。首先利用边长为15 mm的黑白棋盘格,如图5,对CCD进行标定[13],并且消除图像畸变;为了获得被测实物的高度信息将数字光栅投影到被测实物,由CCD采集经被测实物高度调制的3幅数字光栅相移正弦条纹图如图6、图7和图8;对采集到的3幅数字光栅相移条纹图进行均值滤波和高斯滤波,再由三步相移公式计算得到包含被测实物高度信息灰度条纹图如图9;此时得到的相位图像不是被测实物的真正相位,它的相位被截断在-π和π之间,用解包裹算法展开后可以得到被测实物真正相位[14,15]三维分布图如图10所示;重构被测实物如图11所示,最大高度和最小高度差为19.854 8 mm。精度为99.27%,实验结果显示,文中所述方法正确恢复了被测实物的高度信息。

4 结论

相移算法 篇5

相对于光学全息,数字全息可以精确地分析物体的三维分布和相位信息,尤其对于微小物体,是一种理想的形貌和相位分布测量的方法。近年来,随着计算机运算速度及容量的大幅提高和高分辨力光电成像器件的出现,数字全息技术得到了迅速发展,并广泛用于缺陷检测[1],形貌测量[2],显微[3,4]等领域研究。

数字全息是通过光电成像器件实现全息图记录,而在目前的技术条件下,光电成像器件的分辨力与传统的全息感光材料是无法相比的。这样,成像器件所能记录的信息量受到很大的限制。对全息而言,我们可以通过空间自由度数研究全息系统的信息量。根据M.冯.劳德观点,一个光场的空间自由度正比于物(或像)面积和光学系统空间频带宽度的乘积:

式中:Lx,fxmax和Ly,fymax分别为物(或像)沿x,y方向的宽度和频率,SW称为空间带宽积。对成像器件而言,上式实际上是像面(或物面)图像的带宽积。根据劳-鲁克斯自由度不变定理,感光器件空间带宽积是数字全息所能记录的物体或者最终成像信息量的极限。如何充分利用感光器件带宽积,达到信息记录的极限是本文讨论解决的问题。

数字全息再现的原理是光学再现过程的数字模拟,因而同样存在零级像和共轭像。采用离轴全息的方法,可以使再现像和零级像、共轭像分离,但必将使得本来有限的像素数不能被再现像充分利用,降低了再现像的分辨力;采用同轴全息时,零级像和共轭像将与再现像重叠,产生严重的干扰,而不能实现像的再现。目前消除零级像和共轭像的方法,主要有两种:一是频域滤波法[5,6];二是在数字全息记录过程中加入相移技术[7,8],记录具有相对相移的四幅全息图,然后按照一定方式进行数字叠加后再现,从理论上完全可以消除零级像和共轭像,而且可以使由记录器件所决定的像素数全部用于再现实像,扩大原始像的视场,提高再现像的分辨力。但是本技术一般采用四步相移法,需要记录四幅全息图,增加了装置的复杂性和记录时间。

本文首先从信息论的角度证明了同轴傅里叶变换全息光路可以达到信息记录的极限,继而提出简化的相移数字全息消除零级像和共轭像的技术,仅需要拍摄两次全息图和一次物光强度图,记录过程中仅需一次相移。实验结果证实该方法相移操作简单,避免了多次相移带来的误差,很好地去除了零级像和共轭像,有效提高了数字全息记录的信息量,从而提高了再现像的质量。

1 数字全息信息量与记录光路的关系

使用图1来研究数字全息信息最佳记录条件。由图所示,物光和参考光在成像器件面上任一点处的干涉条纹间距由下式决定:

式中:λ为记录光波波长,α是参考光和物光波在成像器件面上的夹角。为了使得成像器件能够记录到干涉条纹,期望δ能够尽可能的大一些,即期望sinα尽可能小一些。下面以一维情况来讨论,参照图1记录光路,设物点的坐标为(xo,0),在记录平面坐标(ξ,zo)处,在满足菲涅耳衍射近似的情况下,可以证明:

假设Lox和LCCDx分别为矩形的物平面和成像器件在x方向的宽度。则上式最大值是:

可以证明,式(4)随Xr绝对值的增大而增大,当Xr=0时,成为同轴全息,此时,sinα最小:

式中:±LCCDx/2分别表示成像器件的上下边缘。

图2给出了在同轴全息情况下,sinα值在成像器件上下边缘随着参考光zr的变化曲线(LCCDx=6.8 mm,zo=300 mm,Lox=4 mm),虚线表示下边缘,实线表示上边缘。从图中可以看出,随着参考点光源的位置向成像器件移动,对于下边缘,sinα值是先变小,然后再变大;而对于上边缘,sinα值一直在变大。当zr=0时,上下边缘对应的sinα值一致,此时成像器件面上的干涉条纹密度达到最小,在满足抽样定理的前提下,对成像器件分辨力的要求最低。显然这正是同轴傅里叶变换全息光路。当zr=0时,式(5)简化为

对应的干涉条纹间隔为

根据抽样定理,干涉条纹的间隔至少是成像器件记录单元间隔两倍,条纹才能被纪录。设成像器件记录单元间隔为dCCDx,即要求δ≥2dCCDx。根据式(7)可得:

在成像器件参数和物体大小以及记录波长都确定的情况下,可以改变zo,使成像器件达到最高效的信息记录,但必须满足:

zo越小,记录到的信息越多,最佳记录距离为zo=dCCDxLox/λ,代入被测量物体的横向分辨力[9]:

得:Lox/domin=LCCDx/dCCDx,可以表示为:Loxfoxmax=LCCDxfCCDxmax,由式(1)可以看出,同轴傅里叶变换全息记录光路满足信息传递空间带宽积不变原理,成像器件分辨力被用到极限。

2 一次相移同轴无透镜傅里叶数字全息

假设从物面出射的光波分布为u(x,y),根据衍射菲涅耳近似,在成像器件接收面的复振幅分布应该为

参考球面波的复振幅分布表示为

则成像器件记录的全息图光强分布可以写为

通过相移装置使参考光的相位改变π/2,相应的成像器件面上全息图的强度分布可写为

式(13)和式(14)中的|O|2和|R|2分别是物光和参考光单独照明成像器件时强度分布,令:

联立式(13)、式(14)和式(15)得:

上式中:,将式(11)的共轭和式(12)带入上式,可以得到:

在傍轴近似的条件下,当参考光是发散球面波的时候,其在感光器件面上的振幅分布是均匀的,即Ar(ξ,η)可以认为是常数[10]。对上式进行傅里叶逆变换得:

显然,上式除了常数和二次位相权重外,我们获得了像面上的物光场分布,这正是无透镜傅里叶变换全息图的本质。为了得到强度像,只要求式(18)的模平方即可。整个过程除了全息图和物光强度的拍摄外,都很方便的利用计算机编程实现。

3 实验与结果分析

图3是本文采用的实验光路,成像器件选用CCD。记录波长为0.632 8µm的He-Ne激光,经分束镜(BS1)分为两路:一路经反射镜(M2)反射,扩束镜(BE2)扩束后照射到物体上。在物体前放置一块毛玻璃,使之产生漫射光均匀照明物体。物体经过显微物镜(MO)在直角分束棱镜(BS2)左侧平面处成一清晰放大的实像,放大倍数为10;另一路经过扩束准直后形成平行光,经反射镜(M1)反射,由透镜2在直角分束棱镜上侧面会聚形成一点光源,这样点光源和物体放大的实像相对CCD面可视为在同一平面上。参考球面波经过直角分束棱镜(BS2)反射,与物光波在CCD面处形成干涉全息图,被CCD接收并存储于计算机。

实验过程中,采用四分之一波片作为相移器(PSD),使参考光的相位改变π/2。首先使波片光轴与激光偏振方向一致,拍摄相移前的数字全息图,然后旋转波片,使其光轴转过90o,参考光的相位将改变π/2,拍摄相移后的数字全息图。然后通过控制快门(Shutter),单独记录物光的强度分布图。实验所用CCD像素尺寸为8.6µm×8.3µm,有效像素数为752×582,光敏面积为6.4 mm×4.8 mm。为了充分利用CCD的感光像素有效的记录信息量,应根据(9)式设计记录距离。

图4为3#分辨力板第9单元作为显微物体(条纹周期为50µm,透光宽度为25µm),通过移动CCD,改变物体到CCD的距离,得到在同轴记录光路下去除零级像和共轭像后的再现像。可以看出,随着记录距离的减小,再现像占据的像素增大,再现像的分辨力得到明显的改善。但是我们发现,当zo小到一定值时,再现像虽然增大,但边缘将变暗。其原因主要是由于再现像的强度分布受到sinc函数的调制,强度沿像场的中心向边缘逐渐减弱[9]。

图5为条纹周期为3.3µm的一正弦光栅(记录距离为60 mm)相移前局部放大全息图(a)和再现像(b),显微物体放大后的横向尺寸为4 mm,由式(9)知同轴记录光路下最小记录距离为zomin≈54.4 mm。此时,因为记录距离已经接近最小距离54.4 mm,且物体的分辨力比较高,物光基本可以充满CCD光敏面积,因而再现像基本占据了所有像素。

4 结论

本文详细研究了一次相移数字全息显微技术,通过同轴傅里叶变换全息记录光路验证并分析了消除零级像和共轭像的效果。结果表明:两步相移技术可以有效地消除零级像和共轭像,简化了实验步骤,避免了多次相移带来的误差。同时,将本技术和无透镜傅里叶变换全息光路相结合,使得记录距离最小,可以充分利用成像器件的带宽,有效提高了数字全息记录的信息量,再现像可以充分利用像素数,有利于提高再现像的分辨力。我们进一步的工作是,在消除零级像和共轭像的基础上,研究有效地消除再现像噪音的方法,更好地提高数字全息显微再现像的质量。

参考文献

[1]ZHANG Fu-cai,Jesus D R VALERA,Ichirou YAMAGUCHI,et al.Vibration analysis by Phase Shifting Digital Holography[J].Opt.Rev(S1340-6000),2004,11(5):297-299.

[2]Francisco PALACIOS,Jorge RICARDO,Daniel PALACIOS,et al.3D image reconstruction of transparent microscopic objects using digital holography[J].Opt.Commun(S0030-4018),2005,248:41-50.

[3]GUSTAFSSON M,SEBESTA M,Bengtsson B,et al.High-resolution digital transmission microscopy—a Fourier holography approach[J].Optics and Lasers in Engineering(S0143-8166),2004,41(3):553-563.

[4]Francisco Palacios,Daniel Palacios,Guillermo Palacios,et al.Methods of Fourier optics in digital holographic microscopy[J].Opt.Commun(S0030-4018),2008,281(4):550-558.

[5]Cuche E,Marquet P,Dpeursinge C.Spatial filtering for zero-order and twin-image elimination in digital off-axis holography[J].Appl.Opt(S0003-6935),2000,39(23):1024-1026.

[6]刘诚,李良钰,李银柱,等.无直透光和共轭像的数字全息[J].光学学报,2002,22(4):427-431.LIU Cheng,LI Liang-yu,LI Yin-zhu,et al.Digital holography free of Zero-order diffraction and conjugate images[J].Acta Optica Sinica,2002,22(4):427-431.

[7]Takaki Y,Kawai H,Ohzu H.Hybrid holographic microscopy free of conjugate and zero-order images[J].Appl.Opt(S0003-6935),1999,38(23):4990-4996.

[8]Yamaguchi I,ZHANG Tong.Phase-shifting digital holography[J].Opt.Lett(S0146-9592),1997,22(16):1268-1270.

[9]刘诚,朱健强.数字全息形貌测量的基本特性分析[J].强激光与粒子束,2002,14(3):328-330.LIU Cheng,ZHU Jian-qiang.Basic characters of digital holographic profiling[J].High Power Laser and Particle Beams,2002,14(3):328-330.

相移算法 篇6

全息技术最早于1948年由英国科学家Dennis Ga bor提出[1], 其基本原理是:利用干涉原理, 将物体发出的特定光波以干涉的形式记录下来, 使物光波前的全部信息都存储在记录介质中, 即“全息图”。当用光波照射全息图时, 由于衍射原理能够重现出原始物光波, 形成与原始物体逼真的三维像。其中, 数字全息技术用CCD等光敏电子成像器件作为记录介质, 采用数值计算的方法再现。因其操作简便、误差小等优点, 数字全息技术在近几年得到快速发展。1997年Yamaguchi又提出了将数字全息技术和相移干涉术相结合的相移数字全息技术[2], 在记录过程中改变参考光的相位, 即在参考光路上加入相位调制器件, 实现参考光相位的改变, 从而记录多幅全息图, 再通过数字计算得到物体的再现图。能够很好地消除共轭像和零级项对再现像的影响, 大大提高了再现象的质量, 为相移数字全息技术的广泛应用奠定了基础。

在对图像的研究和应用中, 人们往往只对图像中某些部分感兴趣, 这些部分通常被称为前景, 余下部分被称为背景。图像分割就是将所需要的前景分离出来。图像分割技术在全息图像处理中起着至关重要的作用[3]。为了更好的研究基于参考光估计的两步相移数字全息算法的可靠性和使用性, 因而对全息再现图的分割研究具有重要意义。

1基于参考光估计的两步相移数字全息原理

相移数字全息技术虽然能较好的消除共轭像和零级项, 但它存在的问题是要记录多幅全息图。记录次数多、曝光次数多, 误差就会较大[4,5]。在此, 提出了基于参考光估计的两步相移数字全息术。该方法的全息记录方法只需记录两张相移量分别是0, π 2的数字全息图。设在全息图平面即CCD平面上的物光波的复振幅分布为uH (xH, yH) , 参考光波为ur1= Arexp (i ⋅ 0) 和ur2= Arexp (i ⋅π/2) 。则两张全息图的强度分布为:

式中I0 (xH, yH) = Ar2+|uH|2。这里把参考光看成是常数计算, I0是零级项。

由全息图的成像原理可得Ar2的大小应该在0~ max (IH1) 之间, 构造一个二维相关系数CC (2D correla- tion coefficient) 如下:

式中abs (·) 代表取绝对值;代表二维相关运算;ET是由单幅加密全息图IH1进行正确反演得到的重建像;E是用两步全息方法所获得的重建像, 其中A2rc是参考光强Ar2的假设值, 通过扫描0~max (IH1) 之间的数值得到, 而不是由提前拍摄的参考光强图得到的。

基于上面的CC相关系数理论, 可通过计算机仿真得到Arc值, 也就是拍摄全息图时的真正参考光振幅值Ar。由此, 就可以求得原始物体的清晰再现像。如图2所示。

2水平集算法

水平集方法是由Sethian和Osher于1988年提出[6], 简单的说来, 水平集就是把低维的一些计算上升到更高一维, 把N维的描述看成是N+1维的一个水平[7]。

水平集方法是将水平闭合曲线隐含的表达为连续函数曲面 Φ (x, y, t) 的一个具有相同函数值的同值曲线。 通常将目标曲线隐含表示在零水平集函数{Φ (x, y, t) = 0} 中。设用于演化的平面闭合曲线为C (p, t) = (x (p, t) , y (p, t) ) , p为任意的参数化变量, t为时间。设曲线的内向法向量为N, 曲率为k, 则曲线沿其法向量方向的演化可以用下面的偏微分方程表示:由φ (C (t) , t) =0, 全微分可得:其中内向法向量整理可得:∂∇

上式为水平集曲线演化的方程。

水平集函数在迭代的过程中可能发生退化, 使它不再保持符号距离函数, 因此必须进行重新初始化操作, 以保证水平集函数接近一个符号距离函数, 从而保证数值解法的稳定性。

标准的重新初始化方法是通过解以下的Hamilton-Jacob方程实现的:

然而, 在演化过程中周期性地对水平集函数进行校正, 即重新初始化为符号距离函数, 这一操作计算量非常大, 为解决该问题, Li, Xu和Fox (简称LXF模型) 等将距离约束信息加入到主动轮廓模型的水平集能量泛函中[8], 这就无需对水平集函数重新初始化就可以驱使水平集函数接近一个符号距离函数, 其方法有很大的优势, 可以节省很多计算时间。

在建立模型前首先引入如下的边缘检测函数来驱使零水平集向物体边界靠拢:

该模型由其内部能量项的约束, 在水平集演化过程中就无需对其重新初始化, 在数值计算上也得到简化。 且迭代步长也允许取的较大, 就加快了曲线演化。

3图像分割方法

水平集算法的核心是解决水平集进化的问题。通过分析比较目前已有的水平集模型, 发现主要存在有两大问题:

(1) 绝大部分水平集模型均涉及到需要在进化过程中对水平集函数重新初始化的问题;

(2) 许多模型都需要借助原始图像的梯度信息来保证水平集进化如期收敛。

采用LXF模型能够较好的解决上述问题, 但是其计算量太大。为了加快计算速度, 本文采用的方法是:先利用某些算法对图像进行粗分割, 得到初始轮廓, 然后以粗分结果作为水平集的初始轮廓, 这样就可以大大降低迭代次数。

4实验结果和分析

4.1粗分割实验过程

粗分割实验如图2所示, 过程如下:

(1) 将原图像转换成灰度图, 并调节亮度。

(2) 屏蔽上半部分, 消除共轭像。

(3) 进行均值滤波。

(4) 先进行膨胀处理, 再腐蚀, 消除微小斑点和空洞。

(5) 利用Sobel算子进行边缘提取。

4.2利用水平集算法进行细分割

以初处理轮廓图作为水平集演化的初始轮廓, 进行基于水平集的细分割。图片大小为534 × 436, 实验采用Matlab 7.6完成。实验结果如图3所示。

以粗分割结果为初始轮廓, 进行水平集迭代处理, 最后迭代次数为144次, 所需时长为1 242.04 s。

5结语

本文分析了全息再现图的图像特点和基于水平集算法分割的算子, 提出以粗分割和水平集细分割相结合的分割方法。实验结果证明:这种分割方法对于全息再现图的处理有较好的效果。

摘要：图像分割是从图像分析到图像处理的关键步骤, 因此对该全息再现图的分割具有重要意义。传统的相移数字全息技术记录次数多、记录时的曝光次数较多, 因而随机误差大。项目组提出了一种基于参考光估计的两步相移数字全息术, 只需记录两幅相移全息图就可实现原物光波的重建。选取了水平集算法进行分割处理, 鉴于水平集算法特点, 为了提高分割精确度和效率, 提出了一种粗分割与细分割相结合的方法。实验证明, 该方法对于两步相移全息再现图的分割具有较好效果。

关键词：相移数字全息技术,全息图,图像分割,水平集,粗分割

参考文献

[1]GABOR D.A new microscope principle[J].Nature, 1948, 161:777-778.

[2]YAMAGUCHI I, ZHANG T.Phase-shifting digital holography[J].Optics Letters, 1997, 22 (16) :1268-1270.

[3]章毓晋.图像处理与分析 (图像工程上册) [M].北京:清华大学出版社, 1999.

[4]熊秉衡, 李俊昌.全息干涉计量:原理和方法[M].北京:科学出版社, 2009.

[5]钟丽云.数字全息的基本问题分析及实现方法研究[D].天津:天津大学, 2004.

[6]OSHER S, SETHIAN J A.Fronts propagating with curvaturedependent speed:algorithms based on Hamilton-Jacobi formulations[J].Journal of Comput Phys, 1988, 79:12-49.

[7]JAVIDI Bahram, OKANO Fumio, SON Jung-young.Three-dimensional imaging, visualization, and display[M].Germany:Springer, 2010.

相移算法 篇7

能源危机问题极大地推动了智能电网、分布式发供电系统的发展。其中,逆变器并联技术作为连接分布式电源与微电网的一个关键性的接口技术,越来越得到国内外学术界与工业界的重视[1,2]。

如何抑制环流从而使各逆变器单元达到功率平均分配或是按一定比例分配,是逆变器并联技术中的一个难点。目前主要通过2类方案来实现逆变单元间的功率分配:一类方案通过检测并互通各并联逆变器单元的输出电流、功率等信息而达到并联的效果[3,4,5];另一类方案不需要逆变单元间进行信息交互,每台逆变器只需检测自己的输出电压电流信息,并采用基于电压幅值—频率下垂的控制方法(PQ下垂法)[6,7,8]即可实现整个逆变并联系统的均流。

由于控制算法的日益更新,其复杂程度也越来越大,数字化的控制方案越来越成为主流。上述2类控制方案都需要对逆变器电压电流检测的离散采样值,通过采用一定的数字算法实时地计算逆变器的输出功率。功率检测及计算的精度和延迟时间等参数会影响到并联系统的稳定性[9]。为了提高逆变并联系统的稳定性和可靠性,需要对数字化功率计算方法的稳态、动态以及运算时间和空间消耗等性能提出更高的要求。

本文首先对逆变器并联技术中现有的数字化功率计算方法做了简要的介绍与总结。在分析了现有方法的优缺点的基础上,提出了改进的基于滑动平均的电压电流相移的数字化功率计算方法。该方法与现有方法的仿真结果对比表明,其具备良好的综合性能。最后通过2台2 kW单相逆变器无互联线并联实验,证实了本文提出的功率计算方法的可行性和有效性。

1 现有数字化功率计算方法介绍

这里只对现有数字化功率计算方法进行最简单的介绍,详细的原理分析及公式推导参见附录A。

1.1 傅里叶分解法

傅里叶分解法[9]的理论基础是周期函数的傅里叶级数理论。通过将周期函数的功率表达式分解为相应的傅里叶级数,然后采用快速傅里叶变换等算法对其系数进行计算,从而最终实现功率的计算。由于计算过程中需要提供相应的正弦表和余弦表,所以傅里叶分解法又称为双表法。

1.1.1 传统傅里叶分解法

传统傅里叶分解法一个工频周期只计算一次功率值,功率计算的延时较大,容易导致逆变器并联系统的不稳定。

1.1.2 基于滑动窗的傅里叶分解法

引入滑动窗的思想,在每个采样周期都计算一次傅里叶系数和功率值,可使零阶保持环节的周期降为仅仅只有一个采样周期。

利用傅里叶分解法计算功率,优点是对高频噪声、直流偏置以及谐波分量不敏感,能够准确地算出基波功率。但计算一次功率值所需的运算次数及数据存储空间较多,占用了过多的系统资源。

1.2 电压电流相移法

该类方法计算无功功率时需要对电流采样序列进行移相操作,故称为电压电流相移法。低通滤波环节一般有周期平均和数字低通滤波器2种选择。

1.2.1 周期平均法[10]

采用周期平均的滤波方法,程序编写简单,运算时间短,但需额外N/4个数据空间来存储前1/4个周期的电流采样信号。该方法存在与传统傅里叶分解法一样的缺点,即在每个工频周期仅计算一次功率值,对功率变化的动态捕捉速度较慢。

1.2.2 数字低通滤波器法[6]

采用数字低通滤波方法,需要特别关注滤波器带宽参数的折中选择。为了获得平稳的功率,要求数字低通滤波器在工频及更高频处的增益尽可能低。但同时为了使并联系统更加稳定可靠,又应尽可能增加低通滤波器的带宽。虽然增加滤波器的阶数可以在提高带宽的同时降低高频增益,但高阶的数字滤波器会增加程序的内存和运算时间的消耗。

1.3 二阶广义积分法

二阶广义积分器(second order generalized integrator,SOGI)在指定频率处可实现极大的增益,能够滤除正弦波形中的杂波,同时具备输入输出同步无延时的优点。SOGI近来被成功应用于并网逆变器系统的锁相环中[11],还被应用在逆变器并联系统中进行功率的实时计算[12]。

利用SOGI滤波器的2个正交输出,相比于电压电流相移法可以节省N/4个采样点的存储空间。因为其在工频处滤波前后的相移为0°,故而跟踪功率动态变化的速度很快。但由于计算一次功率值需要进行4次SOGI滤波运算,降低了功率计算的速度。而且由于SOGI的余弦输出不能滤除直流偏置,当电压电流采样信号中具有直流偏置时,会导致功率计算不够准确。可以通过取周期平均的方法来消除直流偏置量,但会进一步增加运算时间。

1.4 三要素法

三要素法[13]的速度非常快,功率计算的滞后只有仅仅一个采样周期,同时不需要额外的存储空间,占用数控程序的时间和空间的资源都很小。但实际的测量信号中往往存在噪声,哪怕是幅值很细微的一点扰动,都可能导致三要素法的计算结果存在较大误差。另外,该方法忽略了正弦信号的第4个要素——直流偏置量。当电压电流存在直流偏置时,该方法的计算结果也不准确。

2 基于滑动平均的电压电流相移法

2.1 滑动平均思想的引入

对于逆变器并联系统来说,输出电压在一般情况下变化是较小的,但输出电流则可能会随着负载的突变、其他逆变单元的并入或退出等原因瞬间产生较大幅度的变化。传统相移法中,利用当前的电压采样信号与1/4个周期前的电流采样信号来计算功率值,会导致不能快速反映系统中实际无功功率的变化情况。

将电压信号超前1/4个周期后得:

$v^{\prime }(t)=\sqrt{2}V\sin \left(\omega t+\frac{\text{π}}{2}\right)(1)$

同时保持电流不变,此时的瞬时有功功率为:

$\begin{array}{l}{p}^{″}(t)=v^{\prime }(t)i(t)=-V{\rm I}\sin \theta +\\ V{\rm I}\cos \left(2\omega t-\theta -\frac{\text{π}}{2}\right)(2)\end{array}$

将式(2)中的直流分量取反即为无功功率。相应的离散表达式为:

$p^{″}(k)=v^{\prime }(k)i(k)=v\left(k-\frac{{\rm N}}{4}\right)i(k)(3)$

相比其他方法,采用周期平均的电压电流相移法的算法简单,计算时间很短,所需额外的N/4个数据存储空间对于数字处理芯片来说也不构成很大的负担。唯一较大的缺点就是一个工频周期计算一次功率值,这对于逆变器并联系统来说是很难接受的。若能借鉴基于滑动窗的傅里叶分解方法中在滑动窗内进行累加运算的思想,将滑动平均的概念应用于电压电流相移法中,则能得到不错的效果。

构建一时间长度为工频周期(即点数为N)的滑动窗口,并在此窗口中对瞬时功率序列进行平均操作,即从瞬时功率序列的第k+1-N个点一直累加到第k个点,然后除以N,即可得到第k个点的有功功率和无功功率序列。对应的数学表达式为:

${\rm P}(k)=\frac{1}{{\rm N}}{\displaystyle \sum _{n=k+1-{\rm N}}^{k}p}(n)=\frac{1}{{\rm N}}{\displaystyle \sum _{n=k+1-{\rm N}}^{k}v}(n)i(n)(4)$

$\begin{array}{l}Q(k)=-\frac{1}{{\rm N}}{\displaystyle \sum _{n=k+1-{\rm N}}^k{p}^{″}}(n)=\\ -\frac{1}{{\rm N}}{\displaystyle \sum _{n=k+1-{\rm N}}^{k}v}(n-\frac{{\rm N}}{4})i(n)(5)\end{array}$

采用滑动平均后,可在每个采样周期计算一次功率值。相对于传统的电压电流相移法,减小了功率计算延迟时间,改善了逆变器并联系统的动态特性。

2.2 滑动平均法的本质分析

本质上,滑动平均就是一种低通滤波运算,但其与一般的数字滤波器又存在着区别。

根据有功功率计算式(4),可得到Z域内的有功功率Z变换与瞬时功率Z变换的离散传递函数HP(Z)为:

$\begin{array}{l}{\rm H}{}_{{\rm P}}({\rm Z})=\frac{{\rm P}({\rm Z})}{p({\rm Z})}=\frac{1}{{\rm N}}(1+{\rm Z}{}^{-1}+{\rm Z}{}^{-2}+\cdots +\\ {\rm Z}{}^{-({\rm N}-1)})=\frac{1}{{\rm N}}\left(\frac{1-{\rm Z}{}^{-{\rm N}}}{1-{\rm Z}{}^{-1}}\right)(6)\end{array}$

其频率响应为:

${\rm H}{}_{{\rm P}}(\text{e}{}^{\text{j}\omega {\rm T}{}_{\text{s}}})=\frac{1}{{\rm N}}\left(\frac{1-\text{e}{}^{-\text{j}{\rm N}\omega {\rm T}{}_{\text{s}}}}{1-\text{e}{}^{-\text{j}\omega {\rm T}{}_{\text{s}}}}\right)(7)$

式中:Ts为采样周期。

可得频率特性如图1所示。图中:Ts=50 μs,N=400。

如图1所示,在ω=140 rad/s处,传递函数的增益约等于0.707,可知该参数条件下的滑动平均的带宽为140 rad/s。图1同时还给出了带宽都为140 rad/s的一阶及二阶巴特沃兹数字滤波器的频率响应曲线。通过对比可以看出:3种滤波器的传递函数在直流处的增益都为1,在带宽频率内也比较接近;同为巴特沃兹滤波器,二阶的高频增益明显要小于一阶;同巴特沃兹滤波器不同,随着频率的增加,滑动平均的增益呈现以工频为周期的一种变化,且在工频的整数倍处增益都为0。这一特点能够很好地滤除瞬时功率中由谐波带来的整数倍工频分量,该增益特性是一般数字滤波器所不具备的。

3 仿真及对比分析

根据上节对基于滑动平均的电压电流相移法的分析,将其与基于二阶巴特沃兹数字滤波的相移法、基于滑动窗的傅里叶分解法、基于SOGI的方法和三要素法进行对比仿真。

仿真参数如下:工频周期T=0.02 s,采样周期Ts=50 μs,即N=400;巴特沃兹滤波器为二阶,带宽取140 rad/s;SOGI的带宽增益常数设为ξ=0.707。下面分别从稳态精度、噪声抑制、谐波抑制、直流抑制、动态特性等方面作细致的比较分析。这里只提供了有功功率的计算结果,完整的仿真结果参见附录B。

3.1 理想情况

理想情况下,电压电流皆为纯正弦信号,无直流分量,取

v(t)=311sin 314t (8)

$i(t)=10\sin \left(314t+\frac{\text{π}}{4}\right)(9)$

图2为相应的仿真结果。从图2可以看出,除了巴特沃兹滤波的功率计算结果存在振荡,其他方法在准确性和稳定性方面的表现都很好。尤其是三要素法几乎是瞬间就得出了准确的功率计算结果,并且能保持稳定。SOGI法的计算速度也非常快,大约只用了1/2个工频周期即得出了准确结果。巴特沃兹法、滑动傅里叶分解法以及滑动平均这3种方法的功率计算时延相差不多。理想情况下,从功率计算的效果来讲,三要素法性能最佳。

3.2 谐波抑制情况

逆变器并联系统中,由于器件的非线性、死区效应以及非线性负载等因素的影响,电压电流信号往往带有一定量的谐波分量。考虑谐波后,相应的功率计算仿真结果如图3所示。此时的仿真条件变为:

$v(t)=311\sin \omega t+10\sin \left(2\omega t+\frac{\text{π}}{4}\right)+$

$5\sin \left(3\omega t+\frac{\text{π}}{8}\right)(10)$

从图3可以看出,各方法的计算速度同理想情况下一样,但稳态时的功率计算结果却发生了各不相同的变化。滑动傅里叶分解法和滑动平均法的计算结果还是依旧平稳,但其余方法的计算结果中都带有了不同程度的谐波频率的振荡。

3.3 噪声抑制情况

不论采用何种检测方法,电压电流的采样信号中总会含有一定程度的高频噪声。为了对比分析各种方法对噪声的抑制能力,仿真条件改为在电流信号中加入峰—峰值不大于0.1 A的噪声,即

式中:rand(t)为取值在区间(0,1)上的随机数。

仿真结果如图4所示。在电流信号中加入噪声后,三要素法的功率计算结果出现了发散的现象,而且幅度相当大。而其余各种方法除了SOGI法稍有一点点毛刺之外,都能较好地抑制噪声。

3.4 直流偏置抑制情况

由于逆变器的正负电压不对称、采样调理电路的参数随温度漂移等原因,采样信号中也可能会存在直流偏置。信号中的直流分量对各种功率计算方法的影响,可以通过相应的仿真来比较。仿真条件为5 ms时在电流信号中加入2 A的直流偏置量,同时电压不变。

·广告·

此时的电压和电流分别表示为:

仿真结果如图5所示。在电流信号中加入直流偏置的瞬间,三要素法的计算结果中存在很大的尖刺,SOGI法也在一定程度上存在该过冲现象。随后的计算结果中,三要素法和SOGI法都出现了较大幅值的工频分量。巴特沃兹滤波法的结果中也存在类似现象。只有滑动傅里叶分解法和滑动平均法能够完全抑制电流信号中直流偏置的影响。

3.5 信号突变情况

逆变器并联系统运行中,当发生负载突变、其他逆变单元并入或切出等情况时,电流可能会出现瞬间突变的现象。为了模拟该类情况下各种方法的功率计算结果的变化,仿真条件设为2.5 ms时将电流幅值瞬间变为原来的2倍,保持电压不变,即

仿真结果如图6所示。由于电压电流信号在突变前后都是理想的正弦波,故各种方法的功率计算结果都很平稳。

在电流信号突变的瞬间,三要素法仍然能够以最短的时间捕捉到功率的变化,SOGI法紧随其后,但同时这2种方法存在较大的过冲。其余方法跟踪这一功率突变的时间相对较长,大致为一个工频周期。

3.6 功率计算方法的总结

根据仿真结果,结合各方法所需要的运算时间和数据空间情况,各种功率计算方法的总结如表1所示。

注:“m乘n加”表示采用该方法计算一次功率值需要执行m次乘法和n次加法运算。

综合考虑各种方法的稳态和动态性能,尽管滑动傅里叶分解法和滑动平均法这2种方法占用的运算时间和存储空间相对较多,但功率计算效果较为出色。相对于滑动傅里叶分解法,本文提出的基于滑动平均的电压电流相移法由于需要更少的运算次数和存储空间,更加适合于逆变器并联系统中实时功率计算的应用。

4 实验结果

为了进一步验证文中所提出方法的有效性,搭建了2台2 kW单相逆变器无互联线并联系统。数字控制芯片采用TI公司的型号为TMS320F28335的信号处理器。并联控制方法为PQ下垂法,其中的功率计算环节采用本文提出的基于滑动平均的电压电流相移法。输出空载电压有效值为221 V,空载频率为50.1 Hz,采样频率及开关频率为输出电压频率的400倍。

并联系统运行时的稳态均流结果如图7所示。可以看出,采用滑动平均的功率计算方法后,系统的稳态均流效果较理想。

图8为逆变器2稳定运行时逆变器1并入瞬间的动态均流实验结果。并入后经过大约一个工频周期,系统进入稳定均流状态。基于滑动平均的电压电流相移的功率计算方法的可行性和有效性得到验证。

5 结语

本文简要综述了逆变器并联系统中的数字化功率计算方法。在分析总结现有方法优缺点的基础上,借鉴滑动窗口的思想,提出了基于滑动平均的电压电流相移的数字化功率计算方法。从稳态、动态性能以及运算的时间和空间消耗等方面综合考虑,与现有方法的仿真结果对比验证了该方法良好的综合性能。同时逆变器无互联线并联实验系统的稳态和动态的均流效果证实了文中提出的功率计算方法的可行性和有效性。