雨量记录图纸的数字化方法、装置与流程

本发明涉及一种雨量记录图纸的数字化方法、装置，利用计算机自动识别雨量记录图纸所记录的雨量数据，属于水文气象、数字化图像处理、图像识别以及数据库相交叉的技术领域。

背景技术：

在计算机没有普及之前，在水利、气象、生物医学、石油地质等领域，都是通过特定仪器或者人工将观测的有效数据信息绘制在纸质工程图纸上。这些纸质水文资料上的数据均是比较珍贵的历史资料，记录流量、水位、降水量、蒸发量等，其中很多的数据会被经常用来进一步分析与处理，对这些有效信息的利用与共享将为我国水利工程建设提供科学的理论与依据，将会大大促进我国水利事业的发展。

然而这类水文纸质资料并不是完全保存完整，随着时间，会受到存放条件、温度、环境等影响，导致如纸质磨坏破损、时间过久老化、字迹变淡模糊、纸质资料字迹渗透等污染，这些都会对历史数据的安全以及价值带来影响。因此对这类纸质资料进行数字化，将纸质水文资料中的重要的数据转移到计算机中进行分类存储和处理是很有必要的。

由于水文纸质资料种类繁多以及数据多样化，对这些资料进行保存到计算机的过程中，一般是以人工为主，选择代表性的节点及数据，指定插值或曲线拟合的方法，计算出某一时间的水位信息等数据。但是曲线形态各异，拟合或插值并不能通用。另一种方式是通过人工录入数据并建立大型数据库，由于资料种类繁多数据量大，工作量巨大且易出错，数据需要经过严格的审查与校核才能入库，并需要提供各种不同的检索方式，而且并不保证能够在最短时间准确找到需要的数据。

在种类繁多的水文纸质资料中，雨量记录图纸又具有其特殊性。每张雨量图都会记录1天以上的降雨量，其中没有降雨的天对应的记录曲线是一条水平直线。而对应有降雨的曲线记录当天8时到第二天8时的降雨过程线，它从某一较低的 刻度开始，如果此24小时内降雨超过了量杯(满量杯为10mm)，就将量杯的雨水倒出继续计量，因此曲线会形成若干个峰，峰值在9-10mm之间。最后一次量杯里的降雨一直上涨，但只要量杯未满，则对应曲线上的量杯未满记录线。计算该日的降雨量，需要求得各峰值以及量杯未满记录值的和，最终得到降雨量。对于雨量记录纸，有以下几个特点：由于需要记录24小时的降雨过程，因此记录纸幅宽较大，相机拍摄得到的图像像素多，容易产生光照不均匀现象，表现在同一根曲线不同部分颜色不一致；由于降雨过程复杂，在连续下暴雨时，峰比较密集，表现在雨量线陡降陡升；并且由于一张纸上记录了多天降雨信息，记录线之间既可能互相平行，也可能互相交错，传统图像处理方法不易跟踪分析每根线的过程变化。

技术实现要素：

本发明所要解决的技术问题在于克服现有技术不足，提供一种雨量记录图纸的数字化方法，能够从雨量记录图纸图像中快速而准确地自动提取出所记录的降雨量数据，便于后续的数据查询分析。

本发明具体采用以下技术方案解决上述技术问题：

雨量记录图纸的数字化方法，利用计算机自动识别雨量记录图纸所记录的雨量数据，包括以下步骤：

步骤A、从雨量记录图纸的数字图像中分别提取出二值化的网格图像、雨量线图像；

步骤B、对网格图像、雨量线图像进行倾斜校正；

步骤C、利用行列投影的方法从倾斜校正后的网格图像中确定网格区域的边界，并将倾斜校正后的网格图像、雨量线图像中的边界外区域设置为二值化背景；

步骤D、将雨量线图像中的直线设置为二值化背景；

步骤E、将雨量线图像对应到网格图像，求出雨量线上各个峰值点以及最右侧峰值点右侧延长线的最大值点，最后对以上各点所对应的雨量值求和，得到降雨量。

优选地，使用以下方法从雨量记录图纸的数字图像中提取出二值化的雨量线图像：

步骤1、对雨量记录图纸的数字图像进行高斯金字塔重采样处理，得到所述数字 图像的低分辨率图像；

步骤2、对所述低分辨率图像进行初步二值化处理，得到初步二值化的雨量线图像；

步骤3、将雨量记录图纸的数字图像、低分辨率图像、初步二值化的雨量线图像分别按照同样的划分方式划分为从左到右的N个子图像，N为大于1的整数；

步骤4、找出初步二值化的雨量线图像的第i个子图像中每个目标像素点在低分辨率图像的第i个子图像中的对应像素点，i＝1,2,…,N，并在CIE Lab颜色空间中计算出这些像素点与其均值点(L_avg，a_avg，b_avg)之间的最大空间距离dTh；

步骤5、遍历雨量记录图纸的数字图像的第i个子图像中的所有像素点，将满足以下条件的像素点设置为二值化目标，其余像素点设置为二值化背景：在CIE Lab颜色空间中与均值点(L_avg,a_avg,b_avg)之间的距离小于dTh，且在HSI颜色空间中的H分量大于hTh，hTh为大于或等于180的预设H分量阈值；从而得到最终的二值化的雨量线图像。

进一步地，所述对网格图像、雨量线图像进行倾斜校正，具体如下：从初步二值化的雨量线图像中找出最长的直线，然后利用霍夫变换方法检测出该直线的倾斜角度，最后按照检测出的倾斜角度对网格图像、雨量线图像进行旋转校正。

优选地，使用改进的Niblack算法对低分辨率图像进行初步二值化处理，具体如下：首先将所述低分辨率图像灰度化；然后利用以下公式确定各像素点的二值化阈值：

$T(x,y)=m(x,y)\·\[1+k\·(\frac{s(x,y)}{R}-1)\]$

其中T(x,y)为点(x,y)的二值化阈值，m(x,y)是点(x,y)周围预设大小模板内点的灰度均值；s(x,y)表示模板内像素点灰度值的标准差；参数k的取值范围为[0.4，0.6]，R为128；

如点(x,y)处的灰度值小于T(x,y)，则设置其为背景，否则，设置为目标。

根据相同的发明思路还可以得到以下技术方案：

雨量记录图纸的数字化装置，包括：

网格图像提取单元，用于从雨量记录图纸的数字图像中提取出二值化的网格图像；雨量线图像提取单元，用于从雨量记录图纸的数字图像中提取出二值化的雨量线图像；

倾斜校正单元，用于对网格图像提取单元和雨量线图像提取单元所输出的网格图像、雨量线图像进行倾斜校正；

网格裁剪单元，用于利用行列投影的方法从倾斜校正后的网格图像中确定网格区域的边界，并将倾斜校正后的网格图像、雨量线图像中的边界外区域设置为二值化背景；

直线去除单元，用于将雨量线图像中的直线设置为二值化背景；

降雨量计算单元，用于将雨量线图像对应到网格图像，求出雨量线上各个峰值点以及最右侧峰值点右侧延长线的最大值点，最后对以上各点所对应的雨量值求和，得到降雨量。

优选地，所述雨量线图像提取单元具体包括：

高斯金字塔重采样模块，用于对雨量记录图纸的数字图像进行高斯金字塔重采样处理，得到所述数字图像的低分辨率图像；

初步二值化处理模块，用于对所述低分辨率图像进行初步二值化处理，得到初步二值化的雨量线图像；

图像划分模块，用于将雨量记录图纸的数字图像、低分辨率图像、初步二值化的雨量线图像分别按照同样的划分方式划分为从左到右的N个子图像，N为大于1的整数；

阈值获取模块，用于找出初步二值化的雨量线图像的第i个子图像中每个目标像素点在低分辨率图像的第i个子图像中的对应像素点，i＝1,2,…,N，并在CIE Lab颜色空间中计算出这些像素点与其均值点(L_avg,a_avg,b_avg)之间的最大空间距离dTh；

雨量线图像输出模块，用于遍历雨量记录图纸的数字图像的第i个子图像中的所有像素点，将满足以下条件的像素点设置为二值化目标，其余像素点设置为二值化背景：在CIE Lab颜色空间中与均值点(L_avg,a_avg,b_avg)之间的距离小于dTh，且在HSI颜色空间中的H分量大于hTh，hTh为大于或等于180的预设H分量阈值；从而得到最终的二值化的雨量线图像。

进一步地，所述倾斜校正单元使用以下方法对网格图像提取单元和雨量线图像提取单元所输出的网格图像、雨量线图像进行倾斜校正：从初步二值化的雨量线图像中利用霍夫变换方法检测出最长直线的倾斜角度，最后按照检测出的倾斜角度对网格图像、雨量线图像进行旋转校正。

优选地，所述初步二值化处理模块使用改进的Niblack算法对低分辨率图像进行初步二值化处理，具体如下：首先将所述低分辨率图像灰度化；然后利用以下公式确定各像素点的二值化阈值：

$T(x,y)=m(x,y)\·\[1+k\·(\frac{s(x,y)}{R}-1)\]$

其中T(x,y)为点(x,y)的二值化阈值，m(x,y)是点(x,y)周围预设大小模板内点的灰度均值；s(x,y)表示模板内像素点灰度值的标准差；参数k的取值范围为[0.4，0.6]，R为128；

如点(x,y)处的灰度值小于T(x,y)，则设置其为背景，否则，设置为目标。

进一步地，该装置还包括用于获取雨量记录图纸的数字图像的图像获取单元，以及用于存储降雨量计算单元所输出降雨量数据的数据存储单元。

相比现有技术，本发明具有以下有益效果：

本发明能够从雨量记录图纸图像中快速而准确地自动提取出所记录的降雨量数据，便于后续的数据查询分析，大幅提高了雨量数据数字化的效率，降低了人力和时间成本；

本发明在进行雨量线图像提取时，先利用高斯金字塔模型获取低分辨率图像，然后对低分辨率图像进行初步的二值化处理，最后结合雨量记录图纸图像CIE Lab颜色空间和HSI颜色空间特点，对初步的二值化处理结果进行进一步优化，有效降低了计算量，提高了提取速度，且获得的结果更加准确；

本发明以初步二值化的雨量线图像中的最长直线作为参考来进行霍夫变换倾斜校正，算法更简单，准确度更好。

附图说明

图1为实际雨量记录图纸；

图2为具体实施方式中本发明数字化方法的流程示意图；

图3为初步二值化的雨量线图像；

图4为二次优化后的雨量线图像；

图5为提取出的二值化的网格图像；

图6为倾斜校正后的雨量线图像；

图7为倾斜校正后的网格图像；

图8为裁剪后的雨量线图像；

图9为裁剪后的网格图像；

图10为去除直线后的雨量线图像。

具体实施方式

下面结合附图对本发明的技术方案进行详细说明：

图1显示了一份实际雨量记录图纸，如图所示，每张雨量图都会记录1天以上的降雨量，其中没有降雨的天对应的记录曲线是一条水平直线。而对应有降雨的曲线记录当天8时到第二天8时的降雨过程线，它从某一较低的刻度开始，如果此24小时内降雨超过了量杯，就将量杯的水倒出继续计量，因此曲线会形成若干个峰，峰值在9-10mm之间。最后一次量杯里的降雨一直上涨，但只要量杯未满，则对应曲线上的量杯未满记录线。计算该日的降雨量，需要求得各峰值以及量杯未满记录值的和，最终得到降雨量。

本发明首先利用数字图像处理技术从雨量记录图纸的数字图像中提取网格图像和雨量线图像，然后对所提取出的图像进行校正、去噪等处理，然后根据雨量线峰值及延伸线极值计算出雨量值。图2显示了本发明方法的基本流程，具体包括以下步骤：

步骤A、从雨量记录图纸的数字图像中分别提取出二值化的网格图像、雨量线图像。

1、雨量线提取：

由于雨量记录图纸比较大，用数码相机拍摄的雨量图分辨率都很高，通常达到4000*3000的像素级别。如直接对原图(通常为标准的RGB位图)进行处理，显然计算量大，提取速度慢，严重影响雨量记录图纸数字化效率。为解决这一问 题，本发明在提取雨量线时并不直接利用雨量记录原图I₀进行处理，而是通过构造图像高斯金字塔模型来获取较低分辨率图像，即对雨量记录图纸的数字图像进行高斯金字塔重采样处理，然后对低分辨率图像进行处理。具体所采用的金字塔级数可根据原图分辨率确定，例如，对于4000*3000尺度的原图，可采用三级高斯金字塔模型，从而得到低分辨率图像I₂。图像高斯金字塔为现有成熟技术，为节省篇幅，此处不再赘述，具体细节可参考文献(Rafale C Gonzalez,Richard E Woods著.数字图像处理[M].阮秋琦,阮宇智等译.北京:电子工业出版社,2004.)。

对低分辨率图像I₂进行初步的二值化处理，得到初步二值化的雨量线图像。可采用现有的Niblack算法、灰度期望值法、Otsu方法等进行初步二值化分割，本发明优选采用改进的Niblack算法(J.Sauvola,M.Pietikainen.Adaptive document image binarization.Pattern Recognition,33(2000),pp.225–236)，具体如下：

先将低分辨率图像I₂进行灰度化，然后利用以下改进的Niblack公式进行二值化处理：

$T(x,y)=m(x,y)\·\[1+k\·(\frac{s(x,y)}{R}-1)\]$

其中T(x,y)为点(x,y)的阈值，m(x,y)是点(x,y)周围预设大小模板(通常在9*9～15*15之间选取，优选11*11)内点的灰度均值；s(x,y)表示模板内像素点灰度值的标准差；参数k的取值范围为[0.4,0.6]，优选为0.5；由于我们将彩色图像转换成8位灰度图像，所以R可设为128。当点(x,y)处的灰度值小于T(x,y)，则其二值图像的值设为0，否则设为255，从而得到二值图像I₂^B，由此可提取出雨量曲线的大部分。图3显示了采用改进Niblack算法所得到初步二值化雨量线图像的一个实例。

初步提取的雨量曲线由于受到背景干扰很大，所以得到的二值图像含有很多噪点和多余信息，同时如果雨量曲线较细，那么在模板中的白色像素点比例较大，使得阈值偏大，由此得到的曲线断点较多，因此本发明对初步二值图像进行进一步优化，以此来提取更完整的雨量线图像。本发明所采用的优化方法具体如下：

①由于雨量记录图纸长度较大，在拍摄时容易产生光照不均，为了减少光照对雨量曲线的影响，首先将图像I₀、I₂、I₂^B按照同样的划分方式划分为从左到右 的多个子图像，本实施例中将图像均分为左、中、右三个部分；

②以左三分之一为例，将图像I₂在CIE Lab颜色空间中的三个分量归一化；

③将二值图像映射到CIE Lab彩色图像I₂中，统计求得与二值图像中的黑色点(即初步获取的雨量线上的点)所对应的点的L，a，b分量的均值L_avg,a_avg,b_avg；

④计算CIE Lab彩色图像I₂中各点到(L_avg,a_avg,b_avg)的欧几里得距离dis，确定能提取出曲线的阈值dTh，本发明采用初步获取的雨量线上的点到(L_avg,a_avg,b_avg)的最大距离作为dTh；

⑤经大量分析观察发现，原始图像HSI颜色模型中的雨量曲线H分量值通常都大于180，因此本发明进一步以该值或稍大一些的值作为HSI颜色空间的H分量阈值hTh，本发明优选hTh的值为200；

⑥遍历原始图像I₀中的所有点，将满足以下条件的像素点设置为二值化目标，其余像素点设置为二值化背景：在CIE Lab颜色空间中与均值点(L_avg,a_avg,b_avg)之间的欧几里得距离小于dTh，且在HSI颜色空间中的H分量大于hTh；亦即如果距离dis<dTh并且H>hTh，则将该像素点的值设为0，否则设为255；

⑦按照上述方法重复实现中三分之一，右三分之一部分，得到二次优化的雨量曲线。图3的初步雨量线图像采用上述方法优化后如图4所示。

2、网格图像提取

由于网格线的提取需要提取所有的横向刻度线以及纵向刻度线，因此都是在原图I₀上进行处理，本发明优选采用中国发明专利申请(公开号为103971367A，申请日为2014/4/28)中所公开的网格图像提取方法，具体步骤如下：

(1)利用Sobel梯度算子对灰度化的图像进行水平、垂直方向的梯度化，获得两个方向上的梯度图以便获得目标线条的边缘信息；

(2)在水平和垂直方向的梯度图中，分别进行处理(此处以水平方向为例，垂直方向的处理方式相同)：

①为了减少光照不均的影响，对水平梯度图分成左中右三块，在对应小块 内，结合曲线提取的曲线所在位置，剔除掉对应曲线的梯度值的影响，计算块内属于网格线条的梯度均值Th；

②利用梯度均值Th对水平梯度化图像进行二值化，在水平梯度梯度化方向上，另加入像素点的水平梯度比垂直梯度大的条件，使得水平梯度化更加精确，根据此二值化结果初步获得属于橙黄色线条的位置信息；

③根据②中的获得黄色线条位置信息，获得水平方向上对应位置上的像素点在CIE Lab空间的a和b分量大小，并统计这些位置点的a、b均值

④利用③中得到的如果像素点的a、b分量满足则为目标，否则为背景，对原图像进行分割。

⑤获得水平方向二值化图；同样方法获得垂直方向的二值化图。

(3)对水平和垂直方向的二值化图进行或运算合并，得到整个图像二值化图；

(4)对于(4)中获得的二值化图分别进行水平方向1×5的腐蚀并保存结果，和垂直方向上5×1的腐蚀保存结果，去除掉更多的离散点；

(5)对得到水平和垂直两个方向上的二值化图进行合并，得到最终的网格图像。图5显示了提取出的网格图像。

步骤B、对网格图像、雨量线图像进行倾斜校正。

在获取雨量记录图纸的图像过程中，由于其幅面较大，往往会导致图像出现不同程度的倾斜，因此有必要对图像进行倾斜校正。现有的倾斜校正通常是采用霍夫变换的方法，求得图像偏转的角度，然后旋转校正。考虑到雨量线图像中会存在粗长的直线(无降雨的雨量曲线)，便于识别处理，因此本发明从初步二值化的雨量线图像中找出最长的直线，然后利用霍夫(Hough)变换方法检测出该直线的倾斜角度，最后按照检测出的倾斜角度对网格图像、雨量线图像进行旋转校正；其具体过程如下：

(1)在较小尺度的二值图像I₂^B上遍历黑色像素点(i,j)，利用|i*cosα+j*sinα|求得(i,j)到原点的距离，并将其保存到累加器数组中；

(2)找到累加器数组的最大值，即二值图像I₂^B中最长的直线，求得对应偏转角度α；

(3)根据角度α对原图顺时针旋转，转换公式：i1＝i*cosα+j*sinα,j1＝i*cosα-j*sinα,其中(i,j)和(i1,j1)分别为原图和旋转后图像的像素坐标，得到角度倾斜校正 后的雨量线图像和网格图像，分别如图6、图7所示。

步骤C、利用行列投影的方法从倾斜校正后的网格图像中确定网格区域的边界，并将倾斜校正后的网格图像、雨量线图像中的边界外区域设置为二值化背景。

旋转后的雨量线图像和网格图像包含了很多记录图纸以外部分的信息，因此需要将网格部分裁剪出来，避免其他信息的干扰。其具体方法如下：

(1)对网格图像行、列投影，统计各行(列)的黑色像素点个数；

(2)根据行(列)投影的统计设定阈值，此处行阈值取200，列阈值取400；

(3)根据设定的行列阈值找到网格线上下左右边界；

(4)将网格边界外的部分全部置为白色，得到裁剪后的图像；裁剪后的雨量线图像、网格图像分别如图8、图9所示。

步骤D、将雨量线图像中的直线设置为二值化背景。

通常雨量记录图中都包含若干天的降雨信息，而且其中前几天是没有降雨的，其雨量记录就是一根直线，因此可以将雨量线图像中的直线抹去，从而避免对后续有降雨的雨量记录识别产生干扰。具体可采用如下方法：

(1)对雨量线图像在网格的范围内从网格左端开始的450列做水平投影，统计黑色像素点的个数；

(2)设定阈值Th,当某行的水平投影统计值大于Th则认定直线在此行上，并且第一次出现大于Th的行号保留为有雨天的记录；

(3)在统计值大于Th的行号上下20行内将所有像素点的灰度值设为255，得到抹去直线后的雨量线图像，如图10所示。

步骤E、将雨量线图像对应到网格图像，求出雨量线上各个峰值点以及最右侧峰值点右侧延长线的最大值点，最后对以上各点所对应的雨量值求和，得到降雨量。

利用以上各步骤所得到的二值化雨量线图像、网格图像，即可求取出雨量值。雨量值包括两个部分：一是雨量曲线各峰值点所对应的记录值，二是量杯未满记录线(即最右侧峰值点右侧延长线的最大值点)的值。求取雨量值的方法具体如下：

(1)计算网格线每个刻度的宽高

①对网格线左上角300*300区域做行列投影，统计各行(列)黑色像素点个数；

②设定阈值Th，找到网格的横向和竖向刻度线；

③计算横(竖)向刻度线间像素的平均值，记作网格的高(宽)度。

(2)寻找曲线峰

①把刻度9-10的区间按网格宽度分成若干小块，分别统计每个小块内曲线点的个数；

②设定阈值Th，当步骤①统计值大于Th时认定该区域内存在曲线的峰；

③在步骤②得到的包含曲线峰的网格块中找到最高点，并判断该点右边25个像素内是否有比它还高的点，如果没有，则认定该点是一个峰值点，将它坐标保存；

(3)寻找量杯未满记录线的最高点

对最后一个峰值点往右的区域作处理，寻找这片区域中量杯未满记录线的最高点，通过统计该点左右一个网格宽度内黑色像素点的个数判断此点是曲线上的点还是孤立的噪点。如果是曲线上的点，则将该点保存。

(4)计算雨量值

①在网格图像中，极值点(x,y)向下一个网格高度范围内，对y左右各100个像素的区域做水平投影，统计每行黑色像素点的个数，超过140的最小的横坐标认定为(x,y)的下刻度线；

②与步骤①类似，从(x,y)往上一直到网格的顶端的范围内，对y左右各100个像素的区域做水平投影，得到(x,y)上端的各条刻度线；

③取最靠近(x,y)的上下刻度线，计算(x,y)到下刻度线的高度h和上下刻度线间的高度H，计算h/H。

④统计(x,y)上端刻度线的条数i，计算(x,y)所对应的值(10-i)+h/H；

⑤利用上述步骤计算出每个极值点的值，并将量杯未满纪录线上的最高值找出，加上前面的若干峰值，最终得到降雨量。

根据本发明方法可得到本发明雨量记录图纸的数字化装置，包括：

网格图像提取单元，用于从雨量记录图纸的数字图像中提取出二值化的网格图像； 雨量线图像提取单元，用于从雨量记录图纸的数字图像中提取出二值化的雨量线图像；

倾斜校正单元，用于对网格图像提取单元和雨量线图像提取单元所输出的网格图像、雨量线图像进行倾斜校正；

网格裁剪单元，用于利用行列投影的方法从倾斜校正后的网格图像中确定网格区域的边界，并将倾斜校正后的网格图像、雨量线图像中的边界外区域设置为二值化背景；

直线去除单元，用于将雨量线图像中的直线设置为二值化背景；

降雨量计算单元，用于将雨量线图像对应到网格图像，求出雨量线上各个峰值点以及最右侧峰值点右侧延长线的最大值点，最后对以上各点所对应的雨量值求和，得到降雨量。

本发明的数字化装置可使用预先通过扫描仪、数码相机等设备获取的雨量记录图纸的图像作为输入，并将降雨量计算单元得到的雨量数据存入外部数据库，以便后续的分析查询；本发明也可以进一步将图像获取设备和/或数据存储设备整合进来，实现全自动化的雨量数据提取。

实验结果表明，本发明能够对采用数码相机拍摄的雨量记录图纸自动地进行数字化分析，对普通的雨量记录图具有优秀的识别能力，获得的降雨量的值准确度也很高。