非均匀步长分步傅里叶算法的改进

韩星星,赵丽华

(西安交通大学城市学院物理教学部,陕西 西安 710018)

1 引 言

分步傅里叶算法(SSM)是一种数值算法,被广泛用于模拟激光在湍流介质中的传输和散射。此算法最早由Martin等人于1988年用于大气湍流中的光传输中[1],并限定了算法的适用范围。1990年,Martin针对算法中存在的边界反射问题,提出了吸收边界法对算法进行了改进,扩展了其应用范围[2]。2000年,Belmonte对算法中各参数取值问题进行了深入地研究,明确了参数取值范围[3]。2005年,Belmonte 把分步算法扩展应用到了双程传输问题中,研究了漫射目标回波平均强度随着目标高度的变化规律[4]。2006年,Xiao Xifeng等人把分步算法扩展到部分相干光束的传输模拟中,研究了湍流中部分相干光的闪烁问题[5-6]。2009年,钱仙妹等人研究了非均匀湍流中的波束传输模拟算法,比较了三种备用的间距选择方案[7],结果显示过大的相位屏间距会显著增大误差。之后分步算法虽多有应用,但对参数取值问题的讨论并没有太大进展[8-14],算法中的边界反射问题始终限制了相位屏间距的取值范围,并造成较大的误差。

本文针对边界反射问题,对湍流大气中激光波束斜程传输的仿真算法进行了改进,在吸收边界法的基础上,引入滤波函数法,极大拓展了仿真步长的取值范围。在此基础上,模拟计算了激光沿斜程路径上行及下行传输过程中相干长度及闪烁指数的变化规律,并与理论值进行了比较。

2 分步傅里叶算法

激光在湍流大气中的传输过程可以用分步傅里叶算法来模拟,模拟过程如图1所示。

图1 分步傅里叶算法原理图

步进公式[16]为:

(1)

等效相位屏Si(ρ′)由谱反演法生成,离散形式的相位屏生成公式为:

(2)

(3)

其中,r0是大气相干长度:

(4)

大气湍流结构常数随着离地面高度的分布为:

(5)

模拟默认参数如下:D=0.6 m,N=512,λ=1.064 μm。

3 分步傅里叶算法步长取值要求

分步傅里叶算法把连续分布的湍流大气的扰动等效为一系列随机相位薄屏的扰动,这是一种近似行为,其适用范围是有要求的。根据Martin和Flatte的研究结果[15],要保证相位屏近似成立,需:

(6)

从物理意义上讲,分步算法是把光场分解为一系列不同方向的平面波,各平面波在真空中传输后在下一个相位屏前重新合成新的光场。方向不平行于传播方向的光场在真空传播过程中会超出相位屏范围,但是在分步傅里叶算法中,相位屏中的能量是守恒的,这意味着向外发散的能量会被边界反射回来,从而对屏上的场带来误差。

为了解决这个问题,传统方法是给相位屏边界设置一个吸收层,到达边界一定距离的场被逐渐衰减。设吸收层的厚度为Da,一般Da的取值不超过相位屏宽度的1/8。对于距离边界距离为xa的网格点,衰减系数为[15]:

Sa=[1-cos(πxa/Da)]/2

(7)

图2是利用分步傅里叶算法模拟一个准直高斯波束在真空中传播5000 m后形成的光强分布曲线,其中,M是路径上的相位屏数目。在无吸收边界的情况下,光强分布曲线与相位屏数目无关,高频分量反射造成的光强的剧烈起伏充满整个相位屏,此时的光场的模拟显然是无效的。在设置吸收边界的情况下,相位屏数越多,模拟得到的光强分布曲线与理论曲线越吻合。当M=20时,在吸收边界以内,光强模拟值与理论值已经非常接近了。这说明吸收边界可以有效地解决高频分量的边界反射问题。

图2 在真空中利用分步傅里叶算法模拟准直高斯波束传输距离L后的光强分布(附加吸收边界),其中,L=5000 m,W0=0.5 cm

4 滤波函数法

在斜程路径上,湍流主要分布在大气底层,而在高空中湍流很弱。根据公式(6),超过一定高度后,不论路径长度,只需要设置一个相位屏就可以。但图2表明,吸收边界对相位屏间距的拓展是有限的,即使不考虑湍流,步长超过一定值后也会引起较大误差,步长的最大值约为250 m。

这种情况下,现有文献中的解决方法是增加相位屏数目[4],但这会降低模拟效率。为了提高模拟效率,本文提出滤波函数法。

当激光传播极远的距离后,光轴附近的场会趋向于平面波,这意味着在屏间传输过程中高频分量被衰减,低频分量才能到达下一个相位屏,即相位屏间的自由传输相当于一个低通滤波器。而传统分步算法中各空间频率分量的强度是不变的,这显然不符合现实。因此需要给算法中引入一个低通滤波函数。滤波函数形式参照边界吸收函数:

(8)

其中,截止频率为:

κh=kD/Δz

(9)

图3显示了使用滤波法模拟得到的准直激光在传播不同距离后的光强分布情况,在模拟过程中保持M=1。在相位屏中间区域,模拟结果与理论值吻合得很好,但在靠近边界处,模拟结果出现较大的误差。在不同传播距离条件下,模拟结果有效的区间都大致相同,约为-D/4

图3 在真空中利用分步傅里叶算法

5 大气湍流中斜程路径上相位屏的设置方式

在斜程路径上,湍流沿路径分布不均匀,高度越高,湍流越弱,相位屏之间的间距也应该越大。但考虑到滤波函数法对有效区域缩减严重,不宜多用,因此综合吸收边界法和滤波函数法,图4给出了一种斜程路径上的相位屏设置方式,此路径倾斜角(即路径方向与水平面的夹角) θ=10°,地面处湍流结构常数C0=1×10-13m-2/3,图4中的曲线是湍流强度随高度的分布,水平线代表了相位屏的位置。在z<1750 m这一段路径上,相位屏均匀设置且间距250 m,相位屏附加吸收边界。在z=1750 m处设置相位屏,同时附加吸收边界和滤波函数,其后不论路径长度不再设置相位屏。

图4 斜程路径上的相位屏设置方式

图5给出了准直高斯波束(W0=3 cm)在湍流大气中沿图4所示的倾斜路径上行、下行传输距离L=5 km的光斑样本。

图5 高斯波束在湍流中沿斜程路径传输的光斑样本

由图5可见,在上行路径上,光束充满了整个相位屏(除了边界吸收区域内),且形成较大的散斑。在下行路径上,光束仍局限在相位屏中心较小区域内,形成的散斑也较小。这是因为对激光有明显影响的湍流基本上只存在于距地面高度1500 m以下,上行传输的激光在初始阶段就受到湍流的影响,相位受到较大的扰动,反映在接收屏上即光束的扩展较大。在之后的更长的路径上,激光相当于在真空中传输,散斑场在真空中传输越远,散斑尺寸就越大。而下行波束则相反,在初始路径上不受扰动,波束基本保持完全相干,在接近地面时光场才受到湍流的影响,形成了尺寸很小的散斑。

6 光斑统计特征计算验证

6.1 相干长度

在大气湍流中,模拟得到的每个光斑样本都是随机的,但这些光斑符合特定的统计特性。为了验证本文提出的滤波函数法在湍流环境下的有效性,我们在模拟得到100个样本的基础上,统计计算了湍流中光场的散斑场尺寸和强度起伏方差,即相干长度和光强闪烁指数,并与理论值作对比。

要计算光场的相干长度,首先要计算光场的相干度函数[14]:

(10)

其中,I(m,n)是对所有光斑样本上在同一位置(m,n)上对光强求平均得到。

对所有的γ(iΔx,jΔx)求平均得到γT(iΔx,jΔx),沿幅角求平均可得相干度函数γ(ρ),利用γ(ρ0)=e-1得到接收场相干长度ρ0。

图6是准直激光波束激光沿斜程路径传输时的场相干长度。可见模拟值与理论值吻合得很好,说明本文提出的滤波函数及相应的相位屏设置方式是可行的。在大部分情况下,模拟值比理论值略小,尤其在上行路径上,这时的波束发散比较强烈,虽然附加了吸收边界和滤波函数,及仍然不能完全消除部分高频成分的边界反射。

图6 湍流大气中激光波束斜程传输时相干长度模拟结果

在传输路径较短时,激光经过同样距离的上行或下行传输后,相干长度差别很小,这是因为此时路径全部位于湍流范围中。随着路径长度的增大,上行路径与下行路径的差别越来越大,激光相干长度差别也越大。在下行路径上,相干长度随着路径长度的增大而逐渐减小,最终稳定在1 cm左右,所以在大多数下行传输的实际应用中,如卫星对地面的激光通信中,不论卫星高度如何,地面接受到的激光散斑大小是固定的。相干长度逐渐减小。随着路径长度的增大,上行路径和下行路径上波束特性的差别越来越大。在上行路径上,相干长度存在一个最小值,在越过最小值位置后,路径越长,激光散斑越大。这和图5中得到的信息是一致的。

6.2 闪烁指数

光斑的闪烁指数是另一个很重要的统计量,由于光斑不同位置处,闪烁指数是不同的,本文只统计计算了光斑的中心闪烁指数。理论上只能利用相位屏中心点处光强来统计,但实际上为了节省计算量,光斑中心一个圆域内的所有点都可以作为有效点。假设光斑的有效半径为We,有效区域半径取为Ne=0.2We/Δx。

对于每个光斑样本,都得到对应的闪烁指数[14]:

(11)

对所有光斑的闪烁指数求平均得:

(12)

图7是准直激光波束激光沿斜程路径传输时的闪烁指数。理论显示对于同一条路径,上行传输和下行传输时的闪烁指数是相同的。模拟得到的闪烁指数比理论值略大,而下行路径上的闪烁指数又比上行路径上的大,但总体来说三条曲线之间的误差在合理范围内,本文提出的算法是可行的。

图7 湍流大气中激光波束斜程传输时闪烁指数模拟结果

7 结 论

本文主要研究了斜程路径上相位屏的设置问题,结果显示吸收边界法可以有效提高分步算法的效率(更少的相位屏)和精度(更小的误差),但步长存在上限。而滤波函数法彻底消除了对步长的上限(在真空中),但是相位屏有效区域缩小较严重(只有整个相位屏面积的1/4)。两种方法相结合用于大气湍流中激光波束斜程传输的仿真,减小了相位屏的设置,提高了仿真效率。在此基础上,本文对斜程湍流中激光光斑的相干长度及闪烁指数进行了统计计算,结果与Rytov理论给出的结果比较吻合,再一次证明了本文提出的算法的有效性。

参考文献:

[1] Martin J M,Flatte S M.Intensity images and statistics from numerical simulation of wave propagation in 3-D random media[J].Appl.Opt.,1988,27(11):2111-2126.

[2] Martin J M,Flatte S M.Simulation of point-source scintillation through three-dimensional random media[J].J.Opt.Soc.Am.A.,1990,7(5):838-847.

[3] Belmonte A.Feasibility study for the simulation of beam propagation:consideration of coherent lidar performance[J].Appl.Opt.,2000,39(30):5426-5445.

[4] Belmonte A.Coherent return turbulent fluctuations in ground lidar systems profiling along slant paths[J].Optics Express,2005,13(23):9598-9604.

[5] Xiao X,Voelz D.Wave optics simulation of pseudo-partially coherent beam propagation through turbulence:application to laser communications[J].USA:SPIE-The International Society for Optical Engineering,2006:63040-63041.

[6] Xiao X,Voelz D.Wave optics simulation approach for partial spatially coherent beams[J].Optics Express,2006,14(16):6986-6992.

[7] QIAN Xianmei,ZHU Wenyue,RAO Ruizhong.Phase screen distribution for simulating laser propagationalong an inhomogeneous atmospheric path[J].Acta Physica Sinica,2009,58(9):6633-6639.(in Chinese)

钱仙妹,朱文越,饶瑞中.湍流路径上光传播数值模拟的相位屏分布[J].物理学报,2009,58(9):6633-6639.

[8] Xiao X,Voelz D.Wave optics simulation of partially coherent and partially polarized beam propagation in turbulence[J].Proc.SPIE,2009,7464：74640T.

[9] Qian X M,Zhu W Y,Rao R Z.Numerical investigation on propagation effects of pseudo-partially coherent Gaussian Schell-model beams in atmospheric turbulence[J].Optics Express,2009,17(5):3782-3791.

[10] QIAN Xianmei,ZHU Wenyue,RAO Ruizhong.The aperture averaging effect of scintillation of pseudo-partially Gaussian-Schell model beampropagation in turbulent atmosphere[J].Acta Physica Sinica,2013,42(4):044203.(in Chinese)

钱仙妹,朱文越,饶瑞中.伪部分相干高斯-谢尔模型光束在湍流大气中传播的闪烁孔径平滑效应[J].物理学报,2013,42(4):044203.

[11] QIAN Xianmei,RAO Ruizhong.Spatial distribution of Gaussian-beam scintillation in atmosphere by numerical simulation[J].Chinese Journal of Quantum Electronics,2006,23(3):320-324.(in Chinese)

钱仙妹,饶瑞中.间分布的数值模拟研究[J].量子电子学报,2006,23(3):320-324.

[12] QIAN Xianmei,ZHU Wenyue,RAO Ruizhong.Numerical simulation of turbulent effects of laser propagation alonga ground-space slant atmospheric path[J].Infrared and Laser Engineering,2008,37(5):777-782.(in Chinese)

钱仙妹,朱文越,饶瑞中.地空激光大气斜程传输湍流效应的数值模拟分析[J].红外与激光工程,2008,37(5):777-782.

[13] QIAN Xianmei,ZHU Wenyue,HUANG Yinbo,et al.Selection of computing parameters in numerical simulation of laser beam propagationin turbulent atmosphere[J].Acta Photonica Sinica,2008,37(10):1986-1991.(in Chinese)

钱仙妹,朱文越,黄印博,等.激光湍流大气传输数值模拟中计算参量的选取[J].光子学报,2008,37(10):1986-1991.

[14] QIAN Xianmei,ZHU Wenyue,RAO Ruizhong.Simulation of the focal shift of a focused gaussian-beam propagating in turbulent atmosphere[J].Journal of Atmospheric and Environmental Optics,2006,1(2):85-88.(in Chinese)

钱仙妹,朱文越,饶瑞中.湍流大气中聚焦高斯光束焦移的数值模拟分析[J].大气与环境光学学报,2006,1(2):85-88.

[15] Flatte S M,Gerber J S.Irradiance-variance behavior by numerical simulation forplane-wave and spherical-wave optical propagation through strong turbulence[J].J OptSoc Am A,2000,17(6):1092-1097.

[16] WANG Liguo.Characteristics of reflected wave from targets illuminated by laser beam in turbulent atmosphere[D].Xi′an:Xidian University,2014:169.(in Chinese)

王利国.湍流大气中激光波束目标回波特性[D].西安：西安电子科技大学,2014:169.