稀疏矩阵的雷达 dBZ 数据格式方案

稀疏矩阵的雷达 dBZ 数据格式方案雷达 dBZ 数据是我们制作降雨短临预报的重要数据资源。由于该类数据属于栅格数据且时空分辨率较高若以常规气象的栅格点存储格式进行存储则存储压力很大。目前我们主要是用 PNG 图片的格式对 dBZ 矩阵数据进行存储由于 PNG 是一种无损压缩的图片格式且压缩比很高因此可以在不损失数据精度的情况下极大地减轻存储压力。但是使用 PNG 格式存储有一个缺点就是使用时需要手动做时间轴堆叠很不方便。最近我们尝试利用稀疏矩阵的原理设计了一套既节约存储空间又能实现高效随机读取的矩阵数据格式方案smdbz。dBZ 数据的等级特征按照国内中央气象台的制图标准dBZ 的值可以按照 5 的间隔从 5 到 75 分为 14 个等级。再加上 0 值本身dBZ 一共可以分为15个等级。由此我们可以用 0 到 14 这 15 个整数来表示 dBZ 的等级。根据计算机的存储原理一个字节Byte包含 8 个二进制位Bit以无符号整型uint8为例1个字节可以存储一个范围从 0 到 255 的整数2^8256个值)。但是由于 dBZ 的等级其实最大也只到 14所以用一个字节存储会显得比较浪费事实上它仅需要半个字节即4个位2^416即可存下所有等级数字。但由于字节是计算机最小可寻址的内存单位传统编程语言中通常不提供仅包含“半个字节”的数据类型因此我们需要使用位运算来实现使用半个字节的存储效果具体方法后续讲解。dBZ 拼图矩阵的稀疏特征以windy官网我们来看这张全球的雷达拼图可以看到对于整张图的面积来说存在雷达回波的区域面积仅占很小一部分。根据我们统计数据来说绝大部分时间雷达 dBZ 的面积占整张图面积的占比在 5% 以下也就是说超过 95% 的区域的数据其实都是 0这也就符合了稀疏矩阵的特征。而现在在计算机领域有很多成熟的方案来对稀疏矩阵进行高效的压缩存储因此我们可以尝试使用稀疏矩阵的压缩原理来构建 dBZ 矩阵。稀疏矩阵压缩原理稀疏矩阵压缩的核心思路是仅存非零数据。人们基于这个思路发展出很多不同的稀疏矩阵的存储方案。比较常见的包括压缩稀疏行CSR、压缩稀疏列CSC、坐标COO。在这里我们将使用 CSR 方案来实现 smdbz 的稀疏矩阵构建其他方案如果感兴趣可以自行查找。CSR 通过三个数组来存储稀疏矩阵。val(下文出现的value与val为同一变量数组存储所有非零元素的值colind 数组存储这些非零元素所在的列索引rowptr 数组存储每一行的起始索引。如下图所示图片来源Elafrou, Athena Goumas, Georgios Koziris, Nectarios. (2015). A lightweight optimization selection method for Sparse Matrix-Vector Multiplication.图中所示的 val 是将矩阵中所有零值剔除以后把剩余的非零数值进行平铺排列的一维数组。rowptr 存储的是新行位置点以上图矩阵为例该矩阵中每一行的首个非零数据分别是 [7.5, 6.8, 2.4, 9.7, 5.8, 6.6]。这几个数对应在 val 数组中的位置索引分别是 [0, 4, 7, 10, 12, 14]而 rowptr 就是这个索引坐标列表再加一个 val 的总长度 16即 [0, 4, 7, 10, 12, 14, 16]。需要注意的是 rowptr 数组的长度与整个矩阵的行数有关即 rowptr 的长度永远比矩阵行数多1colind 则是 val 中每一个元素在其所在行的列索引colind 列表的长度与 val 的长度保持一致。当我们拿到一个稀疏矩阵只要提取出 val、colind和rowptr 这三个参数就可以构建一个 CSR 结构的压缩稀疏矩阵。基于 scipy 的 CSR 实现详细文章关注微信公众号《关于气象遥感》稀疏矩阵实现雷达存储详情