双波长激光雷达反演技术

双波长激光雷达反演技术

John F. Potter

提出一种用于反演激光雷达回波的数据分析技术并用模拟数据进行测试。这种技术需要同时有两个频率上的雷达回波，且满足以下假设：()1在激光光束路径上后向散射与消光的比值不随空间发生改变()2在激光光束路径上，两个不同频率脉冲的消光因子之比也不随空间发生改变。当雷达脉冲路径上大气分子的散射可被忽略，且气溶胶是由同一种粒子组成时，这些假设可以成立。模拟数据对应的是激光雷达探测1km内光学厚度为1的均匀分布的气溶胶的探测数据。通过分析能确定的物理量有透过率T，两个频率的消光因子之比k及两个频率上的消光廓线。这些量中的误差严重依赖于数据中的噪声水平。当用100个激光脉冲平均来减少噪声时，T和k的有效误差分别为1.93和1.54%，消光廓线中的最大误差为6%。附录将会对考虑大气分子三者的情况进行合理的扩展。

β

Ⅰ.介绍()R

α之类的大通过数学方法反演激光雷达回波以获得例如后向散射系数和消光系数()R

气光学特性系数廓线是激光雷达研究中长期存在并亟待解决的问题。如果能够有独立探测雷达脉冲路径上某些点的消光系数或整条路径上的总透过率的探测方法那我们就能够得到反演数据，但是通常我们并没有这样的探测方法。

在没有辅助数据的情况下，为了得到结果许多方法被开发出来。它们都涉及到对大气的物理性质的假设。在时下最流行的斜率法中通过假设α是在一定间隔内围绕雷达脉冲路

α。Klett已经详细地讨论过这种方法及径上距离为R的某一点的常数以估计这一点的()R

其相关的方法。Spinhirn e等人提出了一些稍微不同的方法。在这种方法中雷达回波通过对不同仰角的探测而获得，例如斜径，且这些反演是在基于假设大气在水平方向上是一致的情况下做出的。

上述的方法都依赖于一些假设假设，即假设气溶胶有某一特定密度分布。由于气溶胶的分布可变性很强，因此对不使用这些假设的方法的研究是一热点问题。

本文要介绍的就是这样一种方法，它需要两个频率上的模拟雷达回波，并且基于以下假设：

()a在激光光束路径上后向散射与消光的比值在空间上平均分布。

()b 在激光光束路径上，两个不同频率脉冲的消光因子之比也在空间上平均分布。 当雷达脉冲路径上大气分子的散射可被忽略，且气溶胶是由同一种粒子组成时，这些假设可以成立。这意味着在任意一点气溶胶粒子都有相同的形状，尺度及折射指数。这种方法也能被用于以下方面，例如：把雷达射线分为许多不同的部分，每一部分对应的气溶胶种类都不同，每一部分都可以独立地进行反演。然而，这样一来，结果的精确度便会下降。 本文给出了关于数学公式的推导，并且测试了一种简单的大气气溶胶模式下的模拟数据。附录将会对考虑大气分子三者的情况进行合理的扩展。

Ⅱ.方法

激光雷达方程可以写为：

式中()P R 是在时间t 接收到的瞬时功率，0P 是时间0t 时的发射功率，c 是光速，τ是脉冲宽度，β总的大气后向散射因子，R 是距离，A 是有效接受面积，0T 是光学透过率，α是大气总的消光因子。

通过假设()a ，我们可以得到:

式中S 是一个未知常数，激光雷达方程可以写为。

式中

我们假设我们有两个波长上的雷达回波，我们把其中消光因子较大的波长称为L ，它对应的物理量会带上一个下标L 来表示。同样的我们把其中消光因子较大的波长称为S ，它对应的物理量会带上一个下标S 来表示。

对于波长L ，方程()3可写成：

通过假设()b ，我们对于所有的R 可以得到：

对于波长S ，方程()3可写成：

通过方程()6和()8，我们可以得到：

其中i R 和j R 是任意两个距离。通过求k 的值我们可以得到：

式中k 的下标i 和j 对应的是距离i R 和j R 。

关于方程()6中透过率的解可以由以下方程得出：

其中

式中0R 和m R 分别是待分析的第一点最后一点。L T 是波长为L 时从0R 到m R 的透过率。L T 定

义为：

显然01L T ≤≤。

一般的，当L T 的值没有给出时，我们并不能得到方程()()11,12的解。当我们的假设满足时，可由两个波长的数据资料估计L T 的值。如果数据中没有噪声且方程()12中L T 的值是准确的，我们使用方程()11中的得到的()L R α带入方程()10时，对于任意的i R 和j R 我们都能得到ij k 准确值。然而，若是方程()12中使用的L T 的值并不准确，那么不同的点得到的ij k 的值便会不一样。因此，我们可以通过尝试大量不同的L T ，把其中使不同的（i R ，j R ）

对应的ij k 差异最小的L T 值作为L T 的估值?L T 。我们将会在第三部分介绍一种回归方法来得到

这一估值，这种方法同样用来计算k 的估值。当我们得到L T 和k 时，我们就能通过方程()

11和()12得到波长为L 时的估值?()L r α

，同理利用方程()7我们可以得到估值?()s r α。 Ⅲ.L T 和k 的估值

这一部分将会介绍两种方法来估算L T 和k 的值。

A. 直接法

第一步是选择一个试验值X T 作为L T 的值，选值应尽可能接近真实值。对于每一组（,i j R R ）我们能得到一个ij k ，其中i R ＞j R 。这些可以被认为是一具有N 行（ⅰ=1 ... N ），和N 列（J =1，...，N ）的矩阵M 的对角线以上区域。为了简化后续运算，我们可以用c N 行c N 列的格子来存放ij k 的平均值，这些平均值记为IJ k 作为11N N ?矩阵M '中的元素。其中1/c N N N =。在以下的模拟中，N 为500且c N 为25，矩阵对角线以下区域的元素值为零，因为我们只计算j ＞i 时的ij k 。当然，对角线上的元素值也为零。

如果数据中没有噪声且用于计算ij k 的L T 的值是准确的，那么矩阵M 中所有的ij k 和矩阵M '中所有的IJ k 都是准确的值k 。当X T 的值小于正确值L T ，则当I 和J 的值增加时IJ k

的值将会增加。当L T 偏离正确值时将会影响到方程()11和()12的解，所以X T 的值小于L T 时会给我们的计算带来麻烦。当的X T 值大于L T 时也会存在类似的问题，即当I 和J 的值增加时IJ k 的值会减少。此外，大量的样本的计算结果表明， IJ k 的值或多或少沿着一条垂直于对角线的任意直线。

鉴于此问题，适当的表达IJ k 的变化的方式是把IJ k 用()/2I J +回归。当用于计算IJ k 的X L T T =是正确值时，回归线的斜率应为零。当X T ＜L T 时，斜率大于零，当X T ＞L T 时，斜率小于零。通过这种方法我们能不断地选择更加合适的X T 的值并通过插值法得到L T 的估

计值。把这个估值记为?L T ，k 的估计值?k 是?L

T 对应的IJ k 的平均值。 B ．最快的方法

这种方法的计算速度比方法一更快，并能得到相似的结果。方法二中首先得到距离为R 上间隔为c N 的平均值()L L R 和()s L R 。1L L 是()L L R 在第一个c N 点的平均值，2L L 是下一个c N 点的平均值，以此类推。用这种方法我们定义1N 个平均值()1,...,1LM L M N =，类似地我们定义1N 个平均值()1,...,1SM L M N =。接下来和方法一的步骤一样，只不过用LM L 和SM L 替换()L L R 和()s L R ，然后我们能在不计算N N ?矩阵M 的情况下直接得到11N N ?矩阵M '。由于这种方法计算要快得多，我们接下来的计算结果都是用方法二计算得出的。

Ⅳ.计算中的误差

()L R α的计算中，误差主要来源于()L L R 和L T ,()L R α的方差可以近似地表达为：

其中()L L R σ是()L L R 的标准差，()L T R σ是L T 的标准差，2

L L L T C 是变量()L L R 和L T 的

协方差。若L T 的值是由所有R 上的()L L R 和()s L R 计算得出的话，那么2

L L L T C 可以忽略，

因此()L L R 在任意的R 上对L T 的影响都可以忽略不计。

方程()14中的导数项为：

其中

用2

()L R α分裂方程()14并利用方程()()15~17我们可以得到：、

我们可以注意到方程()18右边的两项对应了产生误差根本的两部分。第一项对应的是()L L R 产生的误差，对于一个信号脉冲来说它以随机地沿着光束线或正或负地变化。 方程()18右边第二项对应的是与L T 有关的系统误差()/()(,)/L L L L L d R R F R T dT T αα=-，对于一个特定的脉冲这一误差随着R 系统地变化，并且所有R 对应的误差符号一般是相同的。当对于所有满足条件的的R 和()0,01.L L T R Rm T ≤≤≤≤(),F R T 的值都大于零时，误差的符号与L dT 的符号相反。另

外，若是忽略大气的散射作用，由于(),L F R T 的绝对值会随着R 单调递增，因此误差的绝对值也随着R 单调递增。

Ⅴ. 模拟数据

使用用模拟数据来测试这一方法。L 和S 对应的波长分别是0.53和1.06m μ（使用倍频钕激光）。激光雷达的系统参数已经在Table Ⅰ中给出了。我们假设0250,1250m R m R m ==，且在0R 和m R 之间()L R α和()s R α是常数。波长为L 时总的光学厚度取为1.0，对应在所述间隔内每0.001m 的()L R α的值。k 的取值为0.5。在这一波长内，这是大陆霾类型下气溶胶的合理取值。在在0R 和m R 之间的任意一点上，对于两个波长的在激光光束路径上后向散

射与消光的比值我们定为0.03，后向散射因子为()511310L m sr β---=?和

()5111.510S m sr β---=?。

通过这些参数，方程()1被用于计算由检测器的每个积分周期发射的电子的平均数目，相关的噪声的标准偏差可通过这一电子平均数的平方根计算得出。因为信号中的光电倍增管放大期间会产生额外噪声，所以以这种方式计算的噪声还需乘以噪声系数1.5。由背景辐射和暗电流产生的噪声由于被视为模-数转换产生的噪声而被忽略。噪声进使用一个随机数发生器行模拟。由于在激光雷达运用中对大量的激光脉冲平取均值以削减噪声是一种常用的手段，模拟数据使用了四种处理手段分别对1，4,25,100道激光脉冲取平均值。在一道以上的激光脉冲平均中，脉冲实际上并不产生和平均。相反，由于噪声由因子1/2N 而减小，因此对一道激光脉冲取平均是由噪声的标准差产生的，其中N 是激光脉冲平均次数。当不同激光脉冲之间没有相关性时，1/2N

是数据平均产生的理论预测的降噪因子，这是我们这一模

式中的假设条件。

图1 L波段的相对标准差，曲线上的数字指的是进行平均的脉冲数。

图2 S波段的相对标准差，曲线上的数字指的是进行平均的脉冲数。

在图1,2中分别给出了在不同的平均次数下()L L R 和()S L R 的相对标准差()RSDs 。由于在所用的气溶胶模式下，激光雷达回波功率()P R 随着R 的增加而减少且()RSDs 随着()1/2P R 变化，所以相对误差随着R 增加而增加。

Ⅵ.结果

在上述四个方案中，我们对每一个方案都用第二部分所介绍的方法进行了五十次模拟

和反演并取其平均值。每次反演都得出了一个L T 的估计值?L

T （修正后为0.368）和一个k 的估计值?k （修正后为0.5）。设i ε是L

T 第i 次反演的百分误差()1,...,50i =。记ε的平均值和标准差为ε和S ε。对应的均方根误差记作2E ε ，表示为：

类似的，记i δ为k 的百分误差，δ，S δ和E δ分别对应平均值，标准差和均方根误差。T able Ⅱ中列出了四种方案中的这些量的值，并列出了它们的协方差2/L kT L C kT 。在所有方案的t 测试中他们的平均值在95%的置信水平下和零相比并不显著。因此假设平均值的真值是零时均方根误差能最好地估计标准差的真值。

图3,4 给出了相对误差()/()L L R R αδα和()/()L L R R αδα。其中曲线是由方程()()18,19计算得出的，而点是由每种方案的五十次模拟反演运算中取十个R 对应的均方根误差计算得来的。100脉冲平均的点没有给出来，因为它与理论曲线完全重合。方案一和方案四种点和线的差异情况不同是因为方程()()18,19中只包括一阶条件，因此，在这一错误下并不期望有精确的值。在这它们似乎很快地把误差提升到了大约10%。从图3可以看到，

在我们上文中讨论过的由()L L R 的误差产生的随机误差和由L T 的误差产生的系统误差随R 增长而增长的同时，()L R α的误差也随着R 增长而增长。在图4于中可以看出()S R α的

图3.0.53m λμ=时的总误差的相对标准差曲线是由方程()18计算得出的，而点是由模拟数据计算得出的。曲线上的数字指的是进行平均的脉冲数。

误差有着一样的规律。值得注意的是在图三图四中对于任意的R 误差都近似1/2S N - 。

VII.讨论和结论

我们介绍了一种不需要对大气参数进行独立观测的双频率激光雷达回波反演方法。对一个特定模式的模拟显示这种方法的雷达回波信号对噪声十分敏感。但是通过数据平均削减噪声后后，我们获得了绝佳的探测结果。例如，当我们对100个脉冲进行平均时，参数L T 和k 的均方根误差小于2%，在两个波长上激光光束路径上所有点的消光因子的均方根误差都小于6%。

然而，我们的模式作了很多的近似处理，主要由以下几点：

()a 忽略大气分子的瑞利散射。

()b 在激光光束路径上后向散射与消光的比值不随空间发生改变。

()c 在激光光束路径上，两个不同频率脉冲的消光因子之比也不随空间发生改变。 ()d 在激光光束路径上消光因子是一个常数。

()e 忽略由电子热量，背景辐射，电子放大器和模-数转换产生的误差。

()f 两个波长上的系统增益都能用表达为已知的关于距离的函数。

在所用的这个模式中，消光因子足够大，因此瑞利散射可以被忽略。然而当气溶胶的后向散射和消光因子不远大于大气分子的气溶胶的后向散射和消光因子时，这种方法就会失效。在我们所用的波长中，当波长为1.06m μ时瑞利散射可以被忽略，但是0.53m μ波长在气溶胶含量较低的情况下瑞利散射将不能被忽略。

正如我们在第一部分提到的，当激光光束路径上大气分子的散射可被忽略，且气溶胶是由同一种粒子组成即在任意一点气溶胶粒子都有相同的形状，尺度及折射指数时，条件()b 和()c 就能被满足。一般来说，除了雷达光束是水平照射的情况下以上条件都能被满足。 我们之所以在这个模式中做出假设()d 是因为它是可能的假设中最简单的，并在大多数情况下适用。然而，只要气溶胶密度不会低到使散射截面与对应的瑞利散射截面和信噪比相比过小，那么我们即使在气溶胶密度与本模式相差很大的模式下也能得到与本模式类似的结果。

一般来说，传感器热噪声，背景辐射，模拟器噪声和数字化噪声可以通过选择合适的设备来削减到可以忽略的水平。此外我们还能使用多脉冲平均来消除此误差。

假设()f 考虑到许多激光雷达系统的增益随时间变化这一事实，因此用R 修正方程一中的因子21/R .由于这种方法依赖于不同距离的雷达回波之比，他们的增益必须得是稳定且已知的以使回波能被适当地修正。我们并不需要知道增益的准确值，但是所有波长上不同R 上的相对增益必须得是已知的。

由假设不当造成影响可以用更加复杂的模拟来研究，事实上更多的模拟已经在计划中了。然而，对于这个或其他任何反演方法来说，真正的考验是实际数据的测试。我们希望能在不久的将来尝试实际数据。

图4. 1.06m λμ=时的总误差的相对标准差曲线是由方程()18计算得出的，而点是由模拟数据计算得出的。曲线上的数字指的是进行平均的脉冲数。

在此我要感谢M. Bristow, D. Bundy, J.McElroy, 和E. Richardson 对我的帮助。本研究是由与美国内华达州拉斯维加斯环境保护局68-03-3250和68-03-3245合同资助的。

虽然我们提到这个研究的资金的全部或部分是由与美国环境保护局提供的，它并没有受到美国环境保护局审查，但这并不意味着该机构的审查没有官方的认可。

附录A ：考虑瑞利散射的情况

附录中将讨论本文所介绍的反演技术在瑞利散射不可被忽略时的运用。当瑞利散射截面减弱为波长的第四功率时，瑞利散射对于较短的波长来说会变得更重要。通常情况下瑞利散射在波长较长时可以忽略然而在波长较短时不得不被考虑。但是这个问题在两个波长都需要考虑瑞利散射时反而容易解决得多。我们在附录I 中将处理前一种情况，附录II 中将处理后面一种情况。附录中所提到的方法只是一个建议，并没有通过测试，所以不能保证这是唯一的解决方案。从理论上来看，在噪声足够小的情况下这种方法是有效的。

I.瑞利散射在波长较长时可以忽略而在波长较短时需要考虑

我们假设较短的波长是文章中的波长L ，我们通常用可见或近红外光，因为L 是强吸收波段，而一般来说在这些波段中大气气溶胶的消光随波长增加而减少。

当我们考虑瑞利散射是，方程()6变成：

其中(),L R R α和(),L R R β分别是波长为L 时的瑞利散射和消光系数。其他量则和之前的定义一样。

由方程()8可得：

其中i R 和j R 是任意两个距离值，正如我们之前假设的j i R R ＞。从方程()2A 我们可以得到：

其中，和方程()10一样，我们用给k 加下标i 和j 来说明对应的是哪个距离R 的值计算而来的。可以注意到当对所有波长上的瑞利散射都不能被忽略时，我们可以用方程()3A 来代替方程()10，类似的，通过与方程()10一样的方法推导出来的方程也能用这个方法处理但是并不成功，因为它比方程()3更难处理。

当波长L 包括瑞利散射时方程()11和()12可以用下列公式代替：

其中

这里的L T 可通过方程()13得出。

如果我们知道L S 的值，我们就能使用无瑞利散射的方法得出结果。如果我们不知道L S 的值，那么对L T 和L S 都采用不同的初始假设值来获得这些参数和k 的最优估计值。 II.两个波长都要考虑瑞利散射

在这里我们用方程()1A 来代替方程()6，并用下式代替方程()8：

求比值()()/i j s R s R L L ，我们可以得到：

其中

之后我们就能按照之前介绍的方法进行计算了，只不过要考虑到三个变量L T ，L S 和S S 。对于给定的L T 和L S 我们可以通过解方程()4A 和()5A 得到()L R α。通过这个和S S 的假设值我们就能通过计算得到,,ij ij i A B c 和i d 并解出方程()7A 中的k ，把其标记为ij k 。正如我们在第三部分所说的，把能使ij k 的差异最小的L T ，L S 和S S 的取值即是我们要得到的解。 参考文献

1. R. T. H. Collis and P. P. Russell, in Laser Monitoring of the Atmosphere, E.

D. Hinkely, Ed. (Springer-Verlag, New York,1976), p. 117.

2. F. G. Fernald, B. M. Herman, and J. A. Reagan, "Determination of Aerosol Height Distributions by Lidar," J. Appl. Meteorol. 11,482 (1972).

3. J. D. Klett, "Stable Analytical Inversion Solution for Processing Lidar Returns," Appl. Opt. 20, 211 (1981).

4.J. D. Spinhirne, J. A. Reagan, and B. M. Herman, "Vertical Distribution of Aerosol Extinction Cross Section and Interference of Aerosol Imaginary Index in the Troposphere by Lidar Technique," J. Appl. Meteorol. 19, 426 (1980).

5. P. R. Bevington, Data Reduction and Error Analysis for the Physical Sciences (McGraw-Hill, New York, 1969).

6. R. H. Kingston, Detection of Optical and Infrared Radiation (Springer-Verlag, New York, 1978).

7. R. W. Engstrom, Photomultiplier Handbook (RCA Corp., Electro Optics and Devices, Lancaster, PA 17604, 1980).