激光与光电子学进展, 2019, 56 (24): 240102, 网络出版: 2019-11-26   

一种确定消光系数边界值的新算法 下载: 1143次

Method for Determining Boundary Value of Extinction Coefficient
作者单位
1 大连海事大学信息科学技术学院, 辽宁 大连 116026
2 大连市气象装备保障中心, 辽宁 大连 116001
摘要
在反演大气消光系数过程中,一个值得探讨的问题就是消光系数边界值的选择问题。提出了一种基于横向Steffensen型3阶方法求解大气气溶胶消光系数边界值的新算法,并将该方法应用于仿真信号和实际观测到的回波信号中求解大气消光系数,证明了所提算法的真实可靠性。结果表明,所提算法迭代速度快,经过较少的迭代次数就可收敛到消光系数的边界值,能够较准确地反演大气消光系数。
Abstract
During the inversion of atmospheric extinction coefficient, choosing the boundary value of the extinction coefficient is critical. In this study, a new method based on lateral Steffensen-like third-order method is proposed for determining the boundary value of atmospheric aerosol extinction coefficient. To confirm the method reliability, the proposed method is applied to both simulated and actual measured echo signals. Results show that the proposed method has high iterative speed and converges to the boundary value of extinction coefficient after a few iterations, which can invert atmospheric extinction coefficient more accurately.

1 引言

激光雷达是一种将短激光脉冲以一个或多个波长发射到大气中,通过与大气相互作用后测量回波信号的主动式遥感仪器,可以探测气溶胶消光系数、后向散射系数等参数的时空分布。激光雷达已被证明是提供大气边界层(ABL)结构以及连续观测云和气溶胶属性的强大工具之一,也是大气环境的监测手段之一,可以进行连续监测,具有高时空分辨率[1-3]。沈吉等[4]利用拉曼散射激光雷达探测了平流层气溶胶的消光系数廓线,分析了平流层大气气溶胶消光系数的变化特征。项衍等[5]通过讨论气溶胶光学特性反演方法的流程,分析了影响反演结果准确度的各种不确定性因素。常用的大气消光系数反演算法有Collis斜率法[6]、Klett法[7]和Fernald法[8]。当假设大气为均匀分布时,可以使用Collis斜率法对消光系数进行反演。但实际上均匀分布的大气状况不会出现,所以大多采用Klett法或Fernald法求解消光系数。对于这两种算法,不同时间和不同天气条件下消光系数的边界值是一个不确定的值,故精确反演气溶胶消光系数的关键环节就是边界值的选取。1993年,Kovalev[9]提出了一种适用于双组分散射大气的激光雷达信号迭代反演方法,该方法对原始信号进行了变换,使得单组分大气的激光雷达方程求解成为可能。2010年,Marchant等[10]利用迭代最小二乘法从弹性激光雷达数据中得到了有用的信息,该方法能处理任意数量的通道,并在参考点以上和低信噪比的情况下生成稳定的解。2012年,熊兴隆等[11-12]分别通过不动点迭代法和改进的牛顿法求解消光系数的边界值。李红旭等[13]利用弦截法得到了消光系数的初始边界值,然后对边界值进行迭代修正。陈涛等[14]将激光雷达距离校正回波信号在一段范围内的积分值和此范围内的大气透过率作为判断依据,利用穷举法反复迭代确定合适的边界值,但该方法的取值步长固定,迭代速度较慢,运算量大。

本文提出了一种基于横向Steffensen型3阶方法反演消光系数边界值的新方法。首先选取参考点,结合激光雷达方程和Fernald反演算法构建关于消光系数边界值与消光系数均值的非线性方程,然后利用一种横向Steffensen型3阶方法进行求解,最终得到较准确的消光系数边界值。利用本文方法对仿真和实测回波信号进行消光系数的反演,结果表明,本文方法可以得到较准确的消光系数。本文方法不需求导,迭代次数较少,收敛速度可提高到3阶。

2 仪器介绍

本文使用的激光雷达为第二代CAMLTM CE 370-2后向散射激光雷达。对于一个低电平雷达信号的精确探测来说,可通过光子计数系统进行信号的采集,雷达的最大垂直分辨率为15 m。表1为激光雷达的主要技术参数指标。

表 1. CE 370-2系统的主要技术参数指标

Table 1. Main technical parameters of CE 370-2 system

ItemParameter
Transmitter typeNd∶YAG
Wavelength /nmLaser power /mWLaser output energy /μJ53250-1008-20
Repetition frequency /kHzPulse width /nsTelescope diameter /mm4.7<15200
Range resolution /mTelescope field of view (FOV) /μradDetection mode1555/110Photon counting

查看所有表

3 大气消光系数反演算法

激光雷达发射的532 nm的激光处在可见光范围内,激光束通过大气中的分子和气溶胶时会受到散射和吸收,接收系统收集回波信号并进行处理,以获得所需信息。Fernald反演算法[8]将大气散射分为分子散射和气溶胶散射两部分,因此激光雷达接收到的距离z处的回波信号可由激光雷达方程描述为

P(z)=ECz-2[βa(z)+βm(z)]Ta2zTm2(z),(1)

式中:P(z)为从高度z处接收到的回波信号功率;E为输出能量;C为仪器校准常数;βa(z)和βm(z)分别为高度z处气溶胶和分子的后向散射系数;Ta(z)=exp -0zαa(z)dz为仅考虑气溶胶衰减的大气透过率;Tm(z)=exp -0zαm(z)dz为分子的大气透过率;αa(z)和αm(z)分别为高度z处气溶胶和分子的消光系数。

由激光的雷达方程可知,大气气溶胶消光系数和后向散射系数是不确定的,所以激光雷达比是一个必要的参数[15]。对于气溶胶,Sa(z)=αa(z)a(z)。Sa(z)是一个随高度变化的参数,为了简化操作,一般将其假设为不随高度变化的量,本文令Sa(z)=50 sr。对于大气分子,Sm(z)=αm(z)m(z)=8π/3, 其中Sm(z)为大气空气分子消光后向散射比[16]

根据Fernald反演算法,低于边界高度zc处的气溶胶消光系数为

αa(z)=-SaSm·αm(z)+X(z)·exp2SaSm-1zzcαm(z)dzX(zc)αa(zc)+SaSm·αm(zc)+2zzcX(z)·exp2SaSm-1zzcαm(z')dz'dz,(2)

式中:X(z)=P(zz2为距离平方校正信号;zc为参考高度;αa(zc)为zc处的气溶胶消光系数,即边界值,通常选取无云的清洁大气层所处的位置,但激光在大气中传播时会受到不同程度的衰减,地基激光雷达一般达不到这样的高度,可以在激光雷达最大探测距离范围内选取X(z)m(z)最小值所在的高度作为参考高度[17]

4 消光系数边界值的确定方法

4.1 横向Steffensen型三阶方法求解非线性方程

求解形如f(x)=0,x∈R的非线性方程的常用方法有牛顿法、弦截法和抛物线法。在牛顿法中用差分形式代替导数可得到Steffensen法,Steffensen法避免了复杂的求导运算,同时还有2阶的收敛速度。文献[ 18]提出了一种横向Steffensen型3阶方法,用差分代替求导,可将收敛速度提高到3阶。

xk处将f(xk+1)进行泰勒展开,得到

f(xk+1)=f(xk)+f'(xk)Δk+f(xk)2Δk2++f(p-1)(xk)2Δkp-1+O(Δkp),(3)

式中:Δk=xk+1-xk;f'(xk),f″(xk),…,f(p-1)(xk)分别表示函数f(x)在xk处的一阶、二阶,…,(p-1)阶导数;k表示迭代次数;O( Δkp)表示高阶无穷小。

f(xk)+f'(xk)Δk+f(xk)2Δk2=0,(4)

将(4)式中的 Δk2f(xk)f'(xk)2代替,可得到方程

f(xk)+f'(xk)Δk+f(xk)2f(xk)f'(xk)2=0,(5)

求解(5)式得到新的迭代公式为

xk+1=xk-1+12Lf(xk)f(xk)f'(xk),(6)

其中,

Lf(x)=f(x)f(x)f'x2(7)

f~'(xk)=f[xk-2f(xk)]-4f[xk-f(xk)]+3f(xk)2f(xk),(8)f~(xk)=2f[xk,xk-f(xk),xk-2f(xk)],(9)

f(x)的一阶导数和二阶导数用(8)~(9)式[18]代替,并代入(6)~(7)式即可得到横向Steffensen型3阶方法。这种方法不用计算导数值,不仅加快了收敛速度,还增强了解的稳定性。本文将采用这种横向Steffensen型3阶方法进行消光系数边界值的求解。

4.2 消光系数边界值的确定算法

首先构建关于消光系数边界值的非线性方程[19],然后利用横向Steffensen型3阶方法得到迭代公式,求解此方程就可解得消光系数的边界值。算法流程图如图1所示,其中K为最大迭代次数,ε为允许误差。

图 1. 迭代算法流程图

Fig. 1. Flow chart of iterative algorithm

下载图片 查看所有图片

1) 构建非线性方程。首先在激光雷达最大有效探测距离范围内选取X(z)m(z)最小值所在的高度作为参考点zc,设该处气溶胶消光系数的边界值为x,在zc附近选取一小段范围(zb,zc),该范围内的大气消光系数几乎不发生剧烈变化,令zb=zc-10·Δz,(zb,zc)范围内的消光系数均值为

α-(x)=1ni=1nαa(zi,x),(10)

式中:n=floor zc-zbΔz,floor(·)表示向下取整;n为将(zb,zc)按Δz分成的等间隔数,Δz为距离分辨率;zi=zbz·i。则构建的方程为

f(x)=x-α-(x)=x-1ni=1nα(zi,x)(11)

2) 确定迭代公式。

xk+1=xk-1+12Lf(xk)f(xk)f'(xk)(12)

3) 选取适当的迭代初值。

4) 用(12)式生成迭代序列{xk}。

5) 判断迭代是否停止,判断条件为|xk+1-xk|+|f(xk)|<ε[18]k<K

迭代过程中当xkxk+1差值的绝对值与f(xk)的绝对值之和小于允许误差ε时,迭代停止,将xk+1作为消光系数边界值的近似值。ε与想要达到的精确度有关,ε越小,得到的边界值越精确,若迭代次数达到K时还没有找到满足精度的xk,则迭代失败。文中设置K=1000,ε=10-3

5 实验验证及分析

采用CAMLTM CE 370-2后向散射激光雷达系统参数得到的仿真信号和实际测得的回波信号对本文方法进行验证。分别使用本文算法(算法一)、文献[ 13]中提到的弦截法(算法二)和文献[ 11]中的不动点迭代法(算法三)获得大气气溶胶消光系数的边界值,然后将其代入Fernald反演算法中求解大气消光系数,最后对这三种算法得到的结果进行比较。实验中边界高度zc为在激光雷达最大探测距离范围内选取X(z)m(z)最小值所在的高度。评价指标为AARE= 1mj=1mαj-α'jαj,其中αj为消光系数的理论值,α'j为消光系数的实际反演值,j在[1,m]之间取值,m表示在一段连续距离中以分辨率15 m为间隔划分的个数[20]图2为基于CE 370-2激光雷达仿真的回波信号采用上述三种算法反演得到的消光系数廓线对比结果。

图 2. 仿真信号及反演结果。(a)激光雷达仿真信号;(b)算法一的反演结果;(c)算法二的反演结果;(d)算法三的反演结果

Fig. 2. Simulation signals and inversion results. (a) Lidar simulation signal; (b) inversion result of method 1; (c) inversion result of method 2; (d) inversion result of method 3

下载图片 查看所有图片

图2(b)~(d)中可以看出三种算法的反演结果与理论值较为一致。三种算法设置的最大迭代次数K均为1000,ε=10-3,Sa(z)=50 sr,边界高度zc选取10.0050 km,zb=9.8550 km,边界值为0.00018 km-1。算法一设置的迭代初始值为0.4 km-1,迭代3次得到的消光系数的边界值为0.00032 km-1;算法二设置的迭代初始值为0.4 km-1和0.5 km-1,迭代7次得到的消光系数的边界值为0.00032 km-1;算法三设置的迭代初始值为0.02 km-1,迭代90次得到的消光系数的边界值为0.00039 km-1。算法一和算法二得到的消光系数曲线与理论值之间的AARE均为7.89%,算法三得到的结果与理论值之间的AARE为6.97%。表2为三种反演算法的结果比较。可以看出:虽然算法一与算法二得到的边界值近似相等,但算法一的收敛阶提高到3阶,而算法二的收敛阶仅为1.618,并且算法一的迭代次数比算法二的迭代次数少一半左右;算法一的AARE略低于算法三,但是算法一的迭代次数远小于算法三。综合考虑迭代次数和AARE后,认为本文算法得到的消光系数边界值较准确,迭代次数更少,具有明显优势。

表 2. 三种反演算法的结果比较

Table 2. Inversion results of three inversion algorithms

Algorithmα(zc) /km-1AARE /%Iteration times
10.000327.893
20.000327.897
30.000396.9790

查看所有表

2015年11月10日,当日空气质量为中度污染,空气质量指数为160,PM2.5的质量浓度为120 μg/m3。利用算法一对10日11:16—11:53实际测得的激光雷达回波信号进行反演,得到的消光系数边界值为0.00816 km-1,设置的迭代初值为0.2 km-1,迭代次数为3次;算法二得到的大气气溶胶消光系数边界值为0.00816 km-1,设置的迭代初值为0.2 km-1和0.3 km-1,迭代次数为5次;算法三得到的大气气溶胶消光系数边界值为0.00825 km-1,设置的迭代初值为0.2 km-1,迭代次数为82次。三种算法中设置K=1000,ε=10-3,Sa(z)=50 sr,边界高度zc选为3.4950 km,zb=3.3450 km。将参考文献[ 21]中提及的洁净层法得到的消光系数边界值0.00047 km-1作为理论值。图3(a)为激光雷达回波信号,图3(b)~(d)为三种反演算法的结果比较,消光系数边界值取小数点后5位。从图3(b)~(d)可以看出,三种算法的消光系数廓线与理论值较相似,经计算,算法一和算法二反演得到的消光系数分布曲线和理论值之间的AARE均为21.51%,算法二和理论值之间的AARE为21.77%。三种反演算法的对比如表3所示。通过迭代次数和AARE指标可以看出:算法一的AARE略优于算法三,且算法一的迭代次数远小于算法三;虽然算法一与算法二的AARE近似相等,但算法一的收敛速度提高到3阶,算法二的收敛阶仅为1.618,算法一得到的消光系数更准确且迭代次数较少。

图 3. 2015年11月10日实测信号及反演结果。 (a)激光雷达回波信号;(b)算法一的反演结果;(c)算法二的反演结果;(d)算法三的反演结果

Fig. 3. Measured signals and inversion results on November 10, 2015. (a) Lidar echo signal; (b) inversion result of method 1; (c) inversion result of method 2; (d) inversion result of method 3

下载图片 查看所有图片

表 3. 三种反演算法的结果比较

Table 3. Inversion results of three inversion algorithms

Algorithmα(zc) /km-1AARE /%Iteration times
10.0081621.513
20.0081621.515
30.0082521.7782

查看所有表

表4给出了本文算法(算法一)设置不同初始值时的迭代过程数值结果。

表 4. 数值迭代结果

Table 4. Results of numerical iteration

Iteration timesIteration result /km-1
x0=0.1 km-1x0=1.0 km-1
10.011062993506260.21087944416133
20.008161010205090.02457958037373
30.008160824555750.00819052938697
4-0.00816082455596
5-0.00816082455575

查看所有表

表4可知:初始值为0.1 km-1时,迭代3次得到边界值为0.00816 km-1(取小数点后5位);初始值为1.0 km-1时,迭代5次得到边界值为0.00816 km-1。虽然这两个迭代过程的初始值不同,但迭代结果相等,且都能在很少的迭代次数下收敛到边界值,说明本文算法计算得到的消光系数边界值与设置的迭代初始值无关,最终能得到较准确的消光系数边界值。

6 结论

本文提出了一种基于横向Steffensen型3阶方法反演大气气溶胶消光系数边界值的新算法。利用横向Steffensen型3阶方法求解消光系数边界值的非线性方程,得到消光系数边界值。利用该方法计算得到的消光系数边界值与设置的迭代初始值无关,不需要求导运算,且能在一定条件下将收敛速度提高到3阶。通过对仿真信号和实际观测得到的回波信号进行反演,证实了本文方法的可靠性。将本文算法与利用弦截法和不动点迭代法反演得到的消光系数曲线进行比较,结果表明:本文算法的AARE略优于不动点迭代法,迭代次数远小于不动点迭代法的迭代次数;本文算法与弦截法的AARE近似相等,但本文算法的收敛速度提高到3阶,而弦截法的收敛阶为1.618。本文算法的收敛速度更快,得到的消光系数曲线更准确,更具优势,可以应用到实际回波信号消光系数的反演中。

参考文献

[1] Yan Q, Hua D X, Wang Y F, et al. Observations of the boundary layer structure and aerosol properties over Xi'an using an eye-safe Mie scattering lidar[J]. Journal of Quantitative Spectroscopy and Radiative Transfer, 2013, 122: 97-105.

[2] 梅亮. 沙氏大气激光雷达技术及其研究进展[J]. 激光与光电子学进展, 2018, 55(9): 090004.

    Mei L. Atmospheric Scheimpflug lidar technique and its progress[J]. Laser & Optoelectronics Progress, 2018, 55(9): 090004.

[3] 虞历尧, 曹念文, 沈吉. Rayleigh-Raman-Mie激光雷达探测南京北郊对流层气溶胶光学特性[J]. 激光与光电子学进展, 2018, 55(11): 110102.

    Yu L Y, Cao N W, Shen J. Detection of tropospheric aerosol optical properties by Rayleigh-Raman-Mie lidar in the northern suburb of Nanjing[J]. Laser & Optoelectronics Progress, 2018, 55(11): 110102.

[4] 沈吉, 曹念文. 拉曼散射激光雷达反演平流层气溶胶消光系数廓线[J]. 中国激光, 2018, 45(6): 0609002.

    Shen J, Cao N W. Inversion of stratospheric aerosol extinction coefficient profile by Raman scattering lidar[J]. Chinese Journal of Lasers, 2018, 45(6): 0609002.

[5] 项衍, 刘建国, 张天舒, 等. 激光雷达探测气溶胶光学特性的不确定性因素研究[J]. 激光与光电子学进展, 2018, 55(9): 092801.

    Xiang Y, Liu J G, Zhang T S, et al. Uncertainty factors of aerosol optical properties inversion by lidar[J]. Laser & Optoelectronics Progress, 2018, 55(9): 092801.

[6] Collis R TH, Russell PB. Lidar measurement of particles and gases by elastic backscattering and differential absorption[M] ∥Hinkley E D. Laser monitoring of the atmosphere. Topics in applied physics. Berlin, Heidelberg: Springer, 1976, 14: 71- 151.

[7] Klett J D. Stable analytical inversion solution for processing lidar returns[J]. Applied Optics, 1981, 20(2): 211-220.

[8] Fernald F G. Analysis of atmospheric lidar observations: some comments[J]. Applied Optics, 1984, 23(5): 652-653.

[9] Kovalev V A. Lidar measurement of the vertical aerosol extinction profiles with range-dependent backscatter-to-extinction ratios[J]. Applied Optics, 1993, 32(30): 6053-6065.

[10] Marchant C C, Moon T K, Gunther J H. An iterative least square approach to elastic-lidar retrievals for well-characterized aerosols[J]. IEEE Transactions on Geoscience and Remote Sensing, 2010, 48(5): 2430-2444.

[11] 熊兴隆, 蒋立辉, 冯帅, 等. 基于不动点原理的大气气溶胶消光系数边界值确定方法[J]. 光电子·激光, 2012, 23(2): 303-309.

    Xiong X L, Jiang L H, Feng S, et al. Determination of the boundary value of atmospheric aerosol extinction coefficient based on fixed point principle[J]. Journal of Optoelectronics·Laser, 2012, 23(2): 303-309.

[12] 熊兴隆, 蒋立辉, 冯帅, 等. 改进的牛顿法确定大气消光系数边界值[J]. 红外与激光工程, 2012, 41(7): 1744-1749.

    Xiong X L, Jiang L H, Feng S, et al. Using improved Newton method to determine the boundary value of atmospheric extinction coefficient[J]. Infrared and Laser Engineering, 2012, 41(7): 1744-1749.

[13] 李红旭, 常建华, 朱玲嬿, 等. 基于微脉冲激光雷达的能见度反演算法[J]. 电子测量与仪器学报, 2017, 31(10): 1555-1560.

    Li H X, Chang J H, Zhu L Y, et al. Visibility inversion algorithm based on micro pulse lidar[J]. Journal of Electronic Measurement and Instrumentation, 2017, 31(10): 1555-1560.

[14] 陈涛, 吴德成, 刘博, 等. 低层大气中确定气溶胶后向散射系数边界值的新方法[J]. 光学学报, 2010, 30(6): 1531-1536.

    Chen T, Wu D C, Liu B, et al. A new method for determining aerosol baekscatter coefficient boundary value in the lower troposphere[J]. Acta Optica Sinica, 2010, 30(6): 1531-1536.

[15] 宋秀瑜, 曹念文, 杨思鹏. 探究影响南京地区大气气溶胶光学特性反演的因素[J]. 激光与光电子学进展, 2017, 54(4): 040101.

    Song X Y, Cao N W, Yang S P. Influence factors on atmospheric aerosol optical property inversion in Nanjing[J]. Laser & Optoelectronics Progress, 2017, 54(4): 040101.

[16] Sasano Y. Tropospheric aerosol extinction coefficient profiles derived from scanning lidar measurements over Tsukuba, Japan, from 1990 to 1993[J]. Applied Optics, 1996, 35(24): 4941-4952.

[17] 陈莎莎, 徐青山, 徐赤东, 等. 基于微脉冲激光雷达计算整层大气气溶胶光学厚度[J]. 光学学报, 2017, 37(7): 0701002.

    Chen S S, Xu Q S, Xu C D, et al. Calculation of whole atmospheric aerosol optical depth based on micro-pulse lidar[J]. Acta Optica Sinica, 2017, 37(7): 0701002.

[18] Candela V, Peris R. A class of third order iterative Kurchatov-Steffensen (derivative free) methods for solving nonlinear equations[J]. Applied Mathematics and Computation, 2019, 350: 93-104.

[19] 熊兴隆, 蒋立辉, 冯帅, 等. 非线性方程法确定低空探测机载激光雷达消光系数边界值[J]. 光电子·激光, 2012, 23(7): 1356-1362.

    Xiong X L, Jiang L H, Feng S, et al. Constructing and solving the nonlinear equation of airborne lidar for determining the boundary value of the extinction coefficient for atmospheric aerosol in lower atmosphere[J]. Journal of Optoelectronics·Laser, 2012, 23(7): 1356-1362.

[20] 熊兴隆, 冯帅, 蒋立辉, 等. 一种新的大气消光系数边界值确定方法[J]. 光电子·激光, 2011, 22(11): 1699-1705.

    Xiong X L, Feng S, Jiang L H, et al. A novel method for determining the boundary value of the atmospheric extinction coefficient[J]. Journal of Optoelectronics·Laser, 2011, 22(11): 1699-1705.

[21] 周军, 岳古明, 戚福第, 等. 大气气溶胶光学特性激光雷达探测[J]. 量子电子学报, 1998, 15(2): 140-148.

    Zhou J, Yue G M, Qi F D, et al. Optical properties of aerosol derived from lidar measurements[J]. Chinese Journal of Quantum Electronics, 1998, 15(2): 140-148.

陈晓楠, 毕京平, 王凯欣, 韩冰, 柳云雷. 一种确定消光系数边界值的新算法[J]. 激光与光电子学进展, 2019, 56(24): 240102. Xiaonan Chen, Jingping Bi, Kaixin Wang, Bing Han, Yunlei Liu. Method for Determining Boundary Value of Extinction Coefficient[J]. Laser & Optoelectronics Progress, 2019, 56(24): 240102.

本文已被 2 篇论文引用
被引统计数据来源于中国光学期刊网
引用该论文: TXT   |   EndNote

相关论文

加载中...

关于本站 Cookie 的使用提示

中国光学期刊网使用基于 cookie 的技术来更好地为您提供各项服务,点击此处了解我们的隐私策略。 如您需继续使用本网站,请您授权我们使用本地 cookie 来保存部分信息。
全站搜索
您最值得信赖的光电行业旗舰网络服务平台!