发明名称 基于三维随机球仿真的雪粒径谱分布特征参数优化方法
摘要 基于三维随机球仿真的雪粒径谱分布特征参数优化方法,本发明属于遥感与地理信息系统领域。本发明是要解决现有方法成像过程中将三维降维至二维图像时不可避免地存在信息损失并且提取的参数不够精确的问题,而提供了基于三维随机球仿真的雪粒径谱分布特征参数优化方法。一、获得二值化的积雪切片图像;二、利用二值圆拟合算法对二值化的积雪切片图像进行填充;三、根据步骤二中的填充圆列表ResultList对二维模拟切片图像的二维谱分布参数进行拟合与估计;四、统计二维模拟切片图像上的模粒径。本发明主要应用于卫星遥感对地观测领域。
申请公布号 CN103217367B 申请公布日期 2014.11.19
申请号 CN201310140875.7 申请日期 2013.04.22
申请人 中国科学院东北地理与农业生态研究所 发明人 李晓峰;赵凯;武黎黎
分类号 G01N15/02(2006.01)I 主分类号 G01N15/02(2006.01)I
代理机构 哈尔滨市松花江专利商标事务所 23109 代理人 金永焕
主权项 基于三维随机球仿真的雪粒径谱分布特征参数优化方法,其特征在于基于三维随机球仿真的雪粒径谱分布特征参数优化方法按以下步骤实现:一、二值化的积雪切片图像的获得;二、利用二值圆拟合算法对二值化的积雪切片图像进行填充,得到填充圆列表ResultList;三、根据步骤二中的填充圆列表ResultList对二维模拟切片图像的二维谱分布参数进行拟合与估计;四、统计二维模拟切片图像上的模粒径,并进一步和原始的二维谱分布参数进行误差分析统计输出三维优化雪粒径;其中,所述步骤二中利用二值圆拟合算法对二值化的积雪切片图像进行填充,具体步骤如下:步骤A:设定圆直径的阈值参数a<sub>T</sub>,对于雪粒径处理,该范围指定在0.2mm‑2.5mm;步骤B:在二值图像的白色区域放入尽可能大直径的圆,获得图像中填充圆列表ResultList,进一步得到雪粒径的大小及位置数据;B1:对于图像中的每一个白像素,初始化圆的半径r=0;B2:判断以该像元为圆心、r为半径的圆范围内是否包含值为0的像素或者图像边界,如果包含,r=r+1,继续判断;否则,转步骤B3;B3:判断r是否大于0,如果是,将该圆的数据即包括中心位置和半径存储在填充圆列表List,转入步骤B1,计算其他的白像素;否则,直接转步骤B1;B4:对List中的所有圆按照半径从大到小排序;B5:初始化填充圆列表ResultList;B6:判断List中是否还包含数据,如果是,执行步骤B7;否则,跳入B8输出ResultList;B7:删除List中最大元素x;找到List中与x相交叠的所有圆,放入列表TList之中;将TList的所有元素从List中删除;将x加入到ResultList;B8:输出填充圆列表ResultList;步骤C:在二值图像的白色区域继续利用尽可能大的圆填充剩余部分的白色区域;步骤D:检查是否剩余的白色区域能够利用阈值范围内的圆继续填充,如果能,转到步骤C,否则,结束;其中,所述步骤三中对二维模拟切片图像的二维谱分布参数进行拟合与估计具体步骤如下:雪粒径按照gamma分布进行统计:n(a)=K<sub>1</sub>a<sup>P</sup>exp(‑K<sub>2</sub>a<sup>Q</sup>)        (1)其中,a表示雪粒径;雪粒子的占空比f<sub>tot</sub>可以用下式表示为:<maths num="0001" id="cmaths0001"><math><![CDATA[<mrow><msub><mi>f</mi><mi>tot</mi></msub><mo>=</mo><msubsup><mo>&Integral;</mo><mn>0</mn><mrow><mo>+</mo><mo>&infin;</mo></mrow></msubsup><mfrac><msup><mrow><mn>4</mn><mi>&pi;a</mi></mrow><mn>3</mn></msup><mn>3</mn></mfrac><mi>n</mi><mrow><mo>(</mo><mi>a</mi><mo>)</mo></mrow><mi>da</mi><mo>=</mo><mfrac><mrow><mn>4</mn><mi>&pi;</mi></mrow><mn>3</mn></mfrac><mfrac><msub><mi>K</mi><mn>1</mn></msub><mi>Q</mi></mfrac><msubsup><mi>K</mi><mn>2</mn><mrow><mo>-</mo><mfrac><mrow><mi>P</mi><mo>+</mo><mn>4</mn></mrow><mi>Q</mi></mfrac></mrow></msubsup><mi>&Gamma;</mi><mrow><mo>(</mo><mfrac><mrow><mi>P</mi><mo>+</mo><mn>4</mn></mrow><mi>Q</mi></mfrac><mo>)</mo></mrow><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>2</mn><mo>)</mo></mrow></mrow>]]></math><img file="FDA0000553560080000021.GIF" wi="1395" he="145" /></maths>其中,Γ表示gamma函数,n(a)da表示单位体积内雪粒径置于a和a+da之间的雪粒子数,参数K<sub>2</sub>与模粒径a<sub>m</sub>相关,模粒径a<sub>m</sub>从式(1)对a的微分为零得到,<maths num="0002" id="cmaths0002"><math><![CDATA[<mrow><msub><mi>a</mi><mi>m</mi></msub><mo>=</mo><msup><mrow><mo>(</mo><mfrac><mi>P</mi><mrow><msub><mi>K</mi><mn>2</mn></msub><mi>Q</mi></mrow></mfrac><mo>)</mo></mrow><mfrac><mn>1</mn><mi>Q</mi></mfrac></msup><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>3</mn><mo>)</mo></mrow></mrow>]]></math><img file="FDA0000553560080000022.GIF" wi="1088" he="186" /></maths>步骤A:根据步骤二中的不同半径的红色圆盘的统计列表,根据雪粒径范围进行分组,得到k种不同雪粒径的雪粒子的数密度n<sub>1</sub>(a<sub>1</sub>),n<sub>2</sub>(a<sub>2</sub>),…,n<sub>k</sub>(a<sub>k</sub>),以及模粒径<img file="FDA0000553560080000023.GIF" wi="81" he="78" />根据黑色区域和白色区域所占的比例,计算出占空比<img file="FDA0000553560080000024.GIF" wi="97" he="79" />其中,所述<img file="FDA0000553560080000025.GIF" wi="266" he="79" />A<sub>W</sub>表示白色区域的面积,A<sub>B</sub>表示黑色区域所占面积;步骤B:通过获得的这些参数,利用公式(1)、(2)和(3),得到待优化的方程组如下:<maths num="0003" id="cmaths0003"><math><![CDATA[<mrow><mfenced open='{' close=''><mtable><mtr><mtd><msubsup><mi>a</mi><mi>m</mi><mn>0</mn></msubsup><mo>=</mo><msup><mrow><mo>(</mo><mfrac><mi>P</mi><mrow><msub><mi>K</mi><mn>2</mn></msub><mi>Q</mi></mrow></mfrac><mo>)</mo></mrow><mfrac><mn>1</mn><mi>Q</mi></mfrac></msup></mtd></mtr><mtr><mtd><msubsup><mi>f</mi><mi>tot</mi><mn>0</mn></msubsup><mo>=</mo><mfrac><mrow><mn>4</mn><mi>&pi;</mi></mrow><mn>3</mn></mfrac><mfrac><msub><mi>K</mi><mn>1</mn></msub><mi>Q</mi></mfrac><msubsup><mi>K</mi><mn>2</mn><mrow><mo>-</mo><mfrac><mrow><mi>P</mi><mo>+</mo><mn>4</mn></mrow><mi>Q</mi></mfrac></mrow></msubsup><mi>&Gamma;</mi><mrow><mo>(</mo><mfrac><mrow><mi>P</mi><mo>+</mo><mn>4</mn></mrow><mi>Q</mi></mfrac><mo>)</mo></mrow></mtd></mtr><mtr><mtd><msub><mi>n</mi><mn>1</mn></msub><mrow><mo>(</mo><msub><mi>a</mi><mn>1</mn></msub><mo>)</mo></mrow><mo>=</mo><msub><mi>K</mi><mn>1</mn></msub><msubsup><mi>a</mi><mn>1</mn><mi>P</mi></msubsup><mi>exp</mi><mrow><mo>(</mo><mo>-</mo><msub><mi>K</mi><mn>2</mn></msub><msubsup><mi>a</mi><mn>1</mn><mi>Q</mi></msubsup><mo>)</mo></mrow></mtd></mtr><mtr><mtd><msub><mi>n</mi><mn>2</mn></msub><mrow><mo>(</mo><msub><mi>a</mi><mn>2</mn></msub><mo>)</mo></mrow><mo>=</mo><msub><mi>K</mi><mn>1</mn></msub><msubsup><mi>a</mi><mn>2</mn><mi>P</mi></msubsup><mi>exp</mi><mrow><mo>(</mo><mo>-</mo><msub><mi>K</mi><mn>2</mn></msub><msubsup><mi>a</mi><mn>2</mn><mi>Q</mi></msubsup><mo>)</mo></mrow></mtd></mtr><mtr><mtd><mo>.</mo><mo>.</mo><mo>.</mo></mtd></mtr><mtr><mtd><msub><mi>n</mi><mi>k</mi></msub><mrow><mo>(</mo><msub><mi>a</mi><mi>k</mi></msub><mo>)</mo></mrow><mo>=</mo><msub><mi>K</mi><mn>1</mn></msub><msubsup><mi>a</mi><mi>k</mi><mi>P</mi></msubsup><mi>exp</mi><mrow><mo>(</mo><msub><mrow><mo>-</mo><mi>K</mi></mrow><mn>2</mn></msub><msubsup><mi>q</mi><mi>k</mi><mi>Q</mi></msubsup><mo>)</mo></mrow></mtd></mtr></mtable></mfenced><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>4</mn><mo>)</mo></mrow></mrow>]]></math><img file="FDA0000553560080000026.GIF" wi="1246" he="640" /></maths>步骤C:利用Matlab的非线性方程组求解算法对方程组(4)求最优解;步骤D:获得参数K<sub>1</sub>,K<sub>2</sub>,P和Q,从而进一步得到二维谱分布参数表示为f<sub>tot</sub>,a<sub>m</sub>,P和Q;其中,所述步骤四采用三维随机球仿真算法三维优化雪模粒径,具体步骤如下:步骤A:利用二维谱分布参数计算雪粒径的谱分布参数与输入离散粒径的种类数k,然后在三维单位体积的随机位置模拟不同粒径的粒子;其中,所述谱分布参数为P,Q,f<sub>tot</sub>和a<sub>m</sub>;A1、按照初始输入的二维雪粒径谱分布参数P,Q,f<sub>tot</sub>和a<sub>m</sub>将雪粒径以固定步长离散化为N种,N=2000;A2、利用公式(1)‑(3)计算每种离散化后的雪粒径的数密度及其在所有雪粒径中所占比例,设定阈值T1,如果比例大于阈值,则将该雪粒径记录至表Radius‑list;否则,剔除;A3、选择Radius‑list中最大值和最小值,从而确定粒径范围[a<sub>min</sub>,a<sub>max</sub>];A4、根据输入的雪粒子种类k和[a<sub>min</sub>,a<sub>max</sub>]计算第N种粒径的粒子数密度n(a<sub>k</sub>);其中,N≤k;A5、在单位体积内以随机点为球心,根据每种雪粒子的数密度利用三维随机球模拟雪粒子,得到三维随机球仿真雪粒子分布表3DRBS‑list;步骤B:三维图像切片:模拟切片图像的切割制作过程,从不同的位置随机切割;B1、根据设置的间隔步长ds计算第M个球心与切割平面Z=i·ds的距离D<sub>Mi</sub>;Z&lt;1;B2、比较该距离与该球直径d<sub>M</sub>的大小,记录D<sub>Mi</sub>&lt;d<sub>M</sub>的球以及被平面所切得的圆半径;步骤C:根据步骤B2中的切割圆半径统计二维模拟切片图像上的模粒径<img file="FDA0000553560080000031.GIF" wi="95" he="78" />步骤D:误差分析<img file="FDA0000553560080000032.GIF" wi="256" he="72" />并进一步和原始的输入参数进行误差分析;D1、如果误差不满足精度要求即ε&gt;T,T为预设的阈值,那么输入新的参数进行迭代;D2、否则,输出三维优化后的模粒径。
地址 150081 黑龙江省哈尔滨市南岗区哈平路138号