技术领域
[0001] 本发明属于遥感领域，具体涉及一种航空高光谱影像地表反射率反演方法。
相关背景技术
[0002] 高光谱遥感技术在矿产资源勘探、生态环境监测、军事目标探测识别等领域的广泛应用，充分体现了成像光谱技术获取几何、辐射和光谱多维信息“图谱合一”的优势。在地形起伏的山地区域，航空高光谱影像受大气和地形因素的共同影响，地表反射率的高精度反演是高光谱数据应用的基础。大气纠正的主要目的是减小大气吸收、散射和反射对遥感器接收信号的影响，地形纠正旨在减小地形效应导致的阳坡和阴坡同类地物光谱差异问题。
[0003] 当前，围绕航空高光谱影像地表反射率反演，缺少合适且高效的方法将大气纠正和地形纠正有机结合起来，已有的方法存在明显的过校正和校正参数难以获取的问题。本发明基于辐射传输模型开展大气纠正，其参数获取清晰简便，同时，引入SCS‑C模型进行地形纠正，较好地克服了传统地形纠正模型的过校正问题。
具体实施方式
[0062] 下面结合附图和实施例对本发明作进一步详细说明。
[0063] 本发明提供的一种航空高光谱影像地表反射率反演方法，从大气纠正和地形纠正的理论基础出发，基于航空高光谱数据反演大气水汽含量和气溶胶光学厚度，采用MODTRAN模型计算大气透过率、大气程辐射以及大气半球反照率等参数，基于可见光‑短波红外谱段的大气辐射传输方程和MODTRAN模拟的大气参数，同时基于SCS‑C模型计算地形校正参数，实现耦合大气纠正和地形纠正的航空高光谱地表反射率反演。
[0064] 如图1所示，具体包括以下步骤：
[0065] 步骤1：获取航空高光谱遥感影像数据和对应的数字高程模型DEM数据；
[0066] 步骤2：针对步骤1获得的航空高光谱遥感影像数据进行预处理，将传感器的DN值转换为辐亮度；
[0067] 针对步骤1获得的航空高光谱遥感影像数据进行的预处理包括几何纠正和辐射定标，传感器的DN值转换为辐亮度的计算公式为：
[0068] L＝gain×DN+offset (1)
[0069] 其中，L为辐亮度，DN为波段λ的像元值，gain和offset分别为波段λ的增益和偏置。
[0070] 步骤3：基于步骤2预处理后的航空高光谱遥感影像数数据，采用APDA比值法(大气预处理微分吸收法)反演大气水汽含量；
[0071] 采用APDA算法计算水汽含量PW按照以下两个公式完成：
[0072]
[0073] 其中，RAPDA为APDA通道比值，Latm,m为吸收通道大气程辐射，Latm,r为参考通道大气程辐射，λr表示参考通道中心波长，Lm和Lr分别是校正前吸收和参考通道的辐亮度。LIR([x],[y])|a表示通过点(x,y)在点x＝a处估值的线性回归直线。α、β、γ为拟合参数，与水汽含量、气溶胶类型有关；PW为水汽含量。
[0074] 具体实施时，包括以下步骤：
[0075] 步骤3.1：利用MODTRAN辐射传输模型生成水汽含量，查找表LUT1中参数，查找表中参数包括：水汽含量PW、水汽吸收通道和参考通道的程辐射、APDA通道比值RAPDA；
[0076] 步骤3.2：利用公式(3)拟合LUT1中所有的APDA通道比值RAPDA与对应的水汽含量PW，并存储拟合参数α、β、γ；
[0077] 步骤3.3：设置初始水汽含量PW1，从LUT1中获取水汽吸收通道以及参考通道的程辐射；
[0078] 步骤3.4：利用公式(2)计算像元的APDA通道比值RAPDA，并根据公式(4)计算出像元的水汽含量值PW2；
[0079] 步骤3.5：根据PW2结合LUT1，得到新的水汽吸收通道以及参考通道的程辐射值，替换之前获取的程辐射；
[0080] 步骤3.6：再次利用公式(2)计算APDA通道比值RAPDA，根据公式(4)将APDA比值转化成最终的水汽含量PW3。
[0081] 步骤4：计算气溶胶光学厚度；
[0082] 具体实施时，具体包括以下步骤：
[0083] 步骤4.1：根据NDVI值以及近红外波段表观反射率选择暗像元的位置；
[0084] 步骤4.2：根据水汽含量反演得到的平均水汽值，调用MODTRAN建立能见度查找表LUT2；
[0085] 步骤4.3：计算红、蓝波段的地表反射率，以及能见度查找表LUT2中不同能见度VIS对应的参数，包括：大气程辐射La、大气半球反照率S、辐射校正参数Gb、Gt；
[0086] 步骤4.4：将步骤4.3获取的参数代入大气辐射传输方程，结合查找表LUT2得到满足方程的能见度VIS值，再内插得到红、蓝波段对应的能见度VIS，计算二者的平均值即为该暗像元处的能见度VIS；
[0087] 步骤4.5：根据550nm处能见度和气溶胶光学厚度之间的关系，求得气溶胶光学厚度。
[0088] 步骤5：根据步骤3和步骤4获得的大气水汽含量和气溶胶光学厚度，结合MODTRAN模型和影像的元数据，对可见光‑短波红外谱段的大气辐射传输方程进行求解，得到大气程辐射La、大气透过率Gb、大气几何状况的系数Gt和大气半球反照率S等其他参数，进而得到中间反射率；
[0089] 大气辐射传输方程为：
[0090]
[0091] 其中，ρ为像元的地表反射率，ρe为像元与附近像元的平均地表反射率，S为大气向下的半球反照率，La为大气程辐射，Gb和Gt为依赖于大气的透过率和几何状况的系数。
[0092] 表示太阳辐射经大气入射到地表后又反射直接进入传感器的辐射亮度，表示地表像元的反射受到大气散射的影响而进入传感器的辐射亮度。
[0093] 公式中ρ和ρe不同，其主要原因是来自大气散射所引起的“邻近像元效应”，ρe可通过以下方程估算：
[0094]
[0095] 其中，Le为像元与像元周围空间辐射图像的平均值，可以通过原始影像进行计算。
[0096] 具体实施时，包括以下步骤：
[0097] 步骤5.1：根据得到的平均能见度，调用MODTRAN建立地表反射率查找表LUT3，[0098] 步骤5.2：根据像元的水汽含量PW，结合查找表LUT3空间插值得到该像元的大气校正参数：大气程辐射La、大气半球反照率S、辐射校正参数Gb、Gt；
[0099] 步骤5.3：根据公式(6)和步骤5.2获得的大气校正参数求得像元与附近像元的平均地表反射率ρe。
[0100] 步骤5.4：针对步骤5.1‑5.3获得的参数，根据公式(5)求解得到像元的地表反射率ρ，将此地表反射率称为中间反射率；
[0101] 步骤6：针对步骤1获得的数字高程模型DEM数据进行预处理
[0102] 针对步骤1获得的数字高程模型DEM数据DEM数据进行的预处理包括裁剪、平滑处理等操作；
[0103] 步骤7：基于步骤6处理后的数字高程模型DEM数据，进行地形相关因子的计算，得到坡度、坡向数据，并由坡度和坡向进一步计算得到太阳有效入射角；
[0104] 具体计算方法如下：
[0105]
[0106]
[0107] 其中，θt和Asp分别表示坡度和坡向，fx和fy分别是南北方向和东西方向的高程变化，it为太阳有效入射角，θs为太阳天顶角， 为太阳方位角。
[0108] 步骤8：计算SCS‑C地形纠正模型中的经验参数C值；
[0109] 具体计算方法如下：
[0110] DN＝acosit+b (10)
[0111] 在像元DN值与太阳有效入射角的余弦值之间建立回归方程，其中b为回归方程的截距，a为方程的斜率，得到的经验参数C＝b/a。
[0112] 步骤9：基于步骤5获得的中间反射率，结合SCS‑C模型进行地形纠正计算，得到最后的地表真实反射率。
[0113] 具体计算方法如下：
[0114]
[0115] 其中，ρT为坡面地表反射率，即步骤5中得到的中间反射率，ρH为最终的地表真实反射率。
[0116] 上面结合附图和实施例对本发明作了详细说明，但是本发明并不限于上述实施例，在本领域普通技术人员所具备的知识范围内，还可以在不脱离本发明宗旨的前提下作出各种变化。本发明中未作详细描述的内容均可以采用现有技术。
