高斯混合模型自动阈值法遥感冰川信息提取

高斯混合模型自动阈值法遥感冰川信息提取 Extracting glacier information from remote sensing imageries by automatic threshold method of Gaussian mixture model first-author 王亚利，1983年生， 女， 高级工程师，研究方向为遥感地物提取，冰川变化分析。E-mail： wangyali1983021@163.com 王亚利，1983年生， 女， 高级工程师，研究方向为遥感地物提取，冰川变化分析。E-mail： wangyali1983021@163.com wangyali1983021@163.com corresp 都伟冰，1985年生，男，讲师，研究方向为多源遥感冰川监测。E-mail：dwb@hpu.edu.cn 都伟冰，1985年生，男，讲师，研究方向为多源遥感冰川监测。E-mail：dwb@hpu.edu.cn dwb@hpu.edu.cn 河南理工大学 测绘与国土信息工程学院, 焦作 454000 河南理工大学 测绘与国土信息工程学院, 焦作 454000 Henan Polytechnic University, School of Surveying and Land Information Engineering, Jiaozuo 454000, China Henan Polytechnic University, School of Surveying and Land Information Engineering, Jiaozuo 454000, China 利用阈值法进行遥感地物提取效率高、准确率高，但是在阈值的选取方面，传统的手动选取阈值，自动化低，需反复试操作，且易受主观因素影响。文中通过期望最大算法对局部冰川区域归一化雪覆盖指数建立高斯混合模型，去除区域内的混合像元类，再利用高斯混合模型模拟纯净化后的冰川类、非冰川类的NDSI分布情况，根据改进后的高斯混合模型分布情况，自动计算出区域内的冰川提取阈值。本文对不同海拔的3个区域进行算法实验，然后将新疆哈密的哈尔里克山冰川提取边界与冰川编目数据进行对比验证。研究结果表明，该方法自动计算的冰川提取阈值结果可靠、精度高，在差异较大区域仍较稳定，有一定的应用价值。 Extraction of glacier by remote sensing has important value of scientific research and application. The threshold segmentation method has high efficiency and accuracy in remote sensing feature extraction. However, in terms of threshold selection, traditional manual threshold selection demonstrate low automation, repeated trial operation, and is easily affected by subjective factors. For example, the normalized snow cover index can separate the glacier information from the background information effectively. However, the threshold is almost always found manually, and hence is easily affected by subjective factors. Therefore, developing an effective method is necessary to calculate the glacier extraction threshold automatically. Gaussian Mixture Model (GMM) is established for the Normalized Difference Snow Index (NDSI) of the local glacier region by the expected maximum algorithm and to remove the mixed pixel meta classes in the region. Then, the GMM was used to simulate the NDSI distribution of purified glacier and non-glacier. According to the distribution of the improved Gaussian mixture model, the glacier extraction threshold in the region was calculated automatically.Three regions at different altitudes were experimented. Subsequently, the boundary extraction of Karlik mountain glacier in Hami, Xinjiang was compared with the glacier inventory data, and good extraction results were obtained. The experimental results show that the threshold values calculated by different grids in the same experimental area are slightly different, whereas the threshold values calculated by different grids in various experimental areas are largely different, thereby also indicating that the algorithm is still stable in the regions with large differences. As such, the results of glacier extraction threshold calculated by using this method are reliable, accurate, and stable in different regions. Result also shows that this method solves the problem of the unbalanced proportion of glacier and non-glacier in GMM curve simulation to some extent, and the automatic calculation results of glacier extraction threshold are reliable, accurate, and stable in regions with large differences. This GMM method is available for certain applications of delineating glacier boundary in remote sensing. 遥感影像 高斯混合模型 归一化雪覆盖指数 冰川提取 remote sensing image Gaussian Mixture Model NDSI extracting glacier information 1　引言 1 中国是中低纬度高山冰川占比较高的国家，世界冰川编目WGI（World Glacier Inventory）数据结果显示，中国的冰储量仅次于加拿大、俄罗斯和美国（ Guo等，2015 Guo等，2015 ）。高山冰川是区域气候变化的敏感指示器，也是中国重要的固态淡水资源，中国主要的河流都有冰雪融水补给，尤其是西部干旱区的河流径流量很大程度上依赖冰雪融水（ Bolch等，2012 Bolch等，2012 ； Du等，2014 Du等，2014 ； Scherler等，2011 Scherler等，2011 ； 冯童 等，2015 冯童 等，2015 ； 刘时银 等，2015 刘时银 等，2015 ）。冰川变化是气候变化的产物，气候变冷或变暖引起中国高山冰川的积累和消融，导致冰川面积的增减、末端位置的前进或退缩。如果冰川退缩加剧，冰雪融水进入旺盛期，会形成一系列冰川次生灾害。另外，近20年来全球气候变化加剧，导致冰川退缩，部分小冰川消失（ Haeberli和Hoelzle，1995 Haeberli和Hoelzle，1995 ； 施雅风，2001 施雅风，2001 ）。通过对冰川变化研究，可以分析冰川变化特征和波动趋势，认识冰川变化与环境气候变化的关系。因此，山地冰川变化研究备受政府和学术界的关注。其中，利用遥感手段进行冰川提取的方法具有重要科研、应用价值。 中国高山冰川地处偏远，气候环境相对恶劣且数量众多，交通条件极为不便，仅靠实测手段难以实现高效快速地大范围冰川监测。卫星遥感技术的快速发展，促进了冰川监测技术的提高，遥感手段广泛应用于不同区域、不同时相、不同背景的冰川变化研究中（ 张晓博 等，2018 张晓博 等，2018 ）。目前，针对遥感影像自动、半自动化冰川信息提取技术的研究已经取得了一定的成果，阈值分割法和基于分类的方法是提取冰川信息的主要手段（ 都伟冰 等，2015 都伟冰 等，2015 ； 彦立利和王建，2013 彦立利和王建，2013 ）。分类方法在较小局部区域提取精细冰川边界，效果较好，但算法复杂，尤其是监督分类对样本的依赖较强，不利于冰川信息的快速提取。本文主要研究阈值分割法，比值阈值法利用冰川在可见光波段和近红外波段反射率的比值提取冰川边界（ 李佳，2014 李佳，2014 ），归一化雪覆盖指数也能有效地将冰川信息与背景信息分离，但阈值几乎都是手工寻找，易受到主观因素的影响。最大类间方差法以其发明人Otsu命名，又称Otsu方法，是常用的阈值自动计算方法（ 赵健赟和彭军还，2016 赵健赟和彭军还，2016 ）。Otsu方法和本文方法本质上都是利用图像 DN （Digital Number）值的高斯分布特征自动计算阈值。差别在于，（1）Otsu方法大多用于人物照片、景物照片的美化处理，其目标与背景过度区的混合像元干扰因素不大；（2）而本文方法中混合像元是影响冰川边界精确提取的重要因素。本文在“3.4节”对两种方法进行了对比。 由此，本文根据局部遥感影像冰川与非冰川区域具有高斯混合模型分布特点，构建高斯混合模型自动计算冰川提取阈值。主要工作内容如下，（1）对局部冰川区域归一化雪覆盖指数建立高斯混合模型，计算区域内的混合像元，合理去除混合像元所涉及图像 DN 值。（2）利用高斯混合模型模拟纯净化后的冰川类、非冰川类的NDSI值分布情况，根据模型中的两个正态分布，自动计算出区域内的冰川提取阈值。（3）以不同高度、不同时相、不同区域的3个实验区对本文方法进行分析验证等工作，即天山、祁连山、羌塘高原，另选用新疆东部哈尔里克山脉的冰川提取边界与冰川编目数据进行对比和精度分析，对提取结果进行评价。 2　遥感冰川提取的阈值计算原理 1 2.1　遥感冰川提取方法 2 归一化雪覆盖指数NDSI（Normalized Difference Snow Index）利用冰雪在短波红外SWIR（Short Wave InfraRed）波段有强烈的吸收特性，在可见光绿光G（Green）波段有强的反射特性，能很好地区分冰雪与背景信息，将遥感影像中的G波段和SWIR波段的反射率进行归一化处理，以突出冰川与背景信息的差异，其公式如下： <math display="block" id="M1"><mi mathvariant="normal">N</mi><mi mathvariant="normal">D</mi><mi mathvariant="normal">S</mi><mi mathvariant="normal">I</mi><mo>=</mo><mo stretchy="false">(</mo><msub><mrow><mi>ρ</mi></mrow><mrow><mi mathvariant="normal">G</mi></mrow></msub><mo>-</mo><msub><mrow><mi>ρ</mi></mrow><mrow><mi mathvariant="normal">S</mi><mi mathvariant="normal">W</mi><mi mathvariant="normal">I</mi><mi mathvariant="normal">R</mi></mrow></msub><mo stretchy="false">)</mo><mo>/</mo><mo stretchy="false">(</mo><msub><mrow><mi>ρ</mi></mrow><mrow><mi mathvariant="normal">G</mi></mrow></msub><mo>+</mo><msub><mrow><mi>ρ</mi></mrow><mrow><mi mathvariant="normal">S</mi><mi mathvariant="normal">W</mi><mi mathvariant="normal">I</mi><mi mathvariant="normal">R</mi></mrow></msub><mo stretchy="false">)</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20187706&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20187701&type= 48.76800156 4.57200003 式中， <math id="M2"><msub><mrow><mi>ρ</mi></mrow><mrow><mi mathvariant="normal">G</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20187731&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20187728&type= 3.21733332 4.48733330 、 <math id="M3"><msub><mrow><mi>ρ</mi></mrow><mrow><mi mathvariant="normal">S</mi><mi mathvariant="normal">W</mi><mi mathvariant="normal">I</mi><mi mathvariant="normal">R</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20187738&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20187735&type= 6.60400009 4.48733330 分别为可见光G波段和SWIR波段的反射率。NDSI是常用的植被指数方法在冰川提取中的扩展和延伸。NDSI在冰川遥感提取中较为常用，如利用NDSI进行玉龙喀什河源区冰川边界提取（ 上官冬辉 等，2004 上官冬辉 等，2004 ），或利用NDSI来进行青藏高原东缘季风性海洋冰川的提取（ 宋波 等，2007 宋波 等，2007 ）。 另外，比值图像 RI （Ratio Image）方法也是一种常用的冰川提取方法，利用可见光G波段的反射率 <math id="M4"><msub><mrow><mi>ρ</mi></mrow><mrow><mi mathvariant="normal">G</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20187731&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20187728&type= 3.21733332 4.48733330 和SWIR波段的反射率 <math id="M5"><msub><mrow><mi>ρ</mi></mrow><mrow><mi mathvariant="normal">S</mi><mi mathvariant="normal">W</mi><mi mathvariant="normal">I</mi><mi mathvariant="normal">R</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20187738&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20187735&type= 6.60400009 4.48733330 的比值 <math id="M6"><mi>R</mi><mi>I</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20187750&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20187746&type= 3.55599999 2.62466669 来突显冰雪信息，其中， <math id="M7"><msub><mrow><mi>ρ</mi></mrow><mrow><mi mathvariant="normal">G</mi></mrow></msub><mo>∈</mo><mo stretchy="false">[</mo><mn mathvariant="normal">0,1</mn><mo stretchy="false">]</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20187757&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20187753&type= 16.08666801 4.48733330 ， <math id="M8"><msub><mrow><mi>ρ</mi></mrow><mrow><mi mathvariant="normal">S</mi><mi mathvariant="normal">W</mi><mi mathvariant="normal">I</mi><mi mathvariant="normal">R</mi></mrow></msub><mo>∈</mo><mo stretchy="false">[</mo><mn mathvariant="normal">0,1</mn><mo stretchy="false">]</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20187773&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20187763&type= 19.38866615 4.48733330 ，可得 <math id="M9"><mi mathvariant="normal">N</mi><mi mathvariant="normal">D</mi><mi mathvariant="normal">S</mi><mi mathvariant="normal">I</mi><mo>∈</mo><mo stretchy="false">[</mo><mo>-</mo><mn mathvariant="normal">1,1</mn><mo stretchy="false">]</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20187791&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20187785&type= 22.94466591 3.38666677 ，而 <math id="M10"><mi>R</mi><mi>I</mi><mo>∈</mo><mo stretchy="false">[</mo><mn mathvariant="normal">0</mn><mo>,</mo><mo>+</mo><mi mathvariant="normal">∞</mi><mo stretchy="false">)</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20187802&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20187795&type= 20.65866661 3.47133350 。在进行单一或均质冰川提取时，手动调节 RI 的阈值，能够实现冰川边界精确提取。然而， RI 的值域上界为无穷大，不利于利用模型或算法自动计算阈值，从而将冰川和背景分离（ 都伟冰，2014 都伟冰，2014 ）。NDSI值域有界，并且上下界呈对称分布，有利于采用数学模型实现快速的阈值分割和冰川提取（ 都伟冰，2014 都伟冰，2014 ； 李均力 等，2011 李均力 等，2011 ）。 2.2　基于NDSI的高斯混合模型分析 2 高斯混合模型GMM（Gaussian Mixture Model），以下简称GMM模型，由多个正态分布函数组成，故为高斯混合模型（ Zhu等，2012 Zhu等，2012 ， 2016 2016 ； 康一飞 等，2017 康一飞 等，2017 ）。其概率密度函数公式如下： <math display="block" id="M11"><mi>F</mi><mfenced separators="|"><mrow><mi>x</mi></mrow></mfenced><mo>=</mo><mstyle displaystyle="true"><munderover><mo largeop="true">∑</mo><mrow><mi>j</mi><mo>=</mo><mn mathvariant="normal">1</mn></mrow><mrow><mi>c</mi></mrow></munderover></mstyle><mrow><msub><mrow><mi>Φ</mi></mrow><mrow><mi>j</mi></mrow></msub></mrow><mi>f</mi><mfenced separators="|"><mrow><mi>x</mi><mo>;</mo><msub><mrow><mi>μ</mi></mrow><mrow><mi>i</mi></mrow></msub><mo>,</mo><mtext> </mtext><msub><mrow><mi>σ</mi></mrow><mrow><mi>i</mi></mrow></msub></mrow></mfenced></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20187814&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20187809&type= 37.76133347 8.72066593 式中， x 为向量 <math id="M12"><mi mathvariant="bold-italic">X</mi><mo>=</mo><mfenced open="{" close="}" separators="|"><mrow><msub><mrow><mi>x</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub><mo>,</mo><msub><mrow><mi>x</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msub><mo>,</mo><msub><mrow><mi>x</mi></mrow><mrow><mn mathvariant="normal">3</mn></mrow></msub><mo>,</mo><mo>…</mo><mo>,</mo><msub><mrow><mi>x</mi></mrow><mrow><mi>n</mi></mrow></msub></mrow></mfenced></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20187826&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20187822&type= 32.25799942 5.24933338 的 n 个样本值， X 包含 c 个类别，则第 j 类的权重为 <math id="M13"><msub><mrow><mi>Φ</mi></mrow><mrow><mi>j</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20187839&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20187832&type= 3.64066648 3.80999994 ，函数 f （ x ）服正态分布， <math id="M14"><msub><mrow><mi>μ</mi></mrow><mrow><mi>j</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20187850&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20187846&type= 2.62466669 4.48733330 、 <math id="M15"><msub><mrow><mi>σ</mi></mrow><mrow><mi>j</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20187863&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20187858&type= 3.04800010 3.72533321 为第 j 类的均值和方差。 本文中利用NDSI图像将冰川存在的网格区分为3类像元：非冰川像元、混合像元和冰川像元，即 公式（2） 公式（2） 中的 c =3，像元所属类别用 j 表示。由NDSI的定义可知，冰川像元的NDSI为较大（或是正值），非冰川像元NDSI为较小（或为负值），这里定义混合像元类NDSI分布位于冰川类NDSI和非冰川类NDSI之间。 图1 图1 以网格1为例，可见白色部分为冰川像元，NDSI值较大，红色部分为混合像元，其他部分为非冰川像元。混合像元主要分布在冰川类与非冰川类的分界线处， 图1 图1 中共约15万个像元，其中混合像元约2.5万个。 由 式（1） 式（1） 知，高斯混合模型密度函数里包含有参数集 <math id="M16"><mi mathvariant="bold-italic">θ</mi><mo>=</mo><mfenced open="{" close="}" separators="|"><mrow><mi>Φ</mi><mo>,</mo><mi>μ</mi><mo>,</mo><mi>σ</mi></mrow></mfenced></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20187892&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20187889&type= 19.98133469 4.48733330 ，定义此参数集为 <math id="M17"><mi mathvariant="bold-italic">θ</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188131&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188128&type= 1.86266661 2.70933342 。用极大似然法估计参数，就是在给定NDSI集合的情况下，寻找合适的参数 <math id="M18"><mi mathvariant="bold-italic">θ</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188131&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188128&type= 1.86266661 2.70933342 ，使似然函数最大。研究区域NDSI的集合是 <math id="M19"><mi mathvariant="bold-italic">X</mi><mo>=</mo><mfenced open="{" close="}" separators="|"><mrow><msub><mrow><mi>x</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub><mo>,</mo><msub><mrow><mi>x</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msub><mo>,</mo><msub><mrow><mi>x</mi></mrow><mrow><mn mathvariant="normal">3</mn></mrow></msub><mo>,</mo><mo>…</mo><mo>,</mo><msub><mrow><mi>x</mi></mrow><mrow><mi>n</mi></mrow></msub></mrow></mfenced></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20187929&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20187920&type= 32.25799942 5.24933338 ， <math id="M20"><msub><mrow><mi>x</mi></mrow><mrow><mi>i</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20187999&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20187994&type= 2.37066650 3.72533321 是从集合 <math id="M21"><mi mathvariant="bold-italic">X</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20187953&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20187948&type= 2.45533323 2.62466669 中抽取的，集合中元素相互独立，且符合高斯混合模型，未知类别标记为 <math id="M22"><msub><mrow><mi>z</mi></mrow><mrow><mi>i</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188065&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188062&type= 2.03200006 3.72533321 ， <math id="M23"><msub><mrow><mi>Q</mi></mrow><mrow><mi>i</mi></mrow></msub><mo stretchy="false">(</mo><msub><mrow><mi>z</mi></mrow><mrow><mi>i</mi></mrow></msub><mo>=</mo><mi>j</mi><mo stretchy="false">)</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20187973&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20187968&type= 14.39333248 4.57200003 表示NDSI的集合中 <math id="M24"><msub><mrow><mi>x</mi></mrow><mrow><mi>i</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20187999&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20187994&type= 2.37066650 3.72533321 属于 <math id="M25"><mi>j</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20187991&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20187984&type= 0.84666669 3.47133350 类的概率。对于每个NDSI集合中的 <math id="M26"><msub><mrow><mi>x</mi></mrow><mrow><mi>i</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20187999&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20187994&type= 2.37066650 3.72533321 ， <math id="M27"><msub><mrow><mi>z</mi></mrow><mrow><mi>i</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188065&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188062&type= 2.03200006 3.72533321 满足参数为 <math id="M28"><mi>Φ</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188019&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188012&type= 3.04800010 2.70933342 的多项式分布 <math id="M29"><msub><mrow><mi>Q</mi></mrow><mrow><mi>i</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188029&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188025&type= 3.21733332 4.48733330 ，见 式（4） 式（4） 。 <math display="block" id="M30"><msub><mrow><mi>Q</mi></mrow><mrow><mi>i</mi></mrow></msub><mo stretchy="false">(</mo><msub><mrow><mi>z</mi></mrow><mrow><mi>i</mi></mrow></msub><mo stretchy="false">)</mo><mo>=</mo><mi>P</mi><mfenced separators="|"><mrow><msub><mrow><mi>z</mi></mrow><mrow><mi>i</mi></mrow></msub><mo>=</mo><mi>j</mi><mo stretchy="false">|</mo><msub><mrow><mi>x</mi></mrow><mrow><mi>i</mi></mrow></msub><mo>;</mo><mi>Φ</mi><mo>,</mo><mtext> </mtext><mi>μ</mi><mo>,</mo><mtext> </mtext><mi>σ</mi></mrow></mfenced></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188040&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188036&type= 42.67200089 5.16466665 式中， <math id="M31"><mstyle displaystyle="true"><munderover><mo largeop="true">∑</mo><mrow><msub><mrow><mi>z</mi></mrow><mrow><mi>i</mi></mrow></msub><mo>=</mo><mn mathvariant="normal">1</mn></mrow><mrow><mn mathvariant="normal">3</mn></mrow></munderover></mstyle><mrow><msub><mrow><mi>Q</mi></mrow><mrow><mi>i</mi></mrow></msub><mfenced separators="|"><mrow><msub><mrow><mi>z</mi></mrow><mrow><mi>i</mi></mrow></msub></mrow></mfenced></mrow><mo>=</mo><mn mathvariant="normal">1</mn></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188048&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188043&type= 19.72733307 9.22866726 ， <math id="M32"><msub><mrow><mi>Q</mi></mrow><mrow><mi>i</mi></mrow></msub><mfenced separators="|"><mrow><msub><mrow><mi>z</mi></mrow><mrow><mi>i</mi></mrow></msub></mrow></mfenced><mo>≥</mo><mn mathvariant="normal">0</mn></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188058&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188053&type= 14.73200035 5.16466665 。为了确定每个样本所属类别的概率 <math id="M33"><msub><mrow><mi>z</mi></mrow><mrow><mi>i</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188065&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188062&type= 2.03200006 3.72533321 ，需使得 <math id="M34"><mi>P</mi><mfenced separators="|"><mrow><msub><mrow><mi>x</mi></mrow><mrow><mi>i</mi></mrow></msub><mo>,</mo><msub><mrow><mi>z</mi></mrow><mrow><mi>i</mi></mrow></msub></mrow></mfenced></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188076&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188069&type= 11.93799973 5.16466665 取得最大。可得变量 <math id="M35"><mi>H</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188085&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188082&type= 2.53999996 2.62466669 服从参数为 <math id="M36"><mi>θ</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188138&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188135&type= 1.77800000 2.70933342 的高斯分布的概率，即似然函数计算公式如下： <math display="block" id="M37"><mi>L</mi><mfenced separators="|"><mrow><mi>H</mi><mo stretchy="false">|</mo><mi mathvariant="bold-italic">θ</mi></mrow></mfenced><mo>=</mo><mstyle displaystyle="true"><munderover><mo>∏</mo><mrow><mi>i</mi><mo>=</mo><mn mathvariant="normal">1</mn></mrow><mrow><mi>n</mi></mrow></munderover></mstyle><mrow><mi>P</mi><mfenced separators="|"><mrow><msub><mrow><mi>x</mi></mrow><mrow><mi>i</mi></mrow></msub><mo>;</mo><mi mathvariant="bold-italic">θ</mi></mrow></mfenced></mrow></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188110&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188102&type= 32.34266663 8.29733276 式中，似然函数 <math id="M38"><mi>L</mi><mfenced separators="|"><mrow><mi>H</mi><mo stretchy="false">|</mo><mi mathvariant="bold-italic">θ</mi></mrow></mfenced></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188125&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188118&type= 10.32933331 3.80999994 和参数集 <math id="M39"><mi mathvariant="bold-italic">θ</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188131&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188128&type= 1.86266661 2.70933342 之间存在复杂非线性关系，不能直接最大化似然函数。因此，利用期望最大算法在概率模型中寻找极大似然估计的参数 <math id="M40"><mi>θ</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188138&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188135&type= 1.77800000 2.70933342 。期望最大算法是极大似然估计的常用方法，适合具有隐含未知变量的数据，能恰当估计出最优概率模型参数（ 颜凯，2009 颜凯，2009 ）。 期望最大算法经过两个步骤交替进行计算。根据经验为冰川像元、混合像元、非冰川像元（即隐藏变量）设置初始的参数集为 <math id="M41"><mi mathvariant="bold-italic">θ</mi><mo>=</mo><mfenced open="{" close="}" separators="|"><mrow><msub><mrow><mi>θ</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub><mo>,</mo><msub><mrow><mi>θ</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msub><mo>,</mo><msub><mrow><mi>θ</mi></mrow><mrow><mn mathvariant="normal">3</mn></mrow></msub></mrow></mfenced></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188150&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188144&type= 21.08200073 5.24933338 。（1）根据 <math id="M42"><mi mathvariant="bold-italic">θ</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188186&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188182&type= 1.86266661 2.70933342 求每个NDSI变量 <math id="M43"><msub><mrow><mi>x</mi></mrow><mrow><mi>i</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188192&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188189&type= 2.37066650 3.72533321 所属类别；（2）根据归类结果计算最大似然值 <math id="M44"><mi>L</mi><mfenced separators="|"><mrow><mi>H</mi><mo stretchy="false">|</mo><mi>θ</mi></mrow></mfenced></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188179&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188175&type= 10.24466610 3.80999994 所对应的参数集合 <math id="M45"><mi mathvariant="bold-italic">θ</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188186&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188182&type= 1.86266661 2.70933342 ，反过来重新计算每个NDSI变量 <math id="M46"><msub><mrow><mi>x</mi></mrow><mrow><mi>i</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188192&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188189&type= 2.37066650 3.72533321 所属类别。循环迭代，最终参数集最大似然值 <math id="M47"><mi>L</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188202&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188197&type= 1.86266661 2.62466669 趋于稳定，该算法执行完毕。（ 荆晓刚 等，2017 荆晓刚 等，2017 ； 石雪 等，2017 石雪 等，2017 ； 岳佳和王士同，2006 岳佳和王士同，2006 ； 张青 等，2018 张青 等，2018 ；左飞，2016） 当本文设定的3个正态分布组成的高斯混合模型趋于稳定时，输出GMM曲线。根据模拟出的GMM的每个正态分布，分析混合像元对区分冰川像元和非冰川像元的影响程度，对影响大的混合像元，通过去除混合像元类的NDSI，重新模拟GMM，进行算法改进，得到拟合较好的GMM。NDSI分为3个区间，非冰川像元的NDSI∈［-1， a ］，混合像元的NDSI∈（ a ， b ），冰川类的NDSI∈［ b ，1］， a ， b 分别为混合像元类阈值，但不是非冰川像元和冰川像元的阈值，冰川像元提取阈值T根据去除混合像元模拟出的GMM，构建求交方程式获得。 3　研究方法 1 3.1　实验数据准备和预处理 2 实验对不同时相，不同海拔的3个实验区的冰川进行分析，3个实验区的位置如 图2 图2 所示。 实验区域1、2、3分别选区 表1 表1 的Landsat 5 TM遥感影像，数据空间分辨率30 m，实验区所在影像几乎无云雪，影像质量好，便于精度对比分析。 参数 实验区1 实验区2 实验区3 隶属山系 天山 祁连山 羌塘高原 平均海拔/m 3940 4847 5895 传感器 Landsat 5 TM Landsat 5 TM Landsat 5 TM 时相 2007-08-25 2006-09-18 2000-09-28 行列号 P138R030 P135R033 P143R038 将全球尺度范围进行网格划分，是局部信息精确提取和保持时空分析一致性的有效前提。例如，降雨资料按照分析任务不同，包含全球0.25°、0.5°、1°、5°等不同等级的网格（ Cecil等，2014 Cecil等，2014 ）。臭氧数据包含全球1°×1°、1°×1.25°的网格产品（ 彭晓琳 等，2017 彭晓琳 等，2017 ）。Landsat ARD（Analysis Ready Data）产品，包含全球陆地150 km×150 km、189 km×180 km的网格产品（ Egorov等，2019 Egorov等，2019 ）。 由此，引出本文以0.125°×0.125°网格为局部区域，初步确定遥感影像中冰川的位置，将实验区或遥感影像分割为0.125°×0.125°网格，与冰川编目数据（ Liu等，2012 Liu等，2012 ）叠加分析，从而标记出有冰川存在的网格区域。 通过1°×1°，0.5°×0.5°，0.25°×0.25°，0.125°×0.125°网格划分实验，发现当网格大小为1°×1°时，每个网格中冰川占比过低，当网格大小为0.125°×0.125°时，大部分网格内冰川、非冰川的面积比例较为协调。在网格设计的时候，网格尺度要合理，网格过大，冰川初始占比过小，不能满足阈值需求；网格过小，分区太多，造成算法冗余度高。本文设置的网格大小为0.125°×0.125°，通过冰川编目与网格叠加，一景遥感影像中会存在冰川占比大于30%的网格 。 因此本文选用尺度为0.125°×0.125°大小的网格分割研究区域，如 图3 图3 中所示蓝色网格即为此尺度的网格。将中国西部网格化为0.125°×0.125°的网格，然后与冰川编目数据叠加分析，获取有冰川存在的网格，本文实验选取3个实验区中的两个网格进行实验，如 图2 图2 （a）、 图2 图2 （b）、 图2 图2 （c）所示的红色网格，每个实验区分别选用冰川类与非冰川类相对平衡的网格（冰川占比大于30%），和冰川类占比较小的网格（冰川占比小于30%），初始的冰川占比由冰川编目中该网格内冰川面积除以网格面积得到。 图3 图3 中所选区的网格由左往右依次为网格1、网格2、…、网格6，其中网格1、3、5为冰川占比相对平衡网格，而网格2、4、6为冰川占比相对较少的网格。 1,2,3; 3.2　实验数据分析和技术流程 2 （1）根据第2.1节中NDSI值区分冰川和背景信息的原理，根据局部遥感影像NDSI值的统计直方（图（4）中蓝色直方图），可见3类像元呈现3个正态分布的特点。 图4（a） 图4（a） 、 图4（c） 图4（c） 、 图4（e） 图4（e） 对应网格1、3、5的NDSI分布情况。每幅图种的蓝色曲线构成GMM曲线，由左往右，由非冰川像元、混合像元、冰川像元对应高斯分布曲线组合而成。可见，当冰川占比较平衡时，GMM模拟的3类像元的高斯曲线分布合理，冰川像元、非冰川像元频率高、离散度低，混合像元频率低、离散度高。而在 图4（b） 图4（b） 、 图4（d） 图4（d） 、 图4（f） 图4（f） ，GMM曲线模拟3类像元的蓝色曲线的高斯分布极不合理，非冰川像元高斯曲线和混合像元高斯曲线的均值相近，这样势必造成冰川像元提取阈值偏低或者不稳定。这也一定程度上解释了 表3 表3 中Ostu方法阈值计算偏低的原因。 注： 方框标记的数值，表示改进前后标准差大大降低，正文3.2节已详细叙述。 网格 改进前权重 改进后权重 改进前均值 改进后均值 改进前标准差 改进后标准差 1 0.30 0.37 0.86 0.86 0.03 0.03 2 0.01 0.01 0.75 0.84 0.14 0.03 3 0.31 0.34 0.78 0.78 0.04 0.04 4 0.01 0.01 0.74 0.79 0.15 0.05 5 0.35 0.39 0.84 0.84 0.03 0.03 6 0.07 0.06 0.81 0.87 0.14 0.03 网格 GMM法阈值 冰川面积占比 最大类间方差法阈值 网格 GMM法阈值 冰川面积占比 最大类间方差法阈值 网格 GMM法阈值 冰川面积占比 最大类间方差法阈值 1 0.52 0.010 0.16 10 0.36 0.005 0.24 19 0.43 0.017 0.25 2 0.52 0.047 0.24 11 0.30 0.028 0.26 20 0.54 0.044 0.27 3 0.56 0.021 0.23 12 0.30 0.053 -0.18 21 0.54 0.192 0.29 4 0.54 0.327 0.29 13 0.38 0.003 0.18 22 0.58 0.389 0.29 5 0.58 0.218 0.27 14 0.39 0.075 0.22 23 0.47 0.069 0.29 6 0.52 0.072 0.22 15 0.35 0.354 0.24 24 0.54 0.363 0.54 7 0.49 0.006 0.15 16 0.32 0.300 0.25 25 0.33 0.026 0.33 8 0.52 0.044 0.22 17 0.28 0.011 0.23 26 0.49 0.039 0.49 9 0.55 0.005 -0.18 18 0.36 0.045 0.26 27 0.53 0.046 0.53 与冰川面积占比的相关系数 0.36 — 0.48 与冰川面积占比的相关系数 0.03 — 0.18 与冰川面积占比的相关系数 0.60 — 0.65 总体阈值与冰川面积占比相关系数 0.19 — 0.32 — — — — — — — — （2）由此，需要将方差大、干扰强的混合像元去除，降低冰川像元NDSI方差，才能提高冰川像元阈值计算的稳定性和准确性。 （3）从统计直方图上，混合像元正态分布的特征从直观上并不明显，主要是因为混合像元方差大，在网格1、3、5中混合像元与其他两类像元区分明显（ 图4（a） 图4（a） ）。利用同景影像，例如 <math id="M48"><msub><mrow><mi>Φ</mi></mrow><mrow><mi mathvariant="normal">g</mi><mi mathvariant="normal">l</mi><mi mathvariant="normal">a</mi></mrow></msub><mo>></mo><mn mathvariant="normal">0.3</mn></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188282&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188279&type= 14.81666565 3.80999994 的网格（ 图4（a） 图4（a） ）的非冰川类正态分布方程与混合像元类正态分布方程求交，得到混合像元类的下界值 M a =0.0668，根据混合像元正态分布方程和冰川像元正态分布方程求交，得到混合像元类的上界值 M b =0.7669，同景影像中冰川类、非冰川类占比不平衡时（ 图4（b） 图4（b） ），根据以上计算出的混合像元阈值［0.0668，0.7669］，去除混合像元后，模拟出的仅包含冰川像元、非冰川像元2类的GMM更加合理，降低了冰川像元、非冰川像元NDSI方差，计算出网格1、2冰川提取阈值 T 分别为0.54、0.55。 同理，利用网格3混合像元阈值［ M a ， M b ］（ 图4（c） 图4（c） ），去除同景影像网格4混合像元（ 图4（d） 图4（d） ），大大提高了冰川像元、非冰川像元的纯净化程度。在 图4（f） 图4（f） 中表现最为明显，网格6去除混合像元之前，冰川像元方差极大（ 图4（f） 图4（f） 中右侧蓝色曲线），去除混合像元后GMM曲线模拟冰川像元（见 图4（f） 图4（f） 中右侧红色曲线），该红色曲线冰川像元部分与直方图匹配度极高。从侧面证明了，本文方法去除混合像元后，GMM曲线模拟冰川像元NDSI分布精度得到有效地提高。 改进前的GMM的参数值和改进后的GMM参数值如 表2 表2 中用方框标记的数值可见，去除混合像元后，GMM曲线模拟冰川像元NDSI分布的标准差大大降低，提高了GMM模拟冰川像元NDSI的准确程度，这样才能提高冰川提取阈值的准确性。 实验技术路线如 图5 图5 所示。通过大量实验分析可得，网格内冰川类、非冰川类NDSI的占比差异对GMM的估计影响较大，另外，混合像元类会干扰冰川类和非冰川类GMM模拟，采用如下方法解决。 1,2,3;4,5,6;7; （1）在同一景影像中，将网格内NDSI进行GMM模拟，该GMM中包含3类：非冰川类、混合像元类、冰川类。根据期望最大算法，计算每个包含冰川网格内这3类的NDSI占比，分别为 <math id="M49"><msub><mrow><mi>Φ</mi></mrow><mrow><mi mathvariant="normal">n</mi><mi mathvariant="normal">o</mi><mi mathvariant="normal">n</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188373&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188370&type= 6.01133299 3.80999994 、 <math id="M50"><msub><mrow><mi>Φ</mi></mrow><mrow><mi mathvariant="normal">m</mi><mi mathvariant="normal">i</mi><mi mathvariant="normal">x</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188380&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188376&type= 6.01133299 3.80999994 、 <math id="M51"><msub><mrow><mi>Φ</mi></mrow><mrow><mi mathvariant="normal">g</mi><mi mathvariant="normal">l</mi><mi mathvariant="normal">a</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188385&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188382&type= 5.50333309 3.80999994 。 （2）经过大量实验发现，当 <math id="M52"><msub><mrow><mi>Φ</mi></mrow><mrow><mi mathvariant="normal">g</mi><mi mathvariant="normal">l</mi><mi mathvariant="normal">a</mi></mrow></msub><mo>></mo><mn mathvariant="normal">0.3</mn></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188396&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188394&type= 14.64733315 3.80999994 ，该网格非冰川类和冰川类占比平衡，并且网格的混合像元在 <math id="M53"><msub><mrow><mi>Φ</mi></mrow><mrow><mi mathvariant="normal">g</mi><mi mathvariant="normal">l</mi><mi mathvariant="normal">a</mi></mrow></msub><mo>></mo><mn mathvariant="normal">0.3</mn></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188396&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188394&type= 14.64733315 3.80999994 的时候模拟几乎不失真，从 图3（a） 图3（a） 、 图3（c） 图3（c） 、 图3 图3 （e）看出冰川类占比平衡的网格，算法改进前冰川类和非冰川类、以及混合像元类均能很好的模拟出来。其NDSI分布 <math id="M54"><mi mathvariant="normal">N</mi><mi mathvariant="normal">D</mi><mi mathvariant="normal">S</mi><msub><mrow><mi mathvariant="normal">I</mi></mrow><mrow><mi mathvariant="normal">m</mi><mi mathvariant="normal">i</mi><mi mathvariant="normal">x</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188400&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188398&type= 11.00666618 3.80999994 具有较强的代表性，能够代表该景影像的混合像元分布情况。当 <math id="M55"><msub><mrow><mi>Φ</mi></mrow><mrow><mi mathvariant="normal">g</mi><mi mathvariant="normal">l</mi><mi mathvariant="normal">a</mi></mrow></msub><mo><</mo><mn mathvariant="normal">0.3</mn></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188407&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188403&type= 14.64733315 3.80999994 时，混合像元类和冰川类均模拟失真，冰川类因像元太少，几乎模拟不出来，因此需要根据同景影像中 <math id="M56"><msub><mrow><mi>Φ</mi></mrow><mrow><mi mathvariant="normal">g</mi><mi mathvariant="normal">l</mi><mi mathvariant="normal">a</mi></mrow></msub><mo>></mo><mn mathvariant="normal">0.3</mn></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188414&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188411&type= 14.64733315 3.80999994 的网格的GMM计算得到该景影像混合像元的NDSI分布区间为 <math id="M57"><mi mathvariant="normal">N</mi><mi mathvariant="normal">D</mi><mi mathvariant="normal">S</mi><msub><mrow><mi mathvariant="normal">I</mi></mrow><mrow><mi mathvariant="normal">m</mi><mi mathvariant="normal">i</mi><mi mathvariant="normal">x</mi></mrow></msub><mo>∈</mo><mo stretchy="false">(</mo><mi>a</mi><mo>,</mo><mi>b</mi><mo stretchy="false">)</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188425&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188419&type= 23.28333473 4.57200003 。 （3）利用 <math id="M58"><mi mathvariant="normal">N</mi><mi mathvariant="normal">D</mi><mi mathvariant="normal">S</mi><msub><mrow><mi mathvariant="normal">I</mi></mrow><mrow><mi mathvariant="normal">m</mi><mi mathvariant="normal">i</mi><mi mathvariant="normal">x</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188432&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188427&type= 11.00666618 3.80999994 去除掉该景影像中其他网格的混合像元，从而得到每个网格内仅包含冰川、非冰川两类的NDSI分布。 （4）返回第一步，将该景影像每个网格进行GMM模拟，该GMM模型仅包含非冰川类、冰川类，根据两类的正态分布构建求交方程式，计算出冰川提取的阈值 T 。 3.3　高斯混合模型计算冰川提取阈值 2 首先，选择一景影像中冰川像元与非冰川像元分布较平衡的网格，利用期望最大算法分为3类，得到GMM曲线，计算出混合像元类的NDSI分布 <math id="M59"><mi mathvariant="normal">N</mi><mi mathvariant="normal">D</mi><mi mathvariant="normal">S</mi><msub><mrow><mi mathvariant="normal">I</mi></mrow><mrow><mi mathvariant="normal">m</mi><mi mathvariant="normal">i</mi><mi mathvariant="normal">x</mi></mrow></msub><mo>∈</mo><mfenced open="[" close="]" separators="|"><mrow><msub><mrow><mi>M</mi></mrow><mrow><mi mathvariant="normal">a</mi></mrow></msub><mo>,</mo><msub><mrow><mi>M</mi></mrow><mrow><mi mathvariant="normal">b</mi></mrow></msub></mrow></mfenced></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188438&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188436&type= 28.61733437 5.24933338 ，作为同景影像的混合像元分布区间，在模拟同景影像的GMM时，将这些方差较大的混合像元去除，以便模拟出更加合理的GMM。 然后，设置GMM聚类数分为2类，模拟出仅包含冰川像元和非冰川像元的GMM。根据GMM的两个正态分布函数，求交计算出冰川边界的提取阈值。每个网格的NDSI直方图、GMM叠加后，可以清晰地看出GMM的拟合效果，见 图4 图4 中红色曲线为改进后的GMM曲线，计算的阈值 T 在 图4 图4 中已标记出。根据计算出的阈值对NDSI影像进行冰川提取及后处理，得到冰川边界，再根据水体指数（ Li等，2011 Li等，2011 ； 黎慧和赵阳，2018 黎慧和赵阳，2018 ）去除误分为冰川的冰湖， 图6（a） 图6（a） 、（b）中蓝色边界，分别为网格1、3提取的冰湖边界， 图6 图6 的底图为水体指数图。 1,2;3; 最终提取冰川边界，见 图7（a） 图7（a） 、 图7（c） 图7（c） 中的红色边界。用冰川编目数据作为验证数据，见 图7（b） 图7（b） 、 图7（d） 图7（d） 中的黄色边界。其中 图7（a） 图7（a） 、 图7（c） 图7（c） 的底图为NDSI图像， 图7（b） 图7（b） 、 图7（d） 图7（d） 的底图为Landsat 5 TM数据SWIR波段、近红外波段、红光波段合成的假彩色影像。 1,2,3,4;5; 3.4　实验结果分析 2 （1）本文研究方法在3个不同实验区所有冰川存在的网格进行实验，均得到了较好的提取效果，如 表3 表3 列出了不同实验区的冰川提取阈值。网格1—9属于天山实验区，网格10—18属于祁连山实验区，网格19—27属于羌塘高原实验区，每个实验区对应一景Landsat遥感影像（ 表1 表1 ），同一景影像中不同网格计算的冰川提取阈值差异较小，而时相不同、区域不同、海拔不同的网格，其阈值差异较大。3个实验区中羌塘高原实验区海拔 > 祁连山实验区海拔 > 天山实验区海拔；而对于平均冰川提取阈值，天山实验区阈值 > 羌塘高原实验区阈值 > 祁连山实验区阈值，经分析，羌塘高原实验区地类、地形分布复杂，计算的阈值差异较大。另外季节，雪、环境等也是影响阈值计算的重要因素，文中选取夏季，受云、雪干扰较少的时相，保证阈值计算的可比性。 （2）网格的大小选取的原则，就是一景影像中至少有一个网格，这个网格内冰川占比要均衡。本文选用0.125°×0.125°的网格作为局部区域，根据已有冰川编目快速寻找冰川存在网格。经多次实验，冰川在这个尺度的网格中占比较为平衡。如果网格尺度过大，冰川在网格中的占比过于小，造成冰川数据量过小，会因数据不平衡而导致算法不稳定现象。如果网格尺度过小，会增加工作量，不利于大面积的冰川快速提取。如果是小面积冰川提取，网格尺度，可以适当小些，但不宜过大。另外，如2.1节所述，经过归一化处理后的NDSI影像图，能很好的区分冰川和背景信息，因此，网格中的其他地物，不会对本文提出的方法造成影响。 （3）利用Ostu法与GMM方法计算的冰川提取阈值有差异（ 图8 图8 ）， 表3 表3 中对比Ostu法与GMM法的阈值可见，Ostu法计算的冰川提取阈值普遍偏低，许多非冰川区域被误分为冰川区域（ 图8（b） 图8（b） ），而利用GMM法计算的冰川阈值分割效果较好（ 图8（a） 图8（a） ）。如 图8 图8 中蓝色圆圈所标记的范围，最大类间方差法误提现象较为严重。 1,2;3; （4）目标（冰川像元）和背景（陆地像元）的比例大小，是引起阈值分割算法不稳定的重要因素之一。见 表3 表3 中，GMM法计算阈值与冰川占比进行相关性分析，得到3个试验区相关性系数分别为0.36、0.03、0.60；而Ostu法计算的阈值与冰川占比相关性系数分别为0.48、0.18、0.65。本文的GMM方法由于去除了混合像元干扰，受此因素干扰较小。 表3 表3 对比了不同网格Ostu法与GMM法计算的阈值。总体来看，Ostu法计算阈值偏小，误提现象较为严重。这也验证了在3.2节实验分析的论述，假设不去除混合像元，势必造成冰川像元NDSI的GMM曲线方差大，冰川提取阈值偏小。本文将阈值与冰川面积占比的相关性作为评价算法可靠性的指标。对于某一实验区的网格，当冰川提取阈值与冰川面积占比相关性较小时，说明该方法受目标地物占比影响小，算法稳定性高。从 表3 表3 中总体来看，GMM算法阈值、Ostu法阈值与冰川占比的相关性系数分别为0.19和0.32。可得结论，本文改进GMM算法从一定程度上削弱了目标占背景比例较小的干扰。 （5）以实验区1新疆哈尔里克山脉的冰川为验证实验区，通过改进的GMM提取的冰川边界与冰川编目数据对比来验证提取精度。统计提取的冰川面积和冰川编目的面积，进行对比，8个网格的面积精度最低为94.59%，平均面积精度为96.45%。对8个网格的冰川与非冰川区域进行精度评定，总体精度均达95.42%以上，平均值为97.21%，Kappa系数均在0.91以上，平均值为0.93。因为Landsat TM5数据空间分辨率为30 m，随机抽取各个区域冰川末端拐点，使其符合均匀分布，通过提取坐标值与冰川编目的位置进行计算和比较，各个区域的末端位置精度小于15 m（小于影像分辨率30 m），满足提取精度要求。利用本文改进的GMM算法提取定位观测冰川庙儿沟冰帽（冰川目录编号：5Y822C03）边界，与对应冰川编目进行精度比较，总体面积精度达96.57%，Kappa系数0.94，末端精度15.51 m（小于影像分辨率30 m）。因此得出结论，改进的GMM方法提取冰川边界与冰川编目的边界吻合度较高，有一定的应用价值。 4　结论 1 本文将GMM应用于遥感冰川提取过程中，其主要结论和展望总结如下。 （1）选取3个算法实验区，对不同时相、不同高度的冰川进行高斯混合模型计算阈值，结果表明，同一实验区不同网格计算的阈值是差异较小，而不同实验区不同网格，所计算的冰川的阈值差异稍大，这也说明了该算法在差异较大区域仍具有一定的稳定性。 （2）选取哈密地区为验证实验区，利用GMM阈值计算方法，得到该景影像8个网格的冰川提取阈值差异不大，得到的提取边界与冰川编目数据经对比验证，结果可靠，精度高，说明该算法具有一定的应用价值。 （3）本文通过初始化GMM模型为3个类别（冰川像元、混合像元、非冰川像元），从而确定出每景影像中混合像元类的NDSI分布区间 <math id="M60"><mi mathvariant="normal">N</mi><mi mathvariant="normal">D</mi><mi mathvariant="normal">S</mi><msub><mrow><mi mathvariant="normal">I</mi></mrow><mrow><mi>m</mi><mi>i</mi><mi>x</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188565&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=20188563&type= 11.00666618 3.80999994 ，并将其去除，再利用GMM模型模拟纯净化后的冰川类、非冰川类的NDSI分布情况，最大程度保留了有效的混合像元，并排除了无效混合像元对于阈值计算过程的干扰，这种改进后的GMM模型计算的冰川提取阈值稳定，且准确度高，该方法具有一定的创新性。 （4）本文提出了一种快速有效的GMM冰川阈值计算方法，在一定程度上解决了GMM曲线模拟过程中冰川类、非冰川类占比不平衡的问题，但该算法仍有很大的改进空间，例如将占比较小类的直方图进行放大，或将占比过大类的直方图缩小，或者其他能将不同类别的直方图进行均衡化的算法，从而进一步提高GMM冰川阈值计算方法的稳定性和应用范围。 参考文献（References） Bolch T , Kulkarni A , Kääb A , Huggel C , Paul F , Cogley J G , Frey H , Kargel J S , Fujita K , Scheel M , Bajracharya S and Stoffel M . 2012 . The state and fate of Himalayan glaciers . Science , 336 ( 6079 ): 310 - 314 [ DOI: 10.1126/science.1215828 http://dx.doi.org/10.1126/science.1215828 ] The state and fate of Himalayan glaciers Cecil D J , Buechler D E and Blakeslee R J . 2014 . Gridded lightning climatology from TRMM-LIS and OTD: dataset description . Atmospheric Research , 135 - 136 : 404 - 414 [ DOI: 10.1016/j.atmosres.2012.06.028 http://dx.doi.org/10.1016/j.atmosres.2012.06.028 ] Gridded lightning climatology from TRMM-LIS and OTD: dataset description Du W B . 2014 . Research on glacier extraction based on multisource RS imagery and changes monitoring of eastern Tienshan mountains . Jiaozuo : Henan Polytechnic University , 22 - 25 Research on glacier extraction based on multisource RS imagery and changes monitoring of eastern Tienshan mountains 都伟冰 . 2014 . 冰川信息多源遥感提取及天山东段典型冰川变化监测研究 . 焦作 : 河南理工大学 , 22 - 25 冰川信息多源遥感提取及天山东段典型冰川变化监测研究 Du W B , Li J L , Bao A M and Wang B S . 2014 . Mapping changes in the glaciers of the eastern Tienshan mountains during 1977-2013 using multitemporal remote sensing . Journal of Applied Remote Sensing , 8 ( 1 ): 084683 [ DOI: 10.1117/1.JRS.8.084683 http://dx.doi.org/10.1117/1.JRS.8.084683 ] Mapping changes in the glaciers of the eastern Tienshan mountains during 1977-2013 using multitemporal remote sensing Du W B , Li J L , Bao A M , Wang B S and Wang S T . 2015 . Information extraction method of alpine glaciers with multitemporal and multiangle remote sensing . Acta Geodaetica et Cartographica Sinica , 44 ( 1 ): 59 - 66 Information extraction method of alpine glaciers with multitemporal and multiangle remote sensing 都伟冰 , 李均力 , 包安明 , 王宝山 , 王双亭 . 2015 . 高山冰川多时相多角度遥感信息提取方法 . 测绘学报 , 44 ( 1 ): 59 - 66 高山冰川多时相多角度遥感信息提取方法 Egorov A V , Roy D P , Zhang H K , Li Z B , Yan L , and Huang H Y . 2019 . Landsat 4,5 and 7 1982 to 2017) Analysis Ready Data (ARD) observation coverage over the conterminous United States and implications for terrestrial monitoring . Remote Sensing , 11 ( 4 ): 447 ( DOI: 10.3390/rs11040447 http://dx.doi.org/10.3390/rs11040447 ] Landsat 4,5 and 7 1982 to 2017) Analysis Ready Data (ARD) observation coverage over the conterminous United States and implications for terrestrial monitoring Feng T , Liu S Y , Xu J L , Guo W Q , Wei J F and Zhang Z . 2015 . Glacier change of the Yarkant river basin from 1968 to 2009 derived from the first and second glacier inventories of China . Journal of Glaciology and Geocryology , 37 ( 1 ): 1 - 13 Glacier change of the Yarkant river basin from 1968 to 2009 derived from the first and second glacier inventories of China 冯童 , 刘时银 , 许君利 , 郭万钦 , 魏俊锋 , 张震 . 2015 . 1968-2009年叶尔羌河流域冰川变化——基于第一、二次中国冰川编目数据 . 冰川冻土 , 37 ( 1 ): 1 - 13 1968-2009年叶尔羌河流域冰川变化——基于第一、二次中国冰川编目数据 Guo W Q , Liu S Y , Xu J L , Wu L Z , Shangguan D H , Yao X J , Wei J F , Bao W J , Yu P C , Liu Q and Jiang Z L . 2015 . The second Chinese glacier inventory: data, methods and results . Journal of Glaciology , 61 ( 226 ): 357 - 372 [ DOI: 10.3189/2015jog14j209 http://dx.doi.org/10.3189/2015jog14j209 ] The second Chinese glacier inventory: data, methods and results Haeberli W and Hoelzle M . 1995 . Application of inventory data for estimating characteristics of and regional climate-change effects on mountain glaciers: a pilot study with the European Alps . Annals of Glaciology , 21 : 206 - 212 [ DOI: 10.3189/S0260305500015834 http://dx.doi.org/10.3189/S0260305500015834 ] Application of inventory data for estimating characteristics of and regional climate-change effects on mountain glaciers: a pilot study with the European Alps Jing X G , Ge L G and Sun W . 2017 . A clustering algorithm for sea buoy trajectory based on Gaussian mixture model . Modern Computer , ( 36 ): 3 - 5, 8 A clustering algorithm for sea buoy trajectory based on Gaussian mixture model 荆晓刚 , 葛丽阁 , 孙伟 . 2017 . 一种基于高斯混合模型的海上浮标轨迹聚类算法. 现代计算机, (36): 3 - 5 , 8 [ DOI: 10.3969/j.issn.1007-1423.2018.36.001 http://dx.doi.org/10.3969/j.issn.1007-1423.2018.36.001 ] 一种基于高斯混合模型的海上浮标轨迹聚类算法. 现代计算机, (36): Kang Y F , Pan L , Sun M W , Chen Q and Wang Y . 2017 . Gaussian mixture model based cloud detection for Chinese high resolution satellite imagery . Geomatics and Information Science of Wuhan University , 42 ( 6 ): 782 - 788 Gaussian mixture model based cloud detection for Chinese high resolution satellite imagery 康一飞 , 潘励 , 孙明伟 , 陈奇 , 王越 . 2017 . 基于高斯混合模型法的国产高分辨率卫星影像云检测 . 武汉大学学报(信息科学版) , 42 ( 6 ): 782- 788 [ DOI: 10.13203/j.whugis20140875 http://dx.doi.org/10.13203/j.whugis20140875 ] 基于高斯混合模型法的国产高分辨率卫星影像云检测 Li H and Zhao Y . 2018 . Comparative study on water extraction methods based on tm image . Popular Science , ( 9 ): 195 - 196 Comparative study on water extraction methods based on tm image 黎慧 , 赵阳 . 2018 . 基于TM影像水体提取方法比较研究. 科学大众·科学教育, (9): 195 - 196 [DOI: 10.3969/j.issn.1006-3315.2018.09.182] 基于TM影像水体提取方法比较研究. 科学大众·科学教育, (9): Li J . 2014 . The research of glacier variation and climate change in eastern Tianshan (1990-2011). Lanzhou: Lanzhou University , 11 - 12 . The research of glacier variation and climate change in eastern Tianshan (1990-2011). Lanzhou: 李佳 . 2014 . 东天山冰川与气候变化研究(1990-2011) . 兰州 : 兰州大学 , 11 - 12 东天山冰川与气候变化研究(1990-2011) Li J L , Sheng Y W and Luo J C . 2011 . Automatic extraction of Himalayan glacial lakes with remote sensing . Journal of Remote Sensing , 15 ( 1 ): 29 - 43 Automatic extraction of Himalayan glacial lakes with remote sensing 李均力 , 盛永伟 , 骆剑承 . 2011 . 喜马拉雅山地区冰湖信息的遥感自动化提取 . 遥感学报 , 15 ( 1 ): 29 - 43 喜马拉雅山地区冰湖信息的遥感自动化提取 Liu S Y , Guo W Q , Xu J L , Shang guan D H , Wu L Z , Yao X J , Zhao J D , Liu Q , Jiang Z L , Li P , Wei J F , Bao W J , Yu P C , Ding L F , Li G , Ge C M and Wang Y . 2012 . The second glacier inventory dataset of China (Version 1.0)(2006-2011) . National Tibetan Plateau Data Center [ DOI: 10.3972/glacier.001.2013.db http://dx.doi.org/10.3972/glacier.001.2013.db ] The second glacier inventory dataset of China (Version 1.0)(2006-2011) Liu S Y , Yao X J , Guo W Q , Xu J L , Shangguan D H , Wei J F , Bao W J and Wu L Z . 2015 . The contemporary glaciers in china based on the second Chinese glacier inventory . Acta Geographica Sinica , 70 ( 1 ): 3 - 16 The contemporary glaciers in china based on the second Chinese glacier inventory 刘时银 , 姚晓军 , 郭万钦 , 许君利 , 上官冬辉 , 魏俊峰 , 鲍伟佳 , 吴立宗 . 2015 . 基于第二次冰川编目的中国冰川现状 . 地理学报 , 70 ( 1 ): 3- 16 [ DOI: 10.11821/dlxb201501001 http://dx.doi.org/10.11821/dlxb201501001 ] 基于第二次冰川编目的中国冰川现状 Peng X L , Shen H F , Zeng C , Zhang L P and He Z Y . 2017 . A method for the recovery of aura satellite remote sensing ozone products . Geomatics and Information Science of Wuhan University , 42 ( 6 ): 789 - 796 A method for the recovery of aura satellite remote sensing ozone products 彭晓琳 , 沈焕锋 , 曾超 , 张良培 , 何宗宜 . 2017 . 一种aura卫星遥感臭氧产品的修复方法 . 武汉大学学报(信息科学版) , 42 ( 6 ): 789- 796 [ DOI: 10.13203/j.whugis20150390 http://dx.doi.org/10.13203/j.whugis20150390 ] 一种aura卫星遥感臭氧产品的修复方法 Scherler D , Bookhagen B and Strecker M R . 2011 . Spatially variable response of Himalayan glaciers to climate change affected by debris cover . Nature Geoscience , 4 ( 3 ): 156 - 159 [ DOI: 10.1038/ngeo1068 http://dx.doi.org/10.1038/ngeo1068 ] Spatially variable response of Himalayan glaciers to climate change affected by debris cover Shangguan D H , Liu S Y , Ding Y J , Ding L F and Li G . 2004 . Glacier changes at the head of Yurungkax river in the west Kunlun mountains in the past 32 years . Acta Geographica Sinica , 59 ( 6 ): 855 - 862 Glacier changes at the head of Yurungkax river in the west Kunlun mountains in the past 32 years 上官冬辉 , 刘时银 , 丁永建 , 丁良福 , 李刚 . 2004 . 玉龙喀什河源区32年来冰川变化遥感监测 . 地理学报 , 59 ( 6 ): 855- 862 [ DOI: 10.3321/j.issn:0375-5444.2004.06.007 http://dx.doi.org/10.3321/j.issn:0375-5444.2004.06.007 ] 玉龙喀什河源区32年来冰川变化遥感监测 Shi X , Li Y , Li X L and Zhao Q H . 2017 . Gaussian mixture model with neighbor relationship for image segmentation and simplified solving method . Journal of Image and Graphics , 22 ( 12 ): 1758 - 1768 Gaussian mixture model with neighbor relationship for image segmentation and simplified solving method 石雪 , 李玉 , 李晓丽 , 赵泉华 . 2017 . 融入邻域作用的高斯混合分割模型及简化求解 . 中国图象图形学报 , 22 ( 12 ): 1758- 1768 [ DOI: 10.11834/jig.170255 http://dx.doi.org/10.11834/jig.170255 ] 融入邻域作用的高斯混合分割模型及简化求解 Shi Y F . 2001 . Estimation of the water resources affected by climatic warming and glacier shrinkage before 2050 in west China . Journal of Glaciology and Geocryology , 23 ( 4 ): 333 - 341 Estimation of the water resources affected by climatic warming and glacier shrinkage before 2050 in west China 施雅风 . 2001 . 2050年前气候变暖冰川萎缩对水资源影响情景预估 . 冰川冻土 , 23 ( 4 ): 333- 341 [ DOI: 10.3969/j.issn.1000-0240.2001.04.001 http://dx.doi.org/10.3969/j.issn.1000-0240.2001.04.001 ] 2050年前气候变暖冰川萎缩对水资源影响情景预估 Song B , He Y Q , Pang H X , Lu A G , Zhao J D , Ning B Y , Yuan L L and Zhang Z L . 2007 . Identifying automatically the debris-covered glaciers in China’s monsoonal temperate-glacier regions based on remote sensing and GIS . Journal of Glaciology and Geocryology , 29 ( 3 ): 456 - 462 Identifying automatically the debris-covered glaciers in China’s monsoonal temperate-glacier regions based on remote sensing and GIS 宋波 , 何元庆 , 庞洪喜 , 卢爱刚 , 赵井东 , 宁宝英 , 院玲玲 , 张忠林 . 2007 . 基于遥感和GIS的我国季风海洋型冰川区冰碛物覆盖型冰川边界的自动识别 . 冰川冻土 , 29 ( 3 ): 456- 462 [ DOI: 10.3969/j.issn.1000-0240.2007.03.018 http://dx.doi.org/10.3969/j.issn.1000-0240.2007.03.018 ] 基于遥感和GIS的我国季风海洋型冰川区冰碛物覆盖型冰川边界的自动识别 Yan K . 2009 . Study of speaker recognition algorithm based on Gaussian mixture model . Nanjing : Nanjing University of Science and Technology Study of speaker recognition algorithm based on Gaussian mixture model 颜凯 . 2009 . 基于高斯混合模型的说话人识别算法研究 . 南京 : 南京理工大学 基于高斯混合模型的说话人识别算法研究 Yan L L and Wang J . 2013 . Study of extracting glacier information from remote sensing . Journal of Glaciology and Geocryology , 35 ( 1 ): 110 - 118 Study of extracting glacier information from remote sensing 彦立利 , 王建 . 2013 . 基于遥感的冰川信息提取方法研究进展 . 冰川冻土 , 35 ( 1 ): 110 - 118 基于遥感的冰川信息提取方法研究进展 Yue J and Wang S T . 2006 . Algorithm EM and its initialization in Gaussian-mixture-model based clustering . Control and Automation , 22 ( 33 ): 244 - 246, 302 Algorithm EM and its initialization in Gaussian-mixture-model based clustering 岳佳 , 王士同 . 2006 . 高斯混合模型聚类中EM算法及初始化的研究 . 微计算机信息 , 22 ( 33 ): 244 - 246 , 302 [ DOI: 10.3969/j.issn.1008-0570.2006.33.086 http://dx.doi.org/10.3969/j.issn.1008-0570.2006.33.086 ] 高斯混合模型聚类中EM算法及初始化的研究 Zhang Q , Zou X J , Lin G C , Wang B and Peng H X . 2018 . Precise identification and positioning algorithm for work piece based on Gaussian mixture model . Computer Integrated Manufacturing Systems Precise identification and positioning algorithm for work piece based on Gaussian mixture model 张青 , 邹湘军 , 林桂潮 , 王波 , 彭红星 . 2018 . 基于高斯混合模型的零件精确识别与定位算法 . 计算机集成制造系统 基于高斯混合模型的零件精确识别与定位算法 Zhang X B , Zhao X S , Ge D Q , Liu B , Zhang L , Li M and Wang Y . 2018 . Monitoring displacement of Laohugou glacier No . 12 based on Landsat 8 and TerraSAR-X images . Journal of Remote Sensing , 22 ( 1 ): 153 - 160 Monitoring displacement of Laohugou glacier No 张晓博 , 赵学胜 , 葛大庆 , 刘斌 , 张玲 , 李曼 , 王艳 . 2018 . 利用Landsat8和TerraSAR-X影像研究老虎沟12号冰川运动特征 . 遥感学报 , 22 ( 1 ): 153- 160 [ DOI: 10.11834/jrs.20186313 http://dx.doi.org/10.11834/jrs.20186313 ] 利用Landsat8和TerraSAR-X影像研究老虎沟12号冰川运动特征 Zhao J Y and Peng J H . 2016 . Information extraction technologies of mountainous glacier based on Landsat TM image . Journal of Huazhong Normal University (Natural Sciences) , 50 ( 2 ): 309 - 313 Information extraction technologies of mountainous glacier based on Landsat TM image 赵健赟 , 彭军还 . 2016 . 基于Landsat TM影像的山地冰川信息提取技术研究 . 华中师范大学学报(自然科学版) , 50 ( 2 ): 309- 313 [ DOI: 10.3969/j.issn.1000-1190.2016.02.025 http://dx.doi.org/10.3969/j.issn.1000-1190.2016.02.025 ] 基于Landsat TM影像的山地冰川信息提取技术研究 Zhu L L , Guo X , Ruan H H and Lu J . 2016 . Prediction of mechanical properties in bimodal nanotwinned metals with a composite structure . Composites Science and Technology , 123 : 222 - 231 [ DOI: 10.1016/j.compscitech.2015.12.013 http://dx.doi.org/10.1016/j.compscitech.2015.12.013 ] Prediction of mechanical properties in bimodal nanotwinned metals with a composite structure Zhu L L , Shi S Q , Lu K and Lu J . 2012 . A statistical model for predicting the mechanical properties of nanostructured metals with bimodal grain size distribution . Acta Materialia , 60 ( 16 ): 5762 - 5772 [ DOI: 10.1016/j.actamat.2012.06.059 http://dx.doi.org/10.1016/j.actamat.2012.06.059 ] A statistical model for predicting the mechanical properties of nanostructured metals with bimodal grain size distribution Zuo F . 2016 . R Language Combat-Machine Learning and Data Analysis . Beijing: Publishing House of Electronics Industry , 314 - 325 R Language Combat-Machine Learning and Data Analysis 左飞 . 2016 . R语言实战: 机器学习与数据分析 . 北京: 电子工业出版社 , 314 - 325 R语言实战: 机器学习与数据分析 感谢国家地球系统科学数据平台寒区旱区科学数据中心提供的第二次冰川编目数据。 10.11834/jrs.20219153 1007-4619（2021）07-1434-11 王亚利，都伟冰，王双亭.2021.高斯混合模型自动阈值法遥感冰川信息提取.遥感学报，25（7）： 1434-1444 王亚利，都伟冰，王双亭.2021.高斯混合模型自动阈值法遥感冰川信息提取.遥感学报，25（7）： 1434-1444 Wang Y L，Du W B and Wang S T. 2021. Extracting glacier information from remote sensing imageries by automatic threshold method of Gaussian mixture model. National Remote Sensing Bulletin， 25（7）：1434-1444 Wang Y L，Du W B and Wang S T. 2021. Extracting glacier information from remote sensing imageries by automatic threshold method of Gaussian mixture model. National Remote Sensing Bulletin， 25（7）：1434-1444 国家自然科学基金（编号：41601364，41601392）；河南理工大学博士基金（编号：B2016-11）；智慧中原地理信息技术协同创新中心开放课题项目（编号：2016A002）；河南省科技攻关项目（编号：172102210280）；河南理工大学创新型科研团队资助（编号：T2018-4） Supported by National Natural Science Foundation of China (No. 41601364,41601392); Doctoral Foundation Project of Henan Polytechnic university(No.B2016-11); Intelligent Central Plains Geographic Information Technology Collaborative Innovation Center Open Project（No.2016A002）; Key Technology Research and Development Program of Henan Province (No.172102210280); Innovative Research Team of Henan University of Technology (No. T2018-4) 2019-05-24 2021-07-07 2021-09-18 2019-09-28 遥感学报 7 25