孫國棟 秦來安 張巳龍 何楓 譚逢富 靖旭 侯再紅

1)(中國科學(xué)院安徽光學(xué)精密機(jī)械研究所,中國科學(xué)院大氣光學(xué)重點(diǎn)實(shí)驗(yàn)室,合肥 230031)

2)(中國科學(xué)技術(shù)大學(xué)研究生院,合肥 230026)

(2017年9月11日收到;2017年12月4日收到修改稿)

1 引 言

大氣氣溶膠是指由懸浮分布在大氣中的空氣動力學(xué)直徑在0.001—100μm之間的液體或固體顆粒物組成的穩(wěn)定體系.作為大氣的重要組成部分,氣溶膠能通過散射和吸收太陽輻射參與云的形成,直接或間接地影響地球氣候變化.此外,其對能見度、交通安全和人體健康也具有重要影響[1,2].消光系數(shù)作為表征氣溶膠光學(xué)特性的重要參數(shù),是當(dāng)今大氣科學(xué)和環(huán)境科學(xué)等領(lǐng)域的一個熱點(diǎn)研究內(nèi)容.目前,激光雷達(dá)以較高的時空分辨率等優(yōu)勢成為一種探測大氣氣溶膠消光系數(shù)和后向散射系數(shù)等光學(xué)特性的重要技術(shù)手段[3,4].在使用反演算法(Klett后向積分或Fernald后向積分)確定激光雷達(dá)探測范圍內(nèi)消光系數(shù)的分布時首先需要確定最大有效探測距離的消光系數(shù),即消光系數(shù)的邊界值[5?9].確定消光系數(shù)邊界值的常用方法是在對流層頂附近選擇氣溶膠含量最小的高度作為標(biāo)定點(diǎn),假設(shè)該點(diǎn)氣溶膠的散射比為某一很小的定值[10].但對于有效探測距離達(dá)不到對流層頂?shù)牡涂仗綔y激光雷達(dá)而言,這種方法就會失效.為解決該問題,國內(nèi)外學(xué)者相繼提出了一些迭代和估算的方法.Kovalev[11]提出了一種迭代算法,該算法利用大氣氣溶膠和大氣分子的消光系數(shù)定義一個信號校正因子,通過迭代對信號進(jìn)行校正得到大氣氣溶膠消光系數(shù)邊界值.但是該方法校正因子的修正量不易控制,并且迭代次數(shù)較多.王治華等[12]在處理有云層的激光雷達(dá)回波信號時,將激光雷達(dá)信號在大氣中的傳輸過程分為云層區(qū)和非云層區(qū),忽略回波信號中的云層區(qū)信息,利用最小二乘法對非云層區(qū)的激光雷達(dá)回波信號進(jìn)行擬合得到大氣消光系數(shù)的邊界值.但是該方法忽略了真實(shí)大氣是非均勻分布的,導(dǎo)致計(jì)算結(jié)果出現(xiàn)一定的偏差.陳濤等[13]提出了一種散射比迭代法,該方法在探測范圍內(nèi)尋找到一個氣溶膠含量相對較小的高度層,建立回波信號與邊界值之間的等式關(guān)系,利用窮舉法獲得該高度處的氣溶膠散射比并以此來確定邊界值,獲得了較為準(zhǔn)確的反演結(jié)果.但是該方法是等步長迭代,要獲取一個合適的氣溶膠散射比需要較多的迭代次數(shù),計(jì)算量大.

本文提出了基于Broyden算法[14]確定大氣消光系數(shù)邊界值的新方法,把確定大氣消光系數(shù)邊界值的問題轉(zhuǎn)化為求非線性方程數(shù)值解的問題.根據(jù)大氣消光系數(shù)邊界值與激光雷達(dá)回波信號之間的關(guān)系構(gòu)建了一個非線性方程,采用Broyden算法求解該方程的數(shù)值解得到大氣消光系數(shù)邊界值.這種方法具有迭代次數(shù)少,對初始解依賴程度小的優(yōu)點(diǎn).由于這種方法使用的是沿傳輸路徑積分的雷達(dá)信號,所以適用于真實(shí)大氣非均勻的情況.文中使用自行研制的成像激光雷達(dá)對提出的新方法進(jìn)行了驗(yàn)證,并將直接測量得到的大氣透過率與成像激光雷達(dá)最后獲取的大氣透過率進(jìn)行了對比研究.

2 理論和方法

選用Klett方法[5]反演大氣消光系數(shù)時首先需要確定消光系數(shù)邊界值.Klett假設(shè)大氣后向散射系數(shù)β和消光系數(shù)σ之間存在乘冪關(guān)系,即[5,6]

其中B0是常數(shù),k與激光雷達(dá)波長和觀測的氣溶膠性質(zhì)有關(guān),按照其范圍設(shè)定為常用數(shù)值1.因此單次散射激光雷達(dá)方程可以表示為[15]

其中P(z)為從距離z處接收到的散射信號強(qiáng)度,C為系統(tǒng)常數(shù),S(z)為距離平方矯正信號.

方程(2)的穩(wěn)定解為

其中zm為最大探測距離,Sm為zm處的距離平方矯正信號(Sm=S(zm)),σm為zm處的消光系數(shù)(σm=σ(zm)).從方程(3)中可以看出,消光系數(shù)既與相對信號S?Sm有關(guān),也和消光系數(shù)邊界值σm相關(guān).

利用方程(2)可以得到

消光系數(shù)從最小探測距離到最大探測距離的積分可以表示為消光系數(shù)邊界值的函數(shù):

其中z0表示初始距離,

由方程(5)可知消光系數(shù)的平均值主要依賴于信號積分I的幅度.確定消光系數(shù)邊界值的常用方法是把消光系數(shù)邊界值設(shè)定為與平均消光系數(shù)相等來解不等式方程.即

由方程(5),(7)和(8)可以得到以下方程:

將方程(7)代入方程(9)中,令σm=x,設(shè)

即將消光系數(shù)邊界值的確定問題轉(zhuǎn)化為非線性方程數(shù)值解的求解問題.求解非線性方程最常用的方法是牛頓法[16].它具有簡單的迭代形式和二次收斂速度.對于充分接近解的初值,該方法會迅速收斂到精確解.但由于每次迭代都需要計(jì)算導(dǎo)數(shù),導(dǎo)致該方法運(yùn)算量較大.本文中需要求解的方程(10)本身較為復(fù)雜,因此其導(dǎo)數(shù)難以求解,此時牛頓法就不再適用.為了解決這個問題,目前已經(jīng)提出了很多技術(shù)手段,其中最成功的一類為擬牛頓法[17].擬牛頓法利用函數(shù)值逼近導(dǎo)數(shù)值,從而避免了導(dǎo)數(shù)的計(jì)算.函數(shù)在任意點(diǎn)xr處的一組導(dǎo)數(shù)值可以寫成雅可比矩陣的形式:

圖1 Broyden算法流程圖Fig.1.The flow chart of Broyden algorithm.

擬牛頓法提供了一個更新公式,對每次迭代該公式給出了雅可比矩陣的連續(xù)逼近,Broyden等已經(jīng)證明,在特定情況下該更新公式對逆雅可比矩陣也提供了滿意的近似.Broyden算法作為擬牛頓法的一種具體形式,它的算法結(jié)構(gòu)如下[14]:

1)輸入解的初始近似,設(shè)置計(jì)數(shù)器r為零;

2)計(jì)算或假設(shè)逆雅可比矩陣Br的初始近似;

3)計(jì)算Pr=?Brfr,其中fr=f(xr);

4)確定標(biāo)量參數(shù)t,使其滿足‖f(xr+trPr)‖＜‖fr‖,其中符號‖?‖表示待取向量的模;

5)計(jì)算xr+1=xr+trPr;

6)計(jì)算fr+1=f(xr+1),如果‖fr+1‖＜ε(其中,ε是一個很小的預(yù)先設(shè)定的正數(shù)),則退出;否則執(zhí)行第7步;

7)利用更新公式得到雅可比矩陣的近似

其中yr=fr+1?fr;

8)令r=r+1,返回到第3步.

上述算法的結(jié)構(gòu)流程如圖1所示.

3 實(shí)驗(yàn)與系統(tǒng)

3.1 成像激光雷達(dá)系統(tǒng)

成像激光雷達(dá)系統(tǒng)主要由激光器、望遠(yuǎn)鏡、電荷耦合器(CCD)、濾光片和數(shù)據(jù)采集與控制系統(tǒng)組成,其結(jié)構(gòu)和性能在文獻(xiàn)[18]中有詳細(xì)說明.成像激光雷達(dá)系統(tǒng)的實(shí)驗(yàn)測量原理如圖2所示,其中,D為收發(fā)間距,dθ為像元對應(yīng)的張角,γ為CCD傾角.該系統(tǒng)的工作原理可以簡述為:激光器向大氣中發(fā)射一束激光,由于激光與大氣發(fā)生相互作用產(chǎn)生散射光,散射光中包含了大氣分子和氣溶膠的信息,CCD相機(jī)中的像元記錄了光束上各散射高度的散射光子數(shù),即回波信號強(qiáng)度.

圖2 成像激光雷達(dá)系統(tǒng)Fig.2.System of imaging lidar.

3.2 圖像數(shù)據(jù)處理過程

為獲得隨高度變化的回波信號,需對存儲的圖像文件進(jìn)行處理[18?20].處理步驟如圖3所示.經(jīng)過處理得到的回波信號強(qiáng)度個例如圖4所示.然后將各個位置處的回波信號強(qiáng)度以及其他的參數(shù)值代入非線性方程(10)中,運(yùn)用Broyden算法對消光系數(shù)邊界值進(jìn)行求解.由于得到的消光系數(shù)邊界值的準(zhǔn)確性無法定量給出,需要將其轉(zhuǎn)化為可以對比研究的量.本文利用求得的消光系數(shù)邊界值進(jìn)一步求解得到消光系數(shù)的路徑分布,然后按照路徑積分獲得探測范圍內(nèi)的大氣透過率,以大氣透過率作為一個對比量.

圖3 圖像處理過程Fig.3.The procedure of image processing.

圖4 回波信號強(qiáng)度個例Fig.4.Example of return signal.

3.3 對比實(shí)驗(yàn)系統(tǒng)

本文將沒有經(jīng)過反演算法計(jì)算的直接測量的大氣透過率作為對比實(shí)驗(yàn),如圖5所示.直接大氣透過率測量裝置分為發(fā)射裝置和接收裝置.實(shí)驗(yàn)時,發(fā)射裝置將出射光分成功率相等的兩束激光,一束光用于監(jiān)測出射功率,另外一束光到達(dá)接收裝置.接收裝置測量傳輸?shù)浇邮斩说募す馐芰?通過與發(fā)射端的監(jiān)測光束能量對比,獲得大氣透過率.

圖5 1 km實(shí)驗(yàn)地點(diǎn)地圖Fig.5.Experimental site map of 1 km.

4 結(jié)果與分析

4.1 兩種確定消光系數(shù)邊界值方法的對比

圖6 Broyden方法的初值 (a)21:30;(b)22:30Fig.6.Initial value of the Broyden algorithm:(a)21:30;(b)22:30.

利用Broyden方法求解非線性方程確定消光系數(shù)邊界值時首先需要給定初始值,初始值選定的精確與否也決定了其迭代次數(shù).同時也需要給定迭代的容差,本文中設(shè)定為10?6.初始值的選取采用同一天的21:30和22:30兩個時間點(diǎn)獲取的信號來說明,如圖6所示.其中l(wèi)eft表示以等式(9)左邊的項(xiàng)繪出的曲線,right表示以等式(9)右邊的項(xiàng)繪出的曲線,圖6(a)和圖6(b)中的left和right均有交點(diǎn),該交點(diǎn)就是非線性方程(10)的數(shù)值解.從圖中可以得到該數(shù)值解的大致位置,即Broyden算法的初始值.選取不同的初始值得到最后的結(jié)果需要不同的迭代次數(shù),如表1所列.選取的初值越接近方程的解其對應(yīng)的迭代次數(shù)就會越少,但是只要選擇的初值在合理的范圍內(nèi),最后由Broyden算法都能得到最終解且大小相等.這充分說明采用Broyden算法對上述非線性方程式的求解是有效的.

表1 選擇不同初值的迭代結(jié)果Table 1.Iterative results with different initial values.

本文運(yùn)用Broyden算法和最小二乘法分別求取實(shí)驗(yàn)期間內(nèi)大氣消光系數(shù)的邊界值,如圖7所示.兩種方法測量得到的邊界值在趨勢上是一致的,但是在整個測量時間段的前期和后期均具有一定的差異,尤其是在邊界值較大的情況下兩種方法的實(shí)驗(yàn)結(jié)果差異較大,最小二乘法大于Broyden算法求取的邊界值.

圖7 兩種方法測量的邊界值對比Fig.7.Comparison of boundary values measured by two methods.

4.2 不同方法下獲得的大氣透過率對比結(jié)果

通過上述方法獲得了消光系數(shù)邊界值之后,利用激光雷達(dá)方程的穩(wěn)定解(3)式可以得到探測范圍內(nèi)消光系數(shù)的分布.本文分別以21:30和22:30兩個時間點(diǎn)為例,利用兩種方法測量得到1 km以內(nèi)的消光系數(shù)的水平空間分布,實(shí)驗(yàn)結(jié)果如圖8所示.

圖8 水平消光系數(shù)分布 (a)21:30;(b)22:30Fig.8. Distribution of horizontal extinction coeffi-cients:(a)21:30;(b)22:30.

圖8表明消光系數(shù)的水平分布是不均勻的,因?yàn)閷?shí)驗(yàn)場地并不單一,激光光束經(jīng)過的范圍內(nèi)有草地、水庫、馬路等.對探測時間段內(nèi)所有的數(shù)據(jù)進(jìn)行反演計(jì)算得到了其隨時間和水平空間變化的消光系數(shù)分布.由于求得的消光系數(shù)邊界值的不同,造成了兩種方法測量得到的水平消光系數(shù)也有一定的差異.

用兩種方法測量得到的消光系數(shù)分布依據(jù)大氣透過率計(jì)算(12)式可以計(jì)算得到1 km范圍內(nèi)的大氣透過率T隨時間變化的情況:

將這兩種方法與直接大氣透過率測量裝置(3.3節(jié))測得的大氣透過率進(jìn)行對比,結(jié)果見圖9.Direct,Old,New分別表示直接進(jìn)行測量、利用最小二乘法測量以及本文提出的Broyden算法測量的大氣透過率.從對比圖像來看,Broyden算法與直接大氣透過率測量結(jié)果趨于一致,而利用最小二乘法得到的測量結(jié)果偏小.將直接大氣透過率測量結(jié)果作為橫軸,最小二乘法和Broyden算法測量的大氣透過率數(shù)值作為縱軸進(jìn)行相關(guān)性分析,結(jié)果如圖10所示,兩種處理算法與直接測量結(jié)果之間的相關(guān)系數(shù)分別為0.958和0.968.本文提出的方法測量的結(jié)果與直接大氣透過率的測量結(jié)果更為一致.

圖9 不同方法下測量的水平大氣透過率Fig.9.Horizontal transmittance measured by different methods.

圖10 激光雷達(dá)反演大氣透過率結(jié)果直接測量的大氣透過率之間的相關(guān)性Fig.10.The correlation between the transmittance inverted by lidar and direct transmittance results.

4.3 相對誤差分析

為了更進(jìn)一步說明兩種方法的差距,以直接測量得到的大氣透過率作為參考標(biāo)準(zhǔn),將兩種反演方法得到的大氣透過率與參考標(biāo)準(zhǔn)之間進(jìn)行相對誤差分析,如圖11,藍(lán)線是測量時間段內(nèi)所有的相對誤差,紅線是界限誤差,黑線是平均誤差,兩種反演方法與參考標(biāo)準(zhǔn)之間的相對誤差均在20%之內(nèi).以相對誤差值為7.5%作為界限誤差,Broyden算法得到的大氣透過率與參考標(biāo)準(zhǔn)之間的相對誤差在7.5%以上的數(shù)據(jù)占總數(shù)據(jù)的24%,平均相對誤差為4.66%;最小二乘法得到的大氣透過率與參考標(biāo)準(zhǔn)之間的相對誤差在7.5%以上的數(shù)據(jù)占總數(shù)據(jù)的70%,平均相對誤差為9.10%.充分說明了利用本文提出的方法進(jìn)行測量的結(jié)果具有一定的優(yōu)勢.

圖11 不同反演方法與直接測量得到的透過率之間的相對誤差 (a)Broyden方法與直接方法;(b)最小二乘法與直接方法Fig.11.Relative errors of transmittance between different inverse method and direct method:(a)Broyden method and direct method;(b)least square method and direct method.

5 結(jié) 論

當(dāng)牛頓法不適用于求解激光雷達(dá)反演過程中的非線性方程時,Broyden算法可以作為一種恰當(dāng)?shù)姆椒▉硗瓿蛇@一計(jì)算.對于初值的選取只要是在合理的范圍之內(nèi),都能通過有限的迭代次數(shù)獲得一致且較為準(zhǔn)確的解,即消光系數(shù)邊界值.當(dāng)然選取的初值越接近真實(shí)解迭代速度就會越快.通過對比最小二乘法和本文提出的算法計(jì)算得到的大氣透過率,結(jié)果表明消光系數(shù)邊界值求取的準(zhǔn)確性對于最后大氣透過率結(jié)果的準(zhǔn)確性有明顯影響,并且通過與直接測量得到的大氣透過率之間的相關(guān)性和相對誤差分析,進(jìn)一步驗(yàn)證了本文提出的算法相比于最小二乘法有一定的優(yōu)勢.因此本文提出的利用Broyden算法確定消光系數(shù)邊界值的方法是有效的.目前本文中處理的激光雷達(dá)數(shù)據(jù)只考慮了單次散射而未考慮多次散射的影響,因此在以后的工作中將進(jìn)一步加入多次散射的相應(yīng)分析.

[1]Man S W,Kai Q,Hong L,James R C,Kwon H L 2017Atmos.Environ.154 189

[2]Zhao G Y,Liang M,Li Y Y,Duan Z,Zhu S M,Liang M,Sune S 2017Appl.Opt.56 1506

[3]John E B,Sebastian B,Robert B,Parikh N C 2003Appl.Opt.42 2647

[4]Liang M,Mikkel B 2015Opt.Express23 A1613

[5]James D K 1985Appl.Opt.24 1638

[6]James D K 1981Appl.Opt.20 211

[7]Liang M,Peng G,Yang Y,Zheng K 2017Opt.Express25 A628

[8]Cao N W 2015Optik126 2053

[9]Masap M,Nobuo T 1994Appl.Opt.33 6451

[10]Zhou J,Yue G M,Qi F D 1998Chin.J.Quant.Elect.15 140(in Chinese)[周軍,岳古明,戚福第1998量子電子學(xué)報(bào)15 140]

[11]Kovalev V A 1993Appl.Opt.32 6053

[12]Wang Z H,Wang H B,He J,Zheng Y C,Yang J G,Li Y Q,Zhao X B 2008Laser J.29 36(in Chinese)[王治華,王宏波,何捷,鄭玉臣,楊經(jīng)國,李躍清,趙興炳2008激光雜志29 36]

[13]Chen T,Wu D C,Liu B,Cao K F,Wang Z Z,Bo G Y,Yuan L,Zhou J 2010Acta Opt.Sin.30 1531(in Chinese)[陳濤,吳德成,劉博,曹開法,王珍珠,伯廣宇,袁林,周軍2010光學(xué)學(xué)報(bào)30 1531]

[14]George L,John P(translated by Li J,Ren M M)2016Numerical Methods Using MATLAB(Beijing:China Machine Press)pp116–117(in Chinese)[喬治 L, 約翰P著(李君,任明明譯)2016數(shù)值方法-MATLAB版(北京:機(jī)械工業(yè)出版社)第116—117頁]

[15]Ma X M,Tao Z M,Ma M J,Li C J,Wang Z Z,Liu D,Xie C B,Wang Y J 2014Acta Opt.Sin.34 0201001(in Chinese)[麻曉敏,陶宗明,馬明俊,李成軍,王珍珠,劉東,謝晨波,王英儉2014光學(xué)學(xué)報(bào)34 0201001]

[16]Anne G,Timothy P C(translated by Wu Z J,Wang G Y,Fan H J)2016Numerical Methods(Beijing:China Machine Press)p63(in Chinese)[安妮 G,蒂莫西 P C著(吳兆金,王國英,范紅軍譯)2016數(shù)值方法(北京:機(jī)械工業(yè)出版社)第63頁]

[17]Xiong X L,Jiang L H,Feng S,Zhuang Z B,Zhao J Y 2012Infrar.Laser Eng.41 1744(in Chinese)[熊興隆,蔣立輝,馮帥,莊子波,趙俊媛 2012紅外與激光工程 41 1744]

[18]Sun G D,Qin L A,Cheng Z,Hou Z H 2017Laser Optoelect.Prog.54 090102(in Chinese)[孫國棟,秦來安,程知,侯再紅2017激光與光電子學(xué)進(jìn)展54 090102]

[19]Yang C P 2011M.S.Dissertation(Dalian:Dalian Maritime University)(in Chinese)[楊成鵬 2011碩士學(xué)位論文(大連:大連海事大學(xué))]

[20]John E B,Parikh S,Trevor B K 2007Appl.Opt.46 2922