一种基于模糊聚类预处理云系的并行云迹风反演方法

摘要

本发明公开了一种基于模糊聚类预处理云系的并行云迹风反演方法，包括：选取相邻等时间间隔的三幅云图，对云图图像进行数据提取及预处理，将每幅云图分为多块追踪区块，并且将追踪区块分配到指定的多个并行计算处理单元，在每个独立的计算处理单元通过最大相关系数法对所述追踪区块的风矢进行反演计算，得到所述追踪区块风矢的风速及风向。本发明有效提升了云迹风反演的实际执行效率，对低、中、高空全球风场的预测和预警具有重要的应用价值。

说明书

技术领域

本发明涉及卫星云迹风反演领域，具体涉及一种基于模糊聚类预处理云系的并行云迹风反演方法。

背景技术

在天气和灾害预报的研究中，需要大量的风场资料，全球大气风场的观测起到了非常积极的作用。近年来，随着陆地风场观测站网的不断完善，研究人员获得陆地风场资料变的越来越容易。然而，相对于陆地风场资料来说，海洋、高原、沙漠等地区的观测网就显得非常不足，获取这些地区的风场资料较困难。通过卫星云图资料反演大气风场，无疑对获取全球大气风场资料起到了非常有益的补充作用。

同时，一些陆地或水体(称为地表)由于上空无云或少云而直接“裸露”在卫星云图中。研究人员主要是对云图中的云进行研究，例如云分类、云识别、云导风等，因此实现云区域和地表区域的图像分割是必不可少的。针对卫星云图图像分割的研究工作早已展开，目前使用的方法普遍具有以下缺陷：(1)过度依赖于历年的云图图像资料，需要花费人力物力收集和保存；(2)使用基于阈值的图像分割来实现，属于硬性分割，准确率难以保证；(3)云图是一个十分复杂的图像，不仅包含云和地表，还在云和地表之间存在一些处于“过渡阶段”的云(例如生长和消散的云)；所以，在传统方法中，容易将其强制划分到某一类中，但这是不合理的，对云图图像合理分割也是研究的方向之一。

一般来说，卫星云图时间间隔越短、用于反演的示踪云越多、风矢的计算密度越高，云导风反演出的风场资料质量也会越高，但这都对云导风反演算法的效率提出了挑战。虽然有学者提出了云导风反演算法的简易算法，但也难以满足日趋要求实时性的风场计算的效率要求。因此，除研究较低计算代价的反演算法，充分利用现代计算机的体系结构，研究并行云导风反演算法，也是提升反演算法实际运行效率的有效途径之一。

GPU(Graphic Processing Unit)，图形处理器是显卡的“心脏”，相当于CPU在微机上的作用。GPU具有很高的传输带宽和众多的计算单元，计算能力非常强大。随着GPU的不断更新升级，GPU的计算能力越来出众，这为云迹风的反演提供了一种并行加速思路；GPU在高性能运算方面拥有三个非常显著的优势：第一，数据的并行处理能力强大，NVIDIA公司目前最新的GPU峰值浮点运算能力超过3TFLOPS，这几乎和一个小型CPU集群相当；第二，GPU拥有很高的数据吞吐能力，K20的带块超过200GB/s；第三，GPU具有友好的编程接口，NVIDIA公司开发的CUDA框架支持多种高级编程语言，如C、C++和Fortran语言。

综上，虽然在现有技术中，有文章阐述基于多核CPU的卫星云导风并行反演算法，但由于CPU擅长处理具有复杂逻辑计算步骤的任务，如数据压缩、人工智能、物理仿真等，而GPU擅长处理计算密集型且逻辑较简单的重复性任务，如图像领域相关处理操作，两者在物理结构和数据处理上存在很大的不同，在GPU上应用卫星云迹风反演方法是需要克服很多困难和技术障碍，GPU虽然拥有成百上千个计算单元，但每个计算单元的计算能力薄弱，难以执行复杂的分支预测。GPU对显存的存取模式与CPU显著不同，缺少完善的共享存储逻辑，各级Cache容量较小，这给传统CPU上的复杂算法向GPU并行算法移植带来了理论和技术实现上的困难。如何利用GPU强大的计算能力，对高密度、高实时性的卫星云迹风反演进行并行处理是一个挑战。

发明内容

本发明设计开发了一种基于模糊聚类预处理云系的并行云迹风反演方法，本发明的目的之一是解决串行云迹风反演方法中由于迭代次数较多导致的云迹风反演时间长、实时性差的问题，通过并行反演方法将单次迭代并行化，减少迭代次数，进而显著提高云迹风反演的实时性。

本发明的发明目的之二是解决了在规定的搜索区域内对目标云块搜索时间过长导致的云迹风反演实时性差的问题，通过优化示踪云块的追踪方法，大幅度简化对搜索区域目标云块的搜索时间，进而显著提高云迹风反演的实时性。

本发明的目的之三是解决了现有技术中对需要进行反演的云图图像预处理仅仅是采用图像硬分割方法，进而存在的对图像预处理过程中难以取得一个合适的阈值，而且受环境因素影响，准确率难以保证的问题。

本发明提供的技术方案为：

一种基于模糊聚类预处理云系的并行云迹风反演方法，包括：

将云图图像分为多块追踪区块，并且将追踪区块分配到指定的多个GPU并行计算阵列，在每个独立的计算阵列通过最大相关系数法对所述追踪区块中的目标云块进行匹配后进行反演计算，得到所述追踪区块风矢的风速及风向；

在通过所述最大相关系数法对目标云块的匹配过程包括：在所述追踪区块中选中示踪云块，在包含所述示踪云块的搜索区域内，计算出搜索步长大于1的全部与所述示踪云块像素等同区块的相关系数，在这些相关系数中确定最大相关系数，即为与所述示踪云块相似度最高的区块，以该区块为中心增加像素点的区域内，将步长逐步减半再次搜索出相似度最高的区块，直到搜索出步长为1搜索到的最大相关系数的区块，即为目标云块；以及

其中，通过模糊聚类方法对原始的云图图像进行云地分离后确定进行云迹风反演的云图图像，其包括：

采集云图图像中红外通道的灰度特征数据，对所述灰度特征数据进行处理，获得样本灰度特征数据，将所述样本灰度特征数据进行分类，并使每一类样本灰度特征数据分别对应一种云图类型，进而确定所述云图图像。

优选的是，还包括，选取相邻等时间间隔的三幅云图，分别得出每幅云图风矢的风速及风向，将所述风矢的风速、风向以及根据时序将所述云图中相邻云图的风矢速度差、方向差与阈值进行比较，执行以下原则：

若风矢速度小于指定的最低速度阈值，则剔除此风矢；或

若两风矢速度之差大于指定的最大速度差阈值，则剔除此风矢；或

若两风矢方向之差大于指定的最大方向差阈值，则剔除此风矢；或

若两风矢的周围区域内没有与之方向差小于指定的最大方向差阈值的风矢，则剔除此风矢。

优选的是，对所述云图图像数据进行多通道数据提取，保留红外通道数据。

优选的是，还包括：对所述云图图像进行预处理，其包括对图像进行等经纬度投影变换和将原始圆盘图变换为矩形图。

优选的是，对所述云图图像进行预处理还包括对所述图像进行灰度增强。

优选的是，对所述灰度特征数据进行处理包括：

将云图图像中红外通道的灰度特征数据输入得到像素灰度值矩阵，通过计算在矩阵中主对角线像素灰度均值，通过计算在矩阵中辅对角线像素灰度均值，进而通过M＝(M₁+M₂)/2求出像素灰度均值，即样本灰度特征数据。

优选的是，在所述红外通道内选取第一长波红外、第二长波红外以及短波红外的样本灰度特征数据，将所述样本灰度特征数据建立5个灰度特征数据统计量，对原始云图进行图像分割得到所述云图图像。

优选的是，根据每幅云图像素的宽高以及风矢密度计算确定所述追踪区块个数，其为(Width/ρ)×(Height/ρ)，式中，Width为所述每幅云图像素的宽度，Height为所述每幅云图像素的高度，ρ为风矢密度。

优选的是，在通过所述最大相关系数法对目标云块的匹配过程包括：在所述最大相关系数法中，在所述追踪区块中选中示踪云块A，其为N×N像素图像块，在包含所述示踪云块A的M×M像素图像块搜索区域内，计算搜索步长为m的全部与所述示踪云块A像素等同区块的相关系数，在这些相关系数中确定最大相关系数，即为与所述示踪云块A相似度最高N×N像素图像区块A1，在所述区块A1增加个像素点，然后将步长减为个像素点，在以区块A1为中心的边长的区域中计算各N×N像素图像区块的相关系数，在这些相关系数中确定最大相关系数进而确定下一个与所述示踪云块A相似度最高N×N像素图像区块A2，将步长逐步减半直到搜索出步长为1搜索到的最大相关系数的区块，即为目标云块，通过所述示踪云块以及所述目标云块的位置信息计算得到风矢的风速及风向。

本发明与现有技术相比较所具有的有益效果：

1、将卫星云迹风反演算法中的串联算法改为并联算法，大大的降低了迭代次数，通过这种单次迭代并行化的计算方法，从而使整张云图中的风矢可以在GPU中并发执行，进行大规模数量的风矢运算，在不影响反演出的云导风产品的精度及性能的情况下，大大节省了反演和风矢校验所需的时间，进而大幅度的提高了风矢计算效率，显著提高了云迹风反演的实时性；

2、在本发明中，采用变步长迭代计算的方式，大幅度降低在搜索区域对目标云块的搜索时间，进而大幅度提高追踪区块的风矢计算效率，显著提高云迹风反演的实时性。

3、在本发明中，通过模糊聚类方法运用于云图图像分割预处理能够取得更好的云地分离效果，该方法能够有效进行云图不确定性分割和提升分割准确率，在云地分离后的图像上进行云迹风反演，使得反演效果更加准确，同时在进行模糊聚类预处理之前，对灰度特征数据通过快速中值滤波方法或者对角线滤波方法进行预处理去除了椒盐噪声，调高了图像质量，使目标清晰，排除了对云迹风的计算和识别工作中产生的噪声干扰，进一步的提高了反演效果。

附图说明

图1为CPU和GPU在基于GPU平台的卫星云迹风反演方法中的任务划分示意图。

图2为云图图像的FCM结果图。

图3为基于FCM的云图图像分割流程图。

图4为云迹风反演算法并行后的基本流程图。

图5为流处理器数量对效率的影响图。

图6为基于GPU云导风串行和并行反演算法精度对比1的示意图。

图7为基于GPU云导风串行和并行反演算法精度对比2的示意图。

图8为基于GPU并行云导风算法反演出的云导风产品图。

图9为执行时长随线程块增加的趋势图。

图10为加速比随线程块增加的趋势图。

图11为示踪云块追踪示意图。

图12为风矢速度推导示意图。

图13为风矢方向推导示意图。

图14为步长为8的迭代搜索示意图。

图15为步长为4的迭代搜索示意图。

具体实施方式

下面结合附图对本发明做进一步的详细说明，以令本领域技术人员参照说明书文字能够据以实施。

如图1～图4所示，本发明提供了一种基于模糊聚类预处理云系的并行云迹风反演方法，包括如下步骤：

步骤一、使用相邻等时间间隔的三幅红外通道云图，设定反演参数，包括卫星遥感数据在磁盘上的读取路径、云迹风反演的密度、云迹风反演算法相关参数、GPU并行参数、云迹风校验力度参数、云迹风反演产品在磁盘上的保存路径等；

步骤二、使用无损解压缩算法对从气象卫星荷载下发的地球遥感数据进行解压缩；

步骤三、对卫星红外图像数据提取，即对卫星遥感数据解压缩流程产生的卫星多通道数据进行数据提取，保留红外通道数据，舍弃其他通道数据，并且开辟该云图宽高的存储空间，将提取的红外通道的云图数据复制到该存储空间；

步骤四、对提取后的云图数据进行预处理，包括图像投影变换、图像灰度增强以及通过模糊聚类方法对提取的云图图像进行云地分离后确定云图图像；

步骤五、对执行环境进行初始化，包括对反演参数和GPU设备的执行环境的初始化，其包括：根据所述每幅云图像素的宽高以及风矢密度计算需要反演的风矢个数，在GPU的显存中开辟用于保存云迹风的结构体数组，在GPU的显存中开辟用于保存各风矢使用最大相关系数而产生的相关系数值的数组，在GPU的显存中开辟用于存储3幅红外云图数据的数组，CPU将主存中的三幅云图数据依次复制到GPU显存开辟的三个数组中。

步骤六、进行云迹风并行反演计算，CPU将一系列参数和多任务分配和调度规则载入GPU设备，并触发GPU每个线程对应一个风矢的并行反演计算，将每幅云图分为多块追踪区块，并且将追踪区块分配到指定的多个并行计算处理单元，在每个独立的计算处理单元通过最大相关系数法对所述追踪区块中的目标云块进行匹配后进行反演计算，得到所述追踪区块风矢的风速及风向，并保存到GPU的显卡中；

其中，每幅云图像素的宽度、高度以及风矢密度计算确定所述追踪区块个数，其为(Width/ρ)*(Height/ρ)，式中，Width为所述每幅云图像素的宽度，Height为所述每幅云图像素的高度，ρ为风矢密度。

步骤七、对云迹风进行质量校验，校验的内容为风矢速度和风矢方向，校验的对象是3幅相邻时间序列的云图中的前两张云图反演出的结果和后两张云图反演出的结果，将风矢的风速、风向以及根据时序将三幅云图中相邻云图的风矢速度差、方向差与阈值进行比较，按以下原则执行：

若风矢速度小于指定的最低速度阈值，则剔除此风矢；或

若两风矢速度之差大于指定的最大速度差阈值，则剔除此风矢；或

若两风矢方向之差大于指定的最大方向差阈值，则剔除此风矢；或

若两风矢的周围区域内没有与之方向差小于指定的最大方向差阈值的风矢，则剔除此风矢。

步骤八、对上述云迹风产品进行数据保存，得到最终反演产品。

在另一种实施例中，在步骤一中，反演参数中的遥感数据读取路径和产品保存路径均为绝对路径；反演密度的单位是每1度经纬度网格块所包含的反演风矢个数；反演参数中包括示踪云块和追踪区域行列的像素个数；GPU并行参数为GPU内核函数的执行单元的网格划分规则；云迹风校验力度参数为云迹风质量等级，质量等级越高，要求越严格，反演出风矢的可信度高，但个数越少，反之质量低，个数多。

在另一种实施例中，在步骤二中，卫星遥感数据解压缩算法可采用基于游程编码的解压缩算法、基于字典编码技术的LZW解压缩算法、基于哈夫曼编码原理的解压缩算法或者基于算术编码的解压缩算法。

在另一种实施例中，在步骤四中，对图像进行等经纬度投影变换，变换成墨卡托投影，以便进行风矢反演，将卫星原始圆盘图变换为矩形图，方便GPU任务的划分和执行，以提升反演质量；

对图像进行灰度增强处理可采用线性灰度变换增强算法、分段线性灰度变换增强算法、直方图均衡化增强算法或者局域直方图均衡化增强算法对图像的对比度进行增强；

在另一种实施例中，在步骤四中，通过模糊C均值聚类(FCM)对原始的云图图像进行云地分离后确定所述云图图像，其包括如下步骤：

输入：卫星云图图像

输出：消除了地表区域的云图图像

步骤a、首先对云图进行预处理，处理的方法为使用中值滤波法从而去除云图图像中的噪声和干扰。

步骤b、构建FCM的目标数据集。选取红外云图和可见光云图构建IR-VIS二维光谱特征空间，使用一个二维数组data[256,2]按列分别存储每个像素点在红外、可见光图像中的灰度特征值。

步骤c、计算样本一般特征值。若数据库中没有记录8种地表、云类的灰度特征，则对云图样本库进行分析，统计出8种地表、云类的灰度特征并在数据库中记录。若有，则忽略此步，继续执行。

步骤d、将初始聚类中心设置为样本特征值进行FCM。从数据库中获取8类地表、云类的灰度特征作为FCM的初始聚类中心，并在红外-可见光二维光谱特征空间进行聚类。

步骤e、对云图进行图像分割。利用由步骤(4)得到的聚类分区图对所有像素点进行判断，若其二维灰度值在地表灰度特征区域内，则为地表像素点，标记为黑色；若不在，则为云像素点，保留。

最后得到一张矩形的、灰度增强的以及无效区域提出后的云图。

在另一种实施例中，在步骤六中，CPU将一系列参数和多任务分配和调度规则载入GPU设备，并触发GPU每个线程对应一个风矢的并行反演计算，GPU并行执行所有风矢的反演，每个风矢的反演都在独立的线程中执行，风矢的反演使用最大相关系数法计算目标云块和追踪云图的位置关系，从而确定风矢速度和方向，最后将风矢所在的像素位置转换为实际经纬度位置。

在另一种实施例中，在步骤六中，云导风反演并行算法采用进程－线程的Fork-join模型，主进程读取云图数据和参数并生成网络剖分节点间关系后，通过OpenMP编译指导语句，将n个风矢的计算负载动态分配到m个线程上，利用最大相关系数法匹配目标云块、计算风矢的风速和风向、风矢定高、风矢质量校验等步骤由各线程独立计算，再使用最大相关系数法匹配目标云块，写入全局风矢数组步骤中，各线程共享云图数据内存区域和风矢数组内存区域.所有风矢计算完毕后，通过编译指导语句结束线程，主进程将计算好的风矢数据写入磁盘保存。

在另一种实施例中，在步骤六中，最大相关系数法是指在卫星云图追踪区域内寻找和追踪云块相似度最高的云块，云块之间的相似度则由相关系数作为指标，即在云图追踪区域内寻找和追踪云块计算出的相关系数最大的云块，具体包括以下几个步骤：

(1)将云图1分割为多个追踪区块，依次遍历云每个追踪区块；

(2)在当前GPU线程中分配的风矢反演区块中定位中心位置的像素点区域，记录中心点位置；

(3)计算(2)中区域的灰度方差和标准差；

(4)从云图2中的相同区块的左上角的第一个像素点区域开始遍历，直至遍历结束，返回(1)，开始从上到0，下从左到右遍历，直到遍历结束，计算该相同区块的方差和标准差；

(5)计算(3)和(4)的相关系数，并将结果写入存储矩阵；

(6)找出(5)的存储矩阵中的最大相关系数，记录其对应位置；

(7)根据(2)和(6)的位置坐标，计算该区域块的风矢距离和角度，用风矢的距离除以时间间隔，则可获得风矢速度；将反演出的风矢的角度转化为北偏角的形式即为风矢的方向。

具体的说，在步骤六中，在通过所述最大相关系数法对目标云块的匹配过程包括：在所述最大相关系数法中，在所述追踪区块中选中示踪云块A，其为N×N像素图像块，在包含所述示踪云块A的M×M像素图像块搜索区域内，计算搜索步长为m的全部与所述示踪云块A像素等同区块的相关系数，在这些相关系数中确定最大相关系数，即为与所述示踪云块A相似度最高N×N像素图像区块A1，在所述区块A1增加个像素点，然后将步长减为个像素点，在以区块A1为中心的边长的区域中计算各N×N像素图像区块的相关系数，在这些相关系数中确定最大相关系数进而确定下一个与所述示踪云块A相似度最高N×N像素图像区块A2，将步长逐步减半直到搜索出步长为1搜索到的最大相关系数的区块，即为目标云块，通过所述示踪云块以及所述目标云块的位置信息计算得到风矢的风速及风向。

在另一种实施例中，在所述步骤八中，云迹风产品保存方法为CPU将GPU的显存中的风矢数据复制到计算机主存中，然后输出只外存磁盘中。

在另一种实施例中，本发明反演的数据对象为FY-2E卫星红外一通道的云图数据，该云图的大小为2581*1399像素，即Width为2581像素，Height为1399像素，选取等时间间隔相邻的三张FY-2E卫星红外一通道的云图作为本实施例的数据处理对象。

在另一种实施例中，在本发明中的反演参数接收流程，输入卫星遥感数据读取路径Path1和卫星云迹风产品保存路径Path2；输入反演的密度ρ；输入反演中的示踪云块的边长，即示踪云块大小为P*P；输入云迹风反演等级D。

在另一种实施例中，在本发明中的卫星遥感数据解压缩流程，根据上面提供的云图读取路径Path1，读取到FY-2E卫星的遥感数据，使用无损解压缩算法对FY-2E卫星遥感数据进行解压，解压得到各通道的云图像素数据，温度定标数据，经纬度定标数据等。

在另一种实施例中，在本发明中的卫星红外一图像数据提取流程，从上述解压得到的各通道的云图数据中，分别提取三张云图的卫星红外一云图数据，复制到微机主存开辟的数组Data1、Data2和Data3中。

在另一种实施例中，图像灰度增强处理使用直方均衡化算法对云图数据进行对比度增强。

在另一种实施例中，在本发明中的执行环境初始化流程，初始化CPU和GPU的执行环境和内存数据。根据密度和云图的宽高计算，需要反演的云迹风个数n；CPU调用CUDA接口在GPU的共享内存中开辟n个用于保存云迹风的结构体数组windRecords；CPU调用CUDA接口在GPU的共享内存中开辟n个大小为(P+1)*(P+1)的用于保存各风矢使用最大相关系数而产生的(P+1)*(P+1)个相关系数值的数组；CPU调用CUDA接口在GPU的共享内存中开辟三个用于存储三幅红外一云图数据的数组，数组的大小根据权利要求4中所述，为云图的大小Width*Height；CPU调用CUDA接口将主存中的三幅云图数据依次复制到步骤1中GPU设备的共享内存开辟的三个数组中。

在另一种实施例中，在本发明中的云迹风反演流程，CPU将云图上每个风矢的反演划分到一个任务线程中，将所有线程划分到一个GPU执行网格中，GPU根据线程格的任务分配和调度原则，对所有风矢的反演并发执行，反演使用最大相关系数法确定风矢的位移，进而确定风矢的速度大小和方向。最后将风矢所在的像素位置转换为实际经纬度位置，将反演得到的风矢按个执行线程分配的索引位置写入GPU共享内存中开辟的WindRecords数组中。

在另一种实施例中，在本发明的云迹风产品保存流程，CPU将GPU的共享内存中计算出的质量校验后的云迹风数据复制到宿主微机的主存中，然后CPU将主存中的云迹风数据进行保存，保存路径为上述反演参数接收流程接收的云迹风产品保存路径Path2，保存的格式为：风矢编号、风矢经度、风矢纬度、风矢速度以及风矢方向。

下面再结合具体的实施例对并行反演方法以及最大相关系数做具体的说明。

在另一种实施例中，如图2、图3所示，用基于密度的聚类方法对云图进行云地分离；根据云图光谱数字图像像素点间灰度特征的“相似性”对像素进行划分来达到云区域与地表区域分割的目的。传统的阈值方法属于图像硬分割方法，不仅难以取得一个合适的阈值，而且受环境因素影响，准确率难以保证。该方法能有效进行云图不确定性分割和提升分割准确率，在云地分离后的图像上进行云迹风反演，使得反演效果更佳准确；在卫星数字图像的生成、处理和投影变换等过程中，云图通常都不可避免产生一些椒盐噪声，这些噪声不但降低了图像质量，还使目标模糊，对云迹风的计算和识别工作产生了很大的干扰，因此在卫星云图图像中去除噪声和降低干扰是必要的。

通过中值滤波方法对卫星云图图像中灰度特征数据进行去除噪声和降低干扰处理得到灰度样本特征数据，具体如下：

用方形滤波模版，定义了如下二维中值滤波公式：

Z_xy＝Med{g(i，j)}＝Med{g(x+r，y+s)，(r，s)∈A(x，y)∈I²}

其中，Med表示取中值操作；处理时把P(x，y)的灰度值设置为Z_xy即可，如表1所示。

表1 3×3中值滤波模块

第0列第1列第2列第0行P₀P₁P₂第1行P₃P₄P₅第2行P₆P₇P₈

每计算一次中值对模板内各点的灰度值进行排序，根据行或列把模板内的像素分成三组(每组3个像素)进行中值计算，得到的中值即为样本灰度数据特征，使用快速中值滤波方法能够减少像素之间的比较次数，并且可以通过并行实现，极大地提高图像去噪的处理速度，快速中值滤波方法具体的包括：

输入：3×3窗口内的像素灰度值

输出：9个像素的灰度中值

(1)对模板内的每一行，分别求取最大值、最小值和中值，得到数据组A＝[Max0，Max1，Max2]、B＝[Min0，Min1，Min2]、C＝[Med0，Med1，Med2]，如下所示。

A：Max0＝max[P₀，P₁，P₂]，Max1＝max[P₃，P₄，P₅]，Max2＝max[P₆，P₇，P₈]

B：Min0＝min[P₀，P₁，P₂]，Min1＝min[P₃，P₄，P₅]，Min2＝min[P₆，P₇，P₈]

C：Med0＝med[P₀，P₁，P₂]，Med1＝med[P₃，P₄，P₅]，Med2＝med[P₆，P₇，P₈]

(2)求A中的最小值：Maxmin＝min[Max0，Max1，Max2]；

(3)求C中的中值：Medmed＝med[Med0，Med1，Med2]；

(4)求B中的最大值：Minmax＝max[Min0，Min1，Min2]；

(5)求med[Maxmin，Medmed，Minmax]，即为所求中值。

理论上，我们只需要从云图中分割出云区域和地表区域，也就是说把聚类数目设置为2。但是这样得到的云图分割效果极不理想，存在云区域被误划分为地表的错误，因此有必要增加聚类数目来提高云图图像分割的准确率，在本实施例中把云图图像中物体类型划分成8类，即陆地、水体、积雨云、浓积云、卷云、中云、淡积云和低云，因此本文把FCM聚类的划分数量设置为8，通过这样的选择，能够充分利用气象领域的专业知识，约束并指导聚类过程，提高聚类的准确率。极大地减少循环的次数，从而使算法的效果更显著。

在本实施例中，选取红外(IR)以及水汽(WV)两个通道云图的样本灰度数据特征，建立这8类样本的灰度特征统计量(若干样本的均值)，如表2所示。

表2 云图样本的灰度特征统计

确定聚类数目为8，并把样本灰度特征统计量作为FCM聚类的初始中心，以约束和指导聚类过程，聚类迭代过程中不断调整优化聚类中心和隶属度矩阵，最后得到二维特征空间的聚类结果。

具体FCM算法如下：

输入：聚类数目c＝8，数据数量n；固定值λ，1≤λ≤∞；阈值ε；

输出：隶属度矩阵U'＝[u_ij]_n×n

(1)把地表和典型云类样本特征量的均值赋值给V₀＝[v_ij]_2×8；

(2)令循环次数b＝0，利用步骤(5)中公式为隶属度矩阵U⁽⁰⁾赋予初值；

(3)b++；

(4)计算c个聚类中心V_i^(b)(i＝1,2,...,c)，公式如(1-1)。

$V_i=\left(\underset{j=1}{\overset{n}{\Sigma }}{u}_{ij}^{\lambda }{X}_j\right)/\left(\underset{j=1}{\overset{n}{\Sigma }}{u}_{ij}^{\lambda }\right)---(1-1)$

(5)对于j＝1～n，重新计算U^(b)，并使U^(b+1)＝U^(b)；

①计算I_j和

I_j＝{i|1≤i≤c；d_ij＝||X_j-V_i||＝0}>

$\overline{I_j}=\{1,2,...,c\}-I_j---(1-3)$

②对于数据X_j，更新其隶属度。如果则对所有令u_ij＝0，j＝j+1；否则

$u_{ij}=1/\underset{k=1}{\overset{c}{\Sigma }}{(d_{ij}/d_{kj})}^{2(\lambda -1)}---(1-4)$

(6)比较U^(b)和U^(b+1)。若执行步骤(7)；否则返回(3)；

(7)令U’＝U^(b+1)；结束。

FY-2E气象卫星具有可见光、红外和水汽三个光谱云图，因此云图中每个像素点都具有可见光、红外和水汽三种灰度特征属性；FCM在二维平面空间中既有较直观的效果又有较高的执行效率，在本实施例中，从三种属性选择红外及可见光进行具体分析，使用红外通道及可见光通道构建正交的二维光谱特征空间，使图像上像元与二维特征空间中的点一一映射，即云图像元为R_I(m,n),R_V(m,n)，二维光谱特征空间点为S_IV(R_I,R_V)，通过如上述的具体FCM算法对二维灰度特征空间进行FCM聚类。

对云图图像进行过渡云系及不可辨云系的划分：假设云图中的地表、云类别可细分为c类：(X₁,X₂,...,X_c)，S_IV(A_I,A_V)中每个点的隶属度为U(u₁,u₂,...,u_c)，则可根据下列条件确定云图中每个像素点所属的类别：若u_i＝Max(U)＞0.5，则判定点云图中像素点R_I(m,n)、R_V(m,n)属于类别i；若则该像素点介于两种类别之间，隶属度高于阈值下限并低于阈值上限，可归为过渡云系；若u_i≤0.3，则该像素点特征异常模糊，难以辨别，可归为不可辨云系。上述划分步骤中的阈值上下限可根据实际情况调整。

引入过渡云系和不可辨云系的概念，解决了图像中模糊的、不确定的像素点的归属问题，重新对云图进行FCM聚类，得到聚类效果如图2所示。

从图中可看出，隶属度高于阈值上限(0.5)的像素点在二维特征空间内形成8个明显的聚类区域(A、B、C、D、E、F、G、H)，这些区域的范围反映了表2中8类典型样本的灰度特征变化；X、Y区域分别表示过渡云系和不可辨云系；结合表2中的光谱数据，根据像素点在二维特征空间的分布，可确定A为陆地，B为水体，如图2所示。

设灰度二维特征空间中陆地、水体的聚类区域分别为S_IV(陆地)、S_IV(水体)，云图中像元(x,y)在红外和可见光云图中灰度值分别为P_I(x,y)、P_V(x,y)，图像分割后点(x,y)在红外和可见光云图中的灰度值分别是G_I(x,y)、G_V(x,y)。则结合基于阈值的图像分割思想，可建立使用FCM方法的图像分割的准则，如下所示：

在另一种实施例中，通过对角线滤波方法对卫星云图图像中灰度特征数据进行去除噪声和降低干扰处理得到灰度样本特征数据，具体如下：

云图使用灰度阶表示物体发出的辐射，像元具有的灰度值代表了云或地表表面的亮温值。亮温值越低则物体海拔就越高，如高层云；亮温值越高则物体海拔就越低，如陆地、海洋等。在生成、数字处理和投影变换等过程中，卫星云图通不可避免产生一些噪声，噪声不但降低了图像质量，还使目标模糊，更是给后续的图像分割、特征提取、云迹风反演等研究工作造成了不便，因此在卫星云图图像中去除噪声和降低干扰是必要的。

针对云迹风反演主体的特征，云的整体性尤其是云边界对云迹风反演结果影响较大，因此，我们提出既能弥补云团噪声又能保留清晰云边界的对角线(主辅)均值滤波方法，算法如下：

对角线滤波方法

输入：n×n窗口内的像素灰度值矩阵A＝[a_ij]_n×n

输出：对角线像素的灰度均值M

Step1：求主对角线像素灰度均值

Step2：求辅对角线像素灰度均值

Step3：求M＝(M₁+M₂)/2；

通过将云图图像分类的目的有两个，其一是将云和陆地和水体分离开，保证云迹风算法输入的是纯粹的“云”。其二是将云分为“高云”、“中云”、“低云”三大类，与云迹风反演的“高空风”、“中空风”和“低空风”对应，获得更有针对性的反演结果。

选取长波红外1(11.5-12.5μm)、长波红外2(10.3-11.3μm)和短波红外(3.5-4.0μm)三个通道云图的灰度特征，建立水体、地表、低云、中云和高云这5类样本的灰度特征统计量，如表3所示。

表3 云图样本的灰度特征统计

确定地表、水体和云系总数目为5，把三维样本灰度特征统计量作为聚类的初始中心，聚类迭代过程中不断调整优化聚类中心和隶属度矩阵，最后得到三维特征空间的聚类结果，云系聚类算法如下：

云系聚类算法：

输入：聚类数目c＝5，三维灰度特征矩阵值V₀＝[v_ij]_3×5，数据量n；算法结束阈值ε

输出：隶属度矩阵U'＝[u_ij]_n×n

运算过程：

Step1：以5类三维样本特征量初始值赋值给V₀，对隶属度矩阵U⁽⁰⁾按如下公式赋予初值，b＝0；

$u_{ij}=\underset{k=1}{\overset{c}{\Sigma }}{d}_{ij}/\underset{k=1}{\overset{c}{\Sigma }}{(d_{ij}/d_{kj})}^{3/\lambda },i=1,2,...,c;j=1,2,\mathrm{....}n;$

Step2：b++；计算新的聚类中心V_i^(b)(i＝1,2,...,c)：

$V_i=\left(\underset{j=1}{\overset{n}{\Sigma }}{u}_{ij}^{2\lambda }{X}_j\right)/\left(\underset{j=1}{\overset{n}{\Sigma }}{u}_{ij}^{2\lambda }\right)$

Step3：对于j＝1～n，计算U^(b)，并使U^(b+1)＝U^(b)；

①计算I_j和

I_j＝{i|1≤i≤c；d_ij＝||X_j-V_i||＝0}

$\overline{I_j}=\{1,2,...,c\}-I_j$

②对于数据X_j，更新其隶属度。

若则对所有令u_ij＝0；

否则：

$u_{ij}=\underset{k=1}{\overset{c}{\Sigma }}{d}_{ij}/\underset{k=1}{\overset{c}{\Sigma }}{(d_{ij}/d_{kj})}^{3/\lambda },i=1,2,...,c;j=1,2,\mathrm{....}n;$

Step4：若||U^(b)-U^(b+1)||＜ε，执行步骤step5；否则返回Step3；

Step5：U’＝U^(b+1)；结束

如图4所示，由于云图中各风矢的计算都是完全独立的，可以选择云图中风矢进行粗粒度并行，并行云迹风反演算法的基本流程如下：

(1)解压卫星下发数据；

(2)取三张相隔时间为Dt的卫星红外通道云图，按时间早晚顺序命名为云图1、云图2、云图3；

(3)分别将图像灰度增强算法分别将三张图片进行灰度增强，以便更准确追踪到目标云块；

(4)将云图分为(Width/ρ)*(Height/ρ)块追踪区块，式中，Width为所述每幅云图像素的宽度，Height为所述每幅云图像素的高度，ρ为风矢密度；

(5)给云上的所有区块分配指定的计算处理单元，对n个并行处理单元使用一定的规则进行任务分配。此处的并行处理单元包括MPI的进程，openMP中的线程以及GPU的各流处理器，而各并行处理单元中的任务代表一个完整的风矢反演任务；

(6)启动n个并行处理事务对分配好的处理任务进行并行计算；

(7)数据汇总。将所有并行处理单元反演的风矢汇总在一起，并将有效的风矢数据写入磁盘保存；

(8)运算结束。

在本实施例中，如图5所示，使用英伟达公司的CUDA架构，thread是GPU中流处理器执行的基本单位，而包含多个thread的block是处理器调度的基本单元，多个block组成一个grid。在同一个block中的多个thread能够共享同一存储空间。本发明以一个grid做为一次GPU并行反演任务，执行一次CUDA的内核函数kernel_function<<<M，N>>>(arguments)，其中M是grid的尺寸，即一个grid中包括多少个block；N是block的尺寸，即每个block中有多少个thread，每个thread分配一个风矢的计算任务。

使用GPU并行技术对云导风反演算法进行加速，将MN个风矢计算任务分配到GPU的thread上执行，本发明的实验为简单起见，将线程格申请为一维变量blocksPerGrid，每个线程格中开辟threadsPerBlock个线程，本发明的实验每个block上线程数threadsPerBlock为256，线程格中block的数量

blocksPerGrid＝(iMax+threadsPerBlock-1)/threadsPerBlock

每个风矢反演任务获取要反演的风矢的索引为

windIndex＝blockDim.x*blockIdx.x+threadIdx.x

最后调用

dev_makeCloud<<<blocksPerGrid,threadsPerBlock>>>(arguments…)对各风矢反演任务进行分配和执行，每个block中的线程在一个流处理器中被顺序执行，并将反演出的到开辟的全局风矢数组对应的位置windIndex中。

精度对比实验

(1)当反演参数为示踪云块大小为22×22像素，追踪区域大小为44×44像素，风矢密度为2个/度时，选取了赤道上连续的10个风矢作为比较对象，并行与串行算法反演的风矢精度对比如表4所示。(其中S表示串行，T-64代表一个线程块中分配64个thread，T-128代表一个线程块block中分配128个线程thread，则线程格中的block数为总任务数量WN除以每个block中的thread数量)。

表4 基于GPU的串行和并行云导风反演算法精度对比1

(2)反演参数为示踪云块大小为32×32像素，追踪区域大小为64×64，风矢密度为4个/度，选取了北纬32度连续的10个风矢作为比较对象，并行与串行算法反演的风矢精度对比如表5所示(其中S表示串行，T-64代表一个线程块中分配64个thread，T-128代表一个线程块block中分配128个线程thread，则线程格中的block数为总任务数量WN除以每个block中的thread数量)。

表5 基于GPU的串行和并行云导风反演算法精度对比2

从上面的图6和图7中可以看出，基于GPU的并行云导风反演，串行反演和并行反演的获得了大部分相同、少部分相似的数值结果，由于GPU的浮点运算还有精度误差，因此取得这种结果较为理想，由此可以证明基于GPU的并行方法在数值上是可信和有效的。

基于GPU的并行反演算法反演出的云导风产品图如8所示。

性能对比实验

实验使用宿主机器为4核8线程CPU，使用的GPU为英伟达公司的NvidiaGeForce GT650M，核心频率为925MHz，显存为1G，显存频率为5400MHz，显存带宽为84.6GB/s，流处理器数量768个，理论计算能力为1.42TFLOPs。分别在不同的反演参数和线程块的情况下对卫星云图进行反演，算法相对加速比和并行效率随反演参数和线程块数量的变化如表6和表7所示。

表6 基于GPU的并行云导风反演算法加速比和并行效率1

表7 基于GPU的并行云导风反演算法加速比和并行效率2

通过对图9和图10进行分析可以看出，在反演参数相同的情况下，基本遵循线程块数目越多，加速比提升越大的规律。当线程块<64时并行加速比基本成2倍增加趋势，当大于64时并行加速比和线程块的划分有关，两实验最高加速比达到了112和103，加速效果十分显著。

如图11所示，最大相关系数法(Maximum Cross-correlation Coefficient，MCC)的原理可以表述为在云图中选中一示踪云块A，然后云块A(N×N像素图像块)在包含A的目标搜索区S(M×M像素图像块)中滑动，计算S中的所有追踪云块和示踪云块A的相关系数，最大相关系数对应的追踪云块B即为目标云块。

最大相关系数R(Δx,Δy)的计算公式如(2-1)所示：

$R(\Delta x,\Delta y)=\frac{{\Sigma }_i^N{\Sigma }_j^N[A(i,j)-\bar{A}][B(\Delta x+i,\Delta y+j)-\bar{B}]}{\sqrt{{\Sigma }_i^N{\Sigma }_j^N{[A(i,j)-\bar{A}]}^2}\sqrt{{\Sigma }_i^N{\Sigma }_j^N{[B(\Delta y+i,\Delta x+j)-\bar{B}]}^2}}---(2-1)$

其中，分别为追踪块偏离示踪云块的行数和列数，(i,j)为示踪云块的索引值，A(i,j)与B(i,j)分别为示踪云块与追踪块的像素灰度值，分别为示踪云块和追踪块的平均灰度。

获得最大相关系数后，目标搜索区内的追踪云块B相对于示踪云块A在x,y方向上移动的距离分别为Δx和Δy，设示踪云块A的中心点的经纬度为(x₁,y₁)，追踪云块B的中心点的经纬度为(x₂,y₂)，建立地球模型(如图12所示)可以计算出目标云块B相对于示踪云块A的距离d，用d除以两幅图像的观测时间间隔t便可得到风矢的速度；推导过程如下：

风矢的速度

$v=\frac{d}{t}---(2-2)$

其中A、B两点间的距离(RE为地球半径)

d＝RE×γ (2-3)

A，B两点和球心形成的球心角可以用空间向量数量积计算

$cos\left(\gamma \right)=\frac{\overrightarrow{OA}·\overrightarrow{OB}}{\left|\overrightarrow{OA}\right|·\left|\overrightarrow{OB}\right|}---(2-4)$

设向量为(a₁,b₁,c₁)，向量为(a₂,b₂,c₂)，则

a₁＝RE×cos(x₁)×cos(y₁)>

b₁＝RE×sin(x₁)×cos(y₁)>

c₁＝RE×sin(y₁)>

同理

a₂＝RE×cos(x₂)×cos(y₂)>

b₂＝RE×sin(x₂)×cos(y₂)>

c₂＝RE×sin(y₂)>

将公式(2-5)～(2-10)入(2-4)可得球心角γ的大小，然后将求出的γ代入(2-3)可求出A、B两点间的距离d，最后将d代入(2-2)便可求出风矢的速度υ。

风矢的方向需借助球面三角形余弦定理来计算推导。如图13所示，平面AED与地球球面ABC相切于点A，A、B分别为风矢的起点和终点。

对三角形ODE运用三角形余弦定理有

$cos\left(a\right)=\frac{{\mathrm{OE}}^2+{\mathrm{OD}}^2+{\mathrm{DE}}^2}{2OE·OD}---(2-10)$

由于三角形OAE和三角形OAD都是直角三角形，因此有

OE²＝OA²+AE²>

OD²＝OA²+AD²>

将(2-12)(2-13)代入(2-11)并运用三角形余弦定理有

cos(α)＝cos(β)×cos(γ)+sin(β)×sin(γ)×cos(θ) (2-13)

其中

α＝x₂-x₁>

β＝y₂-y₁>

将(2-15)(2-16)和上面求出的γ代入(2-14)，利用反正弦函数即可求出θ的大小。如果A点在赤道以北地区，风矢方向角为γ，若x₂＜x₁，则风矢方向角为2π-θ；若A点在赤道以南地区，则风矢方向角为π-θ，若x₂＜x₁，则风矢方向角为π+θ。

风矢的高度可以通过风矢所在位置的红外一和水汽通道的亮温值来指定。

风矢计算完毕后需要进行风速和方向校验，以便剔除初步计算结果中风速和方向偏差较大的风矢。定义三张连续时间观察的云图为C₁、C₂、C₃。设利用C₁和C₂算出的风矢CV₁的风速和风向分别为υ₁和θ₁，利用C₂和C₃算出风矢CV₂的风速和风向分别为υ₂和θ₂，那么两风矢的风向差Δθ＝|θ₁-θ₂|，风速的相对差Δυ＝|υ₁-υ₂|，当风矢的风向差Δθ或风速差Δυ大于给定的阈值时，则剔除此风矢。

在本实施例中，最大相关系数法使用的场景是找出要反演的风矢所在像素区块中心32*32像素区域相似度最大的32*32区域，相似度与相关系数的关系是单调递增的，所以使用此算法可以找出目标区块的位置，从而计算出风矢的相关数值，示意图如图11所示。

在该算法中，计算两个区块之间的相关系数时间较少，相关系数数组的计算耗时较长，假如区块的大小为N*N，区块所在的追踪区域为M*M，如每个相关系数都计算，则需要计算(M-N+1)*(M-N+1)次相关系数，时间为计算一次相关系数时间的一千多倍，因此，本发明将对这个计算量进行优化。

在本实施例中，使用“变步长迭代计算”，变步长搜索法是指在追踪区域S中，示踪云块的移动步长由原来的始终为1的策略改为变步长移动的策略，搜索的初始步长采用较大值，之后逐步缩小步长同时缩小搜索范围直到步长为1，每种步长都找出当前步长情况下的最佳目标匹配云块，步长为1时便找到了近似全局最佳的匹配云块B。

在本实施例中，假定初始步长为8，则有以下推导过程：

第一遍移动步长为8，计算出(M/8)*(M/8)个相关系数值，找出这些相关系数值中最大的值，从而找到相似度最高的N*N区块A1，方法如图14所示。

找到A1之后，为避免误差，将A1块边长增加step个像素点，然后将步长减少为8/2即4个像素点，即step为4，在以A1为中心的N+step边长的区域中对相关系数最大的N*N区域进行搜索，分别计算各N*N区块和B区块的相关系数，计算出的系数矩阵大小为(N+step)/4*(N+step)/4，即(N+step)²/16，这就代表迭代搜索这些步数就能找到step为4情况下的最大相关系数区块，搜索算法如图15所示。

按照此算法，step逐步减半的方式递减，最后当step步长为1的时候搜索到的最大相关系数区块Block即为与中心区块B相似度最高的区块，从而根据两区块的位置信息计算出风矢的速度和角度。

示踪云追踪的原方法找到最佳相似区块需迭代搜索的次数我们计算可知，

StepCount1＝(M-N+1)²

算法优化后找到最佳相似区块需迭代搜索的次数

StepCount2＝(M/step)²+step/2+step/4+…+1

在步长修改为M为64、N为32、step为8的情况下，原方法需要1089步搜索，优化后的方法需要71步搜索，大大简化搜索需要的时间，理论上可以加速15倍以上，并且后期的精度优化可以针对变步长后再次搜索的范围做与M大小相关的调整。

尽管本发明的实施方案已公开如上，但其并不仅仅限于说明书和实施方式中所列运用，它完全可以被适用于各种适合本发明的领域，对于熟悉本领域的人员而言，可容易地实现另外的修改，因此在不背离权利要求及等同范围所限定的一般概念下，本发明并不限于特定的细节和这里示出与描述的图例。