1.4 现代光学理论
1.4.1 光学传递函数
在望远镜分辨率的讨论中已经介绍了口径场点分布函数的概念。现代光学理论中的另一个重要概念就是调制传递函数(modulation transfer function)。实际上调制传递函数是光学传递函数中模的部分,光学传递函数的另一部分是相位传递函数(phase transfer function),相位传递函数是它的相位部分。如果相位传递函数为零,则光学传递函数就等于调制传递函数。
传递函数一般用来描述线性时不变系统(linear time-invariant system,LTI)的特性。所谓线性是指系统的输入和输出量满足叠加定理。也就是说如果某一个输入量是两个子输入量的和,那么系统在这个输入量的作用下所产生的输出量也将是该系统在这两个子输入量的作用下所分别获得的两个输出量的和。所谓时不变是指如果输入量存在一个时间延迟,所产生的输出量将与原来所产生的输出量完全相同,并且它也同样在时间上有一个相同的延迟。在这样线性的、时间不变的系统中,如果输入量是一个时间极短、同时具有单位面积的脉冲函数,那么它的输出则是系统的脉冲响应函数。系统的传递函数是脉冲响应函数的拉普拉斯变换。系统的输出量可以用脉冲响应函数与输入量函数在时间域的卷积来获得。在拉普拉斯空间,输出量就是传递函数和输入量的乘积。
对一个稳态线性时不变系统,如果它的输入或输出量均可以用复数周期矢量,即具有频率、振幅和相位的相位子(phasor)来表示,那么这个系统的传递函数就可以用输出和输入量的傅立叶变换的比值表示。傅立叶变换是拉普拉斯变换的一个特殊情况,这时原来拉普拉斯变换中的实数部分消失,只保留了虚数部分。也就是去除了它的时间变化部分,而只保留了系统稳态响应。这时传递函数常常被称为系统频率响应。系统频率响应包括两个频率域的分量,一是振幅,另一个是相位。
光学系统常常被看作是线性、空间不变的、非相干(incoherent)的系统。所谓空间不变是和时间不变相对应的,不过这时的输入和输出是发生在空间域内的现象。即当物在空间域内的位置发生变化后,它的像的特点并没有发生任何变化,但是像的空间位置却随着物的空间位置的变化而发生相应的变化。所谓非相干则是指从物空间光源的任何两个不同发光点上所发出的光之间并不存在稳定的相位关系,所以它们之间不能产生干涉,而只有从光源上同一个发光点所发出的光才能够产生干涉。在这种条件下,系统的光学传递函数就是输出和输入量的傅立叶变换的比值。而这个傅立叶平面是在空间频率域上,而不是在时间频率域上。和非相干光相对的是相干(coherent)光。所谓相干包括空间相干(spatially coherence)和时间相干(temporal coherence)。如果一个面光源,在这个光源的一个区域内的所有点上发出的光的波阵面均非常精确地同样进行振动,它们的频率和相位完完全全是相同的,那么我们就称这个光源是空间相干的(spatially coherent)。如果光源不但有如上特性,而且还具有非常稳定的频率,以至于经过很长的一段时间以后,它们的相对相位仍然是完全不变的,那么我们就称它们还具有时间相干性。一般说来,相干的波之间可以发生干涉,而非相干的波之间则不会产生干涉。而介于两者之间的则可能在一定程度上产生干涉,这个干涉程度被称为相干度。
和时间频率类似,空间频率是用来描写一个函数在空间尺度范围内产生重复次数的情况。空间频率可以用来描述一个平面物,一个平面光源,或者一个平面像的细节。因为任何物或像均可以看作是由无穷多个空间频率的不同强度分布的集合。如果p0是某一个物或者像上强度变化的周期,那么相应的空间频率ν0就是单位长度和这个周期的比值。空间频率的单位是长度的倒数。在天文上长度常常是角度的正弦值,而在小角度时,正弦值和角度值相等,所以空间频率的单位也可以是角度的倒数。
光学调制传递函数和对比度(contrast)的变化或调制相关。在空间某个频率上,对比度的定义是C(ν0)=(lmax-lmin)/(lmax+lmin)(图1.49(a))。而光学调制传递函数就是像和物在空间频率上输入和输出值的对比度的比值:
图1.49 空间频率和空间频率的对比度
如果光学系统是一个理想系统,则对于所有的空间频率调制传递函数的值MTF(ν)≡1。但对于没有光强放大机制的实际成像的光学系统,则永远有MTF(ν)<1(图1.49(b))。只有当ν→0时,MTF(ν)→1。而当ν趋向于系统的分辨极限时,则MTF(ν)→0。这个分辨极限所对应的空间频率,我们就叫做该系统的截止频率νc。在光学成像的过程中,空间频率高于截止频率νc的所有物的信息将全部消失。因为在像空间不存在任何高于截止频率的频率,所以在很多情况下常常使用归一化的空间频率,即νn=ν/νc。
如果一个系统,它除了会改变空间频率对比度以外,还会使空间频率上的明暗条纹产生横向的位置移动δ(图1.49(c)),则其相位传递函数就不为零,即:
则该系统的光学传递函数为:
在信息和图像的傅立叶变换,即频率谱的表述中,频率谱的振幅和相位有着不同的地位和作用。在很多情况下,信息中最重要的特点并不是由频率谱的振幅值来保存的,而是由它的相位值来保存的。很多例子表明,单纯利用频率谱的振幅,很难获得原来信息的基本特点,而单纯通过频率谱的相位,却可以获得原来信息的很多重要特点。研究表明如果要恢复原来信息的均方根误差值,相位在数字化过程中应该比振幅的位数平均多1.37到2位。从相关性方面来说,相位值和原来信息的相关度十分接近,而振幅值和原来信息的相关度则差距很大。如果原来的信息是一幅两维图像,那么在原来信息的傅立叶变换中,几乎所有事件的位置均反映在它的频率谱的相位部分。具体在图像上,这些事件就是一些线条,一些点,或者其他局部范围的图像特征。这些事件的空间位置上的变化,在频率谱的振幅图中没有任何显示。而在相位上,就表示为一个固定的相位增加量。同时在振幅上,一般图像会在某一个高频率有截止的倾向,而在相位上则重点显示了高频的一些特征,如点、线或者其他的图像轮廓。在声音信息中,如果是非常短的声音,相位的重要性比较小,而对于较长的声音,它的相位中就保存了原来声音中的最主要的信息(Oppenheim and Lim,1981)。
在线性时不变系统中,脉冲响应函数和传递函数是拉普拉斯变换对。在光学系统中,点分布函数和光学传递函数是傅立叶变换对。一个点光源就相当于在时间域内的一个脉冲函数。总的来说,在光学系统中,物、像、点分布函数和光学传递函数之间总存在着下列的关系:
一个与点分布函数相关的参数是点分布函数在一定角半径内的能量(encircled energy)。要求解这个参数,必须首先计算点分布函数的总能量,然后用角半径内的能量和这个总能量来相除。
如果用口径场复数函数Ρ(x,y)来求解光学传递函数,则光学传递函数是口径场函数在口径范围内的自相关积分,即(图1.50)
图1.50 口径场函数自相关积分的坐标系统
式中在X,Y坐标上口径场Ρ的中心与它的共轭口径场Ρ*的中心距为s。它们的中心点的连线与y轴的交角为α。这时的面积分是在相互重合的区域内进行的。在光学传递函数的表达式中k=(s/2)sinα,j=(s/2)cosα。当s=1时,取截止频率的值νn=2s/λ。
为了确定在光学传递函数中是否存在相位传递函数,有下面的一些判断准则:(1)如果在口径场上任何点的相位均为零,则不存在相位传递函数;(2)如果口径场函数是厄密(Hermitian)函数,即它的复数共轭函数等于它变量符号改变了以后的原函数,那么它的相位传递函数将不存在;(3)如果口径场的相位函数和振幅函数都是偶函数,那么一定存在相位传递函数。总结一下现代光学中有关的重要函数的概念,我们可以看出它们之间有下列的关系,即
在上面四种函数中两两互为傅立叶变换对,但是它们之间的差别是十分明显的。以光学传递函数与点分布函数为例,点分布函数是一个实数函数,而光学传递函数或者是实数函数,或者是一个特殊的复数函数,即厄密函数,它的傅立叶变换永远是实数函数。像场分布函数有两个部分,即振幅分布函数和相位分布函数。光学传递函数通常是一闭区间函数,而点分布函数则是一个开区间函数。同样口径场函数也常常是一闭区间函数,而像场分布函数则是一个开区间函数。当用开区间函数来求闭区间函数时,如果仅仅在有限区域积分,则会引入相当的误差。
当我们应用点分布函数来求光学传递函数时,积分过程不可能在无穷边界内进行。这时点分布函数被截止,如果求解的光学传递函数仍然是在ν=0时作归一处理,则在其空间频率上会产生很大的误差。为了减小这一误差,应该用点分布函数所截取的总能量来进行光学传递函数的归一化处理(见图1.51)。当截取半径为3个艾里斑半径时,则截取的总能量为像点实际总能量的0.943倍,这样归一化处理可以使误差大大减小。另外采样频率也会影响光学传递函数在高频部分的误差。一般在每一个艾里斑半径上应该有10至20个采样点,以保证光学传递函数的精度。
图1.51 利用点分布函数求调制传递函数不同归一处理时的误差(虚线是不同的采样距离所引起的误差)(Wetherell,1980)
点分布函数和光学传递函数互为傅立叶变换对。光学传递函数具有对称性,它或为实数函数,这时相位传递函数消失;或为厄密函数,这时调制传递函数是偶函数,而相位传递函数是奇函数。点分布函数总是实数函数。当光学传递函数是厄密函数时,点分布函数不是偶函数,反之则始终为偶函数。只要口径场上的相位W(ρ,θ)=0,则不论口径场边界是如何不规则,所获得的点分布函数就总是一个偶函数。
口径场函数与像场分布函数互为傅立叶变换对。所以可以用像场分布函数来计算口径场函数,这就是射电天线全息检测面板表面误差的方法。在使用这种方法时,他们采用了一种叫奈奎斯特(Nyquist)的采样方法,采样区间比较大,采样频率较低。采样区间可以达到100个主瓣宽度,采样频率大约是半个主瓣宽度。对于圆形口径,调制传递函数具有如下的形式(见图1.52):
图1.52 中心遮挡率为ε的圆形口径的调制传递函数
对于有中心遮挡并且遮挡率为ε的圆形口径,则有:
1.4.2 波像差和调制传递函数
使用调制传递函数的最大优点就是光学系统中总调制传递函数的值是各分部调制传递函数之乘积。即
同样在系统中的,如果有多个互相独立的,影响传递函数的因素,则总调制传递函数等于各因素的调制传递函数之乘积,即
式中的下标d表示口径,f为几何像差,r为随机分布的波差,p为指向误差。用这种方法可以把复杂的波像差分解成相对简单的波形。同样也可以将Td理解为理想光学系统的传递函数,而实际光学系统的传递函数则是这个理想系统传递函数与一系列衰减(degradation)函数的乘积。光学系统的所有几何像差均可以用口径场上波阵面相位差来表示,这样很容易计算它们相应的调制传递函数,并且用图像来表示。图1.53列出了一些几何像差的调制传递函数。图1.54是离焦误差和口径本身所引起的调制传递函数的图形。这时调制传递函数在一些频率上可能取负值,这种现象被称为分辨率畸形(spurious resolution)。分辨率畸形可以看作相位发生180度变化,从而引起正负号的对换,在实际情况下就是空间频率条纹中明暗条纹的位置产生对调。举例来说如果在物空间是五条明条纹夹在四条暗条纹之间,则在像空间就是四条明条纹夹在五条暗条纹之间。
图1.53 一些几何像差的调制传递函数(Wetherell,1980)
图1.54 圆形口径加上离焦相位差后的调制传递函数(Wetherell,1980)
一部分几何像差是旋转对称的,可以只用一条曲线来表示它的调制传递函数。但是另一部分是非旋转对称的,必须用一组曲线来表示它的调制传递函数。图1.55所示是彗差和像散的调制传递函数的曲线。
图1.55 彗差和像散的调制传递函数(Wetherell,1980)
在总波阵面误差中减去所有的初级像差,剩下的就是高空间频率上的随机误差。这种误差在空间频率一般大于五个周期。比如2.4米望远镜在口径上有周期为l毫米形状误差,则这个误差在半径方向上的空间频率则是1200周,这种高频率的空间误差也叫波纹差(ripple)。奥涅耳(O'nell,1963)对这种波纹型误差进行过统计分析,确定的调制传递函数为:
式中k=2π/λ,ωm是空间频率中值上的均方根波阵面误差,C(ν)是归一化后波阵面误差的自相关函数。一般情况下波阵面随机误差的自相关函数是高斯函数,即,νn是归一化的空间频率,l是以镜面半径归一化后的相关长度。l的倒数就是波纹差的空间频率。这种波纹差的传递函数同样被称为衰减函数,其值如图1.56所示。
图1.56 均方根为0.1波长的波纹型误差的衰减函数,l是相关长度(Schroeder,2000)
空间频率更高的随机误差称为微波纹差(micro-ripple)。这时归一化的相关长度l→0,它的传递函数为:
MTFm和MTFh的乘积是所有高频随机误差的传递函数。
另一种波阵面误差是指整个波阵面的随机摆动,这种摆动会引起像的位置变化,而引起随机指向误差。这种误差可以用旋转对称的随机函数表示:
式中σ′是随机函数的标准误差,α为像的平均移动量,均用角度值表示。将σ′归一化有σ=σ′D/λ。所以随机指向误差的传递函数为:
从上式可知,空间频率增大时随机指向误差的传递函数MTFp(νn)的值会减少。随机指向误差的效果是使星像能量分布更均匀,它的高频细节会消失。
在圆口径传递函数的表达式中已经引入了中心遮挡的因素,但仍有必要讨论一下中心遮挡的传递函数(见图1.57)。中心遮挡的传递函数比较复杂,但是在νn≥(1-ε)/2时:
图1.57 中心遮挡的传递函数(Wetherell,1980)
从公式和传递函数的图形中,可以发现中心遮挡的传递函数在高频部分的数值会大于1,这意味着中心遮挡对点分布函数的能量分布和波阵面误差不同,它会提高星像的分辨率。从数学上讲,这时像场的能量分布是两个贝塞尔函数互相比拼的结果,它们中一个是口径的函数,一个是遮挡部分的函数。这样像场能量会从一个亮环转移到另一个亮环,从而使高频部分的细节得到加强。
1.4.3 斯特尔比(Strehl ratio)
斯特尔比是有像差的光学系统中点分布函数最高点能量与相应的理想光学系统点分布函数最高点能量之比。斯特尔比也可以从实际和理想光学系统的调制传递函数中求得。它等于两个调制传递函数曲线下面的面积之比:
式中MTFi(ν)是理想系统的调制传递函数,而MTF(ν)是实际系统的调制传递函数。由于振幅分布函数是点分布函数的平方根,同时是口径场函数的傅立叶变换,所以斯特尔比等于口径场在零空间的频率,即exp[-2πi(νxx+νyy)]=1时的傅立叶变换:
式中W(ρ,θ)是口径场的相位。当相位差很小,斯特尔比可以取上式泰勒级数的前两项:
式中σ是波阵面的均方根相位差。当均方根相位差为λ/14的时候,斯特尔比等于0.8。当波阵面均方根相位差较小(σ<λ/π),就是弱像差的情况,这时点分布函数基本上仍然是在衍射斑的区域内。不过在它的中心点,能量最大值会下降,而像点的角直径会增大(Dalrymple,2002):
式中θD是口径衍射极限时的像斑张角,。利用实际像点的张角值,则有:
当受到大气扰动,波阵面误差大(σ≥λ/π)的时候,光学系统是处在强像差区域内,这时像的中心区就变得十分模糊。由于能量重新分布,这时不存在艾里斑的暗条纹,像斑是由散射能量和噪声形成的。这时的系统调制传递函数等于:
式中MTFi是理想系统的调制传递函数,lz是光程上的相关长度(第2.4.1节),R是和口径直径相关的长度。它的点分布函数是一个高斯函数:
式中r是像平面上的半径,F是系统焦距。它的含有p%能量的像斑角直径为:
含有50%能量的像斑角直径θ50%=3.33σ/lz。这时斯特尔比St≈lz/σ2。波阵面相关长度对斯特尔比有很大影响,作为极限的情况,如果相关长度等于口径尺寸,那么像面将会倾斜。
1.4.4 星像的空间频率
为了简单易懂,本节中部分公式和讨论均假设口径和星像是一维函数。实际上,所有望远镜的口径和像均是二维平面函数,不过二维公式的推导和一维的原理完全相同。读者可以很容易地将一维的公式推广到二维的情况中去。
由于望远镜的夫朗和弗辐射场A(sinφ)与望远镜的口径场P(xλ)互为傅立叶变换对,则在一维情况下有:
式中xλ=x/λ是一维空间频率。因为望远镜口径场是有界的,其夫朗和弗像场A(sinφ)可以用A(φ)来表示,有:
式中aλ=a/λ为口径场的宽度。在光学望远镜中口径场是均匀照明的,在射电望远镜中口径场可能有边缘衰减,因此其远场方向图与光学望远镜中的点分布函数略有不同。
望远镜光学传递函数等于口径场复数函数的自相关(Kraus,1986):
式中OTF(xλ0)是像斑强度的傅立叶变换(图1.58(c)),P(xλ)是口径场分布(图1.58(b))。当位移xλ0大于口径宽度aλ时自相关函数为零。如果天体源的亮度分布为B(φ0),则在空间频率域的像应该是天体源的亮度分布的傅立叶变换与光学传递函数的积:
这个公式表明望远镜的作用实际是一个空间频率滤波器。在天体源中高于截止频率的任何细节均消失,所保存的仅仅是天体源中小于空间截止频率上的光学信息。
图1.58 口径场复数函数的自相关(Kraus,1986)
为了获得天体源中更多的光学信息,可以借助天文干涉仪(见图1.59)。如果在干涉仪中两个口径的距离为sλ,而口径直径为aλ,则这种干涉仪的截止空间频率为:
图1.59 简单的两个口径的干涉仪
空间截止频率越大,保存的信息越多,所以可以分辨出更多的细节。两个口径组成的干涉仪所产生的夫朗和弗辐射场为(见图1.60):
图1.60 两个口径的干涉仪的口径互相关(Kraus,1986)
式中An(φ)是归一化的单口径的夫朗和弗辐射场。两个口径干涉仪的强度分布为(见图1.60(a)):
这个公式表示干涉仪所形成的像斑是一组经过调制以后的条纹(图1.61(c))。调制后的条纹宽度就是干涉仪基线所对应的空间频率sλ的倒数(图1.61(b)),但是它的振幅受到了单个口径点分布函数的调制(图1.61(a))。
图1.61 (a)单个口径;(b)双缝干涉仪和(c)双口径干涉仪的点分布函数(Kraus,1986)
利用干涉仪对某一天体源观察,其像斑应该是干涉仪点分布函数和天体源强度分布的卷积,即(Kraus,1996)
式中α是天体源张角,S0为源的能量密度,φ0为像平面上的角位移。如果天体源张角远小于sλ的倒数,则星像条纹和口径场的衍射强度分布相同。如果天体源张角略小于sλ的倒数,则星像条纹对比度会大大减小。如果天体源张角正好等于sλ的倒数,则星像条纹完全消失(图1.62)。在上一公式中,第一项是常数,而第二项可以表示为:
式中的变量也可以表示成位移Δφ0的余弦函数,有:
式中的V0(sλ)代表星像条纹的振幅又叫做条纹的能见度(visibility),它是基线长度sλ的函数。从上式可以得到下面的重要公式:
从而有:
式中的V0(sλ)exp(j2πsλΔφ0)叫做复数能见度函数。上面公式对于展源也同样适用,因此可以记为:
图1.62 (a)点源,(b)张角小于sλ的倒数的和(c)张角等于sλ的倒数的天体源所形成的干涉条纹(Kraus,1986)
这个公式表示能见度复数函数是天体源强度分布傅立叶变换,反过来也一样,复数能见度函数的反傅立叶变换就是天体源强度分布。
对于对称分布的天体源,条纹位移值为零或者为半个条纹(Δφ0=sλ/2)。
在射电相关干涉仪中,相对于每个基线的复数能见度函数的值都可以一一测量出来,然后通过傅立叶变换去获得天体源的亮度分布。空间频率和望远镜的基线长度成比例,而波长和基线长度的比就是空间角尺度的大小。在空间频率范围内,单口径望远镜在空间频率为零时系统响应是1,频率响应随频率的增加而减小(图1.63)。对于相关干涉仪,频率响应则集中在一个个的带区内,它的响应峰值分别对应于它们单元之间的不同基线的长度。然而当空间频率为零时,它却没有任何响应,所以天体源中大尺度的信息反而丧失了。
图1.63 单镜面(实线)和相关干涉仪(虚线)的空间频率响应(Emerson,2005)
考虑真实的二维分布的情况,能见度函数是在两个方向上尺度的函数。天体源的强度分布的公式为:
式中的u=sλx,v=sλy分别是口径场在x和y方向上的间距,而x=cosα和y=cosβ分别是相对于这两个方向的角度的方向余弦。用这种方法来实现成像的方法在射电天文学中叫口径综合方法。为了获得好的图像,它需要很好的u,v分布。比较之下光学上的口径综合要比射电上困难很多,有关试验还正在开展。
前面讨论的能见度函数是单色光条件下的定义,当频谱有一定宽度时,该函数将受到一个sin(πτΔ(c/λ))/(πτΔ(c/λ))函数的调制,条纹的能见度将迅速下降,式中τ是两信号的时间延迟,Δ(c/λ)是频谱宽度。频谱宽度的影响以及光的空间和时间相干理论将在第五章5.5.1节讨论,十分重要的Weiner-Khinchin和Van Cittert-Zernike理论将在第八章8.3.3节中介绍。
1.4.5 拼合镜面的成像特点
拼合镜面和其他形式的主镜系统将在第二章讨论,本节仅讨论拼合镜面的成像特点。对于一个理想的拼合镜面,它的点分布函数为(Yaitskova,2003):
式中是像面上的位置矢量,是特定六边形镜面单元的局部位置矢量,表示在这个镜面单元内,而表示不在这个镜面单元内,是镜面单元的面积,d是两个相邻单元之间的距离,N是单元编号,λ是波长,f是焦距,PSFseg是单个镜面的点分布函数,而GF是一栅格形函数,它是表示拼合子镜面中心点的周期性栅格函数的傅立叶变换。这个函数GF以及PSFseg,在位置为w=0时总是取最大值1。而子镜面的点分布函数的零点和这个栅格函数的其他最大值相重合,这样只能看到唯一的主亮点。
如果围绕镜面中心的六边形子镜面共有M圈,那么子镜面总数将是N=3M(M+1)+1。则:
式中α=(πd/2λF)wx和β=(πd/2λF)wy是在像平面上的归一化坐标。这个函数在像平面上具有两倍于π/3 角度的对称性。而子镜面的场分布(field pattern)函数为(Yaitskova,2003):
在像平面上两个相互垂直的方向上子镜的点分布函数为:
如果子镜面位置具有随机分布且平均值为零的轴向误差(piston error),子镜面的点分布函数将不受影响,不过栅格函数GF将产生变化,从而引进杂乱的光斑噪声。这时主光斑的角直径可以用子镜面点分布函数PSFseg的半高全宽(FWHM )表示,它和轴向误差的具体数值无关,它约等于(Yaitskova,2003):
如果d= 1.5 m和λ= 0.5μm,这个角度是0.07″。由于斑点噪声的存在,斯特尔比会减少。当子镜瞬间随机轴向相位误差为δj时,它的斯特尔比为:
当子镜随机轴向相位误差是完全独立的,并且具有平均值为零的高斯分布,那么斯特尔比为:
式中φ是口径场波阵面均方根差。这时像点平均亮度和中心点最高亮度的比为(Yaitskova,2003):
当子镜面具有随机倾斜误差(tip-tilt error)时,拼合镜面的成像要复杂些。这时它的点分布函数也可以表示成几项乘积的形式,即各子镜点分布函数的非相关叠加以及子镜面之间的相互干涉。当子镜面之间仅仅存在微弱相干,而子镜面的倾斜误差随机分布且数值大时,子镜面之间的相干可以忽略不计。这时观察到总点分布函数就是各个具有倾斜角的N个子镜面点分布函数的叠加。不过当子镜面数量不断增加以后,它的第二项对总点分布函数的贡献将很快增加。这种情况也可以看作一个栅格函数和一个经过修改以后的子镜面点分布函数的乘积。这个修改以后的子镜面点分布函数是子镜面的振幅分布函数(1.151式)和一个与倾斜相位误差的分布密切相关的Q 函数乘积的平方:
当倾斜相位误差具有高斯分布时,这个Q函数为:
这时斯特尔比和子镜面数量有直接关联,其表达式为(Yaitskova,2003):