CN109214122B - 一种获取像元尺度地表宽波段半球发射率的方法 - Google Patents
一种获取像元尺度地表宽波段半球发射率的方法 Download PDFInfo
- Publication number
- CN109214122B CN109214122B CN201811214646.4A CN201811214646A CN109214122B CN 109214122 B CN109214122 B CN 109214122B CN 201811214646 A CN201811214646 A CN 201811214646A CN 109214122 B CN109214122 B CN 109214122B
- Authority
- CN
- China
- Prior art keywords
- earth
- emissivity
- kernel function
- indicate
- band
- 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
Links
- 238000000034 method Methods 0.000 title claims abstract description 40
- 239000011159 matrix material Substances 0.000 claims description 42
- 230000002457 bidirectional effect Effects 0.000 claims description 32
- 238000005315 distribution function Methods 0.000 claims description 31
- 238000003384 imaging method Methods 0.000 claims description 9
- 238000013075 data extraction Methods 0.000 claims description 4
- 239000000284 extract Substances 0.000 abstract description 2
- 230000006870 function Effects 0.000 description 169
- 241001061225 Arcos Species 0.000 description 4
- 238000000605 extraction Methods 0.000 description 4
- 238000004364 calculation method Methods 0.000 description 3
- 238000005516 engineering process Methods 0.000 description 3
- 230000017105 transposition Effects 0.000 description 3
- 238000004458 analytical method Methods 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 230000005855 radiation Effects 0.000 description 2
- HCHKCACWOHOZIP-UHFFFAOYSA-N Zinc Chemical compound [Zn] HCHKCACWOHOZIP-UHFFFAOYSA-N 0.000 description 1
- 238000007792 addition Methods 0.000 description 1
- 230000015572 biosynthetic process Effects 0.000 description 1
- 239000000470 constituent Substances 0.000 description 1
- 238000007796 conventional method Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 230000003628 erosive effect Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000008689 nuclear function Effects 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 238000003786 synthesis reaction Methods 0.000 description 1
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 1
- 239000011701 zinc Substances 0.000 description 1
- 229910052725 zinc Inorganic materials 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01J—MEASUREMENT OF INTENSITY, VELOCITY, SPECTRAL CONTENT, POLARISATION, PHASE OR PULSE CHARACTERISTICS OF INFRARED, VISIBLE OR ULTRAVIOLET LIGHT; COLORIMETRY; RADIATION PYROMETRY
- G01J5/00—Radiation pyrometry, e.g. infrared or optical thermometry
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01J—MEASUREMENT OF INTENSITY, VELOCITY, SPECTRAL CONTENT, POLARISATION, PHASE OR PULSE CHARACTERISTICS OF INFRARED, VISIBLE OR ULTRAVIOLET LIGHT; COLORIMETRY; RADIATION PYROMETRY
- G01J5/00—Radiation pyrometry, e.g. infrared or optical thermometry
- G01J5/007—Radiation pyrometry, e.g. infrared or optical thermometry for earth observation
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01J—MEASUREMENT OF INTENSITY, VELOCITY, SPECTRAL CONTENT, POLARISATION, PHASE OR PULSE CHARACTERISTICS OF INFRARED, VISIBLE OR ULTRAVIOLET LIGHT; COLORIMETRY; RADIATION PYROMETRY
- G01J5/00—Radiation pyrometry, e.g. infrared or optical thermometry
- G01J2005/0077—Imaging
Landscapes
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- Spectroscopy & Molecular Physics (AREA)
- Geology (AREA)
- Life Sciences & Earth Sciences (AREA)
- General Life Sciences & Earth Sciences (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- Environmental & Geological Engineering (AREA)
- Photometry And Measurement Of Optical Pulse Characteristics (AREA)
Abstract
本发明公开了一种获取像元尺度地表宽波段半球发射率的方法，具体步骤如下：1)提取多时相相关数据，多时相相关数据包括地表窄波段方向发射率数据和观测天顶角数据；2)根据核驱动模型，将地表窄波段方向发射率进行参数化表达；同时，将地表窄波段方向发射率的核函数进行参数化表达；3)对地表窄波段方向发射率的核函数权重系数进行逐像元率定；4)利用步骤2)的参数化表达结果，对地表窄波段半球发射率进行参数化表达；5)获取像元尺度地表宽波段半球发射率。本发明突破现有方法物理意义不明确的限制，摆脱现有方法较多依赖假设条件且计算精度不高的尴尬现状，提高像元尺度地表宽波段半球发射率遥感获取的精度和效率。
Description
技术领域
本发明涉及一种方法，尤其涉及一种获取像元尺度地表宽波段半球发射率的方法，属于信息技术服务领域。
背景技术
地表发射率是地表固有属性，是地表发射的热辐射与同温度下黑体发射热辐射的比值，与地表组成成分、粗糙度、含水量等因素有关。在地表温度反演、地表分类、土壤形成和侵蚀、稀疏植被覆盖与变化估算、岩床制图和资源勘探、能量平衡研究、军事等方面都有重要的意义。然而，受制于遥感载荷自身限制，当前地表发射率仅能获取载荷观测天顶角所对应的地表窄波段方向发射率，这就使得在后续行业应用过程只能忽略地表发射率的方向性，简单假设地表发射率为朗伯表面，将地表方向发射率近似代替地表半球发射率，这就使得像元尺度地表宽波段半球发射率的获取具有不确定性。当地表发射率各向异性时，这种近似会造成较大的误差，难以得到精确的像元尺度地表宽波段半球发射率遥感产品，严重影响了地表半球发射率产品的质量，无法满足各行业领域对于地表半球发射率遥感应用的精度需求。
现有的像元尺度地表宽波段半球发射率的获取，主要是简单用地表方向发射率来近似代替半球发射率。现有方法主要存在两方面的缺陷：第一，没有从物理层面上考虑地表方向发射率和半球发射率的差别，没有有效建立起利用若干角度的地表方向发射率估算地表半球发射率的物理模型；第二，假设地表发射率为朗伯表面，忽略了地表方向发射率的角度效应，对于地表发射率各向异性的下垫面会造成较大的误差，难以得到精确的像元尺度地表宽波段半球发射率遥感产品。
发明内容
为了解决上述技术所存在的不足之处，本发明提供了一种获取像元尺度地表宽波段半球发射率的方法。
为了解决以上技术问题，本发明采用的技术方案是：一种获取像元尺度地表宽波段半球发射率的方法，具体步骤如下：
1)提取多时相相关数据，多时相相关数据包括地表窄波段方向发射率数据和观测天顶角数据；
2)根据核驱动模型，将地表窄波段方向发射率进行参数化表达；同时，采用不同初等函数形式对地表窄波段方向发射率的核函数进行曲线拟合，将地表窄波段方向发射率的核函数进行参数化表达；
3)依据参数化表达后的地表窄波段方向发射率、步骤1)获取的地表窄波段方向发射率数据和观测天顶角数据，对地表窄波段方向发射率的核函数权重系数进行逐像元率定；
4)利用步骤2)的参数化表达结果，对地表窄波段半球发射率进行参数化表达；
5)利用地表宽波段半球发射率与地表窄波段半球发射率间存在的高度相关性、步骤3)以及步骤4)的结果，获取像元尺度地表宽波段半球发射率。
进一步地，步骤1)中地表窄波段方向发射率数据的提取过程为：根据待获取地区的时空范围，下载相应卫星上搭载的成像光谱仪中的地表方向发射率产品，从地表方向发射率产品中提取成像光谱仪的波段号i的原始计数值，并根据公式①将波段的原始计数值转换到对应的地表窄波段方向发射率，公式①如下所示：
其中，表示第i个波段的地表窄波段方向发射率，θv表示观测天顶角，表示第i个波段地表窄波段方向发射率的原始计数值。
进一步地，步骤1)中观测天顶角数据的提取过程为：从地表方向发射率产品中提取的成像光谱仪的观测天顶角的原始计数值，并同时利用公式②将观测天顶角的原始计数值转换到对应的弧度为单位的观测天顶角，公式②如下所示：
其中，θv表示观测天顶角，表示观测天顶角的原始计数值。
进一步地，步骤2)中地表窄波段方向发射率的参数化表达结果如公式⑦所示：
其中，表示第i个波段的地表窄波段方向发射率，θv表示观测天顶角，Ifvol表示地表窄波段方向发射率的体散射核函数，其数值为双向反射分布函数的体散射核函数半球积分值，Ifgeo表示地表窄波段方向发射率的几何光学核函数，其数值为双向反射分布函数的几何光学核函数半球积分值， 是第i波段对应的核函数权重系数，分别代表了各向同性核函数、体散射核函数、几何光学核函数所占的比例权重。
进一步地，地表窄波段方向发射率的体散射核函数Ifvol的参数化表达结果如公式⑧所示：
其中，Ifvol表示地表窄波段方向发射率的体散射核函数，fvol表示双向反射分布函数体散射核函数，θv表示观测天顶角，θs表示太阳天顶角，表示观测方向与太阳方向间的相对方位角，cos和sin分别是余弦和正弦函数符号，∫代表积分符号，d代表微分符号。
进一步地，地表窄波段方向发射率的几何光学核函数Ifgeo的参数化表达结果如公式⑨所示：
其中，Ifgeo表示地表窄波段方向发射率的几何光学核函数，fgeo表示双向反射分布函数几何光学核函数，θv表示观测天顶角，θs表示太阳天顶角，表示观测方向与太阳方向间的相对方位角，cos和sin分别表示余弦和正弦函数符号，∫代表积分符号，d代表微分符号。
进一步地，步骤3)中地表窄波段方向发射率的核函数权重系数逐像元率定过程为：依照公式⑦，利用步骤1)获取的地表窄波段方向发射率以及观测天顶角数据，逐像元构造超定方程组，超定方程组如公式⑩所示：
Ei＝KWi ⑩
其中，Ei表示地表窄波段方向发射率对应的矩阵，K表示地表窄波段方向发射率核函数值矩阵，Wi表示核函数权重系数矩阵；如式子所示Ei、K、Wi分别表示为：
其中，表示第i个波段顺序号第t天的地表窄波段方向发射率，θv,t表示顺序号第t天的观测天顶角；Ifvol,t表示顺序号第t天的地表窄波段方向发射率的体散射核函数值，由公式⑧计算而得；Ifgeo,t表示顺序号第t天的地表窄波段方向发射率的几何光学核函数值，由公式⑨计算而得；表示第i个波段对应的核函数权重系数，分别代表了各向同性核函数、体散射核函数和几何光学核函数所占的比例权重；
根据公式⑩，并利用最小二乘数学法，获得核函数权重系数的唯一解，如式子所示：
其中，表示估算的第i个波段最优核函数权重系数矩阵，K表示地表窄波段方向发射率核函数值矩阵，Ei表示第i个波段地表窄波段方向发射率对应的矩阵，T代表了矩阵的转置符号,上标-1代表了矩阵的求逆符号。
进一步地，步骤4)中地表窄波段半球发射率的参数化表达结果如公式所示：
其中，表示第i个波段的地表窄波段半球发射率，表示第i波段对应的核函数权重系数，分别代表了各向同性核函数、体散射核函数和几何光学核函数所占的比例权重。
本发明突破现有方法物理意义不明确的限制，摆脱现有方法较多依赖假设条件，且计算精度不高而无法有效应用于行业应用的尴尬现状，提高像元尺度地表宽波段半球发射率遥感获取的精度和效率。即本发明摒弃了传统将地表方向发射率对半球发射率的近似，克服了用方向发射率代替地表半球发射率所带来的精度受限弊端，以核驱动模型为基础，考虑了地表发射率的方向性，构建了地表半球发射率的多时相核驱动物理模型。同时，摆脱了半球积分复杂计算的束缚，利用不同初等函数形式对核函数的半球积分值进行了曲线拟合，获取了积分值的近似解析表达式，简化了地表半球发射率的估算步骤，实现了利用多天不同观测角度的数据进行像元尺度地表宽波段半球发射率的遥感获取。
附图说明
图1为本发明的总体技术流程图。
具体实施方式
下面结合附图具体实施方式对本发明作进一步详细的说明。
如图1所示为本发明的总体技术流程图，主要分为以下五个步骤：
1)提取多时相相关数据，多时相相关数据包括地表窄波段方向发射率数据和观测天顶角数据；
2)根据核驱动模型，将地表窄波段方向发射率进行参数化表达；同时，采用不同初等函数形式对地表窄波段方向发射率的核函数进行曲线拟合，将地表窄波段方向发射率的核函数进行参数化表达；
3)依据参数化表达后的地表窄波段方向发射率、步骤1)获取的地表窄波段方向发射率数据和观测天顶角数据，对地表窄波段方向发射率的核函数权重系数进行逐像元率定；
4)利用步骤2)的参数化表达结果，对地表窄波段半球发射率进行参数化表达；
5)利用地表宽波段半球发射率与地表窄波段半球发射率间存在的高度相关性、步骤3)以及步骤4)的结果，获取像元尺度地表宽波段半球发射率。
本发明提供了一种获取像元尺度地表宽波段半球发射率的方法，是利用多时相不同观测角度的地表窄波段方向发射率和观测天顶角作为输入，结合核驱动模型半球积分的参数化数学表达式，采用最小二乘数学方法，确定核函数的权重系数，最终基于构建的地表半球发射率的核驱动物理模型实现像元尺度地表宽波段半球发射率的获取。
下面对本发明的步骤进行具体说明：
1)地表窄波段方向发射率数据和观测天顶角数据的提取；
首先，地表窄波段方向发射率数据的提取过程为：根据待获取地区的时空范围，下载相应卫星上搭载的成像光谱仪中的地表方向发射率产品，从地表方向发射率产品中提取成像光谱仪的波段号i的原始计数值，并根据公式①将波段的原始计数值转换到对应的地表窄波段方向发射率，公式①如下所示：
其中，表示第i个波段的地表窄波段方向发射率，θv表示观测天顶角，表示第i个波段地表窄波段方向发射率的原始计数值。
观测天顶角数据的提取过程为：从地表方向发射率产品中提取的成像光谱仪的观测天顶角的原始计数值，并同时利用公式②将观测天顶角的原始计数值转换到对应的弧度为单位的观测天顶角，公式②如下所示：
其中，θv表示观测天顶角，表示观测天顶角的原始计数值。
2)地表窄波段方向发射率的参数化表达
对于一个处于热平衡状态的物体来说，根据基尔霍夫定律，地表窄波段方向半球发射率可以表示为公式③：
其中，表示第i个波段的地表窄波段方向发射率，表示第i个波段的地表窄波段双向反射分布函数，θv表示观测天顶角，θs表示太阳天顶角，表示观测方向与太阳方向间的相对方位角，∫表示积分符号，d表示微分符号；
根据核驱动模型，地表双向反射分布函数可以进一步表示为不同核函数的线性加权值，如公式④所示：
其中，表示第i个波段的地表窄通道双向反射分布函数，θv表示观测天顶角，θs表示太阳天顶角，表示观测方向与太阳方向间的相对方位角，fvol和fgeo分别表示双向反射分布函数的体散射核函数和几何光学核函数， 是第i波段对应的核函数权重系数，分别代表了各向同性(常数)核函数、体散射核函数和几何光学核函数所占的比例权重；
双向反射分布函数的体散射核函数fvol表示为公式⑤：
其中，fvol表示双向反射分布函数的体散射核函数，θv表示观测天顶角，θs表示太阳天顶角，表示观测方向与太阳方向间的相对方位角，ξ表示相位角(是关于太阳天顶角、观测天顶角和相对方位角的函数)，cos和sin是余弦和正弦函数符号，arcos是反余弦函数符号。
双向反射分布函数的几何光学核函数fgeo表示为公式⑥：
其中，fgeo表示双向反射分布函数的几何光学核函数，θv表示观测天顶角，θs表示太阳天顶角，表示观测方向与太阳方向间的相对方位角，ξ和T均表示相位角(是关于太阳天顶角、观测天顶角和相对方位角的函数)，cos、sin、tan和sec分别是余弦、正弦、正切和正割函数符号，arcos是反余弦函数符号。
结合公式③和④，地表窄波段方向发射率的参数化为公式⑦：
其中，表示第i个波段的地表窄波段方向发射率，θv表示观测天顶角，Ifvol表示地表窄波段方向发射率的体散射核函数，其数值为双向反射分布函数的体散射核函数半球积分值，Ifgeo表示地表窄波段方向发射率的几何光学核函数，其数值为双向反射分布函数的几何光学核函数半球积分值， 是第i波段对应的核函数权重系数，分别代表了各向同性核函数、体散射核函数、几何光学核函数所占的比例权重。
由于双向反射分布函数的核函数的数学表达式也较为复杂，因此无法给出其半球积分值，即地表窄波段方向发射率的核函数的解析表达式。为了便于后续计算，这里采用不同初等函数形式对地表窄波段方向发射率的核函数进行了曲线拟合，将其重新参数化表示更为简单的形式。
地表窄波段方向发射率的体散射核函数Ifvol的参数化为公式⑧：
其中，Ifvol表示地表窄波段方向发射率的体散射核函数，fvol表示双向反射分布函数体散射核函数，θv表示观测天顶角，θs表示太阳天顶角，表示观测方向与太阳方向间的相对方位角，cos和sin分别是余弦和正弦函数符号，∫代表积分符号，d代表微分符号。
地表窄波段方向发射率的几何光学核函数Ifgeo的参数化为公式⑨：
其中，Ifgeo表示地表窄波段方向发射率的几何光学核函数，fgeo表示双向反射分布函数几何光学核函数，θv表示观测天顶角，θs表示太阳天顶角，表示观测方向与太阳方向间的相对方位角，cos和sin分别表示余弦和正弦函数符号，∫代表积分符号，d代表微分符号。
3)核函数权重系数的逐像元率定
依照公式⑦，假设地表方向发射率的变化仅仅是因为观测角度的差异导致的，那么利用步骤一获取的多时相地表窄波段方向发射率以及观测天顶角数据，可以逐像元构造超定方程组，超定方程组如公式⑩所示：
Ei＝KWi ⑩
其中，Ei表示地表窄波段方向发射率对应的矩阵，K表示地表窄波段方向发射率核函数值矩阵，Wi表示核函数权重系数矩阵；如式子所示Ei、K、Wi分别表示为：
其中，表示第i个波段顺序号第t天的地表窄波段方向发射率，θv,t表示顺序号第t天的观测天顶角；Ifvol,t表示顺序号第t天的地表窄波段方向发射率的体散射核函数值，由公式⑧计算而得；Ifgeo,t表示顺序号第t天的地表窄波段方向发射率的几何光学核函数值，由公式⑨计算而得；表示第i个波段对应的核函数权重系数，分别代表了各向同性核函数、体散射核函数和几何光学核函数所占的比例权重；
根据公式⑩，并利用最小二乘数学法，获得核函数权重系数的唯一解，如式子所示：
其中，表示估算的第i个波段最优核函数权重系数矩阵，K表示地表窄波段方向发射率核函数值矩阵，Ei表示第i个波段地表窄波段方向发射率对应的矩阵，T代表了矩阵的转置符号,上标-1代表了矩阵的求逆符号。
根据步骤1)读取的天顶角，并根据公式⑧和⑨计算地表窄波段方向发射率核函数值，并形成地表窄波段方向发射率核函数值矩阵K，同时根据步骤(1读取第i波段的地表窄波段方向发射率，形成各波段地表窄波段方向发射率对应的矩阵Ei，最后根据公式逐像元确定第i波段的最优核函数权重系数矩阵
4)地表窄波段半球发射率的参数化表达
地表窄波段方向发射率和地表窄波段半球发射率间存在积分关系，其数学形式可表示为公式
其中，表示第i个波段的地表窄波段半球发射率，表示第i个波段的地表窄波段方向发射率，θv表示观测天顶角，cos和sin分别是余弦和正弦函数符号，∫表示积分符号，d表示微分符号；
结合公式⑦，地表窄波段半球发射率进一步表示为公式
其中，表示第i个波段的地表窄波段半球发射率，IFvol表示地表窄波段半球发射率的体散射核函数值，其数值为地表窄波段方向发射率的体散射核函数半球积分值；IFgeo表示地表窄波段半球发射率的几何光学核函数值，其数值为地表窄波段方向发射率的几何光学核函数半球积分值；表示第i波段对应的核函数权重系数，分别代表了各向同性(常数)核函数、体散射核函数和几何光学核函数所占的比例权重；
利用公式⑧可获取地表窄波段半球发射率的体散射核函数值，如式子所示：
相应的，利用公式⑨可获取地表窄波段半球发射率的几何光学核函数值，如式子所示：
结合和公式可进一步表示为公式
其中，表示第i个波段的地表窄波段半球发射率，表示第i波段对应的核函数权重系数，分别代表了各向同性核函数、体散射核函数和几何光学核函数所占的比例权重。
5)像元尺度地表宽波段半球发射率的获取
由于地表宽波段半球发射率与地表窄波段半球发射率间存在高度的相关性，结合公式和获得像元尺度地表宽波段半球发射率的计算公式，再结合步骤3)获取的波段i的最优核函数权重系数矩阵就可实现像元尺度地表宽波段半球发射率的获取。
本发明主要解决的技术问题为两点：第一，摒弃了传统将地表方向发射率近似代替半球发射率的简化方法，通过对核驱动模型中不同核函数进行半球积分，挖掘了地表半球发射率与核函数之间的解析联系，构建了地表半球发射率的核驱动物理模型；第二，利用同一区域多天不同观测角度的地表方向发射率，结合最小二乘数学方法，估算了像元尺度对应的最优核系数，结合对核函数半球积分关系的重新参数化，简化了半球积分的计算过程，最终实现了高精度的像元尺度地表宽波段半球发射率的遥感获取。与现有方法相比，本发明所提出的方法物理意义明确，实现过程简单方便，估算精度可靠，可用于大范围地表宽波段半球发射率的遥感获取。
下面结合具体的实施例对本发明作进一步的说明：
1)多时相地表方向发射率和观测天顶角数据的提取
根据待获取地区的时空范围，相应的下载Terra或者Aqua卫星上搭载的中分辨率成像光谱仪MODIS的L3级的地表方向发射率产品MOD11C1或者MYD11C1。该发射率产品的时间跨度为最接近待获取时间的连续16天的多时相产品，其包含了MODIS不同观测角度对特定目标进行观测的一个完整周期。从HDF格式的不同天多时相产品中提取MODIS的29、31和32波段的原始计数值，并根据公式①将原始计数值转换到对应的地表窄波段方向发射率，公式①如下所示：
其中，表示第i个波段的地表窄波段方向发射率，波段号i的取值为29、31或者32；θv表示观测天顶角，表示第i个波段地表窄波段方向发射率的原始计数值。
相应的，从HDF格式的不同天多时相产品中也提取MODIS的观测天顶角，并同时利用公式②将原始计数值转换到对应的弧度为单位的观测天顶角，公式②如下所示：
其中，θv表示观测天顶角，表示观测天顶角的原始计数值。
2)地表窄波段方向发射率的参数化表达
对于一个处于热平衡状态的物体来说，根据基尔霍夫定律，地表窄波段方向半球发射率可以表示为公式③：
其中，表示第i个波段的地表窄波段方向发射率，波段号i的取值为29、31或者32；表示第i个波段的地表窄波段双向反射分布函数，θv表示观测天顶角，θs表示太阳天顶角，表示观测方向与太阳方向间的相对方位角，∫表示积分符号，d表示微分符号；
根据核驱动模型，地表双向反射分布函数可以进一步表示为不同核函数的线性加权值，如公式④所示：
其中，表示第i个波段的地表窄通道双向反射分布函数，波段号i的取值为29、31或者32；θv表示观测天顶角，θs表示太阳天顶角，表示观测方向与太阳方向间的相对方位角，fvol和fgeo分别表示双向反射分布函数的体散射核函数和几何光学核函数，是第i波段对应的核函数权重系数，分别代表了各向同性(常数)核函数、体散射核函数和几何光学核函数所占的比例权重；
双向反射分布函数的体散射核函数fvol表示为公式⑤：
其中，fvol表示双向反射分布函数的体散射核函数，θv表示观测天顶角，θs表示太阳天顶角，表示观测方向与太阳方向间的相对方位角，ξ表示相位角(是关于太阳天顶角、观测天顶角和相对方位角的函数)，cos和sin是余弦和正弦函数符号，arcos是反余弦函数符号。
双向反射分布函数的几何光学核函数fgeo表示为公式⑥：
其中，fgeo表示双向反射分布函数的几何光学核函数，θv表示观测天顶角，θs表示太阳天顶角，表示观测方向与太阳方向间的相对方位角，ξ和T均表示相位角(是关于太阳天顶角、观测天顶角和相对方位角的函数)，cos、sin、tan和sec分别是余弦、正弦、正切和正割函数符号，arcos是反余弦函数符号。
结合公式③和④，地表窄波段方向发射率的参数化为公式⑦：
其中，表示第i个波段的地表窄波段方向发射率，波段号i的取值为29、31或者32；θv表示观测天顶角，Ifvol表示地表窄波段方向发射率的体散射核函数，其数值为双向反射分布函数的体散射核函数半球积分值，Ifgeo表示地表窄波段方向发射率的几何光学核函数，其数值为双向反射分布函数的几何光学核函数半球积分值，是第i波段对应的核函数权重系数，分别代表了各向同性核函数、体散射核函数、几何光学核函数所占的比例权重。
由于双向反射分布函数的核函数的数学表达式也较为复杂，因此无法给出其半球积分值，即地表窄波段方向发射率的核函数的解析表达式。为了便于后续计算，这里采用不同初等函数形式对地表窄波段方向发射率的核函数进行了曲线拟合，将其重新参数化表示更为简单的形式。
地表窄波段方向发射率的体散射核函数Ifvol的参数化为公式⑧：
其中，Ifvol表示地表窄波段方向发射率的体散射核函数，fvol表示双向反射分布函数体散射核函数，θv表示观测天顶角，θs表示太阳天顶角，表示观测方向与太阳方向间的相对方位角，cos和sin分别是余弦和正弦函数符号，∫代表积分符号，d代表微分符号。
地表窄波段方向发射率的几何光学核函数Ifgeo的参数化为公式⑨：
其中，Ifgeo表示地表窄波段方向发射率的几何光学核函数，fgeo表示双向反射分布函数几何光学核函数，θv表示观测天顶角，θs表示太阳天顶角，表示观测方向与太阳方向间的相对方位角，cos和sin分别表示余弦和正弦函数符号，∫代表积分符号，d代表微分符号。
3)核函数权重系数的逐像元率定
依照公式⑦，假设地表方向发射率的变化仅仅是因为观测角度的差异导致的，那么利用步骤一获取的多时相地表窄波段方向发射率以及观测天顶角数据，可以逐像元构造超定方程组，超定方程组如公式⑩所示：
Ei＝KWi ⑩
其中，Ei表示地表窄波段方向发射率对应的矩阵，K表示地表窄波段方向发射率核函数值矩阵，Wi表示核函数权重系数矩阵；如式子所示Ei、K、Wi分别表示为：
其中，表示第i个波段顺序号第t天的地表窄波段方向发射率，波段号i的取值为29、31或者32；θv,t表示顺序号第t天的观测天顶角；Ifvol,t表示顺序号第t天的地表窄波段方向发射率的体散射核函数值，由公式⑧计算而得；Ifgeo,t表示顺序号第t天的地表窄波段方向发射率的几何光学核函数值，由公式⑨计算而得；表示第i个波段对应的核函数权重系数，分别代表了各向同性核函数、体散射核函数和几何光学核函数所占的比例权重；
根据公式⑩，并利用最小二乘数学法，获得核函数权重系数的唯一解，如式子所示：
其中，表示估算的第i个波段最优核函数权重系数矩阵，波段号i的取值为29、31或者32；K表示地表窄波段方向发射率核函数值矩阵，Ei表示第i个波段地表窄波段方向发射率对应的矩阵，T代表了矩阵的转置符号,上标“-1”代表了矩阵的求逆符号。
根据步骤1)读取的天顶角，并根据公式⑧和⑨计算地表窄波段方向发射率核函数值，并形成地表窄波段方向发射率核函数值矩阵K，同时根据步骤1)读取第29、31和32波段的地表窄波段方向发射率，形成各波段地表窄波段方向发射率对应的矩阵E29、E31和E32，最后根据公式逐像元确定第29、31和32波段的最优核函数权重系数矩阵和
4)地表窄波段半球发射率的参数化表达
地表窄波段方向发射率和地表窄波段半球发射率间存在积分关系，其数学形式可表示为公式
其中，表示第i个波段的地表窄波段半球发射率，波段号i的取值为29、31或者32；表示第i个波段的地表窄波段方向发射率，θv表示观测天顶角，cos和sin分别是余弦和正弦函数符号，∫表示积分符号，d表示微分符号；
结合公式⑦，地表窄波段半球发射率进一步表示为公式
其中，表示第i个波段的地表窄波段半球发射率，波段号i的取值为29、31或者32；IFvol表示地表窄波段半球发射率的体散射核函数值，其数值为地表窄波段方向发射率的体散射核函数半球积分值；IFgeo表示地表窄波段半球发射率的几何光学核函数值，其数值为地表窄波段方向发射率的几何光学核函数半球积分值；表示第i波段对应的核函数权重系数，分别代表了各向同性(常数)核函数、体散射核函数和几何光学核函数所占的比例权重；
利用公式⑧可获取地表窄波段半球发射率的体散射核函数值，如式子所示：
相应的，利用公式⑨可获取地表窄波段半球发射率的几何光学核函数值，如式子所示：
结合和公式可进一步表示为公式
其中，表示第i个波段的地表窄波段半球发射率，波段号i的取值为29、31或者32；表示第i波段对应的核函数权重系数，分别代表了各向同性核函数、体散射核函数和几何光学核函数所占的比例权重。
5)像元尺度地表宽波段半球发射率的获取。
利用地表宽波段半球发射率与地表窄波段半球发射率间存在的高度相关性，像元尺度地表宽波段半球发射率如公式所示：
其中，εh是地表宽波段半球发射率，是29波段的地表窄波段半球发射率，是31波段的地表窄波段半球发射率，是32波段的地表窄波段半球发射率。
结合公式和公式可以进一步表示为公式
其中，εh是地表宽波段半球发射率，和是根据公式逐像元确定的第29、31和32波段的最优核函数权重系数矩阵，C29、C31和C32是地表宽波段半球发射率的系数矩阵，如式子所示，其取值为：
至此，利用和以及步骤3)获取的第29、31和32波段的最优核函数权重系数矩阵和就可实现像元尺度地表宽波段半球发射率的获取。
本发明提供了一种获取像元尺度地表宽波段半球发射率的方法，具有以下优势：
第一，克服了地表方向发射率近似代替半球发射率的弊端，利用三种核函数综合表示不同条件下的地表半球发射率，建立起精度可靠的由地表方向发射率估算地表半球发射率的物理模型；
第二，降低了地表发射率角度效应对地表半球发射率的估算的影响，通过多天不同观测角度数据进行核系数的最小二乘稳健估算，同时借助核函数半球积分的数值模拟，获取了半球积分的近似解析表达式，摆脱了半球积分复杂计算的束缚，实现了像元尺度地表宽波段半球发射率遥感产品的精确高效获取。
上述实施方式并非是对本发明的限制，本发明也并不仅限于上述举例，本技术领域的技术人员在本发明的技术方案范围内所做出的变化、改型、添加或替换，也均属于本发明的保护范围。
Claims (8)
1.一种获取像元尺度地表宽波段半球发射率的方法，其特征在于：所述方法的具体步骤如下：
1)提取多时相相关数据，多时相相关数据包括地表窄波段方向发射率数据和观测天顶角数据；
2)根据核驱动模型，将地表窄波段方向发射率进行参数化表达；同时，采用不同初等函数形式对地表窄波段方向发射率的核函数进行曲线拟合，将地表窄波段方向发射率的核函数进行参数化表达；地表窄波段方向发射率的核函数包括各向同性核函数、体散射核函数、几何光学核函数；
3)依据参数化表达后的地表窄波段方向发射率、步骤1)获取的地表窄波段方向发射率数据和观测天顶角数据，对地表窄波段方向发射率的核函数权重系数进行逐像元率定；
4)利用步骤2)的参数化表达结果，对地表窄波段半球发射率进行参数化表达；
5)利用地表宽波段半球发射率与地表窄波段半球发射率间存在的高度相关性、步骤3)以及步骤4)的结果，获取像元尺度地表宽波段半球发射率。
2.根据权利要求1所述的获取像元尺度地表宽波段半球发射率的方法，其特征在于：步骤1)中所述地表窄波段方向发射率数据的提取过程为：根据待获取地区的时空范围，下载相应卫星上搭载的成像光谱仪中的地表方向发射率产品，从地表方向发射率产品中提取成像光谱仪的波段号i的原始计数值，并根据公式①将波段的原始计数值转换到对应的地表窄波段方向发射率，公式①如下所示：
其中，表示第i个波段的地表窄波段方向发射率，θv表示观测天顶角，表示第i个波段地表窄波段方向发射率的原始计数值。
3.根据权利要求2所述的获取像元尺度地表宽波段半球发射率的方法，其特征在于：步骤1)中所述观测天顶角数据的提取过程为：从地表方向发射率产品中提取的成像光谱仪的观测天顶角的原始计数值，并同时利用公式②将观测天顶角的原始计数值转换到对应的弧度为单位的观测天顶角，公式②如下所示：
其中，θv表示观测天顶角，表示观测天顶角的原始计数值。
4.根据权利要求3所述的获取像元尺度地表宽波段半球发射率的方法，其特征在于：步骤2)中所述地表窄波段方向发射率的参数化表达结果如公式⑦所示：
其中，表示第i个波段的地表窄波段方向发射率，θv表示观测天顶角，Ifvol表示地表窄波段方向发射率的体散射核函数，其数值为双向反射分布函数的体散射核函数半球积分值，Ifgeo表示地表窄波段方向发射率的几何光学核函数，其数值为双向反射分布函数的几何光学核函数半球积分值， 是第i波段对应的核函数权重系数，分别代表了各向同性核函数、体散射核函数、几何光学核函数所占的比例权重。
5.根据权利要求4所述的获取像元尺度地表宽波段半球发射率的方法，其特征在于：所述地表窄波段方向发射率的体散射核函数Ifvol的参数化表达结果如公式⑧所示：
其中，Ifvol表示地表窄波段方向发射率的体散射核函数，fvol表示双向反射分布函数体散射核函数，θv表示观测天顶角，θs表示太阳天顶角，表示观测方向与太阳方向间的相对方位角，cos和sin分别是余弦和正弦函数符号，∫代表积分符号，d代表微分符号。
6.根据权利要求5所述的获取像元尺度地表宽波段半球发射率的方法，其特征在于：所述地表窄波段方向发射率的几何光学核函数Ifgeo的参数化表达结果如公式⑨所示：
其中，Ifgeo表示地表窄波段方向发射率的几何光学核函数，fgeo表示双向反射分布函数几何光学核函数，θv表示观测天顶角，θs表示太阳天顶角，表示观测方向与太阳方向间的相对方位角，cos和sin分别表示余弦和正弦函数符号，∫代表积分符号，d代表微分符号。
7.根据权利要求6所述的获取像元尺度地表宽波段半球发射率的方法，其特征在于：步骤3)中所述地表窄波段方向发射率的核函数权重系数逐像元率定过程为：依照地表窄波段方向发射率的参数化表达结果，并利用步骤1)获取的地表窄波段方向发射率以及观测天顶角数据，逐像元构造超定方程组，超定方程组如公式⑩所示：
Ei＝KWi ⑩
其中，Ei表示地表窄波段方向发射率对应的矩阵，K表示地表窄波段方向发射率核函数值矩阵，Wi表示核函数权重系数矩阵；如式子所示Ei、K、Wi分别表示为：
其中，表示第i个波段顺序号第t天的地表窄波段方向发射率，θv,t表示顺序号第t天的观测天顶角；Ifvol,t表示顺序号第t天的地表窄波段方向发射率的体散射核函数值，由公式⑧计算而得；Ifgeo,t表示顺序号第t天的地表窄波段方向发射率的几何光学核函数值，由公式⑨计算而得；表示第i个波段对应的核函数权重系数，分别代表了各向同性核函数、体散射核函数和几何光学核函数所占的比例权重；
根据公式⑩，并利用最小二乘数学法，获得核函数权重系数的唯一解，如式子所示：
其中，表示估算的第i个波段最优核函数权重系数矩阵，K表示地表窄波段方向发射率核函数值矩阵，Ei表示第i个波段地表窄波段方向发射率对应的矩阵，T代表了矩阵的转置符号,上标-1代表了矩阵的求逆符号。
8.根据权利要求7所述的获取像元尺度地表宽波段半球发射率的方法，其特征在于：步骤4)中地表窄波段半球发射率的参数化表达结果如公式所示：
其中，表示第i个波段的地表窄波段半球发射率，表示第i波段对应的核函数权重系数，分别代表了各向同性核函数、体散射核函数和几何光学核函数所占的比例权重。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201811214646.4A CN109214122B (zh) | 2018-10-18 | 2018-10-18 | 一种获取像元尺度地表宽波段半球发射率的方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201811214646.4A CN109214122B (zh) | 2018-10-18 | 2018-10-18 | 一种获取像元尺度地表宽波段半球发射率的方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN109214122A CN109214122A (zh) | 2019-01-15 |
| CN109214122B true CN109214122B (zh) | 2019-06-04 |
Family
ID=64980847
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201811214646.4A Active CN109214122B (zh) | 2018-10-18 | 2018-10-18 | 一种获取像元尺度地表宽波段半球发射率的方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN109214122B (zh) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111157120B (zh) * | 2020-01-10 | 2021-05-04 | 北京航空航天大学 | 一种具有空间连续性的地表温度仿真方法 |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102298150A (zh) * | 2011-05-23 | 2011-12-28 | 北京师范大学 | 全球陆表宽波段发射率反演方法及系统 |
| CN105373694A (zh) * | 2015-08-31 | 2016-03-02 | 北京师范大学 | Modis宽波段发射率和glass宽波段发射率的融合计算方法 |
| CN105930664A (zh) * | 2016-04-26 | 2016-09-07 | 中国科学院地理科学与资源研究所 | 一种从被动微波数据估算瞬时地表发射率的方法 |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| KR101404430B1 (ko) * | 2013-06-11 | 2014-06-10 | 서울시립대학교 산학협력단 | 적외선 영상을 이용한 지표온도감률 추정 방법 |
-
2018
- 2018-10-18 CN CN201811214646.4A patent/CN109214122B/zh active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102298150A (zh) * | 2011-05-23 | 2011-12-28 | 北京师范大学 | 全球陆表宽波段发射率反演方法及系统 |
| CN105373694A (zh) * | 2015-08-31 | 2016-03-02 | 北京师范大学 | Modis宽波段发射率和glass宽波段发射率的融合计算方法 |
| CN105930664A (zh) * | 2016-04-26 | 2016-09-07 | 中国科学院地理科学与资源研究所 | 一种从被动微波数据估算瞬时地表发射率的方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN109214122A (zh) | 2019-01-15 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Antoine et al. | Assessment of uncertainty in the ocean reflectance determined by three satellite ocean color sensors (MERIS, SeaWiFS and MODIS‐A) at an offshore site in the Mediterranean Sea (BOUSSOLE project) | |
| Li et al. | Satellite-derived land surface temperature: Current status and perspectives | |
| Li et al. | Land surface emissivity retrieval from satellite data | |
| Jain et al. | Accuracy assessment of MODIS, NOAA and IRS data in snow cover mapping under Himalayan conditions | |
| Cheng et al. | Estimating the broadband longwave emissivity of global bare soil from the MODIS shortwave albedo product | |
| Latifovic et al. | A comparison of BRDF models for the normalization of satellite optical data to a standard sun-target-sensor geometry | |
| Alken et al. | Swarm equatorial electric field chain: First results | |
| Mironov et al. | Statistical characterization of short wind waves from stereo images of the sea surface | |
| Liu et al. | A long-term dataset of lake surface water temperature over the Tibetan Plateau derived from AVHRR 1981–2015 | |
| Wu et al. | Vertical distributions and relationships of cloud occurrence frequency as observed by MISR, AIRS, MODIS, OMI, CALIPSO, and CloudSat | |
| Li et al. | Cloud detection and classification algorithms for Himawari-8 imager measurements based on deep learning | |
| Jiao et al. | Evaluation of four sky view factor algorithms using digital surface and elevation model data | |
| Wang et al. | Combined pattern matching and feature tracking for Bohai Sea ice drift detection using Gaofen-4 imagery | |
| Zhang et al. | Development of a high spatiotemporal resolution cloud-type classification approach using Himawari-8 and CloudSat | |
| CN109214122B (zh) | 一种获取像元尺度地表宽波段半球发射率的方法 | |
| Hulley et al. | A new methodology for cloud detection and classification with ASTER data | |
| Luo et al. | Surface bidirectional reflectance and albedo properties derived using a land cover–based approach with Moderate Resolution Imaging Spectroradiometer observations | |
| Ma et al. | Regionalization of land surface heat fluxes over the heterogeneous landscape: from the Tibetan Plateau to the Third Pole region | |
| Ma et al. | Application of an LAI inversion algorithm based on the unified model of canopy bidirectional reflectance distribution function to the Heihe River Basin | |
| Oguro et al. | Comparisons of Brightness Temperatures of Landsat‐7/ETM+ and Terra/MODIS around Hotien Oasis in the Taklimakan Desert | |
| Wang et al. | Study on Remote Sensing of Water Depths Based on BP Artificial Neural Network. | |
| Gao et al. | A cloud detection algorithm over land based on the polarized characteristics difference between cloudless and cloud targets | |
| Tian et al. | Improving MODIS aerosol estimates over land with the surface BRDF reflectances using the 3-D discrete cosine transform and RossThick-LiSparse models | |
| Zhao et al. | Surface temperature retrieval from Gaofen-5 observation and its validation | |
| Lu et al. | Numerical simulation-aided MODIS capture of sediment transport for the Bohai Sea in China |
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 | ||
| CB03 | Change of inventor or designer information | ||
| CB03 | Change of inventor or designer information |
Inventor after: Wu Hua Inventor before: Wu Hua Inventor before: Li Zhaoliang Inventor before: Xue Changdi |
|
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| TR01 | Transfer of patent right | ||
| TR01 | Transfer of patent right |
Effective date of registration: 20240202 Address after: Room 101, Building 1, No. 18 Daoyuan Road, High tech Zone, Suzhou City, Jiangsu Province, 215011 Patentee after: SUZHOU ZHONGKETIANQI REMOTE SENSING TECHNOLOGY CO.,LTD. Country or region after: China Address before: 100101 Beijing city Chaoyang District Datun Road No. 11 Patentee before: Institute of Geographic Sciences and Natural Resources Research, CAS Country or region before: China |