大气湍流随机相位屏的数值模拟和验证

2007年3月 Opto-Electronic Engineering March, 2007 文章编号：1003-501X(2007)03-0001-04

大气湍流随机相位屏的数值模拟和验证

王立瑾1，2，李 强1，2，魏宏刚1，廖 胜1，沈忙作1

( 1．中国科学院光电技术研究所 微细加工光学技术国家重点实验室，四川 成都 610209；

2．中国科学院研究生院，北京 100039 )

摘要：本文使用傅立叶变换法对符合Kolmogorov谱的大气随机相位屏进行了数值模拟，采用三层随机相位屏叠加模拟大气湍流。通过比较模拟相位屏的相位结构函数和理论值的符合度对相位屏的统计特性进行了验证，利用长曝光传递函数和相位结构函数计算大气相干长度r0。结果表明，用傅立叶变换法模拟的随机相位屏是正确的，但相位结构函数存在明显的低频空间频率成分不足，采用三层随机相位屏叠加对相位结构函数和理论值的符合度有所改善，实际r0值比设计值偏大。

关键词：大气湍流；随机相位屏；傅立叶变换法；大气相干长度 中图分类号：TP391 文献标识码：A

Numerical simulation and validation of phase screen

distorted by atmospheric turbulence

WANG Li-jin12，LI Qiang12，WEI Hong-gang1，LIAO Sheng1，SHEN Mang-zuo1

，

，

( 1. The State Key Lab. of Optical Technologies for Microfabrication, the Institute of Optics and Electronics,

the Chinese Academy of Sciences, Chengdu 610209, China；

2．Graduate School of the Chinese Academy of Sciences, Beijing 100039, China )

Abstract：Numerical simulation of Phase Screen (PS) distorted by atmospheric turbulence following the Kolmogorov spectrum by using Fourier Transform is implemented. In order to simulate turbulence, three PSs are placed together. Accuracy of PS is determined by comparing the phase structure function and theoretical results, and atmospheric coherent length r0 is calculated through optical transform function of long exposure and phase structure function. The experimental result shows that the simulated PS with Fourier transform is correct, but low spatial frequency components of PS generated by Fourier Transform significantly deviate from theoretical value. Three PSs put together improve the conformity of phase structure function and theoretical value, and the actual r0 is bigger than the theoretical value. Key words：Atmospheric turbulence; Phase screen; Fourier transform; Atmospheric coherent length

引 言

当光波通过地球大气层时，大气湍流扰动光波的波前相位随时间迅速变化，导致天体目标通过大气的成像质量下降。在0.01~0.001s甚至更短的曝光时间内，大气湍流可以假定为“冻结”的[1位屏来模拟大气湍流。目前，已发展了多种方法生成模拟大气湍流效应的随机相位屏[3-5波前，常用的是首先由McGlamery提出的傅立叶变换法[3

，9]

，2]

，采用随机相，数值模拟的

，9-10]

方法基本上可分为两类：一是频率域间接模拟，这种方法根据大气湍流的功率谱密度函数得到模拟的光学

；二是空间域的直接模拟，利用一组正交完备

的基函数来表示波前相位，比较常用的是Zernike多项式，Noll和Roddier对此作了大量深入的研究[4-5]。

收稿日期：2006-07-10；收到修改稿日期：2006-12-10

基金项目：国家863高技术项目

作者简介：王立瑾(1977-)，女(汉族)，四川成都人，硕士，主要从事光电图像处理研究。E-mail: ljwang456@sina.com

2

光电工程 第34卷第3期

在生成随机相位屏后，需要进一步对其统计特性是否与大气湍流的统计特性相符合进行验证，通常是计算相位结构函数，并与理论值比较，此外还可以通过长曝光传递函数和相位结构函数计算大气相干长度r0，与确定的理论模拟值进行比较。

本文首先用傅立叶变换法生成了单层随机相位屏，然后使用三层相位屏叠加的方法来模拟大气湍流的空间统计特性，模拟了点目标经过大气后的短曝光和长曝光像，并使用大气相位结构函数和长曝光光学传递函数对模拟相位屏的空间统计特性进行了验证，在此基础上计算了大气相干长度r0，证明了模拟方法和模拟结果的正确性。

1 随机相位屏的数值模拟

利用傅立叶变换法生成相位屏，首先是生成一复高斯随机数矩阵，然后用符合Kolmogorov大气相位扰动的功率谱函数对其进行滤波，再进行逆傅立叶变换得到[6]。这一过程可以表示为：

φ(x,y)=∑kx

∑k

y

h(kx,ky)FΦ(kx,ky)exp[j(kxx+kyy)]ΔkxΔky (1)

式中 φ(x，y)是大气扰动相位，h(kx，ky)为频率域里的零均值、单位方差的Hermitian复高斯随机变量，FΦ(kx，

ky)是大气折射率起伏引起的符合Kolmogorov谱的功率谱函数：

FΦ(kr)=2πk2ΔzΦ(k,z)ΔkxΔky (2)

式中 Δz为湍流薄层厚度，Φ(k，z)是折射率功率密度函数，这里采用符合Kolmogorov模型的功率谱密度函数[6]。将(1)式离散化，可以表示为

φ(m,n)=

(Nx/2)−1

m'=−(Nx/2)n'=−(Ny/2)

∑

(Ny/2)−1

∑

h(m',n')f(m',n')exp[j2π(

m'mn'n

+)] (3) NxNy

式中 f(m',n')=

2πGxGy

0.00058r0−56fx2+fy2

()−1112 (4)

式中 Gx，Gy表示相位屏的大小，r0表示大气相干长度。

第4页本文设定望远镜的孔径为1.76m，光波长λ=0.7μm，大气相干长度0.1m，望远镜口径与大气相干长度之比为17.6，其中光学系统的入瞳口径为13.656mm，相位屏大小G=90mm，采样点数为4096×4096。在今后进一步要实现的湍流模拟实验装置中，最简单的装置即是把单层相位屏置于望远镜的入瞳处，但是在单层相位屏上相位分布的重复性较高，因此，采用叠加的三层随机相位板来模拟大气湍流，能够实现更多的空间分布组合，保证相位分布的随机性。为了保持望远镜口径与大气相干长度之比不变，三层相位屏的等效大气相干长度

r0=0.776mm，并且三层相位屏是紧邻的。本文中模拟时，r0i是单层相位屏的大气相干长度，r0和r0i有如下关系：

[7]

3

5001000150020002500300035004000

500

1500

2500

3500

r0=∑r

i=1

[−530i

]−35 (5)

当r0=0.776mm时，则每层相位屏大气相干长度

图1 模拟的Kolmogorov相位屏

Fig.1 Simulated Kolmogorov phase screen

图2 三层相位屏的短曝光像

Fig.2 Short exposure image with three phase screens

r01=r02=r03=1.5mm。本文使用

傅立叶变换法模拟符合Kolmogorov谱的随机相位屏，模拟结果如图1(采样点数为4096×4096)。

本文模拟短曝光图像时，目标通过三层随机相位屏经光学系统成像，在非相干光条件下，系统的点扩散函数(PSF)为：

h(x,y)=IFTp(u,v)e−jφ(u,v){}2

(6)

2007年3月 王立瑾 等：大气湍流随机相位屏的数值模拟和验证

3

式中 p(u，v)为望远镜光瞳函数，φ(u，v)为相位屏，IFT{·}表示逆傅立叶变换，u、v是光瞳上的坐标，x、 y是像面上的坐标。当目标为点目标时，h(x，y)即为点目标所成的像(短曝光像)，如图2(采样点数为300×300)。

2 湍流大气模拟结果的验证

大气湍流相位的统计特性可以用相位结构函数来描述，因此可以将结构函数作为验证模拟的相位屏正确与否的一种判断标准。此外，大气相干长度r0表征了大气湍流的强度，因此通过计算r0也可以验证所模拟的随机相位屏的正确性。

在验证过程中，假设目标位于等晕区内，模拟点目标通过三层随机相位屏成像时，通光口径为

13.656mm。因此，可在三个大随机相位屏上分别随机选取大小为13.656mm的区域相互叠加，多次重复，就可以得到多个不同分布的随机相位屏，进一步由公式(6)可以模拟生成多帧短曝光像。 2.1 相位结构函数

相位结构函数的定义为

D(r)=[φ(ρ+r)−φ(r)] (7)

2

式中 ρ和r表示二维坐标。Fried给出了与Kolmogorov谱对应的结构函数的定义式[1]：

D(r)=6.88(r/r0)53 (8)

图3(a)和图3(b)分别给出了等效大气相干长度r0为0.776mm时三层相位

屏和单层相位屏的相位结构函数，图中横坐标r表示距离，纵坐标表示相位结构函数D(r)。从图3可以看出，用傅立叶变换法生成的相位屏，由于傅立叶变换法采样频率的问题，造成低频成分丢失，

0

0 0.51.01.52.02.53.03.54.04.5

r/mm

0

00.51.01.5 2.0 2.5 3.0 3.5 4.04.5

r/mm

图3(a) 三层相位屏的相位结构函数；(b)单层相位屏的相位结构函数

Fig.3(a) Phase structure function with three phase screens, (b)Phase structure function with a single phase screen

相位屏的结构函数与理论比较存在着明显的低频部分不足，而高频部分比较吻合。比较图3(a)和图3(b)可以看出，三层相位屏的相位结构函数与理论曲线更加吻合。 2.2 计算大气相干长度r0

理论长曝光传递函数定义式为[8]

HLE(f)=exp[−3.44(λfdir0)53] (9)

式中 λ是波长，di是光瞳和像面之间的距离，r0是大气相干长度。由式(10)可以看出，当HLE(f1e)=1e时 相应的频率坐标为f，r0可由式(10)计算得到：

r0=(3.44)×λdif (10)

35因此，如果计算得出模拟相位屏的长曝光传递函数HLE(f)，就能根据式(9)得出大气相干长度r0。本文中使用了如下两种方法计算HLE(f)。

方法一：在各态历经假设下，系综平均传递函数与在无限长积分时间极限下的长曝光传递函数相同[1]。因此，在理论上长时间曝光的点扩展函数等于瞬时点扩展函数的平均，根据傅立叶变换的性质，在傅立叶频率域内系统的传递函数为[8]：

E{H(f)}=HDL(f)HLE(f) (11)

式中 E{⋅}表示平均，FT{·}表示傅立叶变换，f是频率域里的二维坐标，h(f)为瞬时点扩展函数(PSF)，

H(f)=FT{h(f)}，HDL(f)是衍射极限传递函数，HLE(f)是大气的长曝光传递函数。那么，根据模拟的

4

光电工程 第34卷第3期

长曝光传递函数，就可由式(10)计算出模拟相位屏的大气相干长度r0。图4(采样点数为300×300)中给出了

1000帧短曝光像平均后得到的长曝光像。

方法二：长曝光传递函数和相位结构函数存在着如下关系[1]：

HLE(f)=exp(−D(λfdi)/2) (12)

将根据式(7)计算得到的相位结构函数代入式(12)，得到HLE(f)函数后，再根据式(10)也可计算出r0。 图5中分别给出了理论长曝光传递函数，以及由方法一和方法二计算得到的模拟长曝光传递函数，图中横坐标表示空间频率，纵坐标表示长曝光传递函数。比较图5中各条长曝光传递函数曲线，可以看出模拟长曝光传递函数与理论值符合较好。表1中给出了理论r0值和模拟值的比较，其中，True Value表示在数值模拟时预设的理论r0值，Theory表示由理论长曝光传递函数根据式(11)计算的r0值，Method 1表示根据方法一计算得到的r0值，Method2表示根据方法二计算得到的r0值。从表1看出，模拟相位屏的r0值与理论值比较接近，但是存在一定的偏差。

表1 两种方法计算出的结果(r0：mm)Table 1 Results for the two methods (r0 is in millimeter) r0True valueTheory Method 1 Method 20.9684 0.776 0.7839 0.9222 从以上验证结果可以看出，用傅立叶变化法模拟所得的随机相位屏的统计特性与理论值存在一定的误差，虽然已经采用三层相位屏叠加和在大相位屏上截取的办法进行改善，但仍存在低频部分不足的问题，并且使用长曝光传递函数计算得到的r0值本身就会偏大(比较Theory和True value可以看出)。

3 结 论

本文使用傅立叶变换法生成了符合Kolmogorov谱的随机相位屏，用三层随机相位屏叠加来模拟大气湍流，三层随机相位屏是紧邻的，光瞳口径小于相位屏大小，因此可得到多个不同分布的模拟相位屏，模拟计算出短曝光像和长曝光像，并通过比较模拟相位屏的相位结构函数和理论值的符合度对相位屏的统计特性进行了验证，还利用长曝光传递函数和相位结构函数计算了大气相干长度r0。结果分析表明，用傅立叶变换法模拟的大气湍流相位屏是正确的，但相位结构函数存在明显的低频空间频率成分不足，在高频部分比较吻合，r0值与设计值相比偏大。 参考文献：

[1] 张逸新，迟泽英. 光波在大气中的传输与成像[M]. 北京：国防工业出版社，1997.

ZHANG Yi-xin，CHI Ze-ying. Transmitting and Imaging through Atmospheric Turbulence[M]. Beijing：National Defence

Industry Press，1997.

[2] 沈忙作. 斑点干涉术在天文学中的应用[J]. 光学工程，1982，(4)：49-62.

SHEN Mang-zuo. The Application of Speckle Interferometry in Astronomy[J]. Opto-Electronic Engineering，1982，(4)：

49-62.

[3] Benjamin L.McGlamery. Computer simulation studies of compensation of turbulence degraded images[J]. SPIE，1976，74：

225-233.

[4] Robert J.Noll. Zernike polynomials and atmospheric turbulence[J]. J.0pt.Soc.Am，1976，66(3)：207-211. (下转第9页)

2007年3月 全 薇 等：高级像差对人眼成像质量和视觉的影响

9

4 讨 论

本文利用Hartmann-Shack波前传感器人眼像差仪测量了15只眼6mm瞳孔Zernike第2级到第10级波前像差，除Zernike第2级波前像差离焦和像散之外，第3级到第6级的rms值大于或接近Marechal衍射极限，是不可忽略的。对所测量的15只眼的MTF和Strehl 比率的分析表明，Zernik前6级像差对视网膜成像质量的影响较大，Zernik前6级像差被矫正之后，人眼的成像质量才可能接近衍射极限。对所测量的15只眼的视锐度和对比敏感度函数的分析表明，Zernik前6级像差对视觉的影响是不可忽略的。更高级的

像差对人眼成像质量和视觉的影响较小，甚至可以忽略。因此，在对人眼高级像差矫正时，矫正到Zernik第6级会得到相当好的视觉。

参考文献：

[1] EI Hage S. G，Berny F. Contribution of the crystalline lens to the spherical aberration of the eye [J]. J. Opt. Soc. Am，1973，

63(3)：205-211.

[2] Millodot M.，Sivak J. Contribution of the cornea and lens to spherical aberration of the eye [J]. Vision Res，1979，19(6)：

685-687.

[3] Sivak J，Kreuzer R.O. Spherical aberration of the crystalline lens [J]. Vision Res，1983，23(1)：59-70.

[4] Tomlinson A，Hemenger R. P，Garriott R. Method for estimating the spherical aberration of the human crystalline lens in vivo

[J]. Ophthal. Vis. Sci，1993，34(6)：621-629.

[5] Ames A，Proctor C.A. Dioptrics of the eye [J]. J. Opt. Soc. Am，1921，5(1)：22-84.

[6] Von Bahr G. Investigations into the spherical and chromatic aberrations of the eye，and their influence on its refraction [J]. Acta

Ophthal，1945，23(1)：1-47.

[7] Koomen M，Tousey R，Scolnik R. The Spherical aberration of the eye [J]. J. Opt. Soc. Am，1949，39(3)：370-376. [8] Roorda A，Bobier W R. Geometrical technique to determine the influence of monochromatic aberrations on retinoscopy [J]. J.

Opt. Soc. Am. A，1996，13(1)：3-11.

[9] N. M. Jansonius，A. C. Kooijman. The effect of spherical and other aberrations upon the modulation transfer of the defocused

human eye [J]. Ophthal. Physiol. Opt，1998，18(6)：504-513.

[10] Junzhong Liang，Gerald Westheimer. Optical performance of human eyes derived from double-pass measurements [J]. J. Opt.

Soc. Am. A，1995，12(7)：1411-1416.

[11] Junzhong Liang，David R. Williams. Aberrations and retinal image quality of the normal human eye [J]. J. Opt. Soc. Am. A，

1997，14(11)：2873-2883.

JJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJ (上接第4页)

[5] Nicolas Roddier. Atmospheric wavefront simulation using Zernike polynomials[J]. Opt. Eng，1990，29(10)：1174-1180. [6] E.M.Johansson，D.T.Gavel. Simulation of stellar speckle imaging[J]. SPIE，1994，2200：372-383.

[7] S.E.Troxel，B.M.Welsh，M.C.Roggemann. Off-axis optical transfer function calculations in an adaptive-optics system by means

of a diffraction calculation for weak index fluctuations[J]. J.0pt.Soc.Am，1994，11：2100-2111. [8] M.C.Roggemann，B.M.Welsh. Imaging through turbulence[M]. New York：CRC Press，1995. [9] 杨连臣，沈忙作，郭永洪. 天体目标斑点成像的模拟[J]. 光子学报，2000，29(12)：1108-1111.

YANG Lian-chen，SHEN Mang-zuo，GUO Yong-hong. The Speckle Imaing Simulation of Space Objects[J]. ACTA

PHOTONICA SINICA，2000，29(12)：1108-1111.

[10] 魏宏刚，沈忙作，廖 胜. 扩展目标克服大气湍流成像的实验室模拟[J]. 光电工程，2005，32(2)：40-44.

WEI Hong-gang，SHEN Mang-zuo，LIAO Sheng. Laboratory Simulation of Extended Object Imaging Overcoming Atmospheric

Turbulence[J]. Opto-Electronic Engineering，2005，32(2)：40-44.