陈浩，华灯鑫，张毅坤，闫庆，李仕春

(1.西安理工大学 机械与精密仪器工程学院，陕西 西安 710048；2.西安理工大学 计算机科学与工程学院，陕西 西安 710048)

目前，区域性环境监测是大气环境监测的重点。激光雷达作为一种大气遥感探测的主动探测工具，其在探测高度、空间分辨率、时间上的连续监测和测量精度等方面具有其他探测手段无法比拟的优点[1]，已经成为对大气、海洋和陆地进行高精度遥感探测的有效工具[2]。然而，目前激光雷达设备多属于科研单位自行研发，数据结构、表示格式、数据量纲等都不统一，导致激光雷达的观测数据在观测高度点、观测时刻和数据数值范围三个方面具有较大差异[3]。观测高度点，观测时刻可能造成比较对象在比较位置(观测高度点和观测时刻)出现数据空缺，而不同的数据数值范围会导致数据没有可比性，给各站点数据之间的共享与分析利用带来困难，极大地阻碍了多区域性大气环境监测工作的开展。现有的数据格式转换方法一般是根据具体的数据处理要求设计出来的，不适用于激光雷达数据多维、多类型、计量单位和数量级都不同的异构特征[4]，难以在确保原有数据的趋势特征不变化情况下的数据统一。为了解决这种困难，本文主要针对大气激光雷达数据结构的观测时刻、垂直观测高度点以及观测数值范围三个方面，通过借鉴变换矩阵的基本思想，结合线性插值法填补数据分析比较相应位置的空缺数据，最后对归一化的数据进行无量纲处理，给出激光雷达数据统一化转换模型，实现数据格式的统一化。

1 大气激光雷达数据相关定义

通常激光雷达主要采用垂直观测的方法对大气进行探测，在某一时刻的观测数据在垂直的不同探测高度点(Height)上是均匀离散数据，观测高度数据是一组以激光雷达距离分辨率为公差的等差数列。由于观测时刻(Time)的不确定，在观测时刻点是一组不均匀的离散的时间序列数据[5]。基于激光雷达数据的时空特性及其二维结构，通常采用如下形式描述激光雷达数据。

设D={D(t1),D(t2),…,D(tn)}为大气激光雷达时间序列数据，其中集合:

D(tj)={D(h1,tj),D(h2,tj),…,D(hm,tj)}

式中D(hi,tj)表示观测高度为hi(i=1,2,…,m)、观测时刻为tj(j=1,2,…,n)处的观测数据。hi表示观测高度数据中第i个高度值，通常i越大，hi越大，同理tj表示观测时刻数据中第j个时刻值，通常j越大，tj越大。

观测时刻、垂直观测高度点以及观测数值范围表示了大气激光雷达数据最明显的结构特征，也是影响大气激光雷达数据相互共享分析的主要因素。观测时刻tj(j=1,2,…,n)主要影响激光雷达二维数据的列向量数据对比，观测高度hi(i=1,2,…,m)影响激光雷达二维数据行向量数据对比，观测数据的数值D(hi,tj)含有不同的量纲和数量级，影响整个观测对象特征数据的相互比较。

为了便于建立大气激光雷达数据格式的统一化模型，这里给出激光雷达数据的变换矩阵相关定义与性质。

定义1设一组大气激光雷达时间序列观测数据的观测时刻为tj(j=1,2,…,n)，观测高度离散序列为hi(i=1,2,…,m)，则大气激光雷达时间序列观测数据D可表示为：

其中:

di,j=D(hi,tj)

(1)

性质1设激光雷达的距离分辨率为r，那么:

hi+1-hi=r，i=1,2,…,m-1

性质2若i

ti

性质3若i

hi

2 大气激光雷达数据无量纲统一化转换模型

通过对大气激光雷达数据二维结构的分析，将大气激光雷达数据表示为矩阵形式，那么主要通过矩阵的计算分别实现大气激光雷达数据格式在观测高度点和观测时刻点上的统一化，借鉴变换矩阵的原理以及应用场景[6]，这里将大气激光雷达数据做如下表示。

设两组大气激光雷达数据分别表示为矩阵Am×n和矩阵Bp×q，m、p为观测高度点数，n、q为观测时刻数。其中间矩阵的观测时刻为观测时刻的集合，即：

TA={tA1,tA2,…,tAn}

TB={tB1,tB2,…,tBq}

则转换后的观测时刻为TA∪TB，观测时刻数为:

t=C(TA∪TB)

式中C(TA∪TB)表示集合TA∪TB的元素个数。其观测点观测高度点分别为：

HA={hA1,hA2,…,hAm}

HB={hB1,hB2,…,hBp}

则转换后的观测高度点为HA∪HB，观测高度点数为:

h=C(HA∪HB)

对于Am×n和Bp×q的统一化数据中的在观测时刻与观测高度点的空缺数据使用内插法来填补数据。考虑到激光雷达信号总体上呈现距离的平方反比衰减趋势[7]，对于激光雷达系统观测距离分辨率高时，观测高度点的空缺数据位置处于两个相邻的观测位置之间，因此，在观测高度点方向上采用效率较高的线性插值法。对于对流层稳定大气，大气处于一个稳定状态，水平方向上的扩散较弱，分层不明显，在一定范围内变化较小[8]，因此，在观测时刻方向上依然采用效率较高的线性插值法。结合内插法，构造线性变换完成激光雷达原始数据矩阵到中间矩阵的变换。

2.1 观测高度点数据转换方法

设激光雷达数据矩阵Am×n在观测高度点方面的变换矩阵为Lh×m，其中:

h=C(HA∪HB)

令:

Lh×m=(li,j)h×m;Yh×n=(yi,j)h×n;Yh×n=Lh×mAm×n

有:

HY=HA∪HB

Lh×m由以下方法得出。

1) 若hi∈HY，hAj∈HA，hi=hAj，i=1,2,…,h，j=1,2,…,m，则：yi,k=aj,k，k=1,2,…,n,有：

yi,k=(li,1a1,k+li,2a2,k+…+li,j-1aj-1,k+

li,jaj,k+…+li,mam,k)/m

(2)

则:li,j=m,li,1=li,2=…=li,j-1=li,j+1=…=li,m=0。

2) 若hi∈HY，hi∉HA，hAj

yi,k=(li,1a1,k+li,2a2,k+…+li,j-1aj-1,k+

li,jaj,k+…+li,mam,k)/m

(3)

根据公式(1)有aj,k=A(hj,tk)。式中A(hj,tk)为矩阵Am×n所表示的大气激光雷达二维时间序列数据在第j(j=1,2,…,m)个观测高度点，观测时刻为tk(k=1,2,…,n)的观测值，亦有:

aj+1,k=A(hj+1,tk)

由线性内插法可得:

(4)

那么:

li,1=li,2=…=li,j-1=li,j+2=…=li,m=0

根据公式(3)和公式(4)有:

(5)

由性质1可得:

li,j=m(hA(j+1)-hi)/rA

li,j+1=m(hi-hAj)/rA

式中rA为矩阵Am×n的距离分辨率。

3) 若j=m-1，hi>hA(j+1)，则：

li,m-1=m(hAm-hi)/rA

li,m=m(hi-hA(m-1))/rA

li,1=li,2=…=li,m-2=0

4) 当j=1，hi

li,1=m(hA2-hi)/rA

li,2=m(hi-hA1)/rA

li,3=li,4=…=li,m=0

2.2 观测时刻数据转换方法

设观测时刻数据变换矩阵为Rn×t，其中：

t=C(TA∪TB)

令：

Rn×t=(ri,j)n×t;Zm×t=(zi,j)m×t;Zm×t=Am×nRn×t

有：

TZ=TA∪TB，Rn×t的计算同2.1节。

1) 若ti∈TZ，ti=tAj∈TA，i=1,2,…,t，j=1,2,…,n，则：

zk,i=ak,j，k=1,2,…,m

rj,i=n

r1,i=r2,i=…=rj-1,i=

rj+1,i=…=rn,i=0

2) 若ti∈TZ，ti∉TA，i=1,2,…,t，且tAj

rj,i=n(tA(j+1)-ti)/(tA(j+1)-tAj)

rj+1,i=n(ti-tAj)/(tA(j+1)-tAj)

3) 若j=n-1，ti>tA(j+1)，则：

rn-1,i=n(tAn-ti)/(tAn-tA(n-1))

rn,i=n(ti-tA(n-1))/(tAn-tA(n-1))

r1,i=r2,i=…=rn-2,i=0

4) 若j=1，ti

r1,i=n(tA2-ti)/(tA2-tA1)

r2,i=n(ti-tA1)/(tA2-tA1)

r2,i=r3,i=…=rn,i=0

这样，构造出变换矩阵Rn×t，从而完成对观测时刻数据的转换。

2.3 无量纲统一化转换模型

大气激光雷达数据必须进行无量纲处理，消除数据范围量纲的不统一带来的数据分析困难。这里采用线性极值法对不同格式数据进行数据归一化处理，线性无量纲处理方法能够很好地保持原始数据的数据趋势，对数据的影响最小[9-10]。具体方法如下。

设X为一组大气激光雷达时间序列数据经过观测高度点以及观测时刻转换的数据，X的矩阵表示为：Xh×t，Xh×t=Lh×mAm×nRn×t，x为Xh×t中的某一个值，则x归一化处理后的数据x′为：

x′=(x-xmin)/(xmax-xmin)

(6)

式中，xmax与xmin分别为X观测数据的最大值与最小值。x′的范围在0～1之间，且各x′值的分布仍与相应原X各值的分布相同[11-12]。

由以上关于观测点高度数据转换、观测时刻数据转换以及数据的无纲量化处理方法可得，对于大气激光雷达数据Am×n和Bp×q，设Uh×t，令：

ui,j=xmin，i=1,2,…,h，j=1,2,…,t

λ=1/(xmax-xmin)

λ为大气激光雷达数据Am×n的归一化因子，则与Am×n格式相同的数据的数据联合统一化模型为：

(7)

3 应用实例分析

为了验证上述激光雷达数据无量纲统一化转换模型公式(7)的可行性，这里选取两组不同格式类型的激光雷达数据A和B进行数据转换，表1和表2分别是两组大气激光雷达数据的部分数据。数据A为2011年4月10日，利用小型米散射激光雷达，在位于西安理工大学校园内对西安地区上空的气溶胶光学特性进行反演得到的气溶胶消光系数。数据B为2013年6月28日，在渭南环境检测局利用微脉冲米散射激光雷达观测气溶胶得到的消光系数。为了便于验证计算，这里只取部分数据进行数据格式转换。

表1 大气激光雷达数据A部分数据

表2 大气激光雷达数据B部分数据

数据A的观测时刻集为:

TA={6,9,11,12,14,15,18}

其数据矩阵列数为7，数据B的观测时刻集为:

TB={5.5,7.5,9.5,11.5,12,13,14,15,16,17}

其数据矩阵列数为10，则转换后的观测时刻为TA∪TB，观测时刻数为:

t=C(TA∪TB)=14

数据A的观测高度点:

HA={0.0015,0.0045,0.0075,…,0.0315}

其数据矩阵行数为11，数据B的观测高度点集为:

HB={0,0.0075,0.0150,…,0.0300}

其数据矩阵行数为5，转换后的观测高度为HA∪HB，观测高度点数为:

h=C(HA∪HB)=14

对于A11×7构造L14×11使得:

Y14×7=L14×11A11×7

HY=HA∪HB

TY=TA

rA=0.003

令:

HY={h1,h2,…,h14}

HA={hj}，j=2,3,4,5,6,8,9,10,11,12,14

hj=hAi

则:

yj,k=ai,k，i=1,2,…,11，k=1,2,…,7

根据公式(2)可得:lj,i=11，lj,1=lj,2=…=lj,i-1=lj,i+1=…=lj,11=0。

若:HY-HA={hj}，j=1,7,13，hAi

当j=1，hj

当j≠1，lj,i=11(hA(i+1)-hj)/0.003，lj,i+1=11(hj-hAi)/0.003，lj,1=lj,2=…=lj,i-1=lj,i+1=…=lj,11= 0。

对于A11×7,根据2.1节与2.2节可得如下观测高度点数据转换矩阵L14×11和R7×14。

由公式(7)可得A11×7类型的无量纲统一化数据转换模型为：

(8)

表3 数据A统一化转换数据A′14×14

表4 数据B统一化转换数据B′14×14

由表3、4可以看出，统一化转换后的数据A和B具有相同的结构，观测时刻与观测高度点相同，数据量纲已经消除，可以进行相互比较分析等。

由公式(7)可知，矩阵转换的复杂度为O(n3)，随着激光雷达数据的观测时刻点和观测高度点的增加，数据转化所消耗的时间呈三次函数增加。

4 结 语

本文针对大气激光雷达数据的异构性，借鉴了变换矩阵、线性内插法与极值法的数据归一化处理方法，从激光雷达数据之间差异较大的几个因素方面，如激光雷达数据的观测高度点、观测时刻点、测量数据值范围等，建立了不同类型的激光雷达数据格式之间的相互转换模型，对激光雷达数据进行统一化处理和无量纲处理。实验表明，对于在数据的观测高度点、观测时刻点、测量数据值范围等方面产生异构性的激光雷达数据，使用本文给出的方法完全能够实现数据格式的统一化与无量纲化，使数据具有公度性。

然而，本文中的数据插补方法仅适用于激发雷达系统距离分辨率较小的情况，如果距离分辨率较大，需要考虑到激光雷达信号指数衰减、距离平方衰减的曲线曲率，需要研究更加准确的插补方法，本文将对此继续深入研究。

参考文献：

[1]毛建东,华灯鑫,何廷尧,等.银川上空大气气溶胶光学特性激光雷达探测研究[J].光谱学与光谱分析, 2010, 30(7):2006-2010.

Mao Jiandong, Hua Dengxin, He Tingyao, et al.Lidar observations of atmospheric aerosol optical properties over Yinchuan area[J].Spectroscopy and Spectral Analysis, 2010, 30(7):2006-2010.

[2]Deleva A D, Grigorov I V, Avramov L, et al.Raman-elastic-backscatter lidar for observations of tropospheric aerosol[C]∥Proceedings of SPIE, Bulgaria,2008.

[3]尹青,何金海,张华.激光雷达在气象和大气环境监测中的应用[J].气象与环境学报,2009,25(5):48-56.

Yin Qing, He Jinhai, Zhang Hua.Application of laser radar in monitoring meteorological and atmospheric environment[J].Journal of Meteorology and Environment,2009,25(5):48-56.

[4]梁欣廉,张继贤,李海涛,等.激光雷达数据特点[J].遥感信息,2005,(3):71-76.

Liang Xinlian, Zhang Jixian,Li Haitao, et al.The characteristics of lidar data[J].Remote Sensing Information,2005,(3):71-76.

[5]陈浩,华灯鑫,张毅坤,等.基于三次样条函数的激光雷达数据可视化插值法[J].仪器仪表学报, 2013,34(4):831-837.

Chen Hao, Hua Dengxin, Zhang Yikun, et al.Interpolation method for lidar data visualization based on cubic spline function[J].Chinese Journal of Scientific Instrument, 2013, 34(4):831-837.

[6]杨卫东,刘玉树.一种不同坐标系之间的变换矩阵的转换方法[J].计算机辅助设计与图形学报,2000, 12(1),53-56.

Yang Weidong, Liu Yushu.A method for converting transformation matrices between different 3D coordinate systems[J].Journal of Computer Aided Design and Computer Graphics, 2000, 12(1), 53-56.

[7]刘增东,刘建国,陆亦怀,等.基于EMD的激光雷达信号去噪方法[J].光电工程,2008,35(6):79-83.

Liu Zengdong, Liu Jianguo, Lu Yihuai, et al.De-noising lidar signal based on EMD method[J].Opto-Electronic Engineering,2008,35(6):79-83.

[8]陈丽华,林万涛,林一骅,等.一类大气非线性动力热力学反应-扩散系统解的稳定性态[J].物理学报, 2012, 61(14):140202-1-5.

Chen Lihua, Lin Wantao, Lin Yihua, et al.Stable behavior of solution for the reaction-diffusion system of atmospheric nonlinear dynamics and thermodynamics[J].Journal of Physics, 2012, 61(14):140202-1-5.

[9]易平涛,张丹宁,郭亚军,等.动态综合评价中的无量纲化方法[J].东北大学学报:自然科学版, 2009, 30(6):889-892.

Yi Pingtao, Zhang Danning,Guo Yajun, et al.Study on dimensionless methods in dynamic comprehensive evaluation[J].Journal of Northeastern University (Natural Science), 2009, 30(6):889-892.

[10]周玲微,雷廷武,武阳.岔巴沟流域次暴雨产流无量纲模型[J].农业工程学报,2010,26(11):54-60.

Zhou Lingwei,Lei Tingwu, Wu Yang.Event-based dimensionless models for runoff of Chabagou watersheds[J].Transactions of the Chinese Society of Agricultural Engineering, 2010, 26(11):54-60.

[11]郭亚军,易平涛.线性无量纲化方法的性质分析[J].统计研究,2008,25(2):93-100.

GuoYajun, YiPingtao.Character analysis of linear dimensionless methods[J].Statistical Research, 2008, 25 (2): 93-100.

[12]江文奇.无量纲化方法对属性权重影响的敏感性和方案保序性[J].系统工程与电子技术,2012, 34(12):2520-2523.

Jiang Wenqi.Sensibility and alternative COP analysis of dimensionless methods on effect of attribute weight[J].Systems Engineering and Electronics, 2012, 34(12):2520 -2523.