高炉燃烧料面合成孔径雷达双聚焦成像算法

陈相兵，倪梓明，陈先中，侯庆文，于 晴

(1.北京科技大学自动化学院，北京 100083;2.北京科技大学工业过程知识自动化教育部重点实验室，北京 100083)

1 引言

高炉内部料面的可视化[1]，对于高炉布料具有指导性意义，有利于布料策略的调整，能够保证炉内物料和煤气流的合理分布，对于提高高炉产能有很大帮助.结合光学成像原理的高温工业内窥镜装置[2]能够拍摄炉内料面图像，采用自适应最佳分数阶多向微分算子检测料形轮廓曲线[3]，可以有效获取料形轮廓.激光检测法[4]通过在高炉内形成激光网格，再用摄像机捕捉斑点图像获取料面信息，但在高炉存在大量粉尘时收发能量会存在衰减.近年来调频连续波(frequency modulated continuous wave，FMCW)雷 达在工业生产和民用自动化领域获得了广泛的应用，也被应用到高炉料面的高度测量和料面成像方面.在高炉炉喉处对称安装两台雷达能够完整扫描出料面的整条径向料形，因此将双摆动雷达应用于高炉料面形状的测量是料面测量发展的趋势.

由于高炉是一个密闭容器，测量环境复杂，存在大量粉尘、颗粒物等随机干扰以及溜槽等固定干扰，会对雷达回波产生较大影响.为了获取有效的测量点距离值，近年来许多专家学者进行了这方面的研究.传统的测距原理是对雷达混频信号做快速傅里叶变换然后寻找频谱峰值，根据峰值谱线号计算出距离值[5].陈先中等人[6]为了弥补距离分辨率不足的缺点，采用改进的FMCW测量原理，提出一种智能时变阈值信号处理[7]方法，解决了因高炉内料面反射系数变化影响雷达测距精确性的问题.侯庆文等人[8]针对低信噪比雷达信号距离宽带模糊的问题，提出了以先验信息为辅助的关键区域线性调频z变换算法和频带能量加权的距离反演方法，提高算法精度，为后续计算提供依据.

高炉中心区域往往情况复杂，煤气流旺盛，熔融状态喷涌态的情况使得雷达回波包含信息复杂，设定两台雷达均能够扫描炉心区域，并采取融合算法[9-10]处理扫描重叠区域的采样点得到可信度高的料面位置.针对数据融合方面的研究，陈泽宗等人[11]提出一种基于多频雷达的数据融合与跟踪算法，通过加权最近邻关联融合多频数据，对于较高频率的数据选择较高的权重，以此来达到融合数据的目的.Zhang等人[12]为提高多幅存在相同场景的不同图像融合时的性能，提出速降单纯形迭代细化算法，联合考虑配准原图位置和优化不同传感器的融合参数，以获取最优的融合效果.

高炉内部的恶劣测量环境会对FMCW雷达信号造成影响，使其接收到的混频信号信噪比不稳定，测量的料面距离准确性不稳定.特别是炉心处强气流，喷涌态的条件下，谱线峰值处并不一定是料面回波，寻峰求距的效果进一步下降.本文对一维单点信号的频谱沿高炉径向扩展到一个扫描周期内连续多点信号的二维频谱图，并与图像处理的方法结合，把能量值转换成0～255范围内的灰度值，进而从图像的角度处理频谱图.首先通过重采样将雷达扫描的坐标非均匀扇形图插值成常规的矩形图，根据每台雷达的扫描范围确定出炉心处双雷达聚焦区域，按照信噪比关联可信度的原则融合该区域，同时对比直接均值融合的方式，本方法对于处理雷达频谱图像更具科学性，最后采取恒虚警率目标检测算法[13-14]突显料面位置.

2 SAR雷达双聚焦料面成像建模

2.1 成像处理模型描述

SAR成像主要是通过分析散射点对FMCW的散射回波，得到对应各步进频点的散射数据，即散射时域分布信息.目标对测量信号响应时刻与其所处的距离对应，以时间为自变量的散射目标响应结果则可以表征目标的一维距离像，即目标沿着距离向的散射强弱分布情况.图1为成像模型示意图，目标由沿z轴排列的一系列散射点组成，收发探头在(0，z0)处.

图1 成像模型示意图Fig.1 Imaging model sketch

则散射回波信号S(k)可表示为

σ(z)是散射强度，k是波数大小，f表示信号频率，c是光速.整理可得

等式说明以原点为中心的目标所在位置的散射波谱.根据合成孔径雷达(synthetic aperture radar，SAR)成像原理，高炉两边雷达相当于“静态SAR”，通过使雷达摆动扫描，将每个位置的“像”合成在一起，得到料面的SAR成像结果.本文所做的算法处理就是在该散射波谱的基础上进行的.高炉内部双雷达扫描料面示意图如图2所示.

图2 双雷达扫描示意图Fig.2 Scanning diagram of dual radars

炉心部位(图2阴影部分)属于关键区域，其料面形状主要有凸起、平面、漏斗形，冶炼状态也包括喷涌态、鼓泡、平静态等等，都会影响料面测量精度;另外炉内的复杂干扰条件及料面形状可能会对目标物产生遮蔽效应.炉心区域气流旺盛，料面状态较复杂且雷达回波信号突变强，加之离雷达较远，因此对炉心区域采取冗余测量再融合的方式可提高关键区域测量可靠性，即根据需要让两台雷达扫描到对方区域一定范围.重叠范围的左边界由右侧雷达扫描至高炉左侧的最远位置确定，右边界由左侧雷达扫描至高炉右侧的最远位置确定.雷达发收一次得到一个频谱，摆动一周能够得到整个扇形区域内的频谱信息;转换到图像中，在炉心处根据料面的散射波频谱，可以分析出雷达的信噪比高低从而给予融合像素点不同的置信度，得到径向完整料面图.高炉内部高粉尘、多颗粒物的状况使得频谱会出现多个波峰，通过恒虚警率检测算法识别料面目标并弱化干扰.本文提出一种适用于雷达扫描频谱图像处理的重采样方法与双聚焦区域融合方法，并将恒虚警率检测应用于料面识别方面.整体处理流程如图3所示:

图3 料面双聚焦成像处理流程Fig.3 Processing flow of dual-focus imaging on burden surface

2.2 SAR图像像素点重采样插值

雷达每次采样得到的时域数据做快速傅里叶变换(fast Fourier transform，FFT)变换得到一维测量频谱，扩展成由一个扫描周期内的频谱数据组成的二维矩阵:

考虑雷达安装位置及扇形摆角，转换成直角坐标系中的雷达扫描各处坐标点，就能从雷达SAR图像中看到料面变化趋势.把扇形能量值图像转换成0～255范围内的灰度图像，扇形图中各像素点的坐标在直角坐标系中是非均匀的，而常规矩形图的各像素点位置都是整数，为了将扇形频谱图像转换成二维矩形图像，采用重采样插值[15-17]算法得到灰度图像.插值示意图如图4所示.

图4 重采样插值方法示意图Fig.4 Schematic diagram of resampling interpolation method

图中:黑色点代表频谱图像中的某一频率处相应角度处的灰度，位置由直角坐标描述;白色点代表重采样插值点，位置由直角坐标描述，该点的值也用灰度值表示.扇形图坐标的实际意义:横坐标代表高炉的径向距离，纵坐标代表料面高度值，以米(m)为单位.由于待插值点周围的扇形图的像素点个数不固定，因此一般的双线性内插法和3次内插法不适用，而最近邻插值法的效果有待考究.本文提出了一种适合于非均匀坐标像素点转换成均匀坐标像素点的重采样插值方法，即通过以待插值点周围每一点为中心的高斯曲面来确定对待插值点灰度的影响.

设φ是输入图像像素点的离散集合，κ为输出图像像素点集合，选取一个从φ到κ的映射

使得输入图像每一位置的灰度值与输出图像相应位置的灰度值一一对应，即可计算出插值点像素大小.对于二维随机向量(X，Y)具有概率密度:

其中:μ1，μ2，σ1，σ2，ρ 均为常数;; ρ为(X，Y)的相关系数.

当ρ=0且σ1=σ2=σ时，二维高斯曲面的公式为

x，y代表像素的坐标，中心位置为原点.高斯曲面是对称曲面，因此选取高斯函数映射关系判断某一点灰度值对待插值点灰度值的影响度.

图5为插值示意图，对于矩形图像中某一待插值点(x，y):

图5 插值示意图Fig.5 Interpolation schematic

分别以黑点位置(xi，yi)为中心建立如图6标准正态分布的高斯曲面，根据(x，y)点距中心点的距离按照高斯面的变化估算待插值点处的灰度大小.

图6 三维高斯曲面侧视图Fig.6 Side view of three-dimensional Gaussian surface

对于一般的n维实数向量空间ℝn，向量的范数可以用来描述长度或者相应的两点之间的距离:设两向量xi，xj∈ℝn，

一般的p范数的距离Lp为

当p=2时，就是欧氏距离

从而可以得到该点处的灰度值大小:

待插值点离高斯曲面中心的距离代表了受影响程度，待插值点在高斯曲面上的值即为相对于(xi，yi)处灰度值的待插值点的灰度大小.依次计算出待插值点相对于周围每一像素点的高斯曲面插值结果，最终得出周围点对待插值点的映射，是周围各像素点对其贡献大小的综合.

2.3 料面双聚焦区域频谱图像融合算法

高炉上对称安装两台雷达，对双雷达扫描频谱图像双聚焦部分的合理融合，能够得到可信度较高的完整径向料面分布情况.明确雷达在高炉上的安装位置以及喇叭天线摆角的情况下能够计算出两幅图像的重叠位置.采用合适的融合算法是得到较好处理效果的关键.本文采用一种改进的基于信噪比大小的渐入渐出图像融合处理方法，用于融合雷达扫描频谱图像，图像融合示意图见图7.

图7 图像融合示意图Fig.7 Schematic diagram of image fusion

渐入渐出的图像融合算法是直接针对重叠区域内像素点灰度值进行处理的，实质就是对两幅图像对应像素线性分配权值.

确定重叠区域左边界x1和右边界x2的计算方式如下:

式中:N是两侧雷达第1个调频周期获取的有效频带范围的最大谱线号;0.0915为FMCW测距原理的固定系数;α1，α2分别是两个喇叭摆动的起始倾角;R，R′如图2 中标示，代表雷达至炉心距离.用f1(x，y)和f2(x，y)表示同一时间段内两台雷达扫描的频谱图像中的所有像素点.待融合区域的灰度值计算方式为

融合后整幅图像每个像素点的灰度值为

f(x，y)表示融合后图像像素点的灰度值;f1(x，y)表示左边位置待融合图像像素点的灰度值;f2(x，y)表示右边位置待融合图像像素点的灰度值;w1，w2是对应左图和右图的权值，并且满足w1+w2=1，0 ≤w1≤1，0 ≤w2≤1.

w1和w2的计算公式为

式中:xi为当前像素点的横坐标;x1为图像重叠部分左边界像素点横坐标;x2为重叠区域右边界像素点横坐标.

设numbers为雷达发收一次混频时域信号的采样点个数，adData[i]为某采样点的能量值，可得能量均值adAvr:

则该次雷达信号的平均偏差adu表示为

雷达初始化增益Gain=125，每次测量后该次的增益自适应调节:

通过增益值调节雷达硬件运算放大电路使信号采样点平均偏差稳定在900～1000并且放大信号能量.

再考虑频谱信噪比大小，信噪比高的信号更能反映料面的真实位置，测量结果的可信度也就更高.因此在渐入渐出算法融合像素点时参考双雷达对于炉心区域相同位置混频信号的信噪比.对于噪声峰值比目标峰值高的低信噪比信号，这种方法能有效抑制噪声对其融合位置的影响.信号的信噪比大小用频谱估算法来评估，需先确定信号频带范围，再估算出信号和噪声的总能量，fL为频带的起始频率，fH为频带的终止频率，信号和噪声频带示意图见图8.

图8 信号和噪声频带示意图Fig.8 Schematic diagram of signal and noise band

信噪比可用下式直接计算:

计算出重叠区域中2个原始像素点所属两个雷达的频谱中分别的信噪比SNR1与SNR2，则f(x1，y1)处的按信噪比分配权值为

f2(x2，y2)处的权值为

最终的融合权值修改为

2.4 SAR图像恒虚警率料面检测降噪算法

雷达扫描得到的SAR图像中包含了各种杂波等干扰信号，为了准确定位到图像中的料面位置，需要区分开干扰信息和料面信息.由于料面目标的散射特性和其他噪声区别较大，基于贝叶斯统计检测理论的恒虚警率目标检测算法(constant false-alarm radar，CFAR)是目前最为实用的方法，有着较强的鲁棒性，在干扰噪声大的SAR图像目标检测领域应用广泛.

假设SAR图像中不同区域内存在独立同分布的高斯噪声，在雷达回波中只有噪声的回波时设为假设H0，在雷达回波中同时有目标和噪声的回波时设为假设H1.设有关料面的回波样本单元为

图9为检测模型.噪声干扰功率的估计源于样本单元中数据的算术平均值，经过线性检波器后各样本值服从瑞利分布，设噪声功率为ω2，当样本集合中没有料面信息时，其概率密度函数为

图9 CFAR检测模型Fig.9 Detection model of CFAR

当样本数据中包含料面信息时，根据回波的信噪比SNR，此时概率密度函数为

其中I0为0阶第一类贝塞尔函数，对于N个样本单元，各自的联合概率密度为

这时，通过设置一个门限值T，若信号回波功率大于T便可认定信号中存在料面信息，假设H1成立.那么检测概率Pd与虚警概率Pfa可以由下式计算可得

根据纽曼-皮尔逊准则，不断调整门限值T，确保虚警概率Pfa保持稳定，使其小于一定的容忍值ϵ的情况下，使检测概率Pd最大化.因此，引入拉格朗日算子λ，使方程

取得最优解.

将式(31)－(32)代入上式可得

不考虑与T无关的λϵ项，为让fλ取得最大值，对于x ≥T，需满足

可得门限判别准则:

不等式左边用式(27)－(28)代入:

根据贝塞尔函数的级数展开式在x较小时近似简化为I0≈1+x/4，并对式(37)取自然对数，便可得

从而化简式(36)的判别准则，得到

其中K为比例系数，与虚警概率有关，即

若使用的是平方检测器，样本单元x服从指数分布，概率密度函数为

则根据极大似然估计法便可得到噪声功率ω2的估计值为

因此，门限值T的概率密度函数为

虚警概率与门限T相关，是一个不确定的量，其数学期望如下:

最后得出门限值T的表达式为

本文中，使用一个矩形滑动窗口，如图10，在SAR图像上进行滑动，来去除杂波.矩形框包含着用来检测的目标点、用来去除背景噪声的背景杂波区和避免所要检测的像素点被误当成噪声的保护区.保护区P的长度p1、宽度p2与背景杂波区B的长度b1、宽度b2可根据实际情况(图像的大小，所检测目标的形状)设置，以达到更好的检测效果.

图10 滑动窗口模型Fig.10 Model of sliding window

计算背景杂波区B中像素点的个数N，同时将保护区P的像素点xi置为0.

根据式(46)计算出检测门限T，最后得出检测点x是否是目标料面点的CFAR判决:

3 实测数据处理结果及验证

由于高炉危险复杂的内部环境，目前在实际运作中无法通过高炉内的真实准确料形来验证雷达测量的精度.因此为证实摆动雷达扫描得到料面的真实性，在南京钢铁厂区搭建实验环境，通过人工搭建不同料形，对比由雷达扫描冷态料面得到的测量结果来证实雷达测量结果的可靠性.某料形的尺寸设计如图11所示(雷达至料面左边缘4.4 m).

图11 料面模型设计图Fig.11 Design drawing of the burden surface model

该料形的实际搭建效果如图12所示.

图12 人工搭建料面形状Fig.12 Shape of artificially constructed surface

经上方雷达扫描的图像如图13所示.

图13 已知料面扫描结果图像Fig.13 Scanning result image of burden surface

对比雷达扫描图像和设计的料形图像可以得出，雷达扫描图像所反映的料面变化与设计图中的各部分尺寸保持一致，料面左边缘高度为4.4 m，由频谱能量带的上下边界可标识冷态下料面位置误差在5%以下.结合不同料形时的所有测量结果，证实摆动雷达对料面测量结果的可信性，能够对正常生产时的高炉内径向料面提供可靠测量.

本文的数据样本取自南京钢铁厂3高炉上双摆动雷达料面监测系统扫描高温、高压、燃烧状态料面的实测数据.取两台雷达同一布料周期扫描得到的时域数据进行分析验证，检验算法的稳定性.单台雷达每个周期大概能扫描50组数据，经过FFT后由一个周期扫描得到的频谱扩展成的二维灰度SAR成像如图14所示(横坐标表示高炉径向尺度).

图14 原始灰度图像Fig.14 Original gray image

参考高炉半径为3.9 m，按照最近邻插值法得到的二维矩形图像效果如图15所示.

图15 最近邻插值效果Fig.15 Nearest neighbor interpolation effect

本文所用重采样插值方法效果如图16所示.

图16 本文重采样算法插值效果Fig.16 Interpolation effect of the resampling algorithm in this paper

插值图像与原始图像各像素点的均方根误差(root mean square error，RMSE)能够定量对比两种方式的插值效果.

由图15－16可知，最近邻插值法获得的图像像素块明显，图像粗糙;结合表1，该方法处理的图像与原始图像的RMSE更大，表明更多的料面信息丢失.本文所用重采样插值法能在尽可能保存原始信息的基础上平滑料面宽带信息(RMSE较小)，使得料面宽带分布较匀称，料面目标更突出.

表1 插值图像与原图像的均方根误差Table 1 Root mean square error of interpolation image and original image

为了验证本文的双聚焦图像融合算法(dual-focus SAR，DfSAR)对实测雷达信号频谱图像的有效性，采用同一布料周期内两台雷达扫描的径向料形数据，并与双雷达扫描重叠区域采用直接均值融合(mean fusion，MF)方法比较效果，且用客观评价指标评价融合质量.图17为经过插值处理后双雷达扫描的待融合左右两幅料面频谱图.

图17 待融合左右位置图像Fig.17 Left and right position images to be fused

均值融合的处理效果如图18所示.

图18 均值融合算法效果Fig.18 The effect of mean fusion algorithms

本文双聚焦融合算法的处理效果如图19所示.

图19 本文融合算法效果Fig.19 The effect of fusion algorithm in this paper

对比图18－19的融合效果，图18融合区域存在明显叠影现象，弱化了关键的料面区域.图19中的重叠区域过渡自然且反映料面连续变化特征，料面成像效果较优.

对于图像融合效果的客观评判，表2分别给出了基于信息熵、互信息、标准差、清晰度的图像质量评价结果，按以上指标评判雷达SAR图像的融合处理结果.

表2 图像融合算法性能评价指标Table 2 Performance evaluation index of image fusion algorithm

表2显示，双聚焦融合处理算法与均值融合处理算法相比，所得到图像拥有更高的信息熵和互信息，表明其所含信息更丰富，包含原图像中的信息量更大，融合效果更好.标准差评价图像反差大小，值越大，可利用的信息越多，融合效果越好.图像清晰度用平均梯度指标度量，反映了融合图像对微小细节反差和纹理变化的表达能力.DfSAR处理结果的清晰度值较大，融合效果相比MF方法的结果更清晰.因此可以说明，本文所采用双聚焦图像融合算法具有明显的优越性.另使用不同时刻的4组燃烧料面实测数据验证所提融合算法对料面成像的有效性，两种算法处理效果对比和融合图像质量评价如图20和表3所示.

图20 4组扫描数据的不同融合算法效果Fig.20 The effect of different fusion algorithms for four groups of scanned data

表3 两种算法的融合图像质量评价Table 3 Fusion image quality evaluation of two algorithms

融合后的图像中依然存在大量的干扰性杂波信号，使得料面显示不甚清晰.因此，对双雷达扫描频谱图融合处理后的料面SAR图像，采用恒虚警率目标检测算法进行处理，在存在随机噪声的SAR图像中，使得正确识别料面信息的概率达到最佳，降噪后的料面图像处理结果如图21所示.

图21 料面SAR图像降噪效果Fig.21 Noise reduction effect of SAR image of burden surface

对比图19和图21，恒虚警率目标检测算法消除了料面的背景噪声和杂波干扰，识别出料面像素点，所得的雷达SAR料面图像中目标更为明显，料面形状更易分辨，显而易见高炉内料面位置在4 m～6 m之间.

4 结论

为实现双摆动雷达对高炉整条径向料形成像的目的，本文结合SAR成像原理，从图像处理的角度对实际料面回波频谱图像进行处理，先将非均匀扇形能量值图像用重采样插值算法转换成常规灰度图像;其次创新性的对炉心复杂区域采取双聚焦测量方式，设计算法对两台雷达扫描重叠区域的像素点进行双聚焦融合，得到中心区域料面模型;处理结果经过性能评价指标确认，互信息改善10.43%，清晰度改善23.91%.最后通过恒虚警率降噪算法识别整条料面位置.刻画出料面中心区域料形特征，并获取高炉全料面清晰的径向扫描图像.