CN108132220B - 林区机载推扫式高光谱影像的brdf归一化校正方法 - Google Patents

林区机载推扫式高光谱影像的brdf归一化校正方法 Download PDF

Info

Publication number
CN108132220B
CN108132220B CN201711429219.3A CN201711429219A CN108132220B CN 108132220 B CN108132220 B CN 108132220B CN 201711429219 A CN201711429219 A CN 201711429219A CN 108132220 B CN108132220 B CN 108132220B
Authority
CN
China
Prior art keywords
pixel
coordinate system
observation
sun
image
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201711429219.3A
Other languages
English (en)
Other versions
CN108132220A (zh
Inventor
荚文
庞勇
鞠洪波
李增元
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Research Institute Of Forest Resource Information Techniques Chinese Academy Of Forestry
Original Assignee
Research Institute Of Forest Resource Information Techniques Chinese Academy Of Forestry
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Research Institute Of Forest Resource Information Techniques Chinese Academy Of Forestry filed Critical Research Institute Of Forest Resource Information Techniques Chinese Academy Of Forestry
Priority to CN201711429219.3A priority Critical patent/CN108132220B/zh
Publication of CN108132220A publication Critical patent/CN108132220A/zh
Application granted granted Critical
Publication of CN108132220B publication Critical patent/CN108132220B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N21/00—Investigating or analysing materials by the use of optical means, i.e. using sub-millimetre waves, infrared, visible or ultraviolet light
    • G01N21/17—Systems in which incident light is modified in accordance with the properties of the material investigated
    • G01N21/25—Colour; Spectral properties, i.e. comparison of effect of material on the light at two or more different wavelengths or wavelength bands
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N21/00—Investigating or analysing materials by the use of optical means, i.e. using sub-millimetre waves, infrared, visible or ultraviolet light
    • G01N21/17—Systems in which incident light is modified in accordance with the properties of the material investigated
    • G01N21/25—Colour; Spectral properties, i.e. comparison of effect of material on the light at two or more different wavelengths or wavelength bands
    • G01N21/31—Investigating relative effect of material at wavelengths characteristic of specific elements or molecules, e.g. atomic absorption spectrometry
    • G01N21/35—Investigating relative effect of material at wavelengths characteristic of specific elements or molecules, e.g. atomic absorption spectrometry using infrared light
    • G01N21/3563—Investigating relative effect of material at wavelengths characteristic of specific elements or molecules, e.g. atomic absorption spectrometry using infrared light for analysing solids; Preparation of samples therefor
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N21/00—Investigating or analysing materials by the use of optical means, i.e. using sub-millimetre waves, infrared, visible or ultraviolet light
    • G01N21/17—Systems in which incident light is modified in accordance with the properties of the material investigated
    • G01N21/25—Colour; Spectral properties, i.e. comparison of effect of material on the light at two or more different wavelengths or wavelength bands
    • G01N21/31—Investigating relative effect of material at wavelengths characteristic of specific elements or molecules, e.g. atomic absorption spectrometry
    • G01N21/35—Investigating relative effect of material at wavelengths characteristic of specific elements or molecules, e.g. atomic absorption spectrometry using infrared light
    • G01N21/359—Investigating relative effect of material at wavelengths characteristic of specific elements or molecules, e.g. atomic absorption spectrometry using infrared light using near infrared light
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N21/00—Investigating or analysing materials by the use of optical means, i.e. using sub-millimetre waves, infrared, visible or ultraviolet light
    • G01N21/17—Systems in which incident light is modified in accordance with the properties of the material investigated
    • G01N21/55—Specular reflectivity
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N2201/00—Features of devices classified in G01N21/00
    • G01N2201/12—Circuits of general importance; Signal processing
    • G01N2201/127—Calibration; base line adjustment; drift compensation

Landscapes

  • Physics & Mathematics (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Chemical & Material Sciences (AREA)
  • Analytical Chemistry (AREA)
  • Biochemistry (AREA)
  • General Health & Medical Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Immunology (AREA)
  • Pathology (AREA)
  • Image Processing (AREA)

Abstract

本发明提供了一种针对起伏地形的林区机载推扫式高光谱影像的BRDF归一化校正方法，利用机载推扫式高光谱设备观测视场和飞行姿态信息以及数据采集时刻的太阳位置计算出影像像元基于平面的太阳‑观测几何，再基于高精度DEM数据提取对应像元的坡度和坡向信息，将像元的基于平面的太阳‑观测几何旋转到真实太阳‑观测几何，对分类后的影像数据提取各地物类型的像元构成多角度观测反射率数据集并依据真实太阳‑观测几何构建BRDF模型。最后采用乘法归一化因子，将影像内多个角度的方向反射率归一化到指定的某个观测‑太阳角度的反射率。本发明能有效校正起伏地形的林区机载推扫式高光谱影像的BRDF效应，对后续影像的定量研究具有重要意义。

Description

林区机载推扫式高光谱影像的BRDF归一化校正方法
技术领域
本发明涉及林业信息化领域，具体而言，涉及一种针对起伏地形的林区机载推扫式高光谱影像的二项反射分布函数(BRDF)归一化校正方法。
背景技术
自20世纪70年代以来，遥感技术的应用巨幅提升了人们对地物的探测能力并有条件对其开展深入的定量分析。地物的各项异性反射是自然界的宏观现象之一，具体表现为其表面反射随着入射辐射方向以及观测角度的改变而发生变化。地物的多角度观测现象会引起辐射失真，制约着对地物固有的反射属性的准确描述。对地物的反射特性展开二向反射分布函数(Bidirectional Reflectance Distribution Function，BRDF)研究，是定量遥感领域中重要的组成部分。
遥感影像中常见的同一地物表现出不同反射率的现象，主要由于宽传感器观测视场角以及观测区域的地形起伏等原因造成的对该地物进行了多角度地观测。观测视场角越大，观测区地形起伏变化越大，影像中地物的二项反射特性越明显，即影像中地物的反射率随着不断变化的观测角度而产生相应变化。星载MODIS遥感影像采用AMBRALS(algorithmfor MODIS bidirectional reflectance anisotropies of the land surface)算法对宽视场角带来的地物二项反射问题进行校正。而对于起伏地形的林区影像，由于坡度、坡向引起的像元太阳-观测几何等输入参数的变化，主要适用于平坦地形BRDF校正的AMBRALS算法已无法满足实际需求。对于起伏地形的林区影像的辐射校正，通常采用基于经验统计模型的方式。然而统计模型的参数是依据单景影像内地物分布情况而确定，模型参数极大依赖于输入影像的统计信息并非基于地物真实的物理辐射特性，因此基于经验统计模型的技术方法极大地制约了模型的应用拓展性。
机载推扫式高光谱影像具有高空间分辨率(＜5m)和高光谱分辨率(＜10nm)的特性，能对地物信息进行更详细和精准地刻画。机载数据的采集平台高度一般高于当地平均海拔500-2000m，远低于星载传感器的采集高度，由于其宽视场角观测导致的影像多角度效应也将更加明显；再加上起伏地形的干扰，极大地增加了准确描述地物反射特性的难度。目前的机载推扫式影像的BRDF校正方法大多利用传感器本身多角度观测的特性，构建多角度数据集，从而提取地物的BRDF特性进行校正。但该方法无法对起伏地形的区域的影像获取像元的真实太阳-观测几何，所以仍局限于平坦地形。因此如何对起伏地形的林区机载遥感图像进行BRDF校正是一个亟待解决的问题。
发明内容
本发明的目的是解决现有的BRDF校正算法暂时无法依据机载推扫式高光谱设备的成像特点以及在起伏地形的林区的复杂成像环境下的辐射校正问题。充分利用机载推扫式高光谱设备的宽视场角多角度观测和地形起伏带来的太阳-观测几何变化等信息，重新计算影像中像元的真实太阳-观测几何，构建新型多角度影像数据集。新型多角度数据集将满足现有的BRDF校正算法，由此可以解决起伏地形的林区机载推扫式高光谱影像的二项反射分布函数(BRDF)校正问题。
为达到上述目的，本发明提供了一种针对起伏地形的林区机载推扫式高光谱影像的二项反射分布函数(BRDF)归一化校正方法，包括以下步骤：
(1)针对几何校正和大气校正后的机载高光谱影像，根据机载推扫式高光谱成像仪设备观测视场和数据采集时飞行姿态以及太阳几何位置信息，计算影像中每个像元基于平面的太阳-观测几何，包括观测方位角、观测天顶角、太阳方位角和太阳天顶角；
(2)利用机载高光谱数据采集范围内的高精度数字高程模型(Digital FlevationModel，DEM)数据计算影像中每个像元的坡度、坡向信息；
(3)根据上述(2)所提供信息将每个像元的观测几何由全局坐标系变换到局部坐标系，获得像元基于起伏地形的真实的太阳-观测几何；
(4)依据坡度分层抽样提取同一树种像元的各个波段反射率以及该像元局部坐标系的太阳-观测几何数据，构建该树种各个波段的BRDF模型，并提取模型参数；
(5)利用乘法归一化因子，对机载高光谱影像中各个像元的不同角度的方向反射率归一化到指定太阳角度以及传感器观测方向的反射率值；
(6)验证并评估BRDF归一化校正结果。
进一步地，在计算影像像元的太阳-观测几何时，应结合机载推扫式高光谱成像仪设备观测视场信息和数据采集时飞行姿态数据。在传感器扫描过程中，原始影像的每个像素都具有特有的观测天顶角与观测方位角，因此需要对影像中的每个像元进行观测天顶角和观测方位角的计算。依据摄影测量共线方程，若将地面坐标系的原点位置平移到传感器扫描线中心位置，此时像框标坐标系坐标到地面空间坐标系坐标的转换过程为：
Figure BDA0001523406120000041
Figure BDA0001523406120000042
式(1)中，(x，y)和(x0，y0)分别为像素点坐标和像主点S在框标坐标系下的坐标；推扫式成像传感器对地表进行逐行的扫描成像，因此在同一行上的不同像素在像空间坐标系中(y-y0)坐标均为0；传感器扫描行宽度范围内的像素数目为nb，像素尺寸大小为p mm，则像主点坐标为
Figure BDA0001523406120000043
f为摄影中心到成像影像的垂距，即焦距。(x-x0，y-y0，-f)为像素点在像空间坐标系的三维坐标。(u，v，w)为像素点在地面空间坐标系中的三维坐标。(2)式中，
Figure BDA0001523406120000044
ω和k为外方位元素，确定了像空间坐标系三轴在地面坐标系中的方向，上述信息由机载POS数据提供，航向方位角为αaz。像素点的观测天顶角θv与方位角
Figure BDA0001523406120000045
计算公式如下式所示：
Figure BDA0001523406120000046
Figure BDA0001523406120000047
进一步地，将影像像元的太阳-观测几何根据其对应的地形坡度、坡向信息，由全局坐标系变换为局部坐标系，以此获得像元真实的太阳-观测几何。所描述的像元的太阳-观测几何仅针对于影像覆盖范围为平坦地形的地区。对于地形起伏的林区，由于每个像元对应的坡度和坡向角的不同，其太阳-观测几何也随之变化。为了获得像元的真实太阳-观测几何，需要将基于地球平面坐标系(全局坐标系)的像元按照其对应的坡度(α)和坡向(β)信息，转换到该像元对应坡面的局部坐标系中。利用太阳及传感器在像元的局部坐标系中的坐标，即可计算出该像元对应坡面的真实太阳入射及观测几何。全局坐标系到局部坐标系的转换方式包括：首先将全局坐标系绕w轴旋转(π/2-β)，再将旋转后的坐标系绕v轴旋转α，即得到基于坡面的局部坐标系。假设(x′，y′，z′)为通过坐标系转换将传感器在全局坐标系中的位置(u′，v′，w′)转换到局部坐标系中的新坐标换，则转换公式可表示为：
Figure BDA0001523406120000051
其中，r为传感器到像元的直线距离，在后续的角度计算中作为共有项将忽略，亦可将r看作单位距离1带入公式(5)计算。
在局部坐标系中，像元的真实观测天顶角θ′v和方位角
Figure BDA0001523406120000052
分别由式(6)、(7)计算获得：
Figure BDA0001523406120000053
Figure BDA0001523406120000061
同理，坡面像元的真实太阳天顶角和方位角亦可通过上述方式计算。
进一步地，高光谱影像可根据像元的真实太阳-观测几何视为新型多角度观测数据集。在假设影像范围区域具有均匀的森林结构，并且不随坡度和坡向变化的前提下，本发明对新型多角度数据集采用半经验线性核驱动二项反射模型对不同树种的BRDF效应进行模拟，核的组合为体散射Ross-Thick核和几何光学Li-Sparse核。BRDF模型及核的具体计算公式如下：
Figure BDA0001523406120000062
其中
Figure BDA0001523406120000063
是二项反射分布函数，它是太阳天顶角θ′s，观测天顶角θ′v，太阳和传感器的相对方位角
Figure BDA0001523406120000064
不同地物类型c以及波长λ的函数；Kvol和Kgeo分别代表体散射Ross-Thick核和几何光学Li-Sparse核，fiso(c，λ)、fvol(c，λ)和fgeo分别代表某一地物类型在某一波段反射率的BRDF模型中各核函数项所对应的系数。fiso(c，λ)、fvol(c，λ)和fgeo分别与地物类型和波段有关，因此在括号中添加了(c，λ)。
进一步地，由于机载推扫式高光谱影像的高空间分辨率的缘故，单条单波段航带影像通常包含几千万至几亿个像元数。为了快速及相对准确地求解BRDF模型参数，本发明依据坡度和坡向信息对影像分层抽样，抽样比例根据实际情况而定。抽样后获得m(m为正整数)个像元子集，再根据地物树种分类数据提取m个像元子集中同一植被类型n(n为正整数)个像元的各个波段的反射率值以及局部坐标系的太阳-观测几何数据，并利用最小二乘法对该地物某一波段反射率的BRDF模型参数进行解算，建模过程如式(9)所示。
Figure BDA0001523406120000071
求解X·B＝Y，则XT·X·B＝XT·Y，可求的B＝(XT·X)-1XT·Y，即求解出某一地物类型在某一波段反射率的BRDF模型中各核函数项所对应的系数fiso，fvol和fgeo。
本发明的针对起伏地形的林区机载推扫式高光谱影像的二项反射分布函数(BRDF)归一化校正方法，通过利用机载推扫式高光谱设备的宽视场角多角度观测和地形起伏带来的太阳-观测几何变化等信息，重新计算影像中像元的真实太阳-观测几何，准确提取地物的BRDF特性，从而对起伏地形的林区的机载推扫式高光谱影像进行辐射校正。本发明方法适应于大多数推扫式影像数据集，逻辑清楚，适应性强，综合考虑了地表实际情况，充分利用因地形起伏带来的多角度观测效应，弥补了原本传感器稀疏采样多角度观测的不足，能较为准确地提取地物的真实的二项反射分布函数(BRDF)特性，对影像的辐射校正以及后期的遥感定量分析具有理论和应用价值。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1为本发明一个实施例的针对起伏地形的林区机载推扫式高光谱影像的二项反射分布函数(BRDF)归一化校正方法流程图；
图2a-图2d为影像像元基于平面的太阳观测几何和真实观测几何对比图，其中，图2a为基于平面的像元观测天顶角图像，图2b为像元真实观测天顶角图像，图2c为基于平面的像元观测方位角图像，图2d为像元真实观测方位角图像。
图3a、图3b分别为典型树种思茅松和刺栲的近红外波段反射率的BRDF特性图；
图4a、图4b分别为BRDF归一化校正影像前后对比图，图4c、图4d分别为图4a、图4b局部细节图；
图5a、图5b分别为BRDF归一化校正前后思茅松光谱曲线对比图，图中A、B、C和D代表图4c、图4d中不同位置的思茅松光谱数据。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有付出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
如图1所示，一种针对起伏地形的林区机载推扫式高光谱影像的二项反射分布函数(BRDF)归一化校正方法，具体实施步骤如下：
步骤一：计算像元的太阳-观测几何
相对于影像，太阳天顶角是指太阳直射光线与地平面垂线间的夹角。太阳方位角是指太阳相对于像元所处的方位，通常正北方向为0°，正东方向为90°。太阳的天顶角、方位角可以利用地理坐标位置以及影像的拍摄时间计算获得。由于机载设备采集单条航带的影像通常在10至20分钟内，太阳的位置变化不大，故可简化利用影像航带中心位置及该位置采集时间对应的固定太阳天顶角和方位角作为该航带影像所有像元对应的太阳几何信息。若单条航带影像采集时间过长，则可按照上述方法分段计算影像像元对应的太阳几何信息。
本实例根据机载推AISA Eagle II扫式高光谱成像传感器的特性，在传感器扫描过程中，原始影像的每个像素都具有特有的观测天顶角与观测方位角，因此需要对影像中的每个像元进行观测天顶角和观测方位角的计算。依据摄影测量共线方程，若将地面坐标系的原点位置平移到传感器扫描线中心位置，此时像框标坐标系坐标到地面空间坐标系坐标的转换过程为：
Figure BDA0001523406120000101
Figure BDA0001523406120000102
式(1)中，(x，y)和(x0，y0)分别为像素点坐标和像主点S在框标坐标系下的坐标；推扫式成像传感器对地表进行逐行的扫描成像，因此在同一行上的不同像素在像空间坐标系中(y-y0)坐标均为0。传感器扫描行宽度范围内的像素数目为1024，像素尺寸大小为0.012mm，则像主点坐标为
Figure BDA0001523406120000103
传感器摄影中心到成像影像的垂距，即焦距f＝18.50mm。(x-x0，y-y0，-f)为像素点在像空间坐标系的三维坐标，本实例中表示为(
Figure BDA0001523406120000104
0，-18.5)。(u，v，w)为像素点在地面空间坐标系中的三维坐标。
(2)式中，
Figure BDA0001523406120000105
ω和k分别为偏角、倾角和旋角，确定了像空间坐标系三轴在地面坐标系中的方向，上述信息由机载POS数据提供，航向方位角为αaz。像素点的观测天顶角θv与方位角
Figure BDA0001523406120000106
计算公式如下式所示：
Figure BDA0001523406120000111
Figure BDA0001523406120000112
步骤二：数据准备
(1)高光谱影像几何校正与大气校正：根据步骤一所述，构建对应原始影像行列号的每个像元的太阳-观测几何查找表影像数据。将原始高光谱影像和查找表数据同时进行几何校正。几何校正后的高光谱影像原始行列号将根据POS数据中的航迹信息和高分辨率DEM数据重行排列，但同样经过几何级校正后的查找表数据依旧可将太阳-观测几何信息对应到几何校正后的高光谱影像中的像元。
依据机载定标参数文件对上述几何校正后的高光谱影像进行辐射定标并利用ATCOR4机载大气校正软件或模型(进行影像大气校正，最终得到经过几何校正和大气校正后的地物反射率高光谱影像数据。
(2)像元的坡度和坡向信息：利用上述航带影像覆盖范围内的高分辨率DEM数据计算像元的坡度α，坡向β信息。对计算所得的坡度和坡向影像进行空间重采样到，获得与高光谱影像一致的空间分辨率。
(3)地物类型分类：针对林区影像特点，利用上述地物反射率高光谱影像数据进行地物分类，主要分为：针叶林，阔叶林，草地，裸土，水体，道路，建筑物等，其中针叶林和阔叶林可进一步细分至树种。由此，获得该段航带影像范围内的地物分类图。
步骤三：计算像元的真实太阳-观测几何
步骤一中所描述的像元的太阳-观测几何仅针对于影像覆盖范围为平坦地形的地区。对于地形起伏的林区，由于每个像元对应的坡度和坡向角的不同，其太阳-观测几何也随之变化。为了获得像元的真实太阳-观测几何，需要将基于地球平面坐标系(全局坐标系)的像元按照其对应的坡度(α)和坡向(β)信息，转换到该像元对应坡面的局部坐标系中。利用太阳及传感器在像元的局部坐标系中的坐标，即可计算出该像元对应坡面的真实太阳入射及观测几何。全局坐标系到局部坐标系的转换方式包括：首先将全局坐标系绕w轴旋转(π/2-β)，再将旋转后的坐标系绕v轴旋转α，即得到基于坡面的局部坐标系。假设(x′，y′，z′)为通过坐标系转换将传感器在全局坐标系中的位置(u′，v′，w′)转换到局部坐标系中的新坐标换，则转换公式可表示为：
Figure BDA0001523406120000121
其中，r为传感器到像元的直线距离，在后续的角度计算中作为共有项将忽略，亦可将r看作单位距离1带入公式(5)计算。
在局部坐标系中，像元的真实观测天顶角θ′v和方位角
Figure BDA0001523406120000122
分别由式(6)、(7)计算获得：
Figure BDA0001523406120000123
Figure BDA0001523406120000124
同理，坡面像元的真实太阳天顶角和方位角亦可通过上述方式计算。
步骤四：提取地物BRDF特性
(1)BRDF模型
高光谱影像可根据像元的真实太阳-观测几何视为新型多角度观测数据集。在假设影像范围区域具有均匀的森林结构和地物结构，并且不随坡度和坡向变化的前提下，本发明对新型多角度数据集采用半经验线性核驱动二项反射模型对不同树种的BRDF效应进行模拟，核的组合为体散射Ross-Thick核和几何光学Li-Sparse核。BRDF模型及核的具体计算公式如下：
Figure BDA0001523406120000131
其中
Figure BDA0001523406120000132
是二项反射分布函数，它是太阳天顶角θ′s，观测天顶角θ′v，太阳和传感器的相对方位角
Figure BDA0001523406120000133
不同地物类型c以及波长λ的函数；
Kvol和Kgeo分别代表体散射Ross-Thick核和几何光学Li-Sparse核，其计算过程如下：
Figure BDA0001523406120000134
Figure BDA0001523406120000135
Figure BDA0001523406120000136
Figure BDA0001523406120000137
Figure BDA0001523406120000138
Figure BDA0001523406120000141
Figure BDA0001523406120000142
Figure BDA0001523406120000143
Figure BDA0001523406120000144
本发明依照MODIS BRDF/Albedo产品算法参数，设置模型中的
Figure BDA0001523406120000145
Figure BDA0001523406120000146
fiso(c，λ)、fvol(c，λ)和fgeo分别代表某一地物类型在某一波段反射率的BRDF模型中各核函数项所对应的系数。
(2)BRDF模型核参数求解
由于机载推扫式高光谱影像的高空间分辨率的缘故，单条单波段航带影像通常包含几千万至几亿个像元数，为了快速及相对准确地求解BRDF模型参数，本发明依据坡度和坡向信息对新型多角度观测数据集分层抽样，获得m(m为正整数)个像元子集，再根据地物树种分类数据提取m个像元子集中同一植被类型的n(n为正整数)个像元的各个波段的反射率值以及局部坐标系的太阳-观测几何数据，并利用最小二乘法对该地物单一波段反射率的BRDF模型参数进行解算，建模过程如式(18所示。
Figure BDA0001523406120000147
Figure BDA0001523406120000151
求解X·B＝Y，则XT·X·B＝XT·Y，可求的B＝(XT·X)-1XT·Y，即求解出某一地物类型在某一波段反射率的BRDF模型中各核函数项所对应的系数fiso，fvol和fgeo。
步骤五：高光谱影像BRDF校正
计算整条航带高光谱影像每个像元在各个波段的各项异性因子(AnisotropyFactor，ANIF)参数，生成影像对应的ANIF图。ANIF可以反映某一波段的任意方向反射率与指定方向反射率之间的关系，其表达式如下：
Figure BDA0001523406120000152
Figure BDA0001523406120000153
是某一地物在某一波段的任意太阳-观测条件下的拟合反射率，
Figure BDA0001523406120000154
是某一地物在某一波段的固定太阳天顶角
Figure BDA0001523406120000155
观测天顶角
Figure BDA0001523406120000156
以及太阳-传感器相对方位
Figure BDA0001523406120000157
条件下的拟合反射率，通常选择天底方向观测，即
Figure BDA0001523406120000158
利用上述ANIF参数对高光谱影像进行逐BRDF归一化校正，校正方式为：
Figure BDA0001523406120000161
其中，ρBRDF_Cor是经过BRDF归一化校正后的高光谱反射率影像，ρimage是步骤二中，经过几何校正和大气校正后的地物反射率高光谱影像数据。
最后，通过目视比较影像以及对比地面实测地物反射率验证算法精度。
以下1以思茅松为例进行案例说明。
在一台配置有Intel(R)Xeon(R)1.70GHz的双核处理器，32GB内存的工作站计算机上进行实施，以云南起伏地形林区的机载推扫式AISA Eagle II高光谱影像数据为例，采用本发明方法，对影像数据进行BRDF归一化校正(图1)。对比未添加地形因素的像元观测几何(图2a、图2c)和像元真实观测几何(图2b、图2d)，可以看出结合地形因素的像元真实观测几何包含更加丰富的地物多角度观测信息，避免了因稀疏的多角度观测而造成BRDF核驱动模型反演时无约束所带来的BRDF形状剧烈变化的问题。图3a、图3b表示影像中典型植被类型在红波段和近红外波段反射率的BRDF形状，可以看出反演结果合理，当太阳天顶角为30°时，在热点方向上地物表现出较高的反射率。图4a、图4b反映了BRDF归一化校正后的影像较相比于原始反射率影像在一定程度上消除了地形因素对植被反射率的影响，图4c、图4d局部细节展示出影像中同一植被类型表现出相似的光谱反射特性。由图5b可得，BRDF归一化校正后影像中，在阴坡和阳坡的思茅松表现出相似的光谱反射特性，改善了地形因素和机载推扫式传感器观测造成的植被光谱反射差异(图5a)。
上述实施例首先根据机载推扫式高光谱成像仪的特点对影像中每个像元计算其基于平面的太阳-观测几何；再根据DEM数据获取像元对应的坡度和坡向的信息，将像元的太阳-观测几何由全局坐标系转换为局部坐标系，即依据该像元的依据坡度和坡向将像元在平面上的太阳-观测几何旋转到坡面上，从而得到各个像元真实的太阳-观测几何，构建新型多角度数据集；通过半经验核驱动模型模型拟合新型多角度数据集提取不同树种在可见光到近红外波段的BRDF特征；最后采用乘法归一化因子，将影像内各像元的多个角度的方向反射率归一化到指定的某个观测-太阳角度的反射率。本发明能有效校正起伏地形的林区机载推扫式高光谱影像的BRDF效应，对后续影像的定量研究具有重要意义。
以上说明对本发明而言只是说明性的，而非限制性的，本领域普通技术人员理解，在不脱离以下所附权利要求所限定的精神和范围的情况下，可做出许多修改、变化或等效，但都将落入本发明的保护范围内。

Claims (1)

1.一种林区机载推扫式高光谱影像的BRDF归一化校正方法，所述林区机载推扫式高光谱影像为针对起伏地形的林区机载推扫式高光谱影像，其特征在于，包括以下步骤：
(1)针对几何校正和大气校正后的机载高光谱影像，根据机载推扫式高光谱成像仪设备观测视场和数据采集时飞行姿态以及太阳几何位置信息,计算影像中每个像元基于平面的太阳-观测几何信息，基于平面的太阳-观测几何信息包括观测方位角、观测天顶角、太阳方位角和太阳天顶角；
在计算影像像元的太阳-观测几何时，应结合机载推扫式高光谱成像仪设备观测视场信息和数据采集时飞行姿态数据；依据摄影测量共线方程，若将地面坐标系的原点位置平移到传感器扫描线中心位置,此时像框标坐标系坐标到地面空间坐标系坐标的转换过程为:
Figure FDA0002898704400000011
Figure FDA0002898704400000012
式(1)中，(x,y)和(x0,y0)分别为像素点坐标和像主点S在框标坐标系下的坐标；推扫式成像传感器对地表进行逐行的扫描成像，因此在同一行上的不同像素在像空间坐标系中(y-y0)坐标均为0；传感器扫描行宽度范围内的像素数目为nb,像素尺寸大小为p mm，则像主点坐标为
Figure FDA0002898704400000013
f为摄影中心到成像影像的垂距，即焦距，(x-x0，y-y0，-f)为像素点在像空间坐标系的三维坐标，(u,v,w)为像素点在地面空间坐标系中的三维坐标，(2)式中，
Figure FDA0002898704400000021
ω和k为外方位元素，确定了像空间坐标系三轴在地面坐标系中的方向,上述信息由机载POS数据提供,航向方位角为αaz，像素点的观测天顶角θv与方位角
Figure FDA0002898704400000022
计算公式如下式所示：
Figure FDA0002898704400000023
Figure FDA0002898704400000024
(2)利用机载高光谱数据采集范围内的高精度数字高程模型数据计算影像中每个像元的坡度、坡向信息；根据像元的坡度、坡向信息将步骤(1)中各个像元的观测几何由全局坐标系变换到局部坐标系，获得像元基于起伏地形的真实的太阳-观测几何信息；
为了获得像元的真实太阳-观测几何，需要将基于地球平面坐标系(全局坐标系)的像元按照其对应的坡度(α)和坡向(β)信息，转换到该像元对应坡面的局部坐标系中，利用太阳及传感器在像元的局部坐标系中的坐标，即可计算出该像元对应坡面的真实太阳入射及观测几何；全局坐标系到局部坐标系的转换方式包括：首先将全局坐标系绕w轴旋转(π/2-β)，再将旋转后的坐标系绕v轴旋转α，即得到基于坡面的局部坐标系；假设(x′，y′，z ′)为通过坐标系转换将传感器在全局坐标系中的位置(u′,v′,w′)转换到局部坐标系中的新坐标换，则转换公式可表示为：
Figure FDA0002898704400000025
Figure FDA0002898704400000031
其中，r为传感器到像元的直线距离，在后续的角度计算中作为共有项将忽略，亦可将r看作单位距离1带入公式(5)计算；
在局部坐标系中，像元的真实观测天顶角θ′v和方位角
Figure FDA0002898704400000034
分别由式(6)、(7)计算获得：
Figure FDA0002898704400000032
Figure FDA0002898704400000033
同理，坡面像元的真实太阳天顶角和方位角亦可通过上述方式计算；
(3)依据坡度分层抽样提取同一树种像元的各个波段反射率以及该像元局部坐标系的太阳-观测几何信息，构建该树种各个波段的BRDF模型，并提取模型参数；利用乘法归一化因子，对机载高光谱影像中各个像元的不同角度的方向反射率归一化到指定太阳角度以及传感器观测方向的反射率值；验证并评估BRDF归一化校正结果；
其中，高光谱影像可根据像元的真实太阳-观测几何视为新型多角度观测数据集；在假设影像范围区域具有均匀的森林结构和地物结构，并且不随坡度和坡向变化的前提下，对新型多角度数据集采用半经验线性核驱动二项反射模型对不同树种的BRDF效应进行模拟，核的组合为体散射Ross-Thick核和几何光学Li-Sparse核；BRDF模型及核的具体计算公式如下：
Figure FDA0002898704400000041
其中
Figure FDA0002898704400000042
是二项反射分布函数，它是太阳天顶角θ′s，观测天顶角θ′v，太阳和传感器的相对方位角
Figure FDA0002898704400000043
不同地物类型c以及波长λ的函数；Kvol和Kgeo分别代表体散射Ross-Thick核和几何光学Li-Sparse核，fiso(c,λ)、fvol(c,λ)和fgeo分别代表某一地物类型在某一波段反射率的BRDF模型中各核函数项所对应的系数；
再者，依据坡度和坡向信息对影像分层抽样，抽样比例根据实际情况而定，抽样后获得m(m为正整数)个像元子集，再根据地物树种分类数据提取m个像元子集中同一植被类型的n(n为正整数)个像元的各个波段的反射率值以及局部坐标系的太阳-观测几何数据，并利用最小二乘法对该地物某一波段反射率的BRDF模型参数进行解算,建模过程如式(9)所示，
Figure FDA0002898704400000044
Figure FDA0002898704400000051
求解X·B＝Y，则XT·X·B＝XT·Y,可求的B＝(XT·X)-1XT·Y，即求解出某一地物类型在某一波段反射率的BRDF模型中各核函数项所对应的系数fiso，fvol和fgeo。
CN201711429219.3A 2017-12-25 2017-12-25 林区机载推扫式高光谱影像的brdf归一化校正方法 Active CN108132220B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201711429219.3A CN108132220B (zh) 2017-12-25 2017-12-25 林区机载推扫式高光谱影像的brdf归一化校正方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201711429219.3A CN108132220B (zh) 2017-12-25 2017-12-25 林区机载推扫式高光谱影像的brdf归一化校正方法

Publications (2)

Publication Number Publication Date
CN108132220A CN108132220A (zh) 2018-06-08
CN108132220B true CN108132220B (zh) 2021-03-05

Family

ID=62393019

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201711429219.3A Active CN108132220B (zh) 2017-12-25 2017-12-25 林区机载推扫式高光谱影像的brdf归一化校正方法

Country Status (1)

Country Link
CN (1) CN108132220B (zh)

Families Citing this family (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109932341B (zh) * 2019-03-11 2021-03-23 北京环境特性研究所 野外环境下典型目标的双向反射分布函数测量方法
CN109974854B (zh) * 2019-03-18 2020-11-20 石河子大学 一种框幅式fpi高光谱图像的辐射校正方法
CN110083176B (zh) * 2019-05-05 2020-07-24 宁夏大学 一种基于无人机载高光谱成像的brdf数据采集系统和方法
CN110702228B (zh) * 2019-09-25 2021-06-25 华东师范大学 一种航空高光谱影像的边缘辐射校正方法
CN113155740B (zh) * 2020-01-07 2023-05-26 国家卫星气象中心(国家空间天气监测预警中心) 一种定标基准场brdf特性分析方法及系统
CN112504998B (zh) * 2020-09-18 2021-12-17 南京大学 一种泡沫材料的大视场太赫兹无损检测方法
CN113008834A (zh) * 2021-02-09 2021-06-22 中国农业大学 基于遥感影像的地表反射率校正方法及装置
CN113029977B (zh) * 2021-03-11 2022-03-15 武汉大学 一种针对宽视场角多光谱传感器的自动交叉辐射定标方法
CN114581784B (zh) * 2022-05-07 2022-08-12 自然资源部第二海洋研究所 一种长时序逐年红树林遥感监测产品的构建方法

Family Cites Families (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102073038B (zh) * 2010-11-29 2013-06-19 上海大学 基于微小地形的遥感影像的地形校正的方法
CN102902883A (zh) * 2012-09-24 2013-01-30 北京师范大学 一种基于多角度测量构建二向性反射分布函数(brdf)原型库的方法
CN103413014A (zh) * 2013-03-11 2013-11-27 北京师范大学 一种基于二向性反射分布函数(brdf)原型反演地表反照率的方法
CN103324827B (zh) * 2013-04-09 2016-02-24 北京师范大学 一种改善业务化核驱动二向性反射分布函数(brdf)模型热点的方法
US10311163B2 (en) * 2014-06-30 2019-06-04 Microsoft Technology Licensing, Llc Non-parametric microfacet factor models for isotropic bidirectional reflectance distribution functions
CN104156567B (zh) * 2014-07-23 2017-05-03 中国科学院遥感与数字地球研究所 一种耦合卫星遥感影像大气校正和地形校正过程的地表反射率获取技术
CN104834814B (zh) * 2015-04-29 2018-04-10 西北师范大学 遥感影像地形标准化方法
WO2017048674A1 (en) * 2015-09-14 2017-03-23 University Of Florida Research Foundation, Inc. Method for measuring bi-directional reflectance distribution function (brdf) and associated device
CN105716583B (zh) * 2016-01-26 2018-03-30 河海大学 一种基于平行摄影的勘探洞地质编录底图生成方法
CN106971156A (zh) * 2017-03-22 2017-07-21 中国地质科学院矿产资源研究所 一种基于面向对象分类的稀土开采区遥感信息提取方法
CN107330473A (zh) * 2017-07-06 2017-11-07 常州市星图测绘科技有限公司 一种基于空间邻接指数的遥感分类野外调查样方抽选方法

Also Published As

Publication number Publication date
CN108132220A (zh) 2018-06-08

Similar Documents

Publication Publication Date Title
CN108132220B (zh) 林区机载推扫式高光谱影像的brdf归一化校正方法
Verger et al. Green area index from an unmanned aerial system over wheat and rapeseed crops
Li et al. An evaluation of the use of atmospheric and BRDF correction to standardize Landsat data
Smith et al. High spatial resolution data acquisition for the geosciences: kite aerial photography
Ni et al. Mapping three-dimensional structures of forest canopy using UAV stereo imagery: Evaluating impacts of forward overlaps and image resolutions with LiDAR data as reference
CN111563962A (zh) 一种基于几何辐射一体化采样的遥感图像仿真方法
CN103761704B (zh) 基于红外遥感数据的图像生成方法和系统
Susaki et al. Validation of MODIS albedo products of paddy fields in Japan
CN110988909A (zh) 基于tls进行高寒脆弱区沙地植被的植被盖度测定方法
Dutta et al. Characterizing vegetation canopy structure using airborne remote sensing data
CN108898070A (zh) 一种基于无人机平台的高光谱遥感提取薇甘菊装置及方法
CN115187481A (zh) 一种机载推扫式高光谱影像辐射扰动校正方法
He et al. Direct estimation of land surface albedo from simultaneous MISR data
CN109377476A (zh) 遥感影像云检测特征参数的动态阈值获取方法及装置
CN112557325A (zh) 一种果树果品品质近地面遥感监测装置和方法
Näsi et al. Optimizing radiometric processing and feature extraction of drone based hyperspectral frame format imagery for estimation of yield quantity and quality of a grass sward
Harris et al. Radiometric homogenisation of aerial images by calibrating with satellite data
Rengarajan et al. Modeling and simulation of deciduous forest canopy and its anisotropic reflectance properties using the Digital Image and Remote Sensing Image Generation (DIRSIG) tool
Hasegawa et al. Seasonal change of bidirectional reflectance distribution function in mature Japanese larch forests and their phenology at the foot of Mt. Yatsugatake, central Japan
Landier et al. Chapter Remote Sensing Studies of Urban Canopies: 3D Radiative Transfer Modeling
Shell et al. A novel BRDF measurement technique with spatial resolution-dependent spectral variance
CN115524763B (zh) 一种多时相高分辨率山地卫星影像地形辐射校正方法
CN117554300B (zh) 山地地表反照率站点观测遥感空间降尺度方法
CN117371333B (zh) 一种基于核驱动模型的地表温度产品角度归一化方法
Zhou et al. Analysis of topographic effects on vegetation indices

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant