CN103217367A - 基于三维随机球仿真的雪粒径谱分布特征参数优化方法 - Google Patents

Abstract

基于三维随机球仿真的雪粒径谱分布特征参数优化方法，本发明属于遥感与地理信息系统领域。本发明是要解决现有方法成像过程中将三维降维至二维图像时不可避免地存在信息损失并且提取的参数不够精确的问题，而提供了基于三维随机球仿真的雪粒径谱分布特征参数优化方法。一、获得二值化的积雪切片图像；二、利用二值圆拟合算法对二值化的积雪切片图像进行填充；三、根据步骤二中的填充圆列表ResultList对二维模拟切片图像的二维谱分布参数进行拟合与估计；四、统计二维模拟切片图像上的模粒径。本发明主要应用于卫星遥感对地观测领域。

Description

基于三维随机球仿真的雪粒径谱分布特征参数优化方法

技术领域

本发明属于遥感与地理信息系统领域，主要应用于卫星遥感对地观测。

背景技术

地表雪盖对于全球水文循环、大气过程以及全球气候变化均具有重要影响，多年来卫星遥感技术作为重要的工具应用于大尺度的雪盖参数监测和反演，而具有全天时和全天候监测能力的微波遥感技术对于雪深和雪水当量的监测更是其他手段所无法替代。微波遥感辐射传输模型是微波遥感反演地表参数的理论基础，而研究表明积雪粒子的粒径是辐射传输模型最敏感的参数之一，但这一参数的测量一直以来具有很大的主观性和不确定性。

实测雪粒径的方法有多种，一个简单办法是将颗粒放到有毫米刻度的盘上，通过比较颗粒尺寸与盘上刻度线的间隔，可估计平均尺寸，这种目测估计的方法具有一定的主观性；另外的测量方法是筛分法，即选择一系列不同筛孔直径的标准筛，按照孔径从小到大依次摞起，固定在振筛机上，选择适当的模式及时长，自动振动即可实现筛分，这种测量方法易于操作，但操作过程中容易破坏雪层的固有结构和粒径，使得测量结果不够准确。为了准确模拟雪层的微波散射和辐射特征需要明确知道粒径的统计分布特征，而以上两种方法不易做到。

雪粒子的显微图像采集与处理的方法以其采集过程准确，测量结果客观，不破坏雪层结构等优势为许多研究者所采用，如立体测量法，图像采集与处理法等，此类方法通过体视显微镜拍摄雪粒子的图像，并通过数字图像处理技术(二值化，圆形膨胀拟合)获得粒径尺寸的范围和分布，而由此获得的粒径所服从的概论分布对于随机分布的雪粒子的微波(散射、辐射)特性的准确计算和模拟尤为重要。

但雪粒子的显微图像采集与处理的方法同时存在以下问题：

首先，由于积雪切片图像获取本身的步骤较多，实验室控制要求比较严格，所以获取的图像十分有限；

其次，获取图像过程中需要对雪层进行切割处理，而图像也是基于切割的断面采集的，我们实际上最终是为了获得雪层中三维粒径的分布，所以由二维图像推测三维特征的过程中本身就存在损失和较大误差，获得的参数未必能够体现雪层中雪粒径分布的三维特征。

发明内容

本发明是要解决现有方法成像过程中将三维降维至二维图像时不可避免地存在信息损失并且提取的参数自然不够精确的问题，而提供了基于三维随机球仿真的雪粒径谱分布特征参数优化方法。

基于三维随机球仿真的雪粒径谱分布特征参数优化方法按以下步骤实现：

一、二值化的积雪切片图像的获得；

二、利用二值圆拟合算法对二值化的积雪切片图像进行填充，得到不同半径的红色圆盘的统计列表；

三、根据步骤二中的不同半径的红色圆盘的统计列表对二维模拟切片图像的二维谱分布参数进行拟合与估计；

四、统计二维模拟切片图像上的模粒径，并进一步和原始的二维谱分布参数进行误差分析统计输出三维优化雪粒径。

发明效果：

(1)尽量从较少的采样图像中获取更多的信息，能够弥补三维到二维时的信息损失，使得提取的雪粒径分布参数更加准确真实，从而避免了需要更多的二维切片图像所带来的累计误差，总体精度可提高90％以上；

(2)利用初始的雪粒径参数进行计算机模拟得到三维随机球拟合图像，模拟切片过程，多次随机产生切片图像，避免了由于切片图像制作过程切割位置不同所带来的误差，使得获得的优化参数获取更加准确，为积雪微波辐射传输模型的研究以及地表积雪参数的准确提取提供了技术保障。

附图说明

图1是本发明的整体流程图；

图2是具体实施方式一中的积雪样本的二值化切片图像；

图3是具体实施方式二中的输入的图像进行圆填充流程图；

图4是具体实施方式三中的二值图像圆拟合算法获得的粒径拟合结果；

图5是具体实施方式三中的二维谱分布参数的估计流程图；

图6是具体实施方式三中的三维切片未必能够恰好切至雪粒子的直径处的演示图；

图7是具体实施方式三中的二维雪粒径分布函数图；其中，a_m＝0.38，P＝9，Q＝2，f_tot＝0.4；

图8是具体实施方式四中的三维随机球仿真(3DRBS)优化算法的实现流程图；

图9是具体实施方式四中的单位体积内(1cm³)三维模拟雪粒径图；

图10是具体实施方式四中的模拟三维雪粒子的在0.5cm处的切片图像；

图11是具体实施方式四中的模拟三维雪粒子的在0.8cm处的切片图像；

图12是具体实施方式四中的迭代的分布曲线变化过程，虚线是初始的分布曲线，实线为迭代曲线。

具体实施方式

具体实施方式一：本实施方式的基于三维随机球仿真的雪粒径谱分布特征参数优化方法按以下步骤实现：

一、二值化的积雪切片图像的获得；

二、利用二值圆拟合算法对二值化的积雪切片图像进行填充，得到不同半径的红色圆盘的统计列表；

三、根据步骤二中的不同半径的红色圆盘的统计列表对二维模拟切片图像的二维谱分布参数进行拟合与估计；

四、统计二维模拟切片图像上的模粒径，并进一步和原始的二维谱分布参数进行误差分析统计输出三维优化雪粒径。

本实施方式效果：

(1)尽量从较少的采样图像中获取更多的信息，能够弥补三维到二维时的信息损失，使得提取的雪粒径分布参数更加准确真实，从而避免了需要更多的二维切片图像所带来的累计误差，总体精度可提高90％以上；

(2)利用初始的雪粒径参数进行计算机模拟得到三维随机球拟合图像，模拟切片过程，多次随机产生切片图像，避免了由于切片图像制作过程切割位置不同所带来的误差，使得获得的优化参数获取更加准确，为积雪微波辐射传输模型的研究以及地表积雪参数的准确提取提供了技术保障。

具体实施方式二：本实施方式与具体实施方式一不同的是：步骤二中利用二值圆拟合算法对二值化的积雪切片图像进行填充具体步骤如下：

步骤A：设定圆直径的阈值参数a_T，对于雪粒径处理，该范围指定在0.2mm-2.5mm；

步骤B：在二值图像的白色区域(雪粒子)放入尽可能大直径的圆，获得图像中的填充圆列表ResultList，进一步得到雪粒径的大小及位置数据；

B1：对于图像中的每一个白像素(对应数组值为1)，初始化圆的半径r＝0；

B2：判断以该像元为圆心、r为半径的圆范围内是否包含值为0的像素或者图像边界，如果包含，r＝r+1，继续判断；否则，转步骤B3；

B3：判断r是否大于0，如果是，将该圆的数据(包括：中心位置和半径)存储在填充圆列表List，转入步骤B1，计算其他的白像素；否则，直接转步骤B1；

B4：对List中的所有圆按照半径从大到小排序；

B5：初始化填充圆列表ResultList；

B6：判断List中是否还包含数据，如果是，执行步骤B8；否则，跳入B8输出ResultList；

B7：删除List中最大元素x；找到List中与x相交叠的所有圆，放入列表TList之中；将TList的所有元素从List中删除；将x加入到ResultList；

B8：输出填充圆列表ResultList；

步骤C：在二值图像的白色区域(雪粒子)继续利用尽可能大的圆填充剩余部分的白色区域；

步骤D：检查是否剩余的白色区域能够利用阈值范围内的圆继续填充，如果能，转到步骤C，否则，结束。

其它步骤及参数与具体实施方式一相同。

具体实施方式三：本实施方式与具体实施方式一或二不同的是：步骤三中对二维模拟切片图像的二维谱分布参数进行拟合与估计具体步骤如下：

雪粒径按照gamma分布进行统计：

n(a)＝K₁a^P exp(-K₂a^Q)  (1)

其中，a表示雪粒径；

雪粒子的占空比f_t0t可以用下式表示为：

$f_{\mathrm{tot}}={\∫}_0^{+\∞}\frac{{4\mathrm{\πa}}^3}{3}n\left(a\right)\mathrm{da}=\frac{4\mathrm{\π}}{4}\frac{K_1}{Q}{K}_2^{-\frac{P+4}{Q}}\mathrm{\Γ}\left(\frac{P+4}{Q}\right)---\left(2\right)$

其中，Γ表示gamma函数，n(a)da表示单位体积内雪粒径置于a和a+da之间的雪粒子数，参数K₂与模粒径a_m相关，模粒径a_m可以从式(1)对a的微分为零得到，

$a_m={\left(\frac{P}{K_{2}Q}\right)}^{\frac{1}{Q}}---\left(3\right)$

步骤A：根据步骤二中的不同半径的红色圆盘的统计列表，根据雪粒径范围进行分组，得到k种不同雪粒径的雪粒子的数密度n₁(a₁)，n₂(a₂)，...，n_k(a_k)，以及模粒径根据黑色区域和白色区域所占的比例，计算出占空比

其中，所述

A_W表示白色区域的面积，A_B表示黑色区域所占面积；

步骤B：通过获得的这些参数，利用公式(1)、(2)和(3)，得到待优化的方程组如下：

$\left\{\begin{array}{c}{a}_m^0={\left(\frac{P}{K_{2}Q}\right)}^{\frac{1}{Q}}\\ f_{\mathrm{tot}}^0=\frac{4\mathrm{\π}}{3}\frac{K_1}{Q}{K}_2^{-\frac{P+4}{Q}}\mathrm{\Γ}\left(\frac{P+4}{Q}\right)\\ n_1\left(a_1\right)=K_1{a}_1^P\exp \left({-K}_2{a}_1^Q\right)\\ n_2\left(a_2\right)=K_1{a}_2^P\exp \left({-K}_2{a}_2^Q\right)\\ \·\·\·\\ n_k\left(a_k\right)=K_1{a}_k^P\exp \left({-K}_2{a}_k^Q\right)\end{array}\right.---\left(4\right)$

步骤C：利用Matlab的非线性方程组求解算法对方程组(4)求最优解；

步骤D：获得参数K₁，K₂，P和Q，从而进一步得到二维谱分布参数表示为f_tot，a_m，P和Q。

其它步骤及参数与具体实施方式一或二相同。

具体实施方式四：本实施方式与具体实施方式一至三之一不同的是：所述步骤四采用三维随机球仿真算法三维优化雪模粒径具体步骤如下步骤A：利用二维分布参数计算雪粒径的谱分布参数与输入离散粒径的种类数k，然后在三维单位体积的随机位置模拟不同粒径的粒子；其中，所述谱分布参数为P，Q，f_t0t和a_m；

A1、按照初始输入的二维雪粒径谱分布参数P，Q，f_t0t和a_m将雪粒径以固定步长离散化为N种，N＝2000；

A2、利用公式(1)-(3)计算每种离散化后的雪粒径的数密度及其在所有雪粒径中所占比例，设定阈值T1，如果比例大于阈值，则将该雪粒径记录至表Radius-list；否则，剔除；

A3、选择Radius-list中最大值和最小值，从而确定粒径范围[a_min，a_max]；

A4、根据输入的雪粒子种类k和[a_min，a_max]计算第N(N＜＝k)种粒径的粒子数密度n(a_k)；

A5、在单位体积内以随机点为球心，根据每种雪粒子的数密度利用三维随机球模拟雪粒子，得到三维随机球仿真雪粒子分布表3DRBS-list；

步骤B：三维图像切片：模拟切片图像的切割制作过程，从不同的位置随机切割；

B1、根据设置的间隔步长ds计算第M个球心与切割平面Z＝i·ds，(Z＜1)的距离D_Mi；

B2、比较该距离与该球直径d_M的大小，记录D_Mi＜d_M的球以及被平面所切得的圆半径；

步骤C：根据步骤B2中的切割圆半径统计二维模拟切片图像上的模粒径

步骤D：误差分析并进一步和原始的输入参数进行误差分析；

D1、如果误差不满足精度要求(ε＞T，T为预设的阈值)，那么输入新的参数进行迭代；

D2、否则，输出三维优化后的模粒径。

其它步骤及参数与具体实施方式一至三之一相同。

图1是，积雪样本的二值化切片图像；

由图4二值图像圆拟合算法获得的粒径拟合结果可知，从图中可以看出我们的圆拟合算法达到了很好的拟合效果。

由图6可知，三维切片未必能够恰好切至雪粒子的直径处的演示图，说明了对三维雪粒子进行处理生成二维图像并获得雪粒径时会带来不可避免的误差；

由图12可知，迭代的分布曲线变化过程，针对权利要求4中的步骤C给出的具体运算实例，得到的优化模粒径为

算法实施过程中选择的阈值T＝0.001cm，这意味着当三维模粒径为0.047cm时，模拟的模粒径和实际测量的模粒径之间的误差可以控制在0.001cm之内。

Claims (4)

1.基于三维随机球仿真的雪粒径谱分布特征参数优化方法，其特征在于基于三维随机球仿真的雪粒径谱分布特征参数优化方法按以下步骤实现：

一、二值化的积雪切片图像的获得；

二、利用二值圆拟合算法对二值化的积雪切片图像进行填充，得到填充圆列表ResultList；

三、根据步骤二中的填充圆列表ResultList对二维模拟切片图像的二维谱分布参数进行拟合与估计；

四、统计二维模拟切片图像上的模粒径，并进一步和原始的二维谱分布参数进行误差分析统计输出三维优化雪粒径。

2.根据权利要求1所述的基于三维随机球仿真的雪粒径谱分布特征参数优化方法，其特征在于步骤二中利用二值圆拟合算法对二值化的积雪切片图像进行填充，具体步骤如下：

步骤A：设定圆直径的阈值参数a_T，对于雪粒径处理，该范围指定在0.2mm-2.5mm；

步骤B：在二值图像的白色区域(雪粒子)放入尽可能大直径的圆，获得图像中填充圆列表ResultList，进一步得到雪粒径的大小及位置数据；

B1：对于图像中的每一个白像素(对应数组值为1)，初始化圆的半径r＝0；

B2：判断以该像元为圆心、r为半径的圆范围内是否包含值为0的像素或者图像边界，如果包含，r＝r+1，继续判断；否则，转步骤B3；

B3：判断r是否大于0，如果是，将该圆的数据(包括：中心位置和半径)存储在填充圆列表List，转入步骤B1，计算其他的白像素；否则，直接转步骤B1；

B4：对List中的所有圆按照半径从大到小排序；

B5：初始化填充圆列表ResultList；

B6：判断List中是否还包含数据，如果是，执行步骤B8；否则，跳入B8输出ResultList；

B7：删除List中最大元素x；找到List中与x相交叠的所有圆，放入列表TList之中；将TList的所有元素从List中删除；将x加入到ResultList；

B8：输出填充圆列表ResultList；

步骤C：在二值图像的白色区域(雪粒子)继续利用尽可能大的圆填充剩余部分的白色区域；

步骤D：检查是否剩余的白色区域能够利用阈值范围内的圆继续填充，如果能，转到步骤C，否则，结束。

3.根据权利要求1所述的基于三维随机球仿真的雪粒径谱分布特征参数优化方法，其特征在于步骤三中对二维模拟切片图像的二维谱分布参数进行拟合与估计具体步骤如下：

雪粒径按照gamma分布进行统计：

n(a)＝K₁a^P exp(-K₂a^Q)  (1)

其中，a表示雪粒径；

雪粒子的占空比f_tot可以用下式表示为：

$f_{\mathrm{tot}}={\∫}_0^{+\∞}\frac{{4\mathrm{\πa}}^3}{3}n\left(a\right)\mathrm{da}=\frac{4\mathrm{\π}}{4}\frac{K_1}{Q}{K}_2^{-\frac{P+4}{Q}}\mathrm{\Γ}\left(\frac{P+4}{Q}\right)---\left(2\right)$

其中，Γ表示gamma函数，n(a)da表示单位体积内雪粒径置于a和a+da之间的雪粒子数，参数K₂与模粒径a_m相关，模粒径a_m可以从式(1)对a的微分为零得到，

$a_m={\left(\frac{P}{K_{2}Q}\right)}^{\frac{1}{Q}}---\left(3\right)$

步骤A：根据步骤二中的不同半径的红色圆盘的统计列表，根据雪粒径范围进行分组，得到k种不同雪粒径的雪粒子的数密度n₁(a₁)，n₂(a₂)，...，n_k(a_k)，以及模粒径

根据黑色区域和白色区域所占的比例，计算出占空比

其中，所述

A_W表示白色区域的面积，A_B表示黑色区域所占面积；

步骤B：通过获得的这些参数，利用公式(1)、(2)和(3)，得到待优化的方程组如下：

$\left\{\begin{array}{c}{a}_m^0={\left(\frac{P}{K_{2}Q}\right)}^{\frac{1}{Q}}\\ f_{\mathrm{tot}}^0=\frac{4\mathrm{\π}}{3}\frac{K_1}{Q}{K}_2^{-\frac{P+4}{Q}}\mathrm{\Γ}\left(\frac{P+4}{Q}\right)\\ n_1\left(a_1\right)=K_1{a}_1^P\exp \left({-K}_2{a}_1^Q\right)\\ n_2\left(a_2\right)=K_1{a}_2^P\exp \left({-K}_2{a}_2^Q\right)\\ \·\·\·\\ n_k\left(a_k\right)=K_1{a}_k^P\exp \left({-K}_2{a}_k^Q\right)\end{array}\right.---\left(4\right)$

步骤C：利用Matlab的非线性方程组求解算法对方程组(4)求最优解；

步骤D：获得参数K₁，K₂，P和Q，从而进一步得到二维谱分布参数表示为f_tot，a_m，P和Q。

4.根据权利要求1所述的基于三维随机球仿真的雪粒径谱分布特征参数优化方法，其特征在于所述步骤四采用三维随机球仿真算法三维优化雪模粒径，具体步骤如下：

步骤A：利用二维分布参数计算雪粒径的谱分布参数与输入离散粒径的种类数k，然后在三维单位体积的随机位置模拟不同粒径的粒子；其中，所述谱分布参数为P，Q，f_tot和a_m；

A1、按照初始输入的二维雪粒径谱分布参数P，Q，f_tot和a_m将雪粒径以固定步长离散化为N种，N＝2000；

A2、利用公式(1)-(3)计算每种离散化后的雪粒径的数密度及其在所有雪粒径中所占比例，设定阈值T1，如果比例大于阈值，则将该雪粒径记录至表Radius-list；否则，剔除；

A3、选择Radius-list中最大值和最小值，从而确定粒径范围[a_min，a_max]；

A4、根据输入的雪粒子种类k和[a_min，a_max]计算第N(N＜＝k)种粒径的粒子数密度n(a_k)；

A5、在单位体积内以随机点为球心，根据每种雪粒子的数密度利用三维随机球模拟雪粒子，得到三维随机球仿真雪粒子分布表3DRBS-list；

步骤B：三维图像切片：模拟切片图像的切割制作过程，从不同的位置随机切割；

B1、根据设置的间隔步长ds计算第M个球心与切割平面Z＝i·ds，(Z＜1)的距离D_Mi；

B2、比较该距离与该球直径d_M的大小，记录D_Mi＜d_M的球以及被平面所切得的圆半径；

步骤C：根据步骤B2中的切割圆半径统计二维模拟切片图像上的模粒径

步骤D：误差分析

并进一步和原始的输入参数进行误差分析；

D1、如果误差不满足精度要求(ε＞T，T为预设的阈值)，那么输入新的参数进行迭代；

D2、否则，输出三维优化后的模粒径。

Images