本发明属于航空领域，涉及一种斜程能见度的测量方法，特别是一种利用大气消光系数反演斜程能见度的方法。
背景技术：
能见度是影响航空交通运输最常见、最重要的气象要素。在低能见度天气条件(如阴雨、雾霾、低云等)下，大气在垂直方向的分布并不均匀，斜程能见度的大小决定着飞机能否起飞与降落。而跑道的水平能见度与飞机降落时在决断高度上的斜程能见度往往差别很大。当存在地表雾时，水平能见度很低而斜程能见度较高，若根据水平能见度判断可能造成机场关闭。当存在低层云时，水平能见度较高而斜程能见度较低，若根据水平能见度判断则会给飞机降落带来危险。因此，做好斜程能见度的测量，对于保障民用航空飞行安全具有重要意义。
现有的机场测量仪器只能提供近地面的水平能见度，且多用目测法对能见度进行评估，误差较大。虽然也有利用飞机实测或热气球观测的方法，但其在降雨气象条件下难以使用。大气透过率仪虽然能够测量斜程能见度，但其需将接收机放置在塔上，且随着探测角度的改变，接收机的放置高度也要做出相应改变。
激光雷达是一种主动式的光学遥感设备，具有很高的时空分辨率和测量灵敏度，已经成为测量大气气溶胶消光系数及斜程能见度的重要技术手段。中国民航大学自主研制的直接探测多普勒测风激光雷达可以提供大气反射功率信号，我们可以使用此信号对消光系数进行反演，进而求解出斜程能见度。
在使用klett算法对激光雷达回波信号进行大气气溶胶消光系数反演时，我们需要确定的未知参量有两个，一个是大气消光系数边界值，另一个是大气后向散射消光对数比(即k值)。由于这两个参量对反演结果都具有重要影响，因此人们一直比较关注其求解问题。之前的相关研究都是基于某些假设条件或采用迭代法分别得到这两个值，而本专利将给出一种同步确定大气消光系数边界值和大气后向散射消光对数比的计算方法，从而反演出大气消光系数，继而得出斜程能见度。
技术实现要素：
本发明的目的在于提供一种反演结果准确度高、迭代次数少、收敛速度高、易于实现的利用大气消光系数反演斜程能见度的方法。
本发明解决其技术问题是通过以下技术方案实现的：
一种利用大气消光系数反演斜程能见度的方法，其特征在于：所述反演方法包括如下步骤：
(1)建立模型：
通过klett法和激光雷达方程来建立各点大气消光系数的模型如下：
式中，σ(r)为大气消光系数；rm为激光雷达探测最远距离点(边界值点)的高度；σ(rm)是rm处的大气消光系数值，即大气消光系数边界值；k为大气后向散射消光对数比；
(2)构造方程组：
首先，设大气消光系数边界值σ(rm)为x1，大气后向散射消光对数比k为x2；
然后，根据多普勒测风激光雷达对数距离校正信号来寻找大气消光系数近似不变的垂直均匀层，再令边界值点rm和标定点rb处的消光系数相等，构造第一个方程如下：
式中，rb为标定点距离；
最后，根据激光雷达探测的大气气溶胶光学厚度和太阳光度计反演的大气气溶胶光学厚度近似相等再构造另一方程如下：
式中，r′为积分区域的高度。
(3)求解方程组：
1)、利用粒子群算法对上述方程组进行求解，同时求解出大气消光系数边界值和大气后向散射消光对数比:
首先，对回波信号进行2000次累加平均和基于经验模式分解的预处理；
然后，由垂直探测的回波确定出标定点rb＝4.62km及边界值点rm＝4.8km；
最后，通过粒子群算法求解出非线性方程组的解，即消光系数边界值σ(rm)和后向散射消光对数比k；
2)、根据适应度函数求解σ(rm)及k的最优值：
根据适应度函数f(x)的收敛变化情况，可知f(x)迅速下降，在经12次迭代后可达到精度要求，即f(x)最趋近于零，通过迭代可得出σ(rm)＝0.026和k＝0.92；
(4)求解斜程能见度：
将上一步求得的大气消光系数边界值以及大气后向散射消光对数比带入到建立好的各点大气消光系数模型中可求解出大气消光系数均值σp，再把大气消光系数均值带到如下方程中即可得到斜程能见度，
式中，所述的vs为斜程能见度；σp为大气消光系数均值；λ为激光波长；q为修正系数，当能见度较好时(v＞50km)，q＝1.6；当能见度适中时(6km＜v＜50km)，q＝1.3；当能见度较差时(v＜6km)，q＝0.585v1/3。
而且，所述σ(rm)及k的最优值还能通过k的取值范围、斜率法预算的边界值和预设的迭代初始值，来加快迭代速度，减少迭代次数。
本发明的优点和有益效果为：
1.本利用大气消光系数反演斜程能见度的方法，能够同时确定大气消光系数边界值和大气后向散射消光对数比，从而能够求解任意角度的斜程能见度，保障民用航空飞行的安全。定量的分析结果表明用该方法确定的消光系数边界值进行斜程能见度反演，得到的结果具有较高准确性，较少的迭代次数和较高的收敛速度。
2.本发明设计科学合理，具有反演结果准确度高、迭代次数少、收敛速度高、易于实现的优点，是一种具有较高创新性的利用大气消光系数反演斜程能见度的方法。
附图说明
图1为本发明在2月10日多普勒测风激光雷达在探测角为90°时所得到的回波信号图(纵坐标表示距离r，横坐标表示距离回波信号功率对数lnp(r))；
图2为本发明在2月10日探测角为90°时根据回波信号用klett法计算得到的消光系数分布图(横坐标表示消光系数，纵坐标表示距离r)；
图3为本发明在2月10日多普勒测风激光雷达在探测角为30°时所得到的回波信号图(纵坐标表示距离r，横坐标表示距离回波信号功率对数lnp(r))；
图4为本发明在2月10日探测角为30°时根据回波信号用klett法计算得到的消光系数分布图(横坐标表示消光系数，纵坐标表示距离r)；
图5为本发明在2月10日多普勒测风激光雷达在探测角为4°时所得到的回波信号图(纵坐标表示距离r，横坐标表示距离回波信号功率对数lnp(r))；
图6为本发明在2月10日探测角为4°时根据回波信号用klett法计算得到的消光系数分布图(横坐标表示消光系数，纵坐标表示距离r)；
图7为本发明适应度函数f(x)的收敛变化示意图；
图8为本发明适应度函数f(x)的更快收敛的变化示意图；
图9为本发明在2月16日多普勒测风激光雷达在探测角为90°时所得到的回波信号图(纵坐标表示距离r，横坐标表示距离回波信号功率对数lnp(r))；
图10为本发明在2月16日探测角为90°时根据回波信号用klett法计算得到的消光系数分布图(横坐标表示消光系数，纵坐标表示距离r)；
图11为本发明在2月16日多普勒测风激光雷达在探测角为30°时所得到的回波信号图(纵坐标表示距离r，横坐标表示距离回波信号功率对数lnp(r))；
图12为本发明在2月16日探测角为30°时根据回波信号用klett法计算得到的消光系数分布图(横坐标表示消光系数，纵坐标表示距离r)；
图13为本发明在2月16日多普勒测风激光雷达在探测角为4°时所得到的回波信号图(纵坐标表示距离r，横坐标表示距离回波信号功率对数lnp(r))；
图14为本发明在2月16日探测角为4°时根据回波信号用klett法计算得到的消光系数分布图(横坐标表示消光系数，纵坐标表示距离r)。
具体实施方式
下面通过具体实施例对本发明作进一步详述，以下实施例只是描述性的，不是限定性的，不能以此限定本发明的保护范围。
为了说明本专利反演斜程能见度方法的可行性与可靠性，我们分别使用2015年2月10日和2月16日两天的多普勒测风激光雷达回波信号进行了验证。式中，图1、图3、图5分别是2月10日多普勒测风激光雷达在90°、30°、4°三个不同角度探测所得的回波信号；图9、图11、图13分别是2月16日多普勒测风激光雷达在90°、30°、4°三个不同角度探测所得的回波信号。
具体实验流程如下：
(1)建立模型：
通过klett法和激光雷达方程来建立各点大气消光系数的模型如下：
(2)构造方程组：
首先，设大气消光系数边界值为x1，大气后向散射消光对数为x2；
然后，根据多普勒测风激光雷达对数距离校正信号来寻找大气消光系数近似不变的垂直均匀层，再令边界值点和标定点处的消光系数相等，构造第一个方程如下：
最后，根据激光雷达探测的大气气溶胶光学厚度和太阳光度计反演的大气气溶胶光学厚度近似相等再构造另一方程如下：
(3)求解方程组：
1)、利用粒子群算法对上述方程组进行求解，可同时求解出大气消光系数边界值和大气后向散射消光对数比：
首先，对回波信号进行2000次累加平均和基于经验模式分解的预处理；然后，由垂直探测的回波(即图1)确定出标定点rb＝4.62km，边界值点rm＝4.8km；最后，通过粒子群算法求解出非线性方程组的解，即消光系数边界值σ(rm)和后向散射消光对数比k；
2)、根据适应度函数求解σ(rm)及k的最优值：
图7给出了适应度函数f(x)的收敛变化情况，从图中可以看出f(x)迅速下降，在经12次迭代后可达到精度要求，即f(x)最趋近于零，通过迭代可得出σ(rm)＝0.026和k＝0.92。还能通过k的取值范围、斜率法预算的边界值和预设的迭代初始值，来加快迭代速度，减少迭代次数，效果如图8所示。
(4)求解斜程能见度：
将上一步求得的两个值带入到建立好的各点大气消光系数模型中，可求解出大气消光系数均值，再把大气消光系数均值带到如下方程中可得斜程能见度。当消光系数平均值为0.152时，可得斜程能见度为8.337km。
图2、图4、图6分别给出了根据回波信号用klett法计算得到的消光系数分布。下表则详细给出了2月10日各种数据的反演结果。
图10，图12，图14分别给出了根据回波信号用klett法计算得到的消光系数分布。下表则详细给出了2月16日各种数据的反演结果。
实验结果表明，用本专利方法确定的大气消光系数边界值和大气后向散射消光对数比对斜程能见度进行反演，其计算结果更加准确，且收敛速度较快、迭代次数较少。
具体反演过程：
(1)斜程能见度方程的构造过程：
1)水平能见度方程
根据koschmieder提出的目标物视程理论，能见度测量公式为：
式中，v为能见度(km)；ε为视觉对比阈值，是人眼将目标从背景中分辨出来的最小亮度对比，世界气象组织(wmo)推荐ε的值为0.02，而国际民航组织(icao)推荐ε的值为0.05，本专利取ε＝0.05；σ为大气消光系数(km-1)，是表征大气透明度的一个特征量。
2)水平能见度方程的修正
koschmieder公式是对人眼感觉的光波(白光)而言的，作为近似可用人眼最敏感的绿光(0.55μm)来代替白光，即用绿光的大气消光系数来代替白光的大气消光系数系数。我们以绿光波长作为标准，对于其他波长下的能见度求解做如下修正：
式中，λ为激光波长(μm)；q为修正系数，当能见度较好时(v＞50km)，q＝1.6；当能见度适中时(6km＜v＜50km)，q＝1.3；当能见度较差时(v＜6km)，q＝0.585v1/3。
3)斜程能见度方程
对于民航机场飞机的起降，我们更关注的是其斜程能见度。在一般大气条件下，斜程路径是非均匀的。要获得非均匀路径上的真实能见度，我们必须知道光路上的亮度和消光系数，而空对地气柱亮度的测量是比较困难的。本专利采用近水平的空对地斜程能见度方程，避免了使用亮度计测量天空亮度。
根据文献(激光探测机场斜视跑道视程的实验研究，邱金桓1988)，当天气条件为云、雨、雾或v＜4km时，若飞机降落下滑角α＝2°40′，对比感阈值ε＝0.05，向下的平均亮度与2°高度角向下的天空亮度之比d＝2.6，则机场跑道斜程能见度的半经验公式为：
式中，vs是斜程能见度(km)，σp是消光系数平均值。引入波长校正因子后(3)式变成：
式中，λ为激光波长(μm)；q为修正系数，当能见度较好时(v＞50km)，q＝1.6；当能见度适中时(6km＜v＜50km)，q＝1.3；当能见度较差时(v＜6km)，q＝0.585v1/3。
(2)模型的建立过程：
1)激光雷达方程的建立
当多普勒测风激光雷达发射的激光在大气中传输时，其接收单元探测距离r处的大气后向散射回波信号可由以下激光雷达方程表示：
式中，p(r)为探测距离r处大气后向散射回波信号的功率(w)；c1为激光雷达系统常数(w·km3·sr)，它与激光初始能量、系统光学效率、望远镜有效接收面积和探测器量子效率等因素有关；β(r)为探测距离r处的大气后向散射系数(km-1·sr-1)；σ(r)为探测距离r处的大气消光系数(km-1)。
2)水平能见度数据处理方法：
通常情况下，在水平方向上大气是均匀分布的。我们假定沿水平路径大气消光系数恒定为σh，对(5)式两边取自然对数，并对距离r求导可得到：
因为大气在水平方向是均匀的，所以大气后向散射系数不会随距离改变而改变，即
将水平大气消光系数σh代入(2)式，即可求得大气的水平能见度。
3)斜程能见度数据处理方法：
在低层大气中，气溶胶的浓度相对较大，气溶胶的后向散射远大于分子的瑞利散射，故我们可以忽略大气分子的散射作用。在此种情况下，大气消光系数主要由气溶胶消光系数构成，klett法假设大气后向散射系数β(r)和消光系数σ(r)之间存在乘幂关系，即有：
β(r)＝c0σk(r)(13)
式中，c0是常数；幂次值k的大小取决于激光雷达波长与大气气溶胶性质，取值范围一般为0.67≤k≤1。
对于距离r小于标定点距离rm处的大气消光系数可由以下公式求得：
式中，s(r)＝ln(p(r)r2)为激光雷达对数距离校正信号；rm是标定点距离；σ(rm)是rm处的大气消光系数值，即大气消光系数边界值。
(3)方程组的构造过程：
1)以klett法来确定边界值方程，即第一个方程
根据激光雷达对数距离校正信号s(r)寻找大气消光系数近似不变的垂直均匀层，即s(r)曲线中斜率相对不变的一段区域(rb，rm)，再令rm处的消光系数等于rb处的消光系数(σ(rm)＝σ(rb))，用公式(11)求大气消光系数边界值有：
在(15)式中，我们设σ(rm)＝x1，k＝x2则有：
2)以多普勒测风激光雷达探测的大气气溶胶光学厚度和太阳光度计反演的大气气溶胶光学厚度近似相等来建立第二个方程
气溶胶光学厚度可表示为气溶胶消光系数沿垂直路径的积分，其定义式为：
式中，τ为气溶胶光学厚度；σα为气溶胶消光系数；(r1，r2)为积分路径。
以激光雷达探测大气气溶胶光学厚度时，我们可以用klett法反演的大气消光系数减去大气分子的消光系数σm(r)来得到气溶胶的消光系数σα(r)，式中，大气分子的消光系数可用标准大气模式中的大气密度廓线及瑞利散射求出，即气溶胶的光学厚度还可为：
以ce318太阳光度计观测数据所反演的气溶胶光学厚度记为τ0，设太阳光度计探测的整层大气光学厚度为soda，则在激光雷达有效探测范围内气溶胶的光学厚度还可为：
太阳光谱区的大气消光来自大气分子和气溶胶的散射和吸收。由于大部分的气溶胶聚集在0～5km的低层大气，所以气溶胶的消光也主要源于这一层。在低层大气中，气溶胶的消光远大于分子的消光，而随着高度的增加大气分子的消光将占据主导地位。鉴于此，本专利认为，当气溶胶主要集中在以5km下时，5km以上气溶胶颗粒可忽略不计，因此有：
当激光雷达有效探测距离大于5km时，则可用(19)式计算激光雷达有效探测距离内的气溶胶光学厚度。再根据激光雷达探测的大气气溶胶光学厚度和太阳光度计反演的大气气溶胶光学厚度近似相等，即τ＝τ0，因此有第二个方程：
(4)同步确定大气消光系数边界值和大气后向散射消光对数比，再用反演得到的大气消光系数计算斜程能见度
联立式(16)和(20)可得非线性方程组，再用粒子群算法可同步确定σ(rm)和k值。
粒子群算法是美国社会心理学家jameskennedy和电气工程师russelleberhart于1955年共同提出的。该算法利用了生物群体模型，具有全局搜索和局部搜索特点，收敛速度很快且迭代次数较少，在求解方程组的问题中得到了广泛应用。本专利就采用该算法同时求解了大气消光系数边界值和大气后向散射消光对数比。
粒子群算法原理介绍：
对粒子的速度和位置进行初始化，随机选取大小为xi、速度为vi的n个点。
针对所求方程设一个适应度函数，若方程组为：
则其适应度函数设为：
f(x)＝f1(x)2+f2(x)2(8)
找到适应度函数后，根据粒子群算法的迭代公式进行运算，其迭代公式是：
这里，w表示惯性系数，可以写成：
iter表示当前迭代次数，含义是当迭代次数较小时，w比较大，利于快速寻找到解，而当迭代次数较大时，w比较小，利于收敛。pid表示粒子的历史最优值，就是每个粒子在历次迭代中的最优值，即使f(x)最趋于零的值。pgd表示粒子的全局最优值，就是总体的最优值。c1和c2是学习因子，r1和r2表示0到1之间的随机数。xid表示粒子值，vid表示收敛速度因子。
通过迭代不断更新历史最优值，当f(x)最趋于零时，选取历史最优值为所需的解。
得到σ(rm)和k值后，代入式(14)计算出消光系数均值σp，再将代σp入式(4)即可求解出斜程能见度。
本专利选取中国民航大学在国家自然基金项目支持下研制的直接探测多普勒测风激光雷达数据和气溶胶监测网(aeronet)的ce-318太阳光度计数据进行实验。
直接探测多普勒测风激光雷达装置在天津滨海新区中国民航大学停机坪，其提供的大气反射功率信号可用于斜程能见度反演的研究。直接探测多普勒测风激光雷达的激光波长为1064nm，探测距离为10km，距离分辨率为30m，二维扫描系统在水平方向可以旋转0°～360°，垂直方向可以旋转0°～90°，并可通过软件或手动控制。当探测水平或斜程(参考角度为飞机着陆时的下滑角2°～6°)能见度时，首先将扫描系统调整至需测量的方位，然后向大气发射若干激光脉冲，再将回波信号送入光子探测计数器进行检测，最后将结果送入计算机就可反演出大气消光系数和光学厚度，用于能见度的计算。
气溶胶监测网(aeronet)是全球布站的气溶胶特性地基观测网，它以法国cimel公司生产的全自动ce318太阳光度为观测仪器，该仪器提供了1020nm、870nm、657nm和440nm共四个波长的气溶胶光学厚度(aod)数据，利用这四个波长的数据可以确定未知波长的气溶胶光学厚度。为了将激光雷达反演的aod与太阳光度计实测的aod进行对比，考虑到激光雷达安装在距北京香河76公里(直线距离)的天津，激光雷达波长是1064nm，本专利选用aeronet北京站点波长为1020nm的气溶胶光学厚度数据。
尽管为说明目的公开了本发明的实施例和附图，但是本领域的技术人员可以理解：在不脱离本发明及所附权利要求的精神和范围内，各种替换、变化和修改都是可能的，因此，本发明的范围不局限于实施例和附图所公开的内容。