梯度傾斜相關測量水平C和橫向風速廓線的理論與仿真研究?

彭哲靖旭 侯再紅 吳毅

1)(中國科學技術大學環境科學與光電子技術學院,合肥 230026)

2)(中國科學院安徽光學精密機械研究所,中國科學院大氣成分與光學重點實驗室,合肥 230031)

彭哲1)2)靖旭2)?侯再紅2)吳毅2)

1)(中國科學技術大學環境科學與光電子技術學院,合肥 230026)

2)(中國科學院安徽光學精密機械研究所,中國科學院大氣成分與光學重點實驗室,合肥 230031)

(2016年12月1日收到;2017年3月8日收到修改稿)

根據Rytov近似以及泰勒湍流凍結假設,推導出以不同距離的前向散射光為信標的水平路徑上梯度傾斜角的相關表達式.基于該表達式,在理論上提出了計算湍流強度與橫向風速的新方法,并通過數值仿真對該方法進行了初步驗證.結果表明,在5%高斯誤差情況下,大氣折射結構常數和風速的計算結果與理論真值在整體變化上具有較好的一致性,線性相關系數分別能達到0.8與0.9.該方法能夠得到不同湍流與風速條件下的湍流強度廓線以及風速廓線,為反演大氣湍流強度以及風速提供了一種新思路.

∶梯度傾斜,相關,大氣湍流,風速

PACS∶42.68.—w,42.68.Bz,42.25.DdDOI∶10.7498/aps.66.104207

1 引 言

光在大氣中傳播會由于大氣湍流的影響而產生隨機相位抖動和光強起伏,嚴重影響了天文觀測以及激光工程的應用.而湍流折射率結構常數C(h)以及風速v(h)則是與大氣湍流相關的重要參數.自適應光學系統能夠實時校正大氣湍流的影響,但是其校正能力受湍流的時間特征頻率的影響,實時獲取湍流特征頻率需要知道大氣風場和湍流強度廓線.因此,實時獲取大氣湍流廓線以及風速廓線對光在大氣中傳播具有重要意義.

目前,測量大氣湍流廓線的方法有探空氣球法[1]、星光法[2,3]和激光雷達法[4?7].探空氣球法操作簡單,但容易受到風的影響,實時性不高.星光法,如SCIDAR(SCIntillation Detection And Ranging),PML(Profiler of Moon Limb)等,有較好的實時性和精度,但容易受天氣的影響.激光雷達法能夠測量各種不同路徑上的湍流廓線,其測量結果的空間分辨率與儀器自身的空間分辨率相關.測量風速廓線的方法有探空氣球、微波測風雷達、基于多普勒效應的風廓線雷達、基于抖動和閃爍的激光雷達[8,9].探空氣球測量簡單,成本低,但測量周期長,實時性不高.微波測風雷達技術成熟,受天氣影響較小,但回波信號在大氣干燥時偏弱.利用多普勒效應測風精度高,一般采用法布里標準具,受溫度影響較大,且其成本高.基于抖動和閃爍的激光雷達方法大多都是測量路徑上的平均風速,想要獲取風速廓線通常需要模型假設.

本文采用激光雷達的方法,通過質心偏移得到梯度(G)傾斜角相關特性,并利用矩陣變換與近似的方法獲取水平湍流強度廓線和風速廓線.該方法是對現有激光雷達方法對C測量的補充,以及對風速隨路徑分布測量的改進.該方法具有較好的實時性,且可以實現不同路徑上的測量.

2 理論推導

如圖1所示,以圖中z軸所示方向為正方向,O為坐標原點.將S(source z=L0)處的激光器產生的激光在不同距離上的前向散射光作為信標光(圖中F處),并在O處用雙孔望遠鏡采用傍軸方式接收該激光的前向散射,通過在兩孔前安裝有一定楔角的楔鏡和在焦平面傾斜安裝電荷耦合器(CCD),激光信號就能在傾斜的焦平面成像為兩條光柱[5,10].

對光柱上對應F(forwad scattering z=L)處的前向散射有

其中Ri,j(t1,t2,L)表示孔徑i(i=1,2)在t1時刻與孔徑j(j=1,2)在t2時刻接收到的質心抖動相關.從S到F段,同一時刻,激光自身的質心抖動以及傳輸路徑上的質心抖動對兩個孔徑而言都是一樣的,故有(L,t)=(L,t).

梯度傾斜(G傾斜)即平均光線方向,與孔徑的平均相位梯度相關[11].實際測量中,G傾斜角可由質心在CCD靶面上的偏移與成像焦距之比求得.故質心抖動相關與G傾斜角相關呈正比關系,比值與焦距有關.由Rytov近似以及泰勒湍流凍結理論,有[9?13]

圖1 測量原理圖Fig.1.Sketch map of measurement.

其中k0為波數,z表示到望遠鏡的距離,κ表示空間波數,Φ(κ,z)表示湍流譜分布,γ=1?z/L表示傳播因子,余弦函數對應相位抖動的衍射因子,J0表示第一類零階貝塞爾函數,τ表示t1與t2的時間差,F(γκ)表示濾波函數.相對于探測孔徑,光源擴展度遠小于衍射極限,因而可不考慮光源的濾波函數,即F(γκ)僅包含孔徑濾波函數.對應的孔徑濾波函數[11,13]可表示為

其中D表示孔徑直徑,J1為第一類一階貝塞爾函數.

而對于Kolmogorov湍流,有

其中B(L,τ)即為求解湍流折射率結構常數以及風速所需相關量,它對應于單個孔徑G傾斜角的自相關與兩個孔徑G傾斜角的互相關之差;WC(z,L,τ)表示對應的湍流折射率結構常數的權重.由(6a)式以及(2a)和(2b)式可知,B(L,τ)不受光源自身質心抖動以及S到F段質心抖動影響.根據B(L,0)以及WC(z,L,0)即可通過合適的算法求解得到湍流強度沿路徑的分布.為了得到風速,可以通過求解B(L,τ)在τ=0時刻的導數.求導可得

其中,G(L)為B(L,τ)在τ=0時刻的導數,WCv(z,L)為湍流折射率結構常數與風速乘積的權重.通過聯立(6a)和(6b)以及(7a)和(7b)式,可求解得到風速隨路徑的分布.而事實上對時間求導通常需要較高的時間分辨率,對于一個200 Hz的CCD,在不考慮外來誤差影響的情況下,通過泰勒展開求取導數,其相對誤差不會超過1.5%.

若以后向散射作為信標光,則后向散射光(從信標光到望遠鏡)和原始激光(從激光器到信標光)會存在一定的相關性,并不能認為是零.若要計算風速與則需要消除該相關性所帶來的影響.可以通過差分的手段消除原始激光帶來的影響,進而進行相關與求導運算得到與風速的相關表達式及權重函數.

3 算法原理

在不考慮盲區影響的情況下,我們將O(z=0處)到S(z=L0處)看作n段,即l1l2···ln,第一段為l0→ l1,第n段為ln?1→ ln,其中l0=0.并認為每一段的風速與都保持不變,即v1v2···vn與表示第i段的風速,表示第i段的湍流折射率結構常數.令

其中B(li,0)表示將li處的前向散射光作為信標光且時間差為時的差分相關量,G(li)對應B(li,τ)在τ=0時的導數.當孔徑間距、孔徑直徑以及激光波長已知,即可計算得到bi,j與gi,j.對于(9a)與(9b)式,直接通過最小二乘法求解會導致很大誤差.因為矩陣本身是病態的,同時誤差的存在會導致最優解并不是我們想要的結果,往往要加以約束.

取接收孔徑間距ρ=0.235 m,孔徑直徑D=0.12 m,激光波長λ=532 nm,對于不同的L,兩權重函數的曲線如圖2所示

圖2 (網刊彩色)權重隨距離z的變化曲線 (a)WC的變化曲線;(b)WCv的變化曲線Fig.2.(color online)The variation of weight function with z(a)variation of WC;(b)variation of WCv.

從圖2可以看出,當z接近L時,即γ接近0時,權重在量級上快速減小.權重的這種快速下降正是導致矩陣病態的原因.由于這部分權重在量級上遠小于其他部分,因而可以近似認為其值為0.基于這一點,并根據每個分段內權重曲線變化特性對上述矩陣進行近似變換可以得到

α和β為與每一段上曲線變化相關的量.

4 數值仿真結果與分析

為便于計算,我們將距離間隔設置為固定值50 m.而實際測量時,空間分辨率會隨著傳播距離的增加逐漸降低,要得到50 m間隔的相關量,可先對實際觀測結果進行降噪擬合,然后對所得曲線進行間隔為50 m的采樣.同時由于上述各種近似都會引入誤差,隨著誤差傳遞,當累積到一定程度時,則需要增大分段間距以降低該誤差帶來的影響.我們取L0=4000 m,此時,在50 m的分段間距下,近似誤差給計算帶來的影響很小.

而實際上誤差總是存在的.若每一幀的光柱抖動都含有方差5%的高斯噪聲,則對計算得到的B以及通過泰勒展開得到的G進行多幀平均,并對其進行小波降噪和多項式擬合,得到結果如圖4所示.

從圖4可以看出計算得到的結果與設定值有著相同的變化趨勢,但在細節變化上有所丟失,這是降噪以及擬合的平滑效果所致,同時也是由于抖動對的不敏感性.的最大對數相對誤差不超過3.4%,平均對數相對誤差不超過0.9%,風速的最大絕對誤差不超過1.82 m/s,平均絕對誤差不超過0.47 m/s.當湍流與風速條件發生變化時,依然有類似的結果,如圖5所示.

圖5(b)給出的風速結果誤差較大,是由于當風速存在方向變化且變化較快時,多項式擬合誤差相對較大.兩組計算結果一定程度上驗證了該方法在理論上的可行性,更深入的研究需要進一步的實驗驗證.

圖3 (網刊彩色)不考慮誤差的情況下設定值與計算值的對比 (a)的對比;(b)風速的對比Fig.3.(color online)Comparison between the set value and the calculated value regardless of errors(a)comparison of ,(b)comparison of wind velocity.

圖4 (網刊彩色)考慮誤差的情況下設定值與計算值的對比 (a)的對比;(b)風速的對比Fig.4.(color online)Comparison between the set value and the calculated value considering errors(a)comparison of ,(b)comparison of wind velocity.

圖5 (網刊彩色)湍流與風速條件發生變化時設定值與計算值的對比 (a)的對比;(b)風速的對比Fig.5.(color online)Comparison between the set value and the calculated value within change of turbulence and wind(a)comparison of ,(b)comparison of wind velocity.

5 結論與討論

本文基于數值仿真提出了一種基于G傾斜角相關測量水平橫向風速和湍流強度的方法.由于以后向散射光作為信標光測量風速系統較為復雜,因而采用前向散射光作為信標光,并根據統計與相關理論,推導出了相關抖動與湍流強度以及風速之間的表達式,最終通過矩陣變換與近似來實現湍流強度與風速的計算.

仿真結果表明,對于水平測量,在不考慮盲區或盲區較小且不考慮光束擴展的情況下,該方法能有效地反演出橫向風速與湍流強度隨路徑的變化,而不僅僅只是路徑上的一個平均量,且不依賴于模型假設.該方法仍需要通過實驗來進行驗證,并通過改進擴展以實現垂直方向的反演.

[1]Vernin J,Munoz-Tunon C 1994 Astron.Astrophys.284 311

[2]Avila R,Cuevas S 2009 Opt.Express 17 10926

[3]Ziad A,Blary F,Borgnino J,Fante?-Caujolle Y,Aristidi E,Martin F,Lantéri H,Douet R,Bondoux E,Mékarnia D 2013 Astron.Astrophys.559 L6

[4]Cheng Z,TanfF,Jing X,He F,Hou Z H 2016 Acta Phys.Sin.65 074205(in Chinese)[程知,譚逢富,靖旭,何楓,侯再紅2016物理學報65 074205]

[5]Jing X,Hou Z H,Wu Y,Qin L A,He F,TanfF 2013 Opt.Lett.38 3445

[6]Cui C L,Huang H H,Mei H P,Zhu W Y,Rao R Z 2013 High Power Laser Part.Beams 25 1091(in Chinese)[崔朝龍,黃宏華,梅海平,朱文越,饒瑞中 2013強激光與粒子束25 1091]

[7]Cheng Z,He F,Jing X,TanfF,Hou Z H 2016 Acta Opt.Sin.36 401004(in Chinese)[程知,何楓,靖旭,譚逢富,侯再紅2016光學學報36 401004]

[8]Tichkule S,Muschinski A 2012 Appl.Opt.51 5272

[9]Schock M,Spillar E J 1998 Opt.Lett.23 150

[10]Jing X,Hou Z H,Qin L A,He F,Wu Y 2011 Infrared Laser Eng.40 1352(in Chinese)[靖旭,侯再紅,秦來安,何峰,吳毅2011紅外與激光工程40 1352]

[11]Sasiela R J 2007 Electromagnetic Wave Propagation in Turbulence:Evaluation and Application of Mellin Transforms(Bellingham:SPIE)pp69–102

[12]Wang T,Clifford S F,Ochs G R 1974 Appl.Opt.13 2602

[13]Rao R Z 2012 Modern Atmospheric Optics(Beijing:Science Press)pp368–433(in Chinese)[饒瑞中 2012現代大氣光學 (北京:科學出版社)第368—433頁]

[14]Cai D M,Ti P P,Jia P,Wang D,Liu J X 2015 Acta Phys.Sin.64 224217(in Chinese)[蔡冬梅,遆培培,賈鵬,王東,劉建霞2015物理學報64 224217]

[15]Guo Y M,Ma X Y,Rao C H 2014 Acta Phys.Sin.63 069502(in Chinese)[郭友明,馬曉燠,饒長輝2014物理學報63 069502]

[16]Huang K T 2014 Ph.D.Dissertation(Hefei:Hefei Institutes of Physical Science,the Chinese Academy of Sciences)(in Chinese)[黃克濤 2006博士學位論文(合肥:中國科學院合肥物質科學研究院)]

PACS∶42.68.—w,42.68.Bz,42.25.DdDOI∶10.7498/aps.66.104207

*Project supported by the National Natural Science Foundation of China(Grant No.41405014).

?Corresponding author.E-mail:xujing@aiofm.ac.cn

Simulation research and theoretical study on measurement of atmospheric optical turbulence and wind profile using the correlation of gradient-tilt?

Peng Zhe1)2)Jing Xu2)?Hou Zai-Hong2)Wu Yi2)
1)(Department of Environmental Science and Optoelectronic Technology,University of Science and Technology of China,Hefei 230026,China)
2)(Key Laboratory of Atmospheric Composition and Optical Radiation,Anhui Institute of Optics and Fine Mechanics,Chinese Academy of Sciences,Hefei 230031,China)

1 December 2016;revised manuscript

8 March 2017)

In this article,a theoretical method based on thefluctuation of gradient tilt(G-tilt)of active light source is proposed to estimate the horizontal profiles of atmospheric optical turbulence()and transverse wind.The G-tilt,related to the average phase gradient,is in the same direction as the average ray direction.And G-tilt angle is considered to be equal to the ratio between the centroid position offset and the focal length.In this method,a theoretical model based on lidar system is set up,in which forward scatter light beams at different distances are taken as beacons.These beacons are detected by a two-aperture telescope.And two light columns,from which we can obtain the information about G-tilt angle,are imaged by these beacons.In order to obtain the turbulence intensity and wind velocity from G-tilt angle with our theoretical model,the differential cross-correlation expressions of G-tilt angle and its derivative are derived in detail.These two expressions are based on the spatial cross-correlation function obtained from Rytov approximation and Taylor’s frozen- flow hypothesis for Kolmogorov turbulence.Simultaneously,path weighting functions ofand wind velocity are derived,and the effects of path weighting functions on the calculation of our method are analyzed.Based on such an analysis,to realize the inversion of turbulence intensity and transverse wind,the matrix transformation algorithm is proposed.We ignore some minimal values of the path weighting functions in our algorithm so that the ill-conditioned matrix is avoided.Besides,numerical simulation is used for preliminarily validating this method.In our simulation,varies randomly between 10?15m?2/3and 10?14m?2/3while wind velocity ranges from ?5 m/s to 10 m/s.The sign of the wind velocity represents the direction of wind.According to the simulation results,the horizontal profiles of atmospheric optical turbulence and transverse wind calculated are consistent with their theoretical values no matter whether there exists Gaussian noise.When the ratio between the standard deviation of Gaussian noise we added and the original signal is 0.2,the maximum relative error of logarithmicis 3.4%,and the correlation coefficient between the calculated results and theoretical values foris 0.8.Besides,the maximum absolute error of wind velocity is 1.82 m/s,and the correlation coefficient between the calculated results and theoretical values for wind velocity is 0.9.Even if the horizontal profiles of atmospheric optical turbulence and transverse wind vary largely,the calculation results of our method remain stable.Therefore,a new idea is provided for measuring atmospheric turbulence and wind.

∶gradient-tilt,correlation,atmospheric turbulence,wind velocity

?國家自然科學基金(批準號:41405014)資助的課題.

?通信作者.E-mail:xujing@aiofm.ac.cn

?2017中國物理學會Chinese Physical Society