并行计算在基于遥感影像中的应用

转自 http://blog.chinabyte.com/a/285247.html  

在众多自然灾害中地震造成的人员伤亡数量和经济损失最大，近两年来中国汶川。从少量的建筑物的结构破坏，到城市中心或者大型国家社会经济基础建筑的完全摧毁。限于当前技术水平还不能实现有效的地震灾害预报预警，只能通过开展灾后救援，安置灾民和灾后重建等手段，减少灾害带来的生命和财产损失。房屋倒塌作为城镇区域受损程度的重要标志，其数量和分布能够直接反映灾区生命财产损失的情况。因此，快速获取受灾地区房屋倒塌信息，使应急救援在地震灾害后72小时内投送到关键位置是最为紧迫的需要。
    
    　　1.项目应用背景
    
    　　1m以内空间分辨率光学影像能够准确反映倒塌房屋尺度的地表细节信息，是获取灾后倒塌房屋信息的主要数据源，在减灾救灾中发挥重要作用。人工目视解译和变化检测是当前通过高分影像提取倒塌房屋的主要方法。但是，该方法的实现效率很难满足灾害应急阶段的快速获取信息需求，同时，由于操作者个体差异等因素，容易造成多人解译结果存在不一致性，给后续信息分析和决策带来困难。近年来，通过结合图像处理和模式识别算法的变化检测方法，逐渐被应用实现灾后倒塌房屋信息的高精度快速提取。然而，随着这些遥感影像数据的空间分辨率的提高，其数据量也成倍增长，一般一个时相的能反映城市建筑损毁细节的高分辨率遥感数据都要在G级以上，甚至达到T级。这对计算机自动处理这些数据的运算能力提出更高要求。针对灾后72小时应急阶段时间就是生命的快速信息提取的需要下，进一步提高计算机处理高分辨率遥感影像的效率显得至关重要。
    
    　　本项目利用玉树地震期间获取的该地区的高分辨率航空遥感影像，基于这些影像研发了“遥感影像房屋倒塌快速识别系统”，这将为以后类似救灾决策提供有利支持。
    
    　　2.应用的业务瓶颈和规模
    
    　　高分辨率遥感影像是一般是指像素的空间分辨率在 10 m以内的遥感影像。与传统的低空间分辨率的卫星影像相比，高分辨率卫星影像具有单幅影像的数据量显著增加、成像光谱波段变窄和地物的几何结构和纹理信息更加明显等特点。高分辨率卫星遥感影像的出现使得在较小的空间尺度上观察地表的细节变化、进行大比例尺遥感制图以及监测人为活动对环境的影响成为可能，具有广阔的应用前景，是目前针对地震灾害空间信息获取的主要数据。由于地震灾害损毁主要是城市建筑，救灾的目标主要是救助居民人身安全，因此对遥感数据的分辨率要求很高，这造成了单幅影像的数据量成倍增加，计算量更是呈几何基数增长，使计算机处理图像数据的负荷急剧加大，使“遥感影像房屋倒塌快速识别系统”需要大量的计算用时，甚至一次处理达到几十小时。这与从数据中快速提取房屋倒塌信息的目标相违背，影响了该系统在分秒必争的地震灾害抢险决策中发挥的作用，很容易导致该软件研发项目不被用户接受而失败。
    
    　　3.并行化方法
    
    　　根据遥感图像数据是一种栅格数据，其处理算法是以行列为单位反复进行反复运算为主的特点，系统研发选用了NVIDIA的CUDA架构来实现GPU-CPU耦合并行计来提高系统的预算效率。在遥感图像处理的GPU加速算法设计中，针对这种架构分割像素或图像块以映射到并行线程进行处理，再经由多级存储器将被处理的像素或图像块重新组合成结果图像。
    
    　　通过CUDA架构将GPU看作可以并行执行非常多个线程的计算设备(compute device)。它作为主CPU的协处理器或者主机(host)来运作，即：在主机上运行的应用程序中数据并行的、计算密集的部分卸载到此设备上。CUDA将计算任务映射为大量的可以并行执行的线程，有硬件动态调度和执行这些线程。更准确地说，对不同数据上独立执行的应用程序部分可以独立放到在此设备上作为许多不同线程执行的函数。要达到这种效果，可以将这样一个函数在设备的指令集合中编译，并将得到的程序(叫做内核, kernel)加载到设备上。Kernel以线程网格(grid)的形式组织，每个线程网格由若干线程块(block)组成，而每个线程块又有若干个(thread)组成。基于遥感影像栅格化的特点，其处理算法非常适合建立并行计算。针对数据栅格空间其中一个栅格的计算可以映射到一个线程，多个线程可以视作对栅格数据的并行计算。
    
    　　4、实施和效果
    
    　　采用CPU(host thread，主线程)在开始计算前读入计算参数和图像栅格数据，映射到各相应的线程并复制到显存，供GPU计算时调用，其余并行操作均在GPU上完成。
    
    　　除了像素的灰度值，高分辨率遥感影像的纹理是描述图像的一个重要特征，基于小波变换的纹理分割方法是一种理想的图像分割算法。系统的地物类型识别运算中采用了Daubechies1小波，设定小波包分解的层数为m，分别进行m层二维小波包分解，对于每个样本图像子块，均得到4m个分解子频带图像。为了简便，决定进行单尺度的小波变换，m为1，共得到4个分解子频带的图像。每个图像子块的4个子频带的绝对值平均能量值作为一个特征向量。然后对待分类图像进行小波变换，提取各像素特征向量。使用K-Mean算法进行遥感图像分类，计算各点与分类初始中心的距离，确定类别并调整聚类中心，进行迭代。迭代结束后根据分类结果对像素设定颜色，最后计算倒塌房屋和为倒塌房屋比例。
    
    　　在GPU计算中的图像每个像素的特征向量，是对线程块中的每个线程计算像素特征向量，采用线程 将图像像素 和小波变换相关数值，由全局内存读到共享内存中，在同一线程块的各线程均完成数据加载后(通过线程块内同步确保完成)，线程 计算 的像素小波变换后4个子频带的绝对值平均能量值。如此同步控制并行完成每一个像素小波变换的计算，最后将更新后的共享内存中的像素特征向量信息写回各自对应的全局内存。
    
    　　系统实验采用玉树地震期间获取的该地区的高分辨率航空遥感影像，如图1。经系统快速处理后得到如图2结果。
    
    　　
    
    　　图1.原始图像(图像来源于中科院对地观测中心)
    
    　　
    
    　　图2.结果图像
    
    　　在处理大尺寸图像时GPU的总计算效率能够达到CPU的19.2倍，GPU核心计算效率能达到CPU的38.3倍。采用CUDA编程接口，为GPU重新设计和编写原有图像处理算法以获得数十倍的效率提升是值得的。同时考虑到GPU计算相对于CPU计算的成本，以相对低廉的GPU能获取对高分辨率遥感图像快速处理的高性能计算能力。