基于轮廓重建的渔船雷达MAP-MRF 快速反演研究

郑建丽，曹建军，刘 洋

（1.中国水产科学研究院渔业机械仪器研究所，上海 200092；2.辽宁师范大学计算机与信息技术学院，辽宁 大连 116081）

0 引 言

船载雷达是放置在船上用来观测船舶周围各种水面目标和即时气象环境的重要观测工具。在船舶航行中遇到雾、霾、雨等现象导致能见距离受到影响，无法凭目视发现目标的条件下，雷达表现出受环境影响小、观测距离远的优势。目前，常用的航海雷达有SPERRY-MK 系 列、JRC-JMA 系 列 和FURUNO-FAR系列等。根据国际电工委员会IEC62388—2013 标准，船载雷达的显示图像应能够进行连续且平滑的雷达图像更新[1]。但目前主流的雷达终端只能满足部分观测需求，还无法达到连续平滑的显示效果，难以观测海上近地大气层一些快速生成和移动的强对流气象目标。从观测对象的角度分析，海洋环境水汽热量充足，下垫面平滑，气象目标形成和移动速度较快。约70%的龙卷风从出现到消失的生命周期为1～10 min，海上冷风的移动速度可达70 km/h 以上[2]。这些快速生成和变化的气象目标对雷达数据处理提出了高清显示和更快速度的要求。从观测平台的角度分析，船载雷达因平台稳定性差和能量系统相对封闭等原因，探测距离一般为40～80 km。相较于陆基雷达通常200～300 km 的探测距离，船载雷达的目标探测距离更短，因而有效应对时间更少。

陆基雷达处理的数据量通常远远高于船载雷达，但传统陆基雷达对于提升数据处理速度的解决方案不适用于船载雷达。通常陆基雷达产品生成一般分为前端和后端两部分，前端计算机用来控制雷达及其机构并收集反射信号，后端计算机用来将雷达信号反演成图像产品。在陆基雷达的反演过程中后端计算机通常为计算能力较强的超级计算机，通过改善硬件的方式提高反演速度。但由于小型船舶制造成本和能量环境等限制因素，应用改进硬件来提高反演速度的方法并不适用于目前的渔船和游艇等小型船舶的情况。

目前，国内海洋渔船一般采用基于（8～12 GHz）波段的小型雷达，通过天线发射电磁波到海面，电磁波与雷达波长相当的毛细波产生Bragg 散射，后向散射回波后被天线接收，形成“海杂波”，接收终端把连续接收的回波形成图像序列，而波长较长的重力波通过对毛细波的调制作用表现在海杂波图像上，从中提取风、浪、流等海洋环境信息[3-5]。本文分析传统雷达图像反演过程和成像步骤，利用对雷达图层封闭轮廓的重建，提出基于轮廓重建的船载雷达MAP-MRF 快速反演方法，解决渔船雷达图像反演过程中大量重复计算缺失数据的重建插值过程。通过统计力学模型对MRF 中格点划分进行精确建模，利用MAP-MRF 方法对基于轮廓重建的雷达图层边界进行估计。将整幅图像的插值重建转化为对图层轮廓的精确分割，实现雷达图像的快速反演，从而提高显示终端的刷新率，实现对海上极端气象目标的及时发现和连续追踪，满足渔民对天气气象的变化情况特别是坏气象实时了解的需求。在证明基于轮廓反演方法有效性的同时，验证应用基于统计力学模型的MAP-MRF 算法对雷达图像进行相位估计的可行性，为进一步从概率模型方向处理雷达图像提供新的分析思路和处理方法。主要内容分为3 个部分：一是通过对渔船雷达数据和反演过程的分析，明确反演方法改进重点和改进内容；二是对PPI 模式下的雷达基数据图像应用MAP-MRF 算法重建封闭轮廓进行精确重建，进而反演成像；三是将试验结果与较有代表性的传统反演方法在成像质量和运行速度上进行分析和比较。

1 研究内容

1.1 雷达图像反演

雷达图像的重建可分为基于时间维度和基于空间维度的两种过程，二者分别对应雷达反演过程中时间分辨率（temporal resolution）和空间分辨率（spatial resolution）。本文主要讨论单帧雷达图像的空间超分辨率重建方法。

雷达的扫描方式是建立在球坐标系中，但其展示和应用通常在笛卡尔坐标系下，因此雷达图像的反演处理中包含了两个坐标系下的数据坐标转换。当坐标系变换后，雷达的辐射结构会造成近雷达端数据采样密集混叠，远雷达端数据采样稀疏缺失的情况。因此，雷达通常会对坐标系变换后的雷达图像进行重建或插值操作来改善以上问题。将雷达图像重建算法分为雷达格点分析法和通用插值算法两类。

雷达格点分析法是在雷达和气象领域常用的数据插值和分析方法，如Barnes 分析法和Cressman 分析法等[6-9]。Barnes 分析法能够把分布不均匀的数据映射到规则网格坐标系的指数权重，在映射到笛卡尔坐标系时通常应用均一Barnes 插值，在插值到球坐标系时一般应用自适应Barnes 插值。重建雷达图像具有平滑、失真较小的优点，此外Barnes 分析法还可用于基于站点序列的插值和纠错。Cressman 分析法用于从基数据向格点空间的映射，通过加权影响半径内基数据对格点值进行计算，采用迭代影响半径的方式优化精度依据影响半径。Cressman 计算速度快且重建准确率高，但当映射空间密度高于基数据密度时算法稳定性较弱，边缘附近易出现极值点。

通用插值算法因应用范围广，在性能和比较分析方面有比较深入的研究[10-14]。其中最常用的算法包括最邻近插值法、双线性插值法、基于小波的插值法，以及近年来被广泛研究的基于机器学习理论的插值法等。其中最邻近插值法是最简单也是最快速的一种方法，计算复杂度最低，对于噪声敏感。双线性插值法是在实际应用中被广泛采用的算法，时间复杂度较低，峰值信噪比、均方差及主观评价等较最邻近插值法有明显改善。基于小波的插值法是利用小波分解后的高频信息预测高分辨率图像，重构高分辨率图像的方法，是应用变换域插值的一种主要方法。机器学习的图像超分辨率重建法属于基于样例的重建，其中比较有代表性的如DONG W 等[15]基于稀疏编码重建，以及SHI W 等[16]对卷积网络的深度训练图像重建等，对于普通图像深度学习方法具有较好的重建效果，但在气象雷达领域面临着探测目标通常不具有稳定的结构特征，重建内容缺少对应训练数据的问题。

在图像重建的研究中，基于边缘信息的重建是比较有代表性的方向。图像边缘表达了更加丰富的图像信息，通过解决重建过程中的边缘模糊效应和锯齿状不连续等问题能够有效改善整幅重建效果[17-18]。因此有效利用图像边缘信息是有效改善重建效果、提升重建效率的一个有效方向。受此启发，提出精确分割阈值图层轮廓实现雷达图像重建的方法。同时，渔船雷达图像是灰度图像，图像中的舰船目标和陆地区域灰度级较高，海面区域灰度值较低，目标与目标之间不会重叠，对此图像重建方法具有一定的自适应性。

1.2 空气运动描述模型

雷达图像描述的轮廓边缘既对应物理特征的等高线，又描述气象目标内部降水粒子的尺度和密度分布。这些粒子悬浮在大气中，其分布和运动特征可以用大气流体运动模型描述。

此类模型主要是指基于纳维—斯托克斯（Navier-Stokes，NS）方程的各类多相流模型，模拟方法以传统计算流体力学方法为主，光滑粒子方法、格子玻尔兹曼方法（Lattice Boltzmann Method，LBM）应用比较广泛[19]。

2004 年，Dankert 从雷达图像序列中反演海面风场，研究了以光流运动估计技术为基础的风场反演方法。此方法不需要校正雷达系统的相位，具有较强的实用价值。然而，由于反演风场的光流方法对图像质量的要求较高，以及风场反演中光流方程的求解方法不够完善等原因，反演效果仍不够理想，有待进一步改进[20]。而格子玻尔兹曼是基于正四边形或正六边形的均匀对称网络，但因流体压缩性和黏性等影响，在计算精度和计算效率上存在不足。

近年来，新的算法在精度和效率上有所提高但求解过程变得复杂，本文应用Ising 模型对MRF-MAP 计算框架中的势函数进行了替换[21]。首先LEE T D 等[22]在20 世纪50 年代证明了Ising 模型和格子玻尔兹曼模型在力学统计领域的等价性，其次Ising 模型更适应本文的二值分割应用，更方便计算。

2 轮廓重建快速反演算法

参照雷达压缩格式完成雷达基数据sweep{SKs}的读取之后可以得到基数据原始图像Sk，基数据原始图像为一张m×n大 小的灰度图，m为雷达有效探测距离内的径向分辨率，n为Sk扫描一周包含的雷达射线数，即Sk的雷达方向角分辨率。对基数据图像Sk进行PPI 模式的快速反演处理，提出反演过程如图1 所示。

图1 基于轮廓重建的快速反演流程Fig.1 Flowchart of rapid retrieval based on contour reconstruction

2.1 雷达图像预处理

2.1.1 图像去噪雷达接收机是一种处理回波信号的设备，可以对信号进行变换、放大等操作，从而提取接收的微小回波信号[23]。而在雷达接收机接收信号的同时，也会收到一些对设备灵敏度产生负面影响的噪声，这些噪声影响了雷达图像的质量，对准确发现和识别气象目标产生不良影响，因此雷达反演过程中通常需要去噪。在过滤噪声的同时也会过滤掉一些图像的细节特征，使图像纹理变得平滑，这对后续轮廓重建的计算量会产生直接影响，因此以去噪平滑并保持明显轮廓结构为标准，对比了均值滤波、小波阈值滤波和自适应均值滤波等方法，并选取了自适应均值滤波作为去噪处理方法[24]。

2.1.2 阈值化提取轮廓

在雷达图像显示中，通常将连续变化的单通道雷达灰度图像阈值转化成伪彩图显示，这在各阈值图层之间形成了完整的封闭轮廓。将0～70 dBZ 级雷达反射图像阈值化为4-bit 的彩色图层Sk{Th0...Th15}, 并提取图层轮廓。

2.1.3 坐标变换

在PPI 模式坐标变换中，基雷达图像在球坐标下转换为笛卡尔坐标，只对各图层轮廓线像素进行计算，坐标变换后的轮廓线将由连续变为不连续。

2.2 基于MAP-MRF 的封闭轮廓重建

图层轮廓线在坐标变换后从连续线变为离散点，为重建封闭轮廓，以轮廓线相邻像素点为对角，将对角矩形范围内像素点进行双线性插值填充，形成随机场。尝试应用MRF 分割算法来找到随机场中两个相邻图层的精确边缘。以下阐述基于Ising 统计力学模型应用MRF-MAP 算法对随机场进行估计。

2.2.1 定义MRF

在随机场中的轮廓线是一组可以用马尔可夫性描述的随机过程，因此用MRF 方法对轮廓线两侧的格点分布进行估计。定义S为一个二维格点位置集，s为这组位置集中的一个格点，s∈S。设x为一个随机变量，x=0～255，将图像X定义为一个随机场，x∈X，xs为格点s对应的像素值，δs记作与格点s相邻的其他格点，并且相邻关系符合对称性，既s∈δγ⇒γ ∈δs。定义子团（clique）C为一组具有这样相邻关系的格点。那么随

在统计热力学中，一般认为能量越低系统越稳定，出现概率越高。

2.2.2 MRF-MAP 计算框架

在完成图层随机场定义后，应用MRF-MAP 框架对图像中真实的x进行估计。假设y为图像观测值，那么MAP 的求解为

在Ising 模型中粒子间的相互作用关系由他们的自旋方向决定，自旋产生了“↑”和“↓”的两种状态，可以描述图像的二值像素。通过Ising 模型的物理概念可知自旋方向相反的粒子对越多，系统的能量越高，越不稳定，处于这种状态的概率越低。模型中系统能量用自旋方向相反的粒子对相交边界的边长描述[27]。

2.2.3 Metropolis 算法

Metropolis 算法为一种固定温度的模拟退火算法，是通过赋予搜索过程一种时变且最终趋于零的概率突跳性，从而可有效避免陷入局部极小并最终趋于全局最优串行结构的优化算法[28]。应用Metropolis 算法在图像中计算能量函数的全局最优值，具体步骤如下。

2.3 阈值图层集成

在完成雷达图层轮廓线的重建后，对轮廓线内区域进行漫水填充赋值。将整个图像设置为阈值颜色，然后对背景进行填充。因为漫水填充法需要对封闭轮廓内种子像素点标记区域进行填充，而雷达阈值图层的前景多为不连续区域，所以采用背景填充。在完成每张雷达图层后，按阈值从小到大顺序将雷达图层集成为雷达反演图。

3 反演试验和分析

3.1 试验数据和仿真环境

硬件环境为Intel Core2 E7400 CPU 和8Gb 内存，软件环境为Windows7 操作系统，VS2012 编译和OpenCV 库。采用美国CASA（Collaborative Adaptive Sensing of the Atmosphere）项目的X 波段多普勒雷达产品为数据[29]。雷达探测半径40～60 km，具有可对近地空间观测、安装运行成本低、可陆基布置也可车船安装的优点。雷达数据压缩格式为“uf”，数据采集地点为某两个雷达站（Lat：34.812 9，Lon：—97.931 3；Lat：35.031 2，Lon：—97.956 7），采 集 时 间 为2009年5 月14 日，扫描模式为仰角2°～4°的PPI 模式。雷达基数据文件包含CZ 和VE 两个字段，即反射率字段和风速字段，采用反射率字段为例进行试验。在反射率字段中，雷达的方向角分辨率为426，径向分辨率为996，雷达基数据矩阵成像如图2 所示。

图2 雷达基数据文件成像Fig.2 Radar base data file imaging

3.2 雷达反演过程

首先对图像进行去噪处理，去噪处理不仅能提高图像质量，还能平滑图像轮廓，去掉在基数据图像中不能形成封闭轮廓的像素点，减少后续处理的计算量。分析3 种常用去噪结果，并应用自适应均值法，去噪效果对比如图3 所示。

图3 雷达基数据图像去噪比较Fig.3 Comparison of radar cardinality according to image denoising

对处理后的基数据图像进行了4-bit 图层的阈值化并提取了图层轮廓。

在向PPI 模式的坐标变换中，不同于其他方法，只对轮廓线像素位置进行计算。在变换后为了将离散的轮廓线像素点连接起来，以两个相邻像素点为对角的矩形区域进行了双线性插值，使图层重新成为封闭区域。

对连续的矩形区域应用MAP-MRF 算法进行精确轮廓线的计算，保证了雷达远端的轮廓线精度。对封闭轮廓区域内进行漫水填充如图4 所示，并将图层进行合并，形成完整的PPI 模式雷达扫描图像。

图4 雷达图层反演结果Fig.4 Retrieval results of radar layers

3.3 试验结果对比分析

将快速雷达图像反演结果在图像质量和处理时间上与其他雷达图像重建方法进行对比，对比方法包括最邻近插值成像、Barnes 分析成像及双线性插值成像，结果如图5 所示。

图5 雷达图像反演结果Fig.5 Radar image retrieval results

在成像质量的对比分析上，通常依据人的主观感受及一些客观评价指标，如峰值信噪比（PSNR）和均方误差（MSE）等。但因为快速反演方法的图像区域大部分为填充区域，所以为了客观评价反演效果，采用直接验证法，以单束雷达基数据作为真值进行分析。在对比分析中，采样了两个雷达站扫描空域的交集区域，以阈值化后的一个雷达站基数据为真值，对另一个雷达站雷达波束上的400 个像素进行了采样统计。分析不同重建方法处理后的图像中像素能够被正确分类的比例，结果如表1 所示[30-32]。在重建方法对比中，准确率最高的是Barnes 分析法，本文方法优于最邻近插值法，在延单一雷达照射波束上的准确率统计与其他两种方法差距不大。

表1 数据重建准确率对比Tab.1 Data reconstruction accuracy comparison

在算法效率对比中，为了便于分析，对雷达基数据图像进行像素级的插值重建处理。基数据图像大小为996 像素×426 像素，重建图像为2 000 像素×2 000像素，重建图像像素数量约为原始图像9 倍。对10 组基数据文件进行了处理，测试样本数据阈值化后图层为8～13 层。在未对代码进行优化的情况下测得运行时间如图6 所示。

从运行时间来看，基于轮廓重建的MAP-MRF 方法能大大缩短雷达图像的反演时间。但由图6 可知，与其他方法不同，本文方法测试样本的运行时间方差较大。主要影响因素是雷达基数据的范围和图层中轮廓线划分的随机场大小。通过试验得出划分阈值后的图层数量对算法处理时间影响较大。轮廓线像素数量影响划分随机场的大小，进而图像需要的迭代时长不同。

图6 文中方法与传统方法运行时间对比Fig.6 Processing time of proposed method compared with traditional methods

4 结束语

海洋捕捞业一直有“起风鱼”之说，即在坏天气系统如台风侵入前，鱼群往往大而集中并浮表层，因为上述天气系统将带来大风，渔民抓住这一时机进行“迎风头”生产，鱼获量常常较大。但是此时作业常常会伴随大风大浪，重大渔船海损事故几乎都是在此气象条件下发生的，因此准确掌握气象变化显得尤为重要。针对渔船类小型船舶上雷达探测显示时间间隔长，应对海上极端气象目标时间短的问题，提出了基于轮廓重建的MAP-MRF 快速雷达图像反演方法。应用插值过程针对轮廓线计算量较少的原理，阐述了Ising 模型对大气流体运动建模的可行性，引入该模型对轮廓线的分布进行描述，替换了传统的格子气模型对介观尺度大气流体运动的描述。通过MAP-MRF 算法进行轮廓线内外区域的分割，以及CASA 试验场真实数据，验证了该方法的插值准确性。与气象和图像插值领域主要方法进行对比，论证了针对轮廓进行快速反演插值的方法具有更快的处理速度，有利于渔船雷达系统达到更高的刷新率，实现对海上极端气象目标的及时发现和连续追踪。通过试验分析，该方法对于雷达图像的反演速度有较大提升，优化效果受雷达探测对象的空间分布情况影响较大。