一种积雪中黑碳浓度的遥感估算方法
【专利摘要】本发明涉及一种积雪中黑碳浓度的遥感估算方法，利用500米空间分辨率的MODIS每日积雪合成产品（MOD10A1和MYD10A1），结合地面站点实测黑碳浓度数据，进行区域尺度黑碳浓度估算。首先对MODIS每日积雪产品进行影像拼接、地图投影转换和数据格式转换，获取研究区域积雪反照率数据；然后对影像进行去云合成处理，提取地面站点位置的积雪反照率，并参照相应文献，建立积雪反照率和黑碳浓度的关系模型，接着对关系模型进行检验和反复优化；最后利用关系模型处理积雪反照率影像，利用研究区域矢量边界数据截取，得到积雪中黑碳浓度的时空分布结果。与传统方法相比，本发明引入遥感技术手段，可高效获取大范围积雪中黑碳浓度分布的时空特征数据。
【专利说明】
一种积雪中黑碳浓度的遥感估算方法
技术领域
[0001]本发明涉及一种积雪中黑碳浓度的遥感估算方法，属于大气环境遥感应用技术领域。【背景技术】
[0002]黑碳(Black Carbon)气溶胶是一种悬浮于空气中的黑色碳质粒颗粒物，是大气细颗粒物的重要组成部分，对大气环境以及人体健康都有着非常重要的影响。积雪中的黑碳颗粒能显著降低积雪反照率，增强对太阳辐射的吸收，加速积雪融化，造成显著的气候效应和水文循环影响。因此进行冰雪中黑碳含量的研究，获取其时空分布特征，有利于增强对气候变化及区域水循环的认识。
[0003]目前国内外对积雪中黑碳含量进行了较多的观测研究，其研究手段主要是在地面建立观测站点，采集雪样并进行相应测量，从而获取该位置积雪中的黑碳含量。这种方法虽然精确度高，但是费工费时，难以做到大范围积雪中黑碳含量的快速获取和更新。而大范围、实时、动态正是遥感技术的特点。考虑到这一特点，发明人考虑基于遥感技术手段，利用 M0DIS积雪数据产品，建立积雪反照率和黑碳浓度的关系模型，从而实现区域尺度黑碳浓度的快速估算。
【发明内容】
[0004]本发明要解决的技术问题是:克服现有技术手段的缺点，建立一种基于遥感光学影像数据和已有地面测站数据的关系模型，利用M0DIS积雪反照率数据，能够较准确、快速地获取积雪中黑碳含量及空间分布特征，达到比较好的效果。
[0005]为了解决上述技术问题，本发明提出的技术方案是:一种积雪中黑碳浓度的遥感估算方法，包括以下步骤:第一步、利用地面站点实测积雪中黑碳浓度数据，根据地面站点地理位置及测量时间， 下载对应时间和区域的M0DIS每日积雪合成产品M0D10A1和MYD10A1;第二步、对MOD IS影像数据产品MOD 10A1和MYD10A1分别进行预处理，包括影像拼接、地图投影转换和格式转换，并提取产品中的积雪反照率数据；第三步、对产品MOD 10A1和产品MYD 10A1获取的积雪反照率数据进行去云合成处理，获得积雪反照率影像；第四步、依据地面站点地理坐标提取积雪反照率影像上的对应位置积雪反照率； 第五步、参照文献资料，分析地面站点积雪中的黑碳浓度情况，选取合适站点，初步建立研究区域积雪反照率与黑碳浓度的关系模型；第六步、进行积雪反照率影像向黑碳浓度空间数据的转换，并对已建立的关系模型及初步得到的黑碳浓度结果，结合样本数据和相关文献资料进行检验评定，反复优化，直到确立最终的关系模型；第七步、通过研究区域矢量边界数据，根据最终的关系模型提取出研究区域积雪中黑碳浓度的时空分布结果。
[0006]本方法利用500米空间分辨率的M0DIS每日积雪合成产品M0D10A1和MYD10A1，结合地面站点实测黑碳浓度数据，进行区域尺度黑碳浓度估算。首先对M0DIS每日积雪产品进行影像拼接、地图投影转换和数据格式转换，获取研究区域积雪反照率数据;然后对影像进行去云合成处理，提取地面站点位置的积雪反照率，并参照相应文献，建立积雪反照率和黑碳浓度的关系模型，接着对关系模型进行检验和反复优化;最后利用关系模型处理积雪反照率影像，利用研究区域矢量边界数据截取，得到积雪中黑碳浓度的时空分布结果。与传统方法相比，本发明引入遥感技术手段，可高效获取大范围积雪中黑碳含量分布的时空特征数据。
[0007]本发明积雪中黑碳浓度的遥感估算方法，还具有如下改进:1、所述第二步中，影像拼接的方法是，根据研究区范围，将同一时间相应的几景影像拼接为一景完整的影像。
[0008]2、所述第三步中，所述去云合成处理的算法如下:a、对于某一坐标的像元，如果产品MOD 10A1对应的反照率为积雪反照率，则把产品 M0D10A1的该坐标像元作为合成后积雪反照率影像的像元；b、对于某一坐标的像元，如果产品M0D10A1对应的反照率为非积雪反照率，而MYD10A1 对应的反照率为积雪反照率，则把产品MYD10A1的该坐标像元作为合成后积雪反照率影像的像元；c、对于某一坐标的像元，如果产品M0D10A1和产品MYD10A1对应的反照率均非积雪反照率，则合成后积雪反照率影像在该坐标的像元为非积雪反照率像元。
[0009]3、所述第四步中，根据地理位置提取积雪反照率，有两种方法可供选择:地面站点较少时，在ENVI软件主窗口中，利用pixel locator提取相应坐标位置的反照率;地面站点较多时，利用ENVI+IDL语言，通过索引到影像上的行列号，批量获取地面站点对应影像上的反照率。
[0010]4、所述第五步中，为了避免所选样本中的异常值对最终的拟合结果造成较大的影响，应参考其他已有的研究成果对地面站点进行筛选。
[0011]5、所述第五步中，采用最小二乘法对积雪反照率和黑碳浓度数据进行线性拟合， 从而建立积雪反照率和黑碳浓度的关系模型。
[0012]本发明利用地面站点积雪中黑碳浓度数据，以及M0DIS积雪反照率数据产品，数据获取比较容易。影像数据的预处理、站点积雪反照率的提取、利用关系模型处理影像等步骤可利用软件MRT、ENVI及MATLAB处理，简单易懂、操作方便，大大地减少了工作量，提高了处理的效率，增强了本发明的灵活性和可用性。针对地面站点采样情况以及该地区黑碳含量的基本状况，对数据进行了筛选，使建立的关系模型更加符合该地区的自然情况。同时结合 M0DIS积雪反照率产品合成算法，本发明可高效、快速、准确地获取较大尺度地区积雪黑碳浓度。
[0013]综上，本发明方法的执行步骤简单易行，效果较好，解决了黑碳含量研究中样本采集、测量过程中费时费力的问题，有利于促进对黑碳含量分布变化、黑碳对积雪消融的影响、乃至黑碳的气候效应等方面的研究。【附图说明】
[0014]下面结合附图对本发明作进一步说明。
[0015]图1是本发明一种积雪中黑碳浓度的遥感估算方法流程图。
[0016]图2是使用MRT软件对M0DIS积雪产品进行预处理示意图。
[0017]图3是使用ENVI软件提取地面站点积雪反照率示意图。
[0018]图4是使用MATLAB软件根据关系模型处理影像示意图。
[0019]图5是建立的积雪反照率和黑碳浓度的关系模型。
[0020]图6是研究区域北疆地区黑碳浓度结果图。【具体实施方式】[0021 ]下面根据附图详细阐述本发明，使本发明的目的和效果变得更加明显。
[0022]如图1所示，为本发明一种积雪中黑碳浓度的遥感估算方法流程图，具体步骤如下:第一步、利用已有文献(Ye等人2012年在ENVIRONMENTAL RESEARCH LETTERS中的论文 “Black carbon in seasonal snow across northern Xinjiang in northwestern China”)给出的地面站点实测积雪中黑碳浓度数据，根据站点地理位置及测量时间，下载对应时间和区域的MOD IS每日积雪合成产品(MOD 10A1和MYD10A1)。
[0023]本步骤中，已有文献给出的地面站点共有32个，研究区为新疆北疆地区，测量时间为2012年1月28日至2月27日。研究选取美国国家冰雪数据中心(Nat1nal Snow and Ice Data Center, NSIDC)提供的Terra/MODIS每日积雪合成产品 M0D10A1 和Aqua/MODIS每日积雪合成产品MYD10A1，空间分辨率均为500m JODI S每日积雪合成产品共有四个数据子集: 积雪、积雪反照率、积雪覆盖比例和质量验证数据。研究区域覆盖的轨道号为h23v04/ h24v04,所用影像共80景。
[0024]第二步、对两种M0DIS影像数据产品分别进行预处理，提取出覆盖研究区的积雪反照率数据。利用MRT软件进行处理，如图2所示，具体包括如下几个方面的内容:a、根据研究区的范围，将同一时间相应的两轨影像拼接为一景完整的影像。
[0025]b、对拼接后的影像数据进行地图投影转换，将原始数据产品的正弦曲线投影 (SIN)重新投影为 WGS 1984/Geographic 系统。
[0026]c、对投影转换后的数据进行格式转换，将原始数据产品的HDF格式转化成容易处理的Geotiff格式。[〇〇27]d、M0DIS积雪产品中有四个数据子集，只提取其中的积雪反照率数据即可。[〇〇28]第三步、研究区域地面站点较少，只有32个。读取M0DIS积雪反照率影像，可以发现不少区域都容易受到云的影响。M0D10A1和MYD10A1的像元编码方式中，积雪反照率为0? 100，云为150，利用MATLAB软件对两种M0DIS反照率产品，运用像元重新处理的方法进行影像的去云合成。具体的算法如下:a、对于某一坐标的像元，如果产品MOD 10A1对应的反照率为积雪反照率，则把产品 M0D10A1的该坐标像元作为合成后积雪反照率影像的像元；b、对于某一坐标的像元，如果产品M0D10A1对应的反照率为非积雪反照率，而MYD10A1对应的反照率为积雪反照率，则把产品MYD10A1的该坐标像元作为合成后积雪反照率影像的像元；c、对于某一坐标的像元，如果产品M0D10A1和产品MYD10A1对应的反照率均非积雪反照率，则合成后积雪反照率影像在该坐标的像元为非积雪反照率像元。[〇〇29]第四步、根据每个地面站点的地理坐标，利用ENVI软件提取其对应积雪反照率影像上面对应的像元值，如图3所示。该步骤用到了两种方法:在ENVI软件主窗口中，利用 Pixe 1 Locator手动提取影像上某一位置的反照率值;利用ENVI+IDL语言，通过索引到影像上的行列号，批量获取地面站点对应影像上的反照率值。由于所用到的站点数较少，两种方法的效率相差不大，但当站点数目很多时，明显第二种方法效率更高。
[0030]第五步、为了避免所选样本中的异常值对最终的拟合结果造成较大的影响，地面站点的筛选参照了其它文献资料关于研究区域的结论，例如Ming等人2009年在 Atmospheric Research中的论文“Black Carbon (BC) in the snow of glaciers in west China and its potential effects on albedos”得到天山地区冰雪中黑碳含量 112 ±27ng g-1;Xu等人2012年在ENVIRONMENTAL RESEARCH LETTERS的论文 “Post-deposit1nal enrichment of black soot in snow-pack and accelerated melting of Tibetan glaciers”得出新疆天山一号冰川冰雪中冬季表层黑碳含量27?31 ng g^1，夏季可达250?500 ng g<。选取积雪表层黑碳浓度<300 ng的地面站点作为合适的样本点，进行关系模型的构建。[0〇31]对样本站点数据序列Y(x)，x=l，2,…，100,以线性函数Y(x)=ax+ b来拟合，按最小二乘法可求得a和b，其中，x为积雪反照率，Y为黑碳浓度，a为趋势项，b为常数项。a值的正负表示黑碳浓度随积雪反照率变化的方向，a绝对值的大小表示变化的快慢程度。统计样本数为26,通过显著性a=0.1的显著性检验。
[0032]第六步、根据初步建立的关系模型，在MATLAB软件里面编写程序，完成反照率影像向黑碳浓度空间数据的转换，如图4所示。对已建立的关系模型及初步得到的黑碳浓度结果，结合样本数据、相关文献资料进行检验评定，反复优化，直到确立最终的关系模型。
[0033]对关系模型检验评定考虑但不仅限于这些方面:所有样本点拟合值的均方差要尽可能的小;拟合的线性函数尽可能逼近所有样本点，相关系数R越接近1越好;整个反照率区间对应的拟合值都要<300 ng g^1;影像处理后黑碳浓度满足该地区的一些基本情况，如纯净积雪区(积雪深度大、人类活动影响小的地区)黑碳浓度低，城镇密集的地方(如天山北坡经济区)黑碳浓度相应要高等。
[0034]由于缺少积雪反照率在20%以下的样本点数据，所以最终确定的关系模型为Y(x) =-2.4813x + 255.85，R2=0.5243，xG [20,100)。最终的关系模型如图5所示。
[0035] 第七步、在ArcGIS软件中，用空间分析(Spatial Analyst)模块下的Extract1n 工具，利用新疆北疆地区的矢量边界截取获得目标区域的黑碳浓度分布。对结果进行重分类，得到最终的结果图，如图6所示。
[0036]除上述实施例外，本发明还可以有其他实施方式。凡采用等同替换或等效变换形成的技术方案，均落在本发明要求的保护范围。
【主权项】
1.一种积雪中黑碳浓度的遥感估算方法，包括以下步骤:第一步、利用地面站点实测积雪中黑碳浓度数据，根据地面站点地理位置及测量时间， 下载对应时间和区域的MODIS每日积雪合成产品M0D10A1和MYD10A1;第二步、对MOD IS影像数据产品MOD 10A1和MYD10A1分别进行预处理，包括影像拼接、地 图投影转换和格式转换，并提取产品中的积雪反照率数据；第三步、对产品M0D10A1和产品MYD10A1获取的积雪反照率数据进行去云合成处理，获 得积雪反照率影像；第四步、依据地面站点地理坐标提取积雪反照率影像上的对应位置积雪反照率；第五步、参照文献资料，分析地面站点积雪中的黑碳浓度情况，选取合适站点，初步建 立研究区域积雪反照率与黑碳浓度的关系模型；第六步、进行积雪反照率影像向黑碳浓度空间数据的转换，并对已建立的关系模型及 初步得到的黑碳浓度结果，结合样本数据和相关文献资料进行检验评定，反复优化，直到确 立最终的关系模型；第七步、通过研究区域矢量边界数据，根据最终的关系模型提取出研究区域积雪中黑 碳浓度的时空分布结果。2.根据权利要求1所述，其特征在于:第二步中，影像拼接的方法是，根据研究区范围， 将同一时间相应的几景影像拼接为一景完整的影像。3.根据权利要求1所述，其特征在于:第三步中，所述去云合成处理的算法如下:a、对于某一坐标的像元，如果产品M 0 D10 A1对应的反照率为积雪反照率，则把产品 M0D10A1的该坐标像元作为合成后积雪反照率影像的像元；b、对于某一坐标的像元，如果产品M0D10A1对应的反照率为非积雪反照率，而MYD10A1 对应的反照率为积雪反照率，则把产品MYD10A1的该坐标像元作为合成后积雪反照率影像 的像元；c、对于某一坐标的像元，如果产品M0D10A1和产品MYD10A1对应的反照率均非积雪反照 率，则合成后积雪反照率影像在该坐标的像元为非积雪反照率像元。4.根据权利要求1所述，其特征在于:第四步中，根据地理位置提取积雪反照率，有两种 方法可供选择:地面站点较少时，在ENVI软件主窗口中，利用pixel locator提取相应坐标 位置的反照率;地面站点较多时，利用ENVI+IDL语言，通过索引到影像上的行列号，批量获 取地面站点对应影像上的反照率。5.根据权利要求1所述，其特征在于:第五步中，为了避免所选样本中的异常值对最终 的拟合结果造成较大的影响，应参考其他已有的研究成果对地面站点进行筛选。6.根据权利要求1所述，其特征在于:第五步中，采用最小二乘法对积雪反照率和黑碳 浓度数据进行线性拟合，从而建立积雪反照率和黑碳浓度的关系模型。
【文档编号】G06F17/50GK106021655SQ201610302401
【公开日】2016年10月12日
【申请日】2016年5月9日
【发明人】柯长青, 马东辉, 蔡宇, 邵珠德, 赵佳曼
【申请人】南京大学