基于Dini展開的高階Hankel變換及其在光束傳輸中的應用*

游開明 林燕玲 王友文 陳列尊 戴志平 陸世專

(湖南衡陽師范學院物理與電子信息科學系,衡陽 421002)

(2013年2月2日收到;2013年3月30日收到修改稿)

1 引言

由于Hankel變換(Hankel transform,HT)在物理學的很多分支(如光學、聲學、電磁學、分子生物學等)有著廣泛的應用[1-4],引起了很多科學工作者的關注,并對它的數值計算算法進行了研究.1977年,Siegman[5]首先對徑向坐標抽樣點用指數形式量化,將HT化為離散形式,用快速傅里葉變換(FFT)進行計算,得出了準快速Hankel變換算法.此后,出現了很多用于物理學各個方面的HT算法,如文獻[6—18].文獻[9—11]把HT看作是一個二變量函數的簡單積分的傅氏變換,先用一般方法計算簡單積分,再利用FFT得到HT.2003年,Markham和Conchello[12]總結了計算HT的6種方法,并把它們用于計算振蕩函數.文獻[13]和[14]分別用正弦和余弦變換進行快速HT及非均勻快速HT等.2007年,Cerjan[15]把變換函數在單位區間上用Zernike多項式表示,再根據Zernike多項式基礎集與Bessel函數展開的對偶關系進行HT,其展開的系數采用Simpson規則或Gauss-Kronrod求積法計算,使HT變換在精度上有所提高.1998年Yu等[16]對徑向坐標按零階Bessel函數零點分布的比例進行離散,將變換函數做Fourier-Bessel級數展開,實現了準離散HT(quasi-discrete HT,QDHT),2004年Guizar-Sicairos和Gutie′rrez-Vega[17]又把它推廣到p階HT(the p th-order quasi-discrete HT,p QDHT),這類方法利用變換矩陣使計算變得簡單,因而速度也較快,由于能量守恒,逆變換精度很高,不依賴輸入函數,能適應反復變換,具有很高的實用價值.遺憾的是基于Fourier-Bessel級數展開HT方法處于條件收斂,不能計算空間域和空間頻率域端點處的變換值,同時其0階HT也不能計算對稱軸上點的變換值,不能直接用于研究傳輸問題中對稱軸上的特性.文獻[18]給出了一種基于Dini級數展開的0階準離散HT算法,由于Dini級數展開比Fourier-Bessel級數展開收斂更快,精度更高,對于0階HT還能計算對稱軸上的值,所以作者把它應用到了解決光束非線性傳輸問題.然而物理學中更一般的傳輸問題是高階的,如高階Bessel光束傳輸問題等[19-21],又基于Dini級數展開0階HT與高階HT畢竟還有一定區別,因此本文開展基于Dini級數展開的高階p(p＞0)準離散HT(the p th-order quasi-discrete Hankel transform based on Dini series expansion,p DQDHT)算法的研究,并給出它在光束傳輸中的應用.

2 算法的推導

很多有限區域上圓對稱物理量分布問題,數值方法上使用HT更有效.歸一化坐標的整數p(p＞0)階Hankel變換和逆變換分別表示為[7]

對于光束傳輸,式中r表示歸一化空間徑向坐標,ρ表示歸一化空間頻率徑向坐標,Jp表示第一類p階Bessel函數,b和β分別表示空間變量和空間頻率變量的截斷半徑,S表示帶寬乘積bβ的2π倍,即

對于第二類齊次邊界條件,第一類p階Bessel函數滿足如下關系[18]

式中 αn(n=1,2,···)是 p 階 Bessel函數導數 J′p(x)的正實根,δnm為克羅內克符號,m=n時為1,否則為0.

這里考慮p＞0的情況,將方程(1),(2)中的變換函數 f(br)和g(βρ)展開為Jp的級數,分別有

式中的系數 fn和gn由方程(4)確定,再結合方程(1)和(2),分別有

上述級數展開稱為Dini級數展開[22].

將歸一化徑向坐標變量r和ρ分別量化為

綜合方程(5)到(8),得到p階離散Hankel變換對為

這里N為數值計算點數.為了減少方程(10)和(11)右邊求和符中數乘次數,分別定義向量G(m)和F(n)為

于是方程(10)和(11)簡化為如下形式

其中Cnm為N×N變換矩陣C的矩陣元,表示為

根據方程(14)—(16),容易看出矩陣C的一些特性：由方程(16)知C是實對稱矩陣,滿足C=CT;將方程(14)代入(15)或反過來代換,發現F(n)和G(m)均為自身變換,說明正、逆變換毫無誤差時還需要C是正交矩陣,滿足CCT=I(I表示單位矩陣),即C的行列式det[C]=±1;由方程(16)還知C是單參數S的函數,C=C(S),S不同,矩陣C中元素的值不同,所以如何選擇S以使C成為正交矩陣對提高該數值算法的精度非常重要.

根據方程(1),(2)和(9)應取S=αN,然而我們通過數值驗證,得知矩陣C并不是很好的正交矩陣.若取S=jN+1(其中 jN+1為Jp(x)的第N+1個零點,滿足αN＜ jN+1＜αN+1),矩陣C非常接近正交矩陣,這是因為當方程(1)和(2)中的r=1和ρ=1時,其中的變換核函數為Jp(S)=0,可以使變換后的g(ρ)和 f(b)的值為0,所以本算法一般的數值計算中取S=jN+1.若要進一步提高精度,可利用矩陣C的性質對 jN+1進行修正來得到S,其修正表達式為

利用矩陣C不僅能使在一次HT中減少N(3N-4)次復數或實數乘法運算,若將C存儲起來,還可避免下一次變換時對數據的重復計算,這對分步傳輸問題非常重要,可以極大地提高程序的執行速度.

這樣p DQDHT算法為：先由方程(13)將 f(br)轉換為F(n),乘以矩陣C后再按方程(12)將G(m)轉換為g(βρ).逆變換時只要將方程(12)和(13)右邊的ρ和b交換即可.

3 算法的測試與驗證

我們用C++語言編程對p DQDHT算法的正確性與精度進行測試,根據C為實對稱矩陣,對數據的存儲采用三角矩陣形式,可以分別節省N(N-1)/2個雙精度型存儲單元的存儲空間,這用C++語言編程是非常容易實現的.對于Bessel函數及其導數的零點很容易先用公式估算出來[23],然后用牛頓迭代法分別解方程Jp(x)=0和 J′p(x)=0.5(Jp-1(x)-Jp+1(x))=0 很快得到精確值.注意,在本文后面的敘述中r和ρ表示實際坐標,前后的關系分別是r=br和ρ=βρ,同時取S=jN+1.

首先我們用連續函數 f(r)=r2exp(-πr2)作為輸入函數對p DQDHT算法進行測試,它在[0,∞)上的2階HT的精確表達式為

以 g(ρ)作為基準,用 g?(ρ)表示 p DQDHT 算法的變換值,在空間頻率域定義絕對誤差為|g-g?|,并分別用Max E表示最大絕對誤差,Min E表示最小絕對誤差,Mean E表示平均誤差|g-g?|/N,取b=β=(S/2π)1/2,2階HT的計算結果如表1所示.

表1 用p DQDHT計算f(r)的誤差隨抽樣點數N的變化

可見,對于連續函數當抽樣點數N很小時p DQDHT算法就有很高的精度,進一步說明Dini展開收斂很快.

圖1表示函數 f(r)的精確值與經過40次2階p DQDHT后的結果比較,計算中空間域截斷半徑b=4,抽樣點數N=100.圖中的實線表示精確值,圓點線表示2階HT的數值結果,它們幾乎完全重合.在PIV 1.5 GHz的CPU,512 MB內存的PC機上執行程序,p DQDHT算法經過40次變換共花費的時間為0.03 s.可見,p DQDHT算法經過多次變換仍保持極高的精度,而且程序執行速度快.

圖1 函數 f(r)的精確值與經過40次p DQDHT后的結果比較,其中p=2,b=4,N=100

其次,我們用top-hat函數作為輸入函數,它的定義是 f(r)=rp[H(r)-H(r-b)],這里H(r)為Heaviside階躍函數,其p階HT的精確表達式為

分別計算經過一次p DQDHT和p QDHT與(19)式隨空間頻率ρ變化絕對誤差的對數,再進行誤差比較,它們的動態比較結果如圖2所示.計算中空間域截斷半徑b=5,抽樣點數N=200,空間頻率域截斷半徑β在20.025附近.圖2是p=3時p DQDHT和p QDHT兩種算法的結果比較,圖中整個空間頻率范圍內 p DQDHT算法得到的曲線(實線)均位于p QDHT算法得到的曲線(圓點線)的下方,表示p DQDHT算法比p QDHT算法具有更高的精度.圖2還表現出基于Dini級數展開的p DQDHT算法與基于Fourier-Bessel級數展開的p QDHT算法的不同特性,前者誤差較大的地方主要出現在低頻端較小的區域內,在其他的地方保持較高的精度,這對再進行一次HT變換到空間域時減少Gibbs現象是非常有益的;后者在低頻端精度較高,隨著頻率的增大,誤差越來越大,特別是在最高頻率處誤差最大,因此p DQDHT算法更具有優勢.

圖2 “top-hat”函數經過一次 p DQDHT和 p QDHT后與精確值絕對誤差的對數隨空間頻率半徑變化的比較,其中p=3,b=5,N=200

4 p DQDHT在光束傳輸中的應用

為了研究p階準離散Hankel變換在光束傳輸中的應用,我們考慮圓對稱4階Bessel光束f(r)=J4(ktr)exp(i4φ+i kzz)通過凸透鏡后的聚焦演變.讓光束入射到半徑R=4 mm,焦距 f=0.5 m,位于平面z=0的凸透鏡上,并使光束在透鏡光軸上對稱分布,波數k與橫向波數kt和縱向波數kz滿足關系 k2=(2π/λ)2=kt2+kz2.

光束通過透鏡后為非傍軸傳輸,數值方法上采用分步傳輸方法,所以運用p DQDHT解決類似于上述光束傳輸問題的計算步驟可為：

1)選取b,p和N,計算程序中所需的Bessel函數及其導數的零點,將Bessel函數導數的零點存儲于一維數組中;

2)計算矩陣C,并使用指針數組按三角形式存儲起來;

4)計算z=0處初始光束 f(r)=J4(ktr)的值,乘以光束從透鏡前表面傍軸傳輸到后表面的作用因子exp[i kr2/(2f)]后按復數形式存儲于一維數組中(在C和C++中偶地址存實部,奇地址存虛部);

5)選取傳輸步數M,計算傳輸步長Δz=z/M,再進行分步傳輸計算：

for j=0,···,M

對 f(r)進行一次p DQDHT得空間頻率域上的g(ρ);

再進行一次p DQDHT得空間域上的 f(r),存儲光強徑向分布數據,end.

計算中取光的波長λ=632.8 nm,kt=19858.32 m-1,徑向截斷半徑b=4 mm,徑向取樣點數N=256,縱向傳輸距離z=0.75 m,傳輸步數M=300.按照這些參數利用p DQDHT算法編程計算,得到會聚Bessel光束光強I(r,z)的傳播演變如圖3所示.

圖3 4階會聚Bessel光束通過透鏡后光強在焦平面附近的傳播演變

從圖3看出,會聚Bessel光束通過凸透鏡后在焦平面附近三處會聚,會聚的峰值并不出現在對稱軸上,而是在以對稱軸為中心半徑不同的三個圓環上,這是高階貝塞爾光束為空心光束的原因.進一步用數值方法求三個會聚峰值的半徑和位置,得到透鏡焦平面前、焦平面上和焦平面后三個峰值光環的半徑分別為r=0.062647 mm,0.996897 mm和0.110658 mm,會聚平面的位置分別為z=0.38 m,0.5 m和0.72 m,它們的截面光強分布如圖4所示.這些結果與FFT算法得到的結果是一致的,特別是在焦平面上的結果還與文獻[20]用幾何方法得到的會聚圓環的半徑r=fkt/kz很好地相符.

圖4 4階Bessel光束通過透鏡后不同會聚處的截面光強分布 (a)z=0.38 m;(b)z=0.5 m;(c)z=0.72 m

5 結論

基于Dini級數展開,導出了高階Hankel變換的離散表達式,稱為p DQDHT算法,能實現復數或實數形式的HT.用向量表示,可把它們分別寫成為輸出向量是變換矩陣與輸入向量的乘積,計算變得更簡單,計算速度更快.討論了變換矩陣C的性質,以及使它成為正交矩陣與兩個截斷半徑乘積S的關系,定性分析了S的取值,得到S=jN+1時C接近于正交矩陣.為了進一步提高該算法的精度,定量給出了對S進行修正的表達式,同時用數值方法給出了其中的可選參數k與抽樣點數N的關系.

在算法的實施上,事先計算出Bessel函數導數的零點進行靜態存儲,計算變換矩陣C采用三角形式存儲,以免下一次變換對數據重復計算,進一步提高了計算速度.算法中的正變換和逆變換雖然用了兩種表達式,實際上它們的區別只是截斷半徑不同,所以正變換和逆變換可以采用相同的程序段,只是在每次變換后將兩個截斷半徑的值進行交換即可,使程序變得簡單.

通過兩個不同的輸入函數對p DQDHT算法的測試和光束分步傳輸的應用實例,結果表明：p DQDHT算法經過多次變換而不丟失精度,特別適合分步傳輸的場合,這是很多Hankel變換算法做不到的;在整個變換區域內,p DQDHT算法都具有很高的精度,特別是在空間頻率域的高頻端精度更高,在空間域近軸端精度更高,這正符合一般實際應用的需要;程序執行的速度與一般的快速算法相當.

我們相信,借助本文p DQDHT算法在光束分步傳輸中的應用方法,p DQDHT算法在解決物理學各個分支類似的問題中必將發揮重要作用.

[1]Du G H 1988 Acta Phys.Sin.37 769(in Chinese)[杜功煥1988物理學報37 769]

[2]Yan CC,Xue G G,Liu C,Gao SM 2005 Acta Phys.Sin.54 3058(in Chinese)[閏長春,薛國剛,劉誠,高淑梅2005物理學報54 3058]

[3]Yang X J,Zang WP,Tian JG,Liu ZB,Zhou WY,Zhang CP,Zhang G Y 2005 Acta Phys.Sin.54 2735(in Chinese)[楊新江,臧維平,田建國,劉智波,周文遠,張春平,張光寅2005物理學報54 2735]

[4]Chen G B,Wang H N,Yao JJ,Han Z Y,Yang SW 2009 Acta Phys.Sin.58 1608(in Chinese)[陳桂波,汪宏年,姚敬金,韓子夜,楊守文2009物理學報58 1608]

[5]Siegman E 1977 Opt.Lett.1 13

[6]Agrawal G P 1981 Opt.Lett.6 171

[7]Magni V,Cerulle G,De Silvestri S 1992 J.Opt.Soc.Am.A 9 2031

[8]Agnesi A,Reali GC,Patrini G,Tomaselli A 1993 J.Opt.Soc.Am.A 10 1872

[9]Ferrari JA 1995 J.Opt.Soc.Am.A 12 1812

[10]Ferrari JA,Perciante D,Dubra A 1999 J.Opt.Soc.Am.A 16 2581

[11]Perciante CD 2004 J.Opt.Soc.Am.A 2 1911

[12]Markham J,Conchello A J2003 J.Opt.Soc.Am.A 20 621

[13]Knockaert L 2000 IEEETrans.Signal Process.48 1695

[14]Liu Q H,Zhang ZQ 1999 Appl.Opt.38 6705

[15]Cerjan C 2007 J.Opt.Soc.Am.A 24 1609

[16]Yu L,Huang M,Chen M,Chen W,Huang W,Zhu Z 1998 Opt.Lett.23 409

[17]Guizar-Sicairos M,Guti′errez-Vega JC 2004 J.Opt.Soc.Am.A 21 53

[18]You K M,Wen SC,Chen L Z,Wang Y W,Hu Y H 2009 Chin.Phys.B 18 3893

[19]Wang Z,Gao C Q,Xin JT 2012 Acta Phys.Sin.61 124209(in Chi-nese)[王錚,高春清,辛璟燾2012物理學報61 124209]

[20]Gutie’rrez-Vega JC,Rodr′?guez-Masegosa R,Cha’vez-Cerda S 2003 Pure Appl.Opt.5 276

[21]Zhang Q A,Wu FT,Zheng W T,Pu JX 2011 Sci.Sin.Phys.Mech.Astron.41 1131(in Chinese)[張前安,吳逢鐵,鄭維濤,蒲繼雄2011中國科學：物理學·力學·天文學41 1131]

[22]Arfken G B,Weber H J 2001 Mathematical Methods for Physicists(California：Harcourt-Academic)

[23]Wolfram S 1991 Mathematica：A System for Doing Mathematics by Computer(Massachusetts：Addison-Wesley)