差分方程及Z变换工程实例

工程实例

机械电子工程 谈卓雅

一、基于时域有限差分方法求解薛定谔方程 在量子力学理论中一维时域薛定谔方程的具体形式如下:

()()()()t x x V x t x m h t t x jh ,,2,222ψψψ+??-=?? (1)

式中,参数h-为普朗克常数；m 为粒子质量；()y x ,ψ定义为状态变量;()x V 为势函数。在求解含时薛定谔方程时,本文将利用时域数值计算方法——时域有限差分(finite difference time domain, FDTD)来差分离散薛定谔方程。其基本思想是利用微分方程的中心差分离散形式建立空间和时间上的迭代。时域有限差分法与传统的量子力学计算方法相比更加直观、计算效率更高、操作性更强。因此,该数值方法在量子力学等其他领域中的应用亟待进一步发展与完善。在运用FDTD 方法时,数值稳定性条件是首要考虑因素之一,它直接影响着数值计算结果的精度和有效性。本文主要是分析用FDTD 方法求解薛定谔方程时所需要满足的稳定性要求,又进一步提出了在不同势能情况下从一维到三维的统一的数值稳定性表达方式。

1 薛定谔方程的基本形式

将式(1)改写为：

()()()()t x x V h

j x t x m h j t t x ,,2,22ψψψ-??=?? （2）

为了便于计算,将()y x ,ψ复函数的实部和虚部分别考虑： ()()()()t x x V h x t x m h j t t x R R I ,1,2,2

2ψψψ-??=?? （3a ） ()()()()t x x V h x

t x m h j t t x I I R ,1,2,22ψψψ+??=?? （3b ）

FDTD 方法在时间上的差分格式如下所示： ()

()()()()I or R t

t n x t n x t t x t n t =??-?+=???+=ξψψψξξξ,,1,,2/12 (4)

及空间上的差分格式为： ()

()()()()()()[]t n k x t n xk t n k x x x t x t n t x k x ?-?+??-?+??=???=?=,1,2,11,222ξ

ξξξψψψψ (5) 同时记()()k t n x k n ψ

ψ???,,将式(4)与式(5)代入式(3)后化简为： ()()()()()()[]

()()k k V h t k k k x t m h k k n I n I n I n I n R n R 2/12/12/12/1211212+++++?+-+-+??-=ψψψψψψ (6a) ()()()()()()[]()()k k V h t k k k x t m h k k n R n R n R n R n I n I 2/12/12/12/1211212+++++?--+-+??+

=ψψψψψψ (6b)

由此,利用时间步步的迭代可以算出每时刻的函数值。然而,在实际操作时需要使空间步长Δx 及时间步长Δt 满足一定的关系以达到数值稳定的目的。

2 数值稳定性分析

FDTD 方法是实现用一组差分方程的解来代替时域薛定谔微分方程的解。只有当差分方程的解是稳定和收敛的,这种代替才有意义。稳定是指寻求一种离散间隔所满足的条件,在此条件下差分方程的数值解与严格解之间的差为有界量即为收敛。因此,FDTD 方法中空间步长Δx 与时间步长Δt 不是相互独立的,以避免数值模拟结果的不稳定性。

首先,将ψ函数表示为一维空间变量和时间变量函数的乘积,即()()()t x y x ?φψ=,,考虑函数随时谐变化：

()()()t j x t x ωφψex p , （7）

其中的稳态解是下面函数对时间进行一阶偏微分的方程解：ωψψj t

=?? （8） 将式(8)进行中心差分离散,则 n n n j t ωψψψ=?--+2

/12/1 （9）

式中，()t n x n

?=,ψψ。当t ?足够小时，定义数值增长因子2/12/1-+==n n

n n q ψψψψ 化简式（9）为： 012=-?-tq j q ω

（10） 则式（10）的解为： 2212??

?????-±?=t t j q ωω （11）

根据冯·诺依曼的数值稳定性要求,当时间步趋于无穷大,即Δt 足够小时,q 必须满足|q|≤1，即 12≤?t

ω （12）

此为方程解满足稳定的要求。

其次,在量子力学中的一维时间变量的函数值()()t h E j t j e e

t /--==ω?代入薛定谔方程

中，简化得到： ()()()()02222=-+??x x V E h m x x φφ （13）

采用平面波本征模思想，一维空间变量的函数值为()()x jk x x -=ex p ?，将其代入下面的二阶精度差分近似

()()()()φφφφφ

2222

222sin 2??

????????????-=??-?+?-?+?=???=x x k x x x k x k x x k x x x k x （14） 则得到 ()()0222sin 222=?--??

????????????φφx k V E h m x x k x （15） 且在量子力学中()h h f fh E ω=??

????∏∏==22，则简化式（15）为 2

222sin 422t h t V x x k m t h x ?=?+??

?????????????ω （16）

结合上述的稳定性要求式（16），可得： 1222sin 422≤?+??

?????????????h

t V x x k m t h x （17） 考虑到12sin 2≤??

?????x k x ，因此得到最终的稳定性条件为： ()22max 2≤?+??h t V x m t

h （18）

在取不同势能V 值时，稳定性条件也发生变化，下面分情况讨论：

当0=V 时，稳定性条件为 ()

222≤??x m t

h （19） 当0≠V 时，稳定性条件为 ()22max 2≤?+??h t V x m t

h （20）

且其中()x k V V ?=m ax max 。

本文所采用的稳定性分析方法在时间、空间离散格式上是独立的,所以稳定性条件可以直接推广到三维：

()()()21112max 222≤?+???????+?+??h

t V z y x m t h （21）

二、 基于Z 变换及模糊加权均值滤波的匀速运动模糊图像恢复

传感器成像过程中,记录介质一次积分时间内拍摄目标和摄像机之间的相对运动会造成图像的模糊,给后续运动图像的处理分析带来困难(比如在捕获)。运动模糊图像常用恢复方法有基于时域的差分恢复法及基于频域的逆滤波、维纳滤波、约束去卷积法等。但这些方法都相对比较复杂,一般需要进行大量的矩阵运算,运算速度较慢,且可能由于噪声、近似取值等原因,使得图像质量进一步降低。文献提出了一种基于Z 变换的运动模糊图像快速恢复算法,该算法能在保证图像恢复质量的同时,大幅度提高图像恢复运算的速度,基本达到实时的要求。但其最大的缺点是算法的“迭代特性”使得前面像素点的噪声分量会在图像恢复过程中不断放大、累加,影响后续恢复的像素点,使复原图像质量大为下降。为此,有必要在图像恢复过程中,实时地对复原像素点进行修正或平滑滤波处理。常见的图像平滑滤波方法有空间域方法如邻域平均法、中值滤波法和频率域方法如低通滤波法等。但低通滤波法适合对整个图像进行处理,而本文需要实时地对复原点

进行滤波处理,因此,较适合应用空间域方法处理;中值滤波法对脉冲噪声有较好的抑制作用,且能较好保留边缘信息,但它的平滑能力较差,这大大限制了其在图像降噪中的应用。均值法对高斯噪声有较好的滤波能力,但由于:(1)在模糊图像恢复过程中产生的噪声点与有效点的差值较大,采用相同权值的均值滤波会造成图像特征边缘的失真。(2)在均值运算中,各个点的权值都一样,当滤波区间中存在奇异点(脉冲噪声)时,奇异点在很大程度上影响滤波效果,且奇异点的存在经均值滤波还会影响、扩散到周围点,故均值滤波对脉冲噪声十分敏感,而恢复过程中的噪声正好符合脉冲噪声的特征。(3)均值滤波没有充分利用邻域内测点间的相关性和测点的位置信息。因此,本文在严格数学推导下,提出一种基于Z 变换及模糊加权均值滤波的匀速运动模糊图像恢复算法,该算法在保证恢复速度的同时,可消除噪声扩散,提高复原图像的质量。

1.基于Z 变换的模糊图像恢复算法

（1）退化原理及过程

在运动检测过程中,假设运动目标的运动方向与像平面始终平行,一帧理想的清晰图像为原图像()y x f ,,()y x f ,)经过一个退化过程()y x h , (如光学系统的点扩散函数PSF)得到退化图像()y x g ,,这一过程可描述为:

()()()y x f y x h y x g ,,,*= (1)

式中,(*)表示卷积,为了简化运动模糊恢复问题,这里暂不考虑噪声的影响。

2、基于Z 变换的退化及恢复模型

（1）退化模型

设图像()y x f ,在像平面内作匀速直线运动,运动方向与X 轴方向平行(对于旋转匀速运

动,可按文献[4]的方法转换成匀速直线运动;对于其它方向的匀速运动,可通过三角公式旋转变换成X 轴的运动),速度为v,曝光时间为T,则成像后可得到模糊图像g(x,y):

()()()?-=T dt y t v x f T

y x g 0,1, （2） 设每帧图像的总曝光时间T 内运动对象移动N 个采样点,则有:N=v ×T,图像经过采样转化为离散图像g(n):

()()()()∑-=-=*=1

01N i i n f N n f n h n g （3）

上式的物理意义为模糊路径上任一像素点模糊后的灰度值是该像素点及其前面N-1个像素点的原灰度值的加权累积,N 为模糊宽度,用像素数来表示,并且有N ≥1。设G(Z)及F(Z)分别为g(n)和f(n)的Z 变换,对式(4)进行Z 变换后得：

()()∑-=-=1

01N i i z f N

z G （4） 由Z 变换的移位性质得: ()()()()()()

1211+---++++=N z z F z z F z z F z F N z G Λ （5） 对式(5)进一步整理得到运动模糊图像的退化模型为:

()()()()()∞≤<--===---=-∑z z N z z N z F z G z H N

N k k 0111110 （6）

退化模型H(Z)及F(Z)本身均存在零点,这将导致系统产生振荡,形成病态系统,所以应选择H(Z)的收敛域,避开其零点的影响,使系统得到满意的解。

（2）恢复模型

在假设条件不变的情况下,图像的恢复是退化的逆过程。式(6)可整理为:

()()()()()z F z z G z z G N z F N --+-?=1

（7）

对式(7)进行Z 逆变换得： ()()()()()N n f n g n g N n f -+--?=1 （8）

从式(8)可以看出,本文的恢复模型只需要进行简单代数运算即可得到复原图像,可实现运动模糊图像的实时恢复,较其它算法具有一定的优越性。

由于式(8)为一个递推表达式,如果像素点f(n-N)存在着恢复误差或噪声(H(Z)及F(Z)由于本身均存在零点,因此容易出现畸变点),则f(n)也会受到影响,且噪声将进一步扩散到

()()N

,

,2,1,0

≤

+Λ,width为模糊图像宽度。因此,有必要在f(n) =

,

'-

width

n

m

m

n

mN

f/

恢复的同时,对其输出值进行估计或者滤波,以减小其对后续输出点的计算误差。

在滤波算法的选择方面,我们考虑到噪声主要表现为局部极个别的畸变点(当然,如果不处理就会影响后续需要恢复的点),可以将图像转换到频域进行处理,将高频部分去掉,再反变换到时域,但针对每个窗口都进行频域处理相对还是较为复杂的,因此,这里考虑在时域对其进行处理。均值滤波算法是常用的时域算法,但容易造成过平滑或图像失真问题,而加权均值的滤波算法在保持原来有效信号的同时,能对畸变点有较好的滤波作用,因此,本文提出在恢复运算进行的同时结合模糊加权均值滤波算子对恢复点进行降噪处理,以改善复原图像的质量。