一种抗大气气溶胶植被指数计算方法转让专利
申请号 : CN201510864360.0
文献号 : CN105527229B
文献日 : 2018-01-02
发明人 : 陈云浩 , 王丹丹
申请人 : 北京师范大学
摘要 :
本发明公开了一种抗大气气溶胶植被指数计算方法，本方法将卫星影像分割成若干个采样区域，并根据采样区域中邻域像元在红光和近红光外波段的表观反射率，来反演得到该采样区域中心像元受气溶胶影响前的植被指数。将所有采样区域中反演出的中心像元的植被指数，融合成相应的图像。通过本发明能够实现无需获得和输入复杂的大气轮廓参数，无需寻找特定条件的暗像元，避免参数空间分布异质性带来的模型误差以及避免待校正影像不存在暗像元而无法进行大气校正的技术效果。
权利要求 :
1.一种抗大气气溶胶植被指数计算方法，其特征在于，包括如下步骤；
步骤1：对卫星影像进行预处理，得到表观反射率；
步骤2：对卫星影像进行大气校正；
步骤3：去除卫星影像中云像元；
步骤4：设置卫星影像采样区域大小；并将卫星影像根据采样区域大小进行划分；所述采样区域根据均质强度进行设置；
步骤5：遍历采样区域，并判断是否遍历完成；
若遍历未完成时，则执行步骤501；
若遍历完成时，则执行步骤502；
步骤501：遍历计算采样区域中每个邻域像元与中心像元表观反射率之间连线的斜率值；并判断是否遍历完成；
若遍历完成时，则执行步骤50101至步骤50102；
若遍历未完成时，则执行步骤501；
步骤50101：对步骤501中是正数的所有斜率值求平均值；
步骤50102：通过步骤50101中的平均值，恢复对应中心像元的植被覆盖指数；并执行步骤5；
所述步骤50102中通过步骤50101中的平均值恢复对应中心像元的植被覆盖指数的计算公式如下，其中NDVIi’表示采样区域中心像元的植被覆盖指数；ki’表示平均值；
步骤502：根据恢复对应中心像元的植被覆盖指数对卫星影像进行图像融合。
2.如权利要求1所述的抗大气气溶胶植被指数计算方法，其特征在于，所述步骤501中对每个邻域像元与中心像元表观反射率之间连线的斜率值的计算公式如下，其中，kij’表示斜率值，Nj和Rj分别表示第j个邻域像元在红光和近红外波段的表观反射率，Ni和Ri分别为中心像元i在红光和近红外波段的表观反射率。
3.如权利要求1所述的抗大气气溶胶植被指数计算方法，其特征在于，所述步骤50101中对步骤501中是正数的所有斜率值求平均值的计算公式如下，其中，ki’表示平均值。
4.如权利要求1所述的抗大气气溶胶植被指数计算方法，其特征在于，所述采样区域根据均质强度进行设置的大小是5个像元*5个像元。
5.如权利要求1-4中任一项所述的抗大气气溶胶植被指数计算方法，其特征在于，所述步骤50101中还包括先对是正数的所有斜率值进行统计方法过滤。
6.如权利要求5所述的抗大气气溶胶植被指数计算方法，其特征在于，所述对是正数的所有斜率值进行统计方法过滤的方法包括，去除最小值、去除最大值以及去除均值两个标准差以外的值。
说明书 :
一种抗大气气溶胶植被指数计算方法
技术领域
[0001] 本发明涉及一种抗大气气溶胶植被指数计算方法,特别是涉及一种基于邻域像元计算抗大气气溶胶植被指数的方法。
背景技术
[0002] 基于遥感的植被覆盖指数(NDVI)指数在监测全球地表植被变化方面起重要作用。太阳辐射到达地表，经反射进入卫星传感器的过程受到大气扰动。其中，气溶胶是主要干扰因素。因此，气溶胶校正成为通过遥感观测得到地表真实反射率的关键环节。大气校正的误差直接影响后续应用包括图像分类，地表覆盖变化以及其他地表参数的反演；例如对NDVI影响。大气成分导致观测NDVI与地表真实NDVI产生误差中以气溶胶散射为主(降低0.04-
0.20个单位)，其次为水汽(降低0.04-0.08个单位)，最后是瑞丽散射(降低0.02-0.04个单位)。随着中国城市化进程加快，气溶胶对气候变化的影响越来越受到重视，也为光学遥感的发展带来挑战。大气气溶胶通过降低红光和近红外信号之间的对比度来降低植被指数值。Pinty and Verstraete指出校正后的NDVI值比观测NDVI值大0.15-0.2个单位(M.M.Verstraete,B.P.(1992).GEMI：a non-linear index to monitor global vegetation from satellites.Vegetation.)。
0.20个单位)，其次为水汽(降低0.04-0.08个单位)，最后是瑞丽散射(降低0.02-0.04个单位)。随着中国城市化进程加快，气溶胶对气候变化的影响越来越受到重视，也为光学遥感的发展带来挑战。大气气溶胶通过降低红光和近红外信号之间的对比度来降低植被指数值。Pinty and Verstraete指出校正后的NDVI值比观测NDVI值大0.15-0.2个单位(M.M.Verstraete,B.P.(1992).GEMI：a non-linear index to monitor global vegetation from satellites.Vegetation.)。
[0003] 目前，有两类方法用来消除气溶胶对植被指数的影响：第一，通过大气校正得到地表真实反射率。该类方法主要包括辐射传输模型方法(RTM)，改进的暗目标方法(IDOS)和经验线性方法(ELM)。(Zhou,J.,Wang,J.,Li,J.,&Hu,D.(2011).Atmospheric correction of PROBA/CHRIS data in an urban environment.International Journal of Remote Sensing,32(9),2591-2604.doi:10.1080/01431161003698443)。目前广泛应用的大气校正方法包括6S、MODTRAN和LOWTRAN。它们均基于辐射传输模型，AOD是其中重要输入参数之一。目前通过遥感手段反演气溶胶浓度的方法主要包括the ocean法、brightness法,contrast reduction法和Densely Dark Vegetation(DDV)方法(Yoram J.Kaufman,A.E.W.,Lorraine A.Remer,Bo-Cai Gao,Rong-Rong Li,and Luke Flynn.(1997).The MODIS 2.1-μm Channel-Correlation with Visible Reflectance for Use in Remote Sensing of Aerosol.IEEE Transactions on Geoscience and Remote Sensing,35(5),1286-1298.)。
第二，定义新的抗大气植被指数。Kaufman和Tanré利用蓝光波段校正气溶胶对红光波段的影响并提出新植被指数(ARVI)(Tanré,Y.J.K.a.D.(1992).Atmospherically resistant vegetation index(ARVI)for EOS-MODIS.IEEE Transactions on Geoscience and Remote Sensing,30(2),261-270.)。之后，既抗土壤背景又抗大气影响的植被指数如SARVI、MNDVI和SARVI2(Huete,A.R.,Liu,H.Q,Batchily,K,and Van Leeuwen,W.(1997).A comparison of vegetation indices over a global set of TM images for EOS-MODIS.Remote Sensing of Environment,59,440-451.)。Karnieli等人提出一种新的抗植被指数AFVI，该指数运用短波红外(1.6μm或2.1μm)代替红光波段。
第二，定义新的抗大气植被指数。Kaufman和Tanré利用蓝光波段校正气溶胶对红光波段的影响并提出新植被指数(ARVI)(Tanré,Y.J.K.a.D.(1992).Atmospherically resistant vegetation index(ARVI)for EOS-MODIS.IEEE Transactions on Geoscience and Remote Sensing,30(2),261-270.)。之后，既抗土壤背景又抗大气影响的植被指数如SARVI、MNDVI和SARVI2(Huete,A.R.,Liu,H.Q,Batchily,K,and Van Leeuwen,W.(1997).A comparison of vegetation indices over a global set of TM images for EOS-MODIS.Remote Sensing of Environment,59,440-451.)。Karnieli等人提出一种新的抗植被指数AFVI，该指数运用短波红外(1.6μm或2.1μm)代替红光波段。
[0004] 目前已有的得到抗气溶胶植被指数的方法存在一些局限性，主要表现在以下几个方面：1、DDV方法很难应用于没有暗像元存在的影响；2、气溶胶空间分布的异质性增加6S模型的复杂性；3、运用第三波段代替受气溶胶影响的波段会有输入参数空间分辨率与产品空间分辨率不一致和波段之间空间分辨率不一致的现象。
发明内容
[0005] 本发明的目的是提供一种抗大气气溶胶植被指数计算方法，以解决需获取以及输入复杂的大气轮廓参数、参数空间分布异质性带来的模型误差以及待校正影像不存在暗像元而无法进行大气校正的技术问题。
[0006] 为实现以上发明目的，本发明提供一种抗大气气溶胶植被指数计算方法，包括如下步骤；
[0007] 步骤1：对卫星影像进行预处理，得到表观反射率；
[0008] 步骤2：对卫星影像进行大气校正；
[0009] 步骤3：去除卫星影像中云像元；
[0010] 步骤4：设置卫星影像采样区域大小；并将卫星影像根据采样区域大小进行划分；所述采样区域根据均质强度进行设置；
[0011] 步骤5：遍历采样区域，并判断是否遍历完成；
[0012] 若遍历未完成时，则执行步骤501；
[0013] 若遍历完成时，则执行步骤502；
[0014] 步骤501：遍历计算采样区域中每个邻域像元与中心像元表观反射率之间连线的斜率值；并判断是否遍历完成；
[0015] 若遍历完成时，则执行步骤50101至步骤50102；
[0016] 若遍历未完成时，则执行步骤501；
[0017] 步骤50101：对步骤501中是正数的所有斜率值求平均值；
[0018] 步骤50102：通过步骤50101中的平均值，恢复对应中心像元的植被覆盖指数；并执行步骤5；
[0019] 步骤502：根据恢复对应中心像元的植被覆盖指数对卫星影像进行图像融合。
[0020] 进一步地，所述步骤501中对每个邻域像元与中心像元表观反射率之间连线的斜率值的计算公式如下，
[0021]
[0022] 其中，kij’表示斜率值，Nj和Rj分别表示第j个邻域像元在红光和近红外波段的表观反射率，Ni和Ri分别为中心像元i在红光和近红外波段的表观反射率。
[0023] 进一步地，所述步骤50101中对步骤501中是正数的所有斜率值求平均值的计算公式如下，
[0024]
[0025] 其中，ki’表示平均值。
[0026] 进一步地，所述步骤50102中通过步骤50101中的平均值恢复对应中心像元的植被覆盖指数的计算公式如下，
[0027]
[0028] 其中NDVIi’表示采样区域中心像元的植被覆盖指数。
[0029] 进一步地，所述采样区域根据均质强度进行设置的大小是5个像元*5个像元。
[0030] 进一步地，所述步骤50101中还包括先对是正数的所有斜率值进行统计方法过滤。
[0031] 进一步地，所述对是正数的所有斜率值进行统计方法过滤的方法包括，去除最小值、去除最大值以及去除均值两个标准差以外的值。
[0032] 与现有技术相比，本发明的有益效果是：
[0033] 1.运用将卫星影像进行采样区域划分，并对采样区域中每个邻域像元与中心像元表观反射率之间连线的斜率值的平均值进行计算，以此平均值计算出植被覆盖指数来对图像进行融合的技术方案，获得无需获取以及输入复杂的大气轮廓参数、避免参数空间分布异质性带来的模型误差以及避免待校正影像不存在暗像元而无法进行大气校正的技术效果；
[0034] 2.运用对斜率值进行统计方法处理的技术方案，获得了避免出现奇异值的像元而增大斜率平均值导致误差的技术效果。
附图说明
[0035] 图1是本发明的背景技术中的流程图；
[0036] 图2是本发明中的方法对气溶胶校正得到的植被覆盖指数与表观植被覆盖指数浓度变化的对比图。
具体实施方式
[0037] 下面结合附图和具体实施例对本发明作进一步说明。
[0038] 实施例1：
[0039] 如图1所示，本发明的抗大气气溶胶植被指数计算方法，包括如下步骤；
[0040] 步骤1：对卫星影像进行预处理，得到表观反射率；
[0041] 步骤2：对卫星影像进行大气校正；
[0042] 步骤3：去除卫星影像中云像元；
[0043] 步骤4：设置卫星影像采样区域大小；并将卫星影像根据采样区域大小进行划分；所述采样区域根据均质强度进行设置；
[0044] 步骤5：遍历采样区域，并判断是否遍历完成；
[0045] 若遍历未完成时，则执行步骤501；
[0046] 若遍历完成时，则执行步骤502；
[0047] 步骤501：遍历计算采样区域中每个邻域像元与中心像元表观反射率之间连线的斜率值；并判断是否遍历完成；
[0048] 若遍历完成时，则执行步骤50101至步骤50102；
[0049] 若遍历未完成时，则执行步骤501；
[0050] 步骤50101：对步骤501中是正数的所有斜率值求平均值；
[0051] 步骤50102：通过步骤50101中的平均值，恢复对应中心像元的植被覆盖指数；并执行步骤5；
[0052] 步骤502：根据恢复对应中心像元的植被覆盖指数对卫星影像进行图像融合；
[0053] 具体来说，当需要对卫星影像进行分析时，先对影像进行辐射定标来获得表观反射率，再运用大气校正软件；例如完整的遥感图像处理平台(ENVI)的Flash模块，对影像中水汽和臭氧进行校正，并对气溶胶模型以及气溶胶反演设置为不处理，并将初始能见度设置为100km。根据数据BQA波段质量文件去除影像中的云像元。设置采样区域大小，对均质性较强的地表；例如：大片森林、草地或农田，为了保证模型稳定性，可设置采样区域越大越好；对于均质性较弱的地表，或虽然内部均质但离散化、破碎化的地表，此时采用区域不宜设置过大。一般为保证模型的稳定性，将采样区域设置为5*5(单位为像元)最好。采样区域设置好后，将卫星影像分割成若干个采样区域。在每个采样区域中计算每个邻域像元与中心像元表观反射率之间连线的斜率值ki’，计算公式如下：
[0054]
[0055] 其中，kij’表示斜率值，Nj和Rj分别表示第j个邻域像元在红光和近红外波段的表观反射率，Ni和Ri分别为中心像元i在红光和近红外波段的表观反射率。
[0056] 将一个采样区域中的所有邻域像元与中心像元表观反射率之间连线的斜率值计算完成后，选取所有斜率值大于零的斜率kj’，计算其平均值，计算公式如下：
[0057]
[0058] 其中ki’表示所有斜率值大于零的斜率kj’的平均值。
[0059] 由于在计算斜率kj’的过程中，可能出现奇异值，对后续平均值的计算产生误差，所以在计算平均值的时候，对斜率值大于零的斜率kj’进行统计方法的处理；例如：在大片森林或者大片耕地区域选择去除最小值法；在斑块破碎化区域选择去除均值两个标准差以外的值的方法。
[0060] 一个采样区域的平均值ki’计算完成后，可以根据ki’对该采用区域的植被覆盖指数进行计算，计算公式如下：
[0061]
[0062] 其中NDVIi’是该采用区域的植被覆盖指数。
[0063] 最后根据所有采样区域的植被覆盖指数对卫星影像进行图像融合。
[0064] 实施例2：
[0065] 以下是利用模拟数据对本发明的验证。
[0066] 如图2，该图是运用本发明中的方法，对卫星影像校正得到的NDVI与表观NDVI随气溶胶浓度变化的对比图。具体利用的实验参数如下：利用ESPA提供的2014年8月19日上午11时中国天津地区Landsat8地表真实反射率产品，运用6S大气校正模型地表真实反射率影像在AOD＝0.05、AOD＝0.3、AOD＝0.5以及AOD＝1.0的情况下对应的表观反射率。
[0067] 通过图2可以看出，气溶胶对表观反射率有很大影响，导致表观NDVI无法真实反映地表植被覆盖情况。本发明很好的消除了气溶胶对NDVI的影响，实现了在复杂气溶胶浓度下对地表植被覆盖变化进行快速评价的效果。
[0068] 实施例3：
[0069] 表格1
[0070]
[0071] 表格2
[0072]
[0073] 如表格1、表格2所示，表格1是伊春、绥化地区运用本发明方法对气溶胶校正的精度验证结果，表格2是北京地区运用本发明方法对气溶胶校正的精度验证结果。
[0074] 图中“region”表示从遥感影像中选择的感兴趣区；“AE(absolute error)”表示绝对误差；“size”表示感兴趣区中像元的个数；“vegetation coverage”表示植被覆盖度等级；“extent of correction”表示气溶胶校正的程度(由表格第二列/第一列求得)。
[0075] 具体利用的实验参数如下：利用2014年8月19日上午11时中国天津地区与2014年8月9日上午10时中国黑龙江省伊春市和绥化市交界处的Lansat8OLI影像和ESPA提供的对应时间和地区的地表真实反射率产品。
[0076] 由于中国天津地区地处京津冀经济圈，气溶胶浓度较其他地区高；黑龙江省伊春市和绥化市境内河谷密布，有大片林地、草地和农田。先利用ENVI软件中Flash大气校正模块对影像进行水汽和臭氧校正，得到校正后的红光和近红外波段。后利用本发明中的方法在此基础上进行气溶胶校正，校正结果与地表真实反射率产品进行对比验证算法精度。
[0077] 验证结果表明，本发明中的方法用于在分子散射校正基础上的森林、农田或草地等均质区之上具有较高的精度。
[0078] 除上述实施例外，本发明还可以有其他实施方式，凡采用等同替换或等效变换形成的技术方案，均落在本发明的保护范围内。