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

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

背景技术：
在天气和灾害预报的研究中，需要大量的风场资料，全球大气风场的观测起到了非常积极的作用。近年来，随着陆地风场观测站网的不断完善，研究人员获得陆地风场资料变的越来越容易。然而，相对于陆地风场资料来说，海洋、高原、沙漠等地区的观测网就显得非常不足，获取这些地区的风场资料较困难。通过卫星云图资料反演大气风场，无疑对获取全球大气风场资料起到了非常有益的补充作用。同时，一些陆地或水体(称为地表)由于上空无云或少云而直接“裸露”在卫星云图中。研究人员主要是对云图中的云进行研究，例如云分类、云识别、云导风等，因此实现云区域和地表区域的图像分割是必不可少的。针对卫星云图图像分割的研究工作早已展开，目前使用的方法普遍具有以下缺陷：(1)过度依赖于历年的云图图像资料，需要花费人力物力收集和保存；(2)使用基于阈值的图像分割来实现，属于硬性分割，准确率难以保证；(3)云图是一个十分复杂的图像，不仅包含云和地表，还在云和地表之间存在一些处于“过渡阶段”的云(例如生长和消散的云)；所以，在传统方法中，容易将其强制划分到某一类中，但这是不合理的，对云图图像合理分割也是研究的方向之一。一般来说，卫星云图时间间隔越短、用于反演的示踪云越多、风矢的计算密度越高，云导风反演出的风场资料质量也会越高，但这都对云导风反演算法的效率提出了挑战。虽然有学者提出了云导风反演算法的简易算法，但也难以满足日趋要求实时性的风场计算的效率要求。因此，除研究较低计算代价的反演算法，充分利用现代计算机的体系结构，研究并行云导风反演算法，也是提升反演算法实际运行效率的有效途径之一。GPU(GraphicProcessingUnit)，图形处理器是显卡的“心脏”，相当于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＝(M1+M2)/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，保存的格式为：风矢编号、风矢经度、风矢纬度、风矢速度以及风矢方向。下面再结...