范有臣,赵洪利,孙华燕,郭惠超,赵延仲
(解放军装备学院,北京101416)
基于距离选通的激光主动成像是近年来激光成像的研究热点,但是其成像目标大部分是静止目标,在成像时要预设距离值才能得到目标图像[1-2],只适用于静止的已知距离的合作目标,在实际应用时受到很大限制,无法完成非合作目标的成像,尤其是运动目标。文中针对非合作的运动目标,在距离选通的基础上,采用红外热像仪跟踪运动目标,同时引导激光成像系统,并加入了测距机的功能,实时采集运动目标的距离信息,而后根据此距离信息修正选通距离,进而实现对运动目标的激光成像。
在对运动目标成像时,由于测距机的误差及目标快速运动时的速度影响,需要增加选通门宽,以实现运动目标始终在选通成像范围之内的目的。由于门宽的增加,后向散射的影响也越来越严重。文中分析了后向散射对图像质量的影响,并提出将提升小波变换引入同态滤波的算法,采用高通滤波对小波变换系数进行处理,降低后向散射的影响,提高图像信噪比。
激光主动成像系统对运动目标成像时,相比于静止目标有两个不同点:(1)运动目标的捕获跟踪。由于运动目标空间位置是不断变化的,激光主动成像系统必须跟随运动目标一起运动才有可能获得图像,但是由于激光主动成像系统的视场角小(一般为几mrad),其自身无法完成对运动目标的跟踪,必须借助其他成像系统完成外部引导才能精确跟踪运动目标。(2)运动目标的距离信息获取。由于运动目标在不断运动中,其距离信息在不断地变化之中,而基于距离选通的激光主动成像系统必须依赖于精确的距离信息,距离信息获取对于能否成像起着至关重要的作用。
从以上两点分析,针对运动目标的捕获跟踪,采用红外热像仪在大视场角内发现跟踪目标,针对运动目标的距离信息获取采用激光测距机实时获取目标距离信息。整个系统结构如图1所示。
图1 系统结构图
整个系统工作流程如下:首先,运动目标在大视场范围出现后,红外热像仪捕捉运动目标,稳定跟踪后触发激光器出光,然后,测距机和接收系统开始工作,同时,测距机将距离信息实时传送至微机,控制ICCD探测器的快门开启时间和选通距离,实现对目标的激光成像。系统实物如图2所示。
装置主要由ICCD成像系统、照明激光系统、接收光学系统、同步控制系统、激光测距系统和计算机系统构成。选通型探测器选用了北方夜视技术有限公司的1XC18/18WHS-G选通型ICCD。激光器为1064/532 nm双波段半导体泵浦的YAG激光器,其单脉冲能量可调,1064nm波长输出时最大为200 mJ,532 nm波长输出时最大为80 mJ,重复频率50 Hz。接收光学系统为透射式天文望远镜,焦距1800 mm、口径250 mm。激光测距系统为1570 nm人眼安全激光器,其测距精度为2 m,光束发散角≥1 mrad,在能见度大于12 km时作用距离≥20 km。
图2 系统实物图
基于距离选通的激光主动成像,其选通成像的优点之一即为抑制后向散射,但是在对运动目标成像时,选通门宽比静止目标大,后向散射影响也更为明显。图3为对运动目标成像时,不同门宽时的成像。目标为一簇氢气球,目标距离500 m,其运动速度约为5 m/s。
图3 不同门宽的后向散射
图3 (a)选通门宽为200 ns,后向散射的光基本被滤除,图3(b)选通门宽为1μs,有部分的后向散射光进入到门宽内。从图中可以看出,后向散射光对图像质量影响严重,需要对其进行处理以便于后续目标分析。
同态滤波算法能增强图像对比度并可以将图像的动态范围进行适当的压缩。用f(x,y)二次函数来表示二维空间区域时,可以引入两个分量来进行表征[3-4]:一个是入射时的光源总量i(x,y),另一个是物体在反射时的反射光总量r(x,y),即:
一般情况下后向散射的部分是变化缓慢的,与图像的照明分量变化一致,对应于低频成分,而目标本身的反射分量则会有强烈的变化,对应于高频成分。如果能够设计出提高高频分量,抑制低频分量的高通滤波器H(u,v)就可以达到降低后向散射影响的结果。H(u,v)在选择时分成以下2个部分:低频部分的H(u,v)<1;高频部分的H(u,v)>1。这样就能够使图像中低频成分得到抑制。同态滤波的整个流程用图4进行表示:
图4 同态滤波相应的流程图
图中的对数运算用ln表示,线性高通滤波器用H(u,v)表示,傅里叶变换由FFT表示,傅里叶的逆变换一般由(FFT)-1进行表示,原始图像的函数表达式为f(x,y),经过处理之后的亮度图像为g(x,y)。
选用理想滤波器时,会出现振铃效应,所以必须使用具有光滑特性的滤波器来避免振铃效应的影响。文中采用改进过的二价巴特沃思滤波器[5],不仅能够消除振铃效应也能够尽量减小图像边缘的模糊程度。二阶巴特沃思滤波器的公式如下:
其中,ρ=[(u-u0)2+(v-v0)2]1/2是频率(u,v)到滤波器中心(u0,v0)的距离。当(u0,v0)=(0,0)时,ρ=[u2+v2]1/2,ρc是截止频率,当ρ/ρc〉〉1时,有rH=r1;当ρ/ρc〈〈1时,有rL=r1-r2。
在原始的同态滤波中,需要通过Fourier变换才能够得到图像的频率特征,但Fourier变换还是有很多不足的地方,例如它不能同时进行时间局部化处理,所以图像的频域和时域不能够同时获得。而小波变换在频域和时域中都显示出了良好的局部性能特征,因此进行同态滤波计算时普遍采用小波分析方法。但是小波分析需要借助于傅里叶变换,在进行卷积运算时多采用浮点数运算,这样会对系统的实时处理能力产生很大的影响。
早在1996年,Sweldens等[6-7]学者就曾经提出了一种小波构造,该方案不需要依赖于傅里叶变换的提升方案(Lifting-Scheme),这就是第二代小波变换(SGWT)。通过这种形式的处理能够避免浮点数的运算,完成整数之间的小波变换,这样就不需要通过Fourier变换,降低了运算量。Dauhechies[8]对此做过相关研究,研究表明Mallat算法能够对小波转换进行有效的提升,这种运算形式能够大大降低传统的小波变换的运算量。
SGWT变换的变换过程一般由以下三个不同的部分构成:
(1)分解。一般而言分解是指把原始波形分成不同的部分,分解后的波形由两个相互垂直的子集构成。这样分解后的每个子集的长度只有原来的子集长度的1/2。其具体步骤是将一个数列分为偶数序列ej-1和奇数序列oj-1,即:
(2)预测。预测的过程是有根据的,它的依据是奇数和偶数序列之间存在的联系,可以通过奇数序列预测出偶数序列或是通过偶数序列预测出奇数序列。预测的具体表达式如下:
式中,oj-1表示实际值,p(ej-1)表示预测值,它们之间的差值用dj-1表示。用预测函数pk来表示p,函数pk可取为ej-1本身:pk(ej-1,k)=ej-1,k=sj,2k或者其他比较复杂的函数。
(3)更新。经过对应的分裂步骤产生的子集和原来的数据并不是完全一致的,所以还需要一个更新过程来保证原始数据的整体性不变。用算子U来表示其更新过程,具体算法如下:
其中,sj-1是sj的低频部分。
小波提升算法的过程是完全可逆的,和正向变换一样,小波提升反向的变换过程为:
图5表示的是一幅图像f(x,y)经过两次提升变换的过程。
图5 提升小波变换示意图
图5 中,f是原始图像,分解后各个不同频率的分量分别用HH,HL,LH,LL来进行表示。
从频谱分析中可以看出,原始信号被划为不同的系数,在频带上的表现就是一个低频系数和多个高频系数,这也是小波分解的最终结果。经过小波的不断分解,假设分解了n级,当j=1,2,…,n时,可以得到LLn和LHj、HLj、HHj,其中反映图像各亮度分布和面貌参数的部分是LLn。一般如果要降低图像边缘的模糊程度,可以使用高通滤波处理,从而衰减一些低频的信息。因为后向散射部分在频率域中多表现为低频,所以可以通过这种方法弱化后向散射的影响。对于低频系数而言,巴特沃思滤波器的公式可以简化为r1-r2,并对低频系数采用线性调节,将HLLn=(r1-r2)(k(x-m)+m)作为增强后的小波系数。其中,x为小波系数,m为小波系数的平均值,k(x-m)+m完成对LLn系数的线性均衡,k为对比度调节因子,其满足1≥k≥0。
算法流程如图6所示。
图6 算法流程图
利用上述方法对拍摄的一组实验图像进行后向散射抑制处理,算法选取参数为r1=4,r2=2,kc=1/8,k=1,对原图像进行3层提升小波分解,选用的是’sym3’小波,图7为实验的结果。
图7 算法处理结果
图7 (a)为原始图;图7(b)为经过传统同态滤波的处理后图像,选择二阶的巴特沃思滤波器,滤波器的截止频率是0.75;图7(c)是经过传统小波同态滤波处理后的图像,其参数选择与提升小波保持一致;图7(d)为提升小波的同态滤波处理后的结果。
三种对烟雾弱化处理的方法在视觉上的差别为图7(c)、图7(d)两图明显好于图7(b)。对于三种方法的处理结果进行客观评价:信息熵反映的是信息量,随着信息量的增加,其信息熵也会越来越大图像越清晰;图像对比度用标准差表示,对比度和标准差成正相关关系,图像的渐变等反映的越清晰。图像层次越多,图像的平均梯度值对图像细节变换与纹理变化反映的越详细,其图像也就层次感越好,愈加清晰。
其中,图像行数用M表示,图像列数用N表示,处理后的图像用g(x,y)表示,像素值为i的像素出现的概率为pi,经过计算处理后得到ga(x,y)2和gb(x,y)2。最后的结果如表1所示。
表1 各图像参数计算结果Tab.1 The result of every image
可以看出处理后的图像的灰度明显地降低了,这说明达到了抑制后向散射的效果,但是处理后图像的标准差、信息熵、平均梯度与边缘强度都比原图小,原因是激光图像不同于自然可见光图像,在激光图像中,照亮部分即被认为是目标,由于后向散射的影响,造成除了目标之外,其他非目标部分也有亮斑,这样,图像的各个客观评价参数将不能反映图像的真实情况,必须寻找一种反映后向散射影响大小的参数。图8为边缘检测的结果,从图中可以看出在原图中后向散射占据了很大一部分,处理后的图像减少了后向散射部分,只保留了目标部分。此时,针对后向散射的特点,采用以下参数进行评价。
后向散射值:
其中,g(x,y)表示整幅边缘检测图像,M,N分别为图像的行数和列数。f(x,y)表示目标所在位置的图像,m,n分别为图像的行数和列数。式中,分子为后向散射在整幅图像中的亮度值总和,分母为目标图像的亮度值总和。这样能够有效检测出后向散射在整幅图像中的比例。经计算后的结果如表2所示。
图8 边缘检测结果(图8(a)、(b)、(c)、(d)分别与图7(a)、(b)、(c)、(d)对应)
表2 各图后向散射值
从表2可以看出,处理后的散射值明显小于原图,达到了抑制后向散射的效果。且文中算法明显优于前两种算法。
文中采用红外热像跟踪运动目标,并利用测距机实时采集运动目标的距离信息,实现了对运动目标的激光主动成像。针对激光运动目标成像时独有的后向散射引起的图像降质问题,将后向散射的噪声看成一种低频信息,提出了基于提升小波的同态滤波处理算法,并针对后向散射提出了一种评价参数。实验结果表明,算法能够有效地去除后向散射的影响,并对原始图像质量有一定改善,为后续的目标分析奠定了基础。
[1] ZHU Xiaopeng,LIU Jiqiao,HE Yan,et al.Range gated imaging lidar at wavelength of 532nm[J].Infrared and Laser Engineering,2012,41(2):358-362.(in Chinese)竹孝鹏,刘继桥,贺岩,等.532nm激光距离选通成像系统[J].红外与激光工程,2012,41(2):358-362.
[2] GUO Huichao,SUN Huayan,DU Lin.Range information calculation method for 3D imaging based on serial images by time-slice technology[J].Infrared and Laser Engineering,2012,41(12):3258-3262.(in Chinese)郭惠超,孙华燕,都琳.利用时间切片序列图像的三维成像距离信息计算方法[J].红外与激光工程,2012,41(12):3258-3262.
[3] SHEN Wenshui,ZHOU Xinzhi.Algorithm for removing thin cloud from remote sensing digital images based on homomorphic filtering[J].High Power Laser and Particles,2010,1:45-48.(in Chinese)沈文水,周新志.基于同态滤波的遥感薄云去除算法[J].强激光与粒子束,2010,1:45-48.
[4] ZHENG Dongmei,SHI Junsheng,SONG Xiaohui,et al.ICT image enhancement by wavelet based homomorphic filtering method[J].Journal of Changchun University of Science and Technology,2007,30(3):44-46.(in Chinese)郑东梅,石俊生,宋晓辉,等.基于小波的同态滤波算法在ICT图像增强中的应用[J].长春理工大学学报,2007,30(3):44-46.
[5]ZHANG Jinquan,YANG Jinhua,WANG Xiaokun.Research on the weakening of smoke in images[J].Journal of Image and Graphics,2010,15(12):1733-1737.(in Chinese)张金泉,杨进华,王晓坤.图像烟雾弱化方法研究[J].中国图象图形学报,2010,15(12):1733-1737.
[6] ZHU Xifang,WU Feng,TAO Chunkan.A new algorithm of cloud removing for optical images based on wavelet threshold theory[J].Acta Photonica Sinica,2009,38(12):3312-3317.(in Chinese)朱锡芳,吴峰,陶纯堪.基于小波阈值理论的光学图像去云处理新算法[J].光子学报,2009,38(12):3312-3317.
[7] Voicu L I,Myler H R,Weeks A R.Practical considerations on color image enhancement using homomorphic filtering[J].Journal of Electronic Imaging,1997,6(1):108-113.
[8] Sweldens W.The lifting scheme:A custom-design construction of biorthogonal wavelets[J].Application Computer Harmonic Analysis,1996,3(2):186-200.