1 引言

激光雷达技术以其高时空分辨率的优势,日益成为研究气溶胶垂直结构和时空演变特征的重要工具。由于发射光束与接收视场不完全重合,激光雷达在探测盲区存在信号丢失、在过渡区存在信号不完全接收的问题,导致对过渡区的气溶胶光学参数反演结果不合理,但近地表气溶胶与人类活动最密切,因此激光雷达必须具备近地表气溶胶的有效探测能力。通过理论计算或实验测量得到激光雷达的几何因子,是获取该区域内气溶胶参数的有效手段之一。

在理论计算研究方面, Stelmaszczyk等^[1]通过分析牛顿望远镜中像的简单模型,推导激光雷达的几何因子;狄慧鸽等^[2]基于数学推导和仿真,给出激光雷达重叠因子表达式,并分析了光束在均匀分布和高斯分布时激光雷达几何因子的变化;汪少林等^[3]引入可完全聚焦填充系数,从几何光学角度推导同轴与离轴激光雷达的几何因子;宋跃辉等^[4]通过分析激光束光轴与望远镜光轴夹角对同轴激光雷达几何因子的影响,提出同轴米氏散射激光雷达准直判据;王威等^[5]提出了基于激光强度分布的激光雷达几何因子计算方法。理论算法虽然可以计算激光雷达几何因子,并指导激光雷达光学系统的优化配置,但它依赖于在实际中很难获得的系统参数,包括激光能量分布、传输特性和光学系统性能等。

在实验测量研究方面,有学者提出了水平测量法、成像法^[6]、振动拉曼-米氏法^[7-10]、纯转动拉曼-米氏法^[11-12],以及星载与地基激光雷达联合测量(即Su Jia的方法)^[13]等。水平测量法需要激光雷达结构具备水平测量能力,并且需要假设水平方向上气溶胶均匀分布,因此在实际应用中具有一定的局限性;成像法对光学系统精确度要求较高且较为复杂;振动拉曼-米氏方法需要假设米氏散射信号接收通道和拉曼信号接收通道的几何因子一致,并且由于拉曼波长与米氏散射波长之间波长差距较大,需要假定气溶胶消光波长指数;纯转动拉曼-米氏方法虽然可以弥补波长指数差异的问题,但同样存在拉曼和米氏通道几何因子一致问题的假设,并且应用范围仅限在纯转动拉曼-米氏激光雷达上。Su Jia的方法,需要利用星载数据反演近地表的消光系数,受制于星载信号在近地面信噪比较低、激光雷达比未知等因素的影响,星载激光雷达(CALIPSO)反演气溶胶消光系数存在一定的误差,致使几何因子反演结果也存在一定误差^[14]。

目前应用最广泛的单一功能米氏散射激光雷达,只能通过水平法确定其几何因子,限制了其针对边界层大气的探测能力。为了解决上述矛盾,本文根据CALIPSO可以完全接收近地表回波信号的特点,联合全接收的CALIPSO探测信号,以及同步观测的、受接收过渡区影响的地基激光雷达信号,提出反演几何因子的新方法。

2 基本原理

根据激光雷达方程,地基激光雷达与CALIPSO的距离校正后的回波信号X_g(r)和X_s(r)可以分别表示为

Xg(r)=c1g(r)β(r)exp(-2τ0r),(1)Xs(r)=c2β(r)exp[-2η(r)τrrc],(2)

式中:r为以地面为原点的垂直高度;c₁,c₂分别为地基激光雷达和CALIPSO的系统常数;g(r)为地基激光雷达的几何因子,由于CALIPSO位于太阳同步极地轨道上,因此在低层大气中几何因子为常数,并包含在系统常数中;β(r)为高度r处的总后向散射系数,包括大气分子后向散射系数与气溶胶后向散射系数;η(r)为多次散射因子。为了减弱粒子多次散射对测量结果的影响,选择大气相对清洁、透过率高的状况测量,此时认为η(r)=1;τ_0r为从地面至测量高度r处的光学厚度;τ_rrc为从测量高度r至CALIPSO校准高度r_c处的光学厚度,其表达式为

τ0r=∫0rα(r)dr,(3)τrrc=∫rrcα(r)dr,(4)

式中α(r)为高度r处的总消光系数。

若不校准地基激光雷达几何因子,则地基激光雷达几何因子影响区域内的距离修正信号可写为

Xg(r)=c3β'(r)exp(-2τ'0r),(5)

式中c₃为地基激光雷达系统不校准几何因子时的系统常数;β'(r)为不校准几何因子得到的总后向散射系数;τ'_0r为从地面至测量高度r处对应的光学厚度。由(1)式和(5)式可以得到校正几何因子前后,后向散射系数β'(r)与β(r)之间的关系为

cβ(r)g(r)=β'(r)exp(-2τ'0r)exp(-2τ0r)。(6)

对于CALIPSO,其探测方向与地基激光雷达相反,距离修正信号X_s(r)可由CALIPSO数据产品Level-1B中获得,它是经过校准和距离修正后的衰减后向散射回波信号。将(1)式和(2)式相乘,星载与地基距离修正信号乘积可表示为

XsrXg(r)=c1c2g(r)β2(r)exp[-2(τ0r+τrrc)],(7)

式中τ_0r+τ_rrc表示地面至CALIPSO校准高度的光学厚度,为常数。设常数c₄=exp[-2(τ_0r+τ_rrc)], (7) 式化简为

XsrXg(r)=c1c2c4g(r)β2(r)。(8)

由(8)式可以看出,地基激光雷达几何因子g(r)可由星载与地基距离修正回波信号表示,但实际上无法得到几何因子校准后的后向散射系数β(r)。将(8)式中地基距离修正信号X_g(r)由不校准几何因子计算的后向散射系数β'(r)表示,同时两边乘以g(r),可以得到

g(r)XsrXg(r)=g(r)c3β'(r)exp(-2τ'0r)c2β(r)exp(-2τrrc)。(9)

根据(6)式中β(r)与β'(r)、τ_0r与τ'_0r之间的关系将(9)式中几何因子转化为由未校正几何因子的后向散射系数β'(r)表示的形式,即

g(r)=c2c3c4β'2rcXsrXgr·exp(-4τ'0r)exp(-4τ0r),(10)

式中g(r)为联合星载与地基激光雷达反演的最终几何因子。但其中exp(-4τ_0r)为未知数,因此首先假设(10)式中后半部的比值为常数,即

exp(-4τ'0r)=c5·exp(-4τ0r),(11)

则几何因子可以近似表示为

g'(r)=kβ'2rXsrXgr,(12)

式中g'(r)为近似几何因子,可以在几何因子为1处选择干净大气,由大气分子后向散射系数得到;k为常数;β'(r)由地基弹性散射回波信号通过Fernald法确定。根据(12)式确定近似几何因子g'(r)后,先校正地基激光雷达距离修正信号得到经g'(r)校正后的后向散射系数β(r),计算从地面至探测高度处的透过率exp(-2τ_0r),再代入(10)式确定最终几何因子g(r)。

3 实验结果与分析

3.1 星载与地基激光雷达系统

CALIPSO有效载荷共包含3部分,周期为16 d,其中正交极化云-气溶胶激光雷达(CALIOP)包含1064 nm和具有偏振探测能力的532 nm通道,能提供气溶胶、云后向散射及退偏比的垂直剖面信息等。地基激光雷达(BIT-Lidar,北京理工大学激光雷达)位于北京市海淀区北京理工大学(39.9° N,116.3° E),是一台纯转动拉曼-米氏散射激光雷达,包含532 nm弹性米氏散射通道和两个高低阶纯转动拉曼通道(正斯托克斯533.7 nm、535.1 nm,反斯托克斯530.3 nm、529.0 nm),并可根据探测要求,在同轴与离轴两种接收模式之间进行切换,CALIOP与BIT-Lidar系统均采用Nd∶YAG激光器,其他相关参数如表1所示。

3.2 结果验证与数据分析

CALIPSO每月经过BIT-Lidar所在地附近4次。为了验证本方法反演几何因子的可行性与准确性,分别在BIT-Lidar离轴与同轴两种情况下,利用BIT-Lidar 实际弹性散射回波信号与CALIPSO的532 nm距离修正信号,反演地基激光雷达几何因子。同时,与Su Jia的方法反演出的几何因子作对比。另外,纯转动拉曼-米氏方法由于不受振动拉曼-米氏方法中波长指数假设的制约,反演激光雷达几何因子较为成熟,可靠性较高。因此,由BIT-Lidar中的米氏信号与纯转动拉曼信号,根据纯转动拉曼-米氏方法^[12]确定系统几何因子,并与本方法反演结果进行比对。

图1为2016年8月25日晚CALIPSO(黑线)与BIT-Lidar(点线)的相对距离修正信号,当晚BIT-Lidar为离轴模式。为了提高星载数据信噪比,在数据处理过程中将200条回波廓线累加平均,同时对星载与地基回波信号进行插值使垂直分辨率统一为15 m,对地基信号进行背景去噪处理。为了验证当晚地基与CALIPSO探测的气溶胶空间分布较为均一,将BIT-Lidar回波信号构造成与CALIPSO探测方向相同的距离修正信号廓线,并与CALIPSO的距离修正信号廓线进行直接对比。图1中的虚线为用BIT-Lidar回波信号构造的与CALIPSO探测方向相同的距离修正信号^[15],可以看出转化后的地基信号与星载信号整体结构趋于一致,且总回波信号与构造的大气分子回波信号(红线)比较接近。这表明当天对流层中的高层无云层和明显的气溶胶层,气溶胶浓度相对较低。但在1.5 km以下,星载信号反映出近地面存在气溶胶边界层结构,而地基回波信号由于受到几何因子的限制,随着高度的降低,强度急剧衰减,二者在近地面偏差较大。

在反演几何因子时,根据(12)式,利用地基弹性散射回波信号,由Fernald法反演后向散射系数,但由于气溶胶雷达比垂直廓线是未知的,需要假设当天的激光雷达比S_a,这将使几何因子的反演结果产生一定的不确定性。图2是利用8月25日星载与地基回波信号,通过假设不同的激光雷达比,得到的几何因子反演结果。可见,激光雷达比对几何因子的影响在过渡区比较明显,且不正确的激光雷达比会使几何因子在充满区出现单调递减或递增等不合理情况。为了解决这个问题,根据Tao等^[16]的方法,结合星载与地基弹性散射回波信号,确定当天最佳的激光雷达比。但近地面信号因未经几何因子校正,且考虑到地基激光雷达在对流层高层的信噪比会下降,在计算不同雷达比下的气溶胶后向散射系数与参考气溶胶后向散射系数廓线差值的平方和时,最低高度应选择在几何因子充满区以上,本文选择计算高度为3~9 km,最终确定当天最佳雷达比为61 sr。

图 1. 2016年8月25日晚相对距离修正信号

Fig. 1. Relative range-corrected signals on the evening of August 25, 2016

下载图片 查看所有图片

图 2. 不同激光雷达比反演的几何因子廓线

Fig. 2. Geometric form factor profiles retrieved by different lidar ratios

下载图片 查看所有图片

图3中点线为计算的近似几何因子g'(r),实线为计算的最终几何因子g(r)。二者在1.8 km以上基本重合,该区域几何因子修正前后透过率的比值近似为常数。1.8 km以下透过率修正前后的几何因子之间只存在微小的改变,表明透过率的修正对几何因子的结果影响较小。由激光雷达方程可以看出,几何因子主要作用在后向散射系数上,透过率是对后向散射系数取负指数值,所以削弱了透射率对几何因子的影响。本文方法反演的离轴几何因子到达充满区后,在1附近有微小波动,这可能是由星载与地基信号的信噪比不一致造成的^[17-18]。

图4中实线为联合CALIPSO与BIT-Lidar反演的离轴模式几何因子廓线,其中将达到充满区以后的几何因子设为常数1,点线为利用纯转动拉曼-米氏法确定的几何因子,两者基本吻合,Su Jia方法反演的几何因子(虚线)达到充满区距离较低。可以看出,BIT-Lidar离轴系统在0.6 km以下的几何因子几乎为0,对应几何因子的盲区,在0.6~2.3 km之间的几何因子急剧增大,在2.3 km附近达到1,以上则为充满区。

图 3. 近似几何因子g'(r)与最终几何因子g(r)

Fig. 3. Approximate geometric form factor profile g'(r) and final geometric form factor profile g(r)

下载图片 查看所有图片

图 4. 三种方法反演的离轴几何因子

Fig. 4. Off-axial geometric form factor profiles retrieved by three different methods

下载图片 查看所有图片

图5为2017年2月17日晚CALIPSO(黑线)、BIT-Lidar(点线)的相对距离修正信号与转化后的地基距离修正信号(虚线),当晚BIT-Lidar为同轴模式。图6中实线表示的是根据本文方法反演的BIT-Lidar同轴几何因子廓线,当天选择的激光雷达比是55 sr,可以看出同轴几何因子在1.0 km附近变为常数。相比于离轴系统,同轴系统盲区范围小,过渡区小,有利于近地面回波信号的接收。充满区最低高度与拉曼-米氏方法反演的几何因子相比约低0.17 km,Su Jia的方法比拉曼-米氏方法约低0.46 km。

图 5. 2017年2月17日晚相对距离修正信号

Fig. 5. Relative range-corrected signals on the evening of February 17, 2017

下载图片 查看所有图片

图 6. 三种方法反演的同轴几何因子

Fig. 6. Coaxial geometric form factor profiles retrieved by three different methods

下载图片 查看所有图片

运用本文方法确定的激光雷达几何因子,分别校正当天的地基激光雷达回波信号,反演气溶胶后向散射系数。图7给出了2016年8月25日晚大气分子后向散射系数(红线),几何因子校准前气溶胶后向散射系数(蓝线),本文方法(黑线)、拉曼-米氏方法(点线)与Su Jia的方法(虚线)反演的离轴几何因子校准后的气溶胶后向散射系数。结果表明,本文方法校准后的气溶胶后向散射系数在近地面相比于校准前有所提高,可以反演出更低区域的气溶胶后向散射系数。离轴模式时,0.6 km以下为探测盲区。0.6~1.0 km区域内,3种方法反演的几何因子对气溶胶后向散射系数的校正均不完全。1 km以上至充满区高度范围内,本文方法与拉曼-米氏方法校正的气溶胶后向散射系数之间平均相对误差为20.8%,与Su Jia方法相比,在几何因子过渡区范围内,平均相对误差提高了25.4%。

图8为2017年2月17日晚几何因子校正激光雷达回波信号前、后反演的气溶胶后向散射系数。可以看出,本文方法与Su Jia的方法校正的气溶胶后向散射系数较拉曼-米氏方法的结果均整体偏小,探测盲区高度降低至0.18 km,本文方法与拉曼-米氏方法校正不完全范围为0.18~0.25 km。0.25 km至充满区范围内,本文方法与拉曼-米氏方法的平均相对误差为18.6%,与Su Jia方法相比,提高了10.4%。

图 7. 2016年8月25日晚几何因子校准前、后气溶胶后向散射系数

Fig. 7. Aerosol backscatter coefficient of geometric form factor before and after correction on the evening of August 25, 2016

下载图片 查看所有图片

图 8. 2017年2月17日晚几何因子校准前、后气溶胶后向散射系数

Fig. 8. Aerosol backscatter coefficient of geometric form factor before and after correction on the evening of February 17, 2017

下载图片 查看所有图片

4 结论

通过联合CALIPSO 532 nm与地基532 nm弹性散射回波信号,可以反演地基激光雷达几何因子。实验结果表明,本方法反演的同轴和离轴几何因子与纯转动拉曼-米氏法的结果吻合,证明了本文所提方法的可行性,且与Su Jia的方法相比,精度有所提高。在离轴模式下,平均相对误差提高了25.4%;在同轴模式下,平均相对误差提高了10.4%。

该方法具备如下优点:与米氏散射激光雷达水平方法获取几何因子相比,避免了几何因子水平测量法涉及的大气水平均匀假设的不利因素;与Su Jia提出的基于星地同步观测方法相比,减小了利用星载数据反演几何因子过渡区内透过率时引入的误差,提高了反演几何因子的精度;与振动拉曼-米氏方法相比,不需要假设气溶胶消光波长指数;与转动拉曼-米氏方法相比,避免了拉曼和米氏通道几何因子一致问题的假设,且适用于应用范围最广泛的米氏散射激光雷达。

相比于其他方法,这种新的方法可以应用于只含有532 nm弹性散射通道的激光雷达反演几何因子,并利用卫星过境时间规律的特点,对系统进行常规标定,为激光雷达近距离数据的精确反演提供了一种可供参考的方法。