CN117554300A - 山地地表反照率站点观测遥感空间降尺度方法 - Google Patents
山地地表反照率站点观测遥感空间降尺度方法 Download PDFInfo
- Publication number
- CN117554300A CN117554300A CN202410035742.1A CN202410035742A CN117554300A CN 117554300 A CN117554300 A CN 117554300A CN 202410035742 A CN202410035742 A CN 202410035742A CN 117554300 A CN117554300 A CN 117554300A
- Authority
- CN
- China
- Prior art keywords
- observation
- radiometer
- radiation
- pixel
- site
- 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.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 34
- 230000005855 radiation Effects 0.000 claims abstract description 162
- 230000004907 flux Effects 0.000 claims abstract description 82
- 238000000354 decomposition reaction Methods 0.000 claims abstract description 20
- 238000004458 analytical method Methods 0.000 claims abstract description 14
- 238000002310 reflectometry Methods 0.000 claims description 50
- 238000012937 correction Methods 0.000 claims description 21
- 230000005540 biological transmission Effects 0.000 claims description 20
- 230000004044 response Effects 0.000 claims description 20
- 238000004088 simulation Methods 0.000 claims description 15
- 238000011144 upstream manufacturing Methods 0.000 claims description 5
- 230000000007 visual effect Effects 0.000 claims description 4
- 238000004364 calculation method Methods 0.000 claims description 3
- 230000009467 reduction Effects 0.000 claims description 3
- 238000012876 topography Methods 0.000 abstract description 10
- 238000012795 verification Methods 0.000 description 8
- 230000000694 effects Effects 0.000 description 5
- 238000011161 development Methods 0.000 description 4
- 238000010586 diagram Methods 0.000 description 4
- 238000007689 inspection Methods 0.000 description 3
- 230000001788 irregular Effects 0.000 description 3
- 230000008569 process Effects 0.000 description 3
- 230000008859 change Effects 0.000 description 2
- 238000002474 experimental method Methods 0.000 description 2
- 238000005259 measurement Methods 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 238000012545 processing Methods 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 230000001360 synchronised effect Effects 0.000 description 2
- 238000012952 Resampling Methods 0.000 description 1
- 238000010521 absorption reaction Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000000701 chemical imaging Methods 0.000 description 1
- 238000000605 extraction Methods 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 238000010606 normalization Methods 0.000 description 1
- 230000003595 spectral effect Effects 0.000 description 1
- 238000001228 spectrum Methods 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
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
-
- 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
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
- G01S19/38—Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system
- G01S19/39—Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system the satellite radio beacon positioning system transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
-
- 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
- G01N2021/1793—Remote sensing
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
- Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
- Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation
Landscapes
- Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Remote Sensing (AREA)
- Radar, Positioning & Navigation (AREA)
- Analytical Chemistry (AREA)
- General Health & Medical Sciences (AREA)
- Biochemistry (AREA)
- Immunology (AREA)
- Pathology (AREA)
- Spectroscopy & Molecular Physics (AREA)
- Chemical & Material Sciences (AREA)
- Life Sciences & Earth Sciences (AREA)
- Health & Medical Sciences (AREA)
- Computer Networks & Wireless Communication (AREA)
- Geophysics And Detection Of Objects (AREA)
- Photometry And Measurement Of Optical Pulse Characteristics (AREA)
Abstract
本发明属于遥感技术领域,涉及山地地表反照率站点观测遥感空间降尺度方法。本发明基于辐射计观测特性,通过视域分析方法提取山地辐射计站点观测足迹,并基于地表反射特性计算像元反射辐射贡献度,分解站点观测上行短波辐射通量,基于局地地形特征分解山地地表站点观测下行短波辐射,得到高分辨率像元尺度地表反照率。本发明提出的基于辐射计观测特性、地表反射特性及局地地形特征的山地地表反照率站点观测遥感空间降尺度方法能够有效实现山地地表反照率辐射计站点观测降尺度,得到较高精度高分辨率像元尺度地表反照率。
Description
技术领域
本发明属于遥感技术领域,具体而言,涉及山地地表反照率站点观测遥感空间降尺度方法。
背景技术
地表反照率是决定地表能量平衡的重要参数,太阳辐射计的站点观测是获得地表反照率的重要手段,也是开展地表反照率遥感产品验证的重要数据源之一。传统地表反照率站点观测主要用于验证低空间分辨率卫星遥感产品精度,在高空间分辨率遥感产品验证中存在突出的尺度不匹配问题。此外,山地地表反照率站点观测范围显著区别于平坦地表,且随着辐射计架设高度的变化,其观测足迹往往大于高分辨率像元空间分辨率。随着中高分辨率地表反照率产品算法发展和产品生产的快速发展,山区高分辨率像元尺度地表反照率参考“真值”数据已成为真实性检验的迫切需求。
在平坦均一地表条件下,辐射计站点观测没有尺度效应,能够直接应用于像元尺度地表反照率验证。然而,真实地表很难满足平坦均一条件,异质地表或山地辐射计站点观测与粗分辨率像元之间存在明显尺度差异。
一方面,由于山地地形起伏,站点观测的足迹范围受地形影响呈现不规则特征;另一方面,不同站点观测架设高度的不同导致站点观测覆盖多个高分辨率遥感像元。山地地表反照率站点观测仅能代表仪器所覆盖范围内的反照率,直接将其应用于遥感产品真实性检验存在误差,验证结果不能代表遥感产品的真实情况。因此在山地复杂地形情况下,特定站点观测条件得到的地表反照率难以直接应用于山地地表反照率遥感反演模型的发展与精度验证。如何从站点观测足迹地表反照率获得准确高分辨率像元尺度地表反照率,支撑高分辨率山地遥感产品反演模型发展与验证成为难点问题。
发明内容
为了解决上述技术问题,本发明基于辐射计观测特性、地表反射特性及局地地形特征,提出了一种山地地表反照率站点观测遥感空间降尺度模型,并利用模型模拟的方式验证了降尺度方法的效果。为解决山地地表反照率尺度问题、高分辨率地表反照率遥感产品验证的问题,本发明提供山地地表反照率站点观测遥感空间降尺度方法。
第一方面,本发明提供了山地地表反照率站点观测遥感空间降尺度方法,包括:
获取通量观测站点的辐射观测数据、卫星遥感反射率数据、数字高程模型数据;
根据辐射计观测特性与数字高程模型数据,利用视域分析方法提取山地辐射计站点观测足迹;
根据辐射计站点观测足迹各像元观测角度,计算辐射计观测余弦响应度,并对辐射计观测短波上行辐照度进行校正,得到校正后的辐射计观测上行短波辐射通量;
利用地形辐射校正中的余弦校正模型的校正系数作为下行短波辐射分解的权重因子;
基于辐射计观测下行辐照度,计算辐射计观测下行短波辐射通量;
根据辐射计观测下行短波辐射通量与下行短波辐射分解的权重因子,确定像元短波下行辐射通量;
利用卫星遥感反射率数据,计算辐射计观测范围内像元反射率在短波波段内的积分;
计算各像元反射率波段积分值占观测范围内所有像元反射率波段积分值总和的比例,确定像元反射辐射贡献度;
根据像元反射辐射贡献度与校正后的辐射计观测上行短波辐射通量,确定各像元上行短波辐射通量;
根据像元尺度上行辐射通量与像元尺度下行辐射通量的比值,得到降尺度后的目标像元地表反照率。
在上述技术方案的基础上,本发明还可以做如下改进。
进一步,站点观测上行短波辐射通量与下行短波辐射通量的比值即为站点观测地表反照率,设为地表反照率,/>为辐射计观测下行太阳短波辐射,/>为辐射计观测上行太阳短波辐射,计算公式为:
。
进一步,卫星遥感反射率数据选择Sentinel-2多光谱高分辨率数据;数字高程模型数据由机载激光雷达数据生成。
进一步,根据辐射计观测特性与数字高程模型数据,利用视域分析方法提取山地辐射计站点观测足迹,包括:
设置视点参数、视场参数,通过求解视线与地形的交点得到视域范围;
设视域范围为,可视点为/>,视域范围为由所有以可视点/>为中心的像元组成的空间范围,观察点为/>,目标点分别为/>,观察点的高程为,目标点的高程为/>,/>与/>视线间的任意地形点为点/>,点/>的高程为/>,观察点/>与目标点/>之间的水平距离为/>,点/>与目标点/>之间的水平距离为/>,不遮挡观察点与目标点之间视线的点/>的最高高程为/>,点/>沿垂直于水平面方向与视线交点的高程为/>,则:
;
。
进一步,根据辐射计站点观测足迹各像元观测角度,计算辐射计观测余弦响应度,并对辐射计观测短波上行辐照度进行校正,包括:辐射计观测余弦响应为对于太阳辐射计的任意观测天顶角,设观测天顶角为,设站点观测足迹辐射计余弦响应度为/>,辐射计观测像元/>中心位置相对垂直方向的角度为/>,观测角度为/>,像元/>在/>观测角度下对应的辐射计余弦响应误差为/>,辐射计观测足迹内高分辨率像元个数为/>,辐射计观测上行短波辐照度为/>,余弦响应校正后上行短波辐照度为/>,则:
;
。
进一步,基于辐射计观测下行短波辐照度,计算辐射计观测下行短波辐射通量,包括:设水平地面等效辐射通量为,倾斜地表观测的辐射通量为/>,坡面局地入射天顶角为/>,太阳天顶角为/>,坡度为/>,太阳方位角为/>,坡向角度为/>,像元/>下行辐射占整个观测足迹下行辐射的比例为/>,分解得到的像元/>下行辐射为/>,辐射计观测下行短波辐照度为/>,像元面积为/>;
水平地面等效辐射通量为:;
坡面局地入射天顶角为:;
像元下行辐射占整个观测足迹下行辐射的比例为:/>;
分解得到的像元下行辐射为:/>。
进一步,根据像元反射辐射贡献度与校正后的辐射计观测上行短波辐射通量,确定各像元上行短波辐射通量,包括:设像元反射率波段积分为/>,波长为/>,/>表示卫星第1个短波波段对应波长,卫星短波波段数量为/>,像元/>地表反射率为/>,像元/>反射率积分占整个观测足迹比例为/>,像元/>由辐射计观测分解所得上行短波辐射通量为/>,余弦响应校正后上行短波辐照度为/>,则:
像元反射率波段积分为:/>;
像元𝑖反射率积分占整个观测范围比例为:;
像元由辐射计观测分解所得短波上行辐射通量为:/>。
进一步,还包括根据三维激光雷达数据构建真实场景,并采用三维辐射传输模型对像元尺度上行辐射通量与像元尺度下行辐射通量进行验证,包括:
模拟10 m分辨率类Sentinel反射率数据,同时利用三维辐射传输模型模拟场景总的上行辐射、下行辐射及10 m分辨率像元尺度上行辐射、像元尺度下行辐射;分解场景总的上行辐射通量,计算得到10 m分辨率像元尺度地表反照率,并对模拟得到10 m分辨率像元尺度地表反照率进行比较,计算均方根误差和相关系数平方对高分辨率像元尺度地表反照率进行验证。
进一步,计算均方根误差和相关系数平方,包括:
设模型模拟10 m分辨率像元尺度参考地表反照率为,降尺度后的目标地表反照率为/>,像元/>模型模拟参考地表反照率均值为/>,像元数量为/>,均方根误差为/>,相关系数平方为/>,则:
;
。
本发明的有益效果是:本发明提出的基于辐射计观测特性、地表反射特性及局地地形特征的山地地表反照率站点观测遥感空间降尺度方法效果良好,能够有效实现山地站点观测降尺度,得到较高精度高分辨率像元尺度地表反照率。
附图说明
图1为本发明实施例提供的山地地表反照率站点观测遥感空间降尺度方法的原理图;
图2为本发明实施例提供的山地地表反照率站点观测遥感空间降尺度方法的具体实施方式的原理框图;
图3为太阳天顶角、坡面局地入射天顶角与坡度的关系示意图;
图4为三维辐射传输模型模拟地表反照率与降尺度地表反照率对比的散点密度图。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚,下面将结合本发明实施例中的附图,对本发明实施例中的技术方案进行清楚、完整地描述,显然,所描述的实施例是本发明一部分实施例,而不是全部的实施例。通常在此处附图中描述和示出的本发明实施例的组件可以以各种不同的配置来布置和设计。
作为一个实施例,如附图1所示,为解决上述技术问题,本实施例提供山地地表反照率站点观测遥感空间降尺度方法,包括:
获取通量观测站点的辐射观测数据、卫星遥感反射率数据、数字高程模型数据;
根据辐射计观测特性与数字高程模型数据,利用视域分析方法提取山地辐射计站点观测足迹;
根据辐射计站点观测足迹各像元观测角度,计算辐射计观测余弦响应度,并对辐射计观测短波上行辐照度进行校正,得到校正后的辐射计观测上行短波辐射通量;
利用地形辐射校正中的余弦校正模型的校正系数作为下行短波辐射分解的权重因子;
基于辐射计观测下行辐照度,计算辐射计观测下行短波辐射通量;
根据辐射计观测下行短波辐射通量与下行短波辐射分解的权重因子,确定像元短波下行辐射通量;
利用卫星遥感反射率数据,计算辐射计观测范围内像元反射率在短波波段内的积分;
计算各像元反射率波段积分值占观测范围内所有像元反射率波段积分值总和的比例,确定像元反射辐射贡献度;
根据像元反射辐射贡献度与校正后的辐射计观测上行短波辐射通量,确定各像元上行短波辐射通量;
根据像元尺度上行辐射通量与像元尺度下行辐射通量的比值,得到降尺度后的目标像元地表反照率。
山地地表反照率站点观测遥感空间降尺度方法的具体实施方式的原理框图如附图2所示。
在实际应用过程中,首先基于辐射计观测特性和高分辨率DEM数据,利用视域分析方法提取山地辐射计站点观测足迹,然后对辐射计余弦响应特性进行校正,基于地表反射特性,对站点观测上行辐射通量进行分解,同时引入地形校正余弦模型,基于局地地形特征分解站点观测短波下行辐射通量,分解得到的像元尺度上下行辐射通量的比值即为地表反照率,由此实现站点观测地表反照率降尺度。
可选的,站点观测上行短波辐射通量与下行短波辐射通量的比值即为站点观测地表反照率,设为地表反照率,/>为辐射计观测下行太阳短波辐射,/>为辐射计观测上行太阳短波辐射,计算公式为:
。
可选的,卫星遥感反射率数据选择Sentinel-2多光谱高分辨率数据;数字高程模型数据由机载激光雷达数据生成。
具体的,获取通量观测站点的辐射观测数据、卫星遥感反射率数据、数字高程模型数据的过程中,优先选择当地正午时刻(约为下午13:00时左右)前后半小时,即12:30、13:00、13:30观测的上行短波辐照度、下行短波辐照度,计算平均值得到当日站点观测的上行短波辐照度、下行行短波辐照度。选择正午时刻一方面可以避免太阳天顶角对地表反照率的影响,另一方面正午时刻与多数卫星过境时刻相近,一定程度减小地表反照率遥感产品真实性检验误差。采用高分辨率卫星遥感地表反射率数据表征观测足迹内的地表反射特征。
卫星遥感反射率数据选择Sentinel-2A/B多光谱高分辨率数据,Sentinel-2卫星是高分辨率多光谱成像卫星,分为Sentinel-2A和Sentinel-2B两颗相同的卫星。Sentinel-2卫星在太阳同步轨道上以180 km的相位相互同步,平均轨道高度为786 km。Sentinel-2卫星各携带一枚多光谱仪器(MSI),包括13个光谱波段,地面分辨率分别有10 m、20 m和60 m。选择晴空过境时刻Sentinel-2A/B数据,将Sentinel反射率数据进行时空谱归一化处理,空间分辨率统一为10 m。
数字高程数据由空-天-地立体观测实验机载激光雷达数据生成,得到的数字高程模型DEM原始数据分辨率为0.1 m。机载激光雷达数据记录了数字地表高程信息,从机载激光雷达数据中提取地面点云,通过普通克里金插值得到高精度DEM数据,获取对应的数字地形信息,能够反映地面的精细地形特征,经过重采样处理后空间分辨率与像元几何位置同Sentinel高分辨率反射率数据保持一致。
辐射计观测范围分析和提取是站点观测降尺度的前提和基础,在平坦地表,基于一定视场角的辐射计观测范围可以较为方便被刻画,但是在复杂山地,由于显著的地形特征(海拔、坡度、坡向)以及临近像元遮挡等因素的影响,山地辐射计观测呈现出不规则的观测足迹,因此根据辐射计观测特性分析其在山地的特定观测足迹。
基于视域分析方法和高分辨率数字高程数据,分析得到了三个站点区域的辐射计观测范围。总体可以看出,三个站点观测范围均呈现出不规则足迹,且均覆盖了不同的地表类型,显著区别于平坦地表。山地局地地形特征和辐射计架设高度显著影响山地辐射计站点观测范围。
具体的,站点通量塔辐射计站点观测范围试验,三个站点中,10 m塔架的坡下方由相对平缓地表构成,该塔观测足迹为扇形,区域内覆盖地表主要为灌丛、林地、草地、道路等地表覆盖。统计得出,该塔辐射计观测足迹面积约为67492.3 m2,对应Sentinel-2卫星10米分辨率像元约674个。对比可以看出,30米通量塔观测足迹为细长扇形型,主要受两侧高大山体遮蔽影响,统计得出,该塔辐射计观测足迹面积约为77716.9 m2。对应Sentinel-2卫星10米分辨率像元约777个。75米塔架在三座塔中架设高度最高,且观测范围最大,观测足迹呈椭圆形并受两侧山地影响,统计得出,该塔辐射计观测足迹面积约为154649 m2,对应Sentinel-2卫星10米分辨率像元约1546个。通过试验,该塔的观测范围与理论观测足迹有很大差异。
综上可以看出,复杂山地辐射计观测足迹显著受到观测塔的架设位置、架设高度和局地地形影响,且观测足迹内地物贡献会随着高度变化而有显著差异。因此,在开展山地地表反射率分析时,需要考虑站点观测足迹。
基于获取的高分辨率数字高程数据,利用GIS视域分析的方法,模拟辐射计在山地区域的观测范围。设置视点参数、视场参数,通过求解视线与地形的交点得到视域范围,一般的,辐射计视场角小于180°(如150°)。
可选的,根据辐射计观测特性与数字高程模型数据,利用视域分析方法提取山地辐射计站点观测足迹,包括:
设置视点参数、视场参数,通过求解视线与地形的交点得到视域范围;
设视域范围为,可视点为/>,视域范围为由所有以可视点/>为中心的像元组成的空间范围,观察点为/>,目标点分别为/>,观察点的高程为,目标点的高程为/>,/>与/>视线间的任意地形点为点/>,点/>的高程为/>,观察点/>与目标点/>之间的水平距离为/>,点/>与目标点/>之间的水平距离为/>,不遮挡观察点与目标点之间视线的点/>的最高高程为/>,点/>沿垂直于水平面方向与视线交点的高程为/>,则:
;
。
“余弦响应”用于描述测量太阳辐照度的传感器特性,指对于太阳辐射计的任意观测天顶角,理想辐射计在观测天顶角/>方向测量的辐照度是其在垂直方向上测量的辐照度与该角度的余弦的乘积,是辐射计观测的基本特性。然而,实际情况下辐射计观测与理想辐射计余弦响应特性并不完全相同,具有余弦响应误差。辐射计观测上行短波辐照度是整个观测足迹内上行辐照度的综合,因此对于辐射计观测足迹内不同的亚足迹高分辨率像元,特定位置的辐射计观测不同位置像元反射辐射时的观测角度不同,其观测值受到辐射计余弦响应特性的影响,因此需要综合辐射计余弦响应特性及余弦响应误差,对辐射计观测上行短波辐照度进行余弦响应校正,以获取垂直观测方向观测足迹内更准确的上行短波辐照度。
可选的,根据辐射计站点观测足迹各像元观测角度,计算辐射计观测余弦响应度,并对辐射计观测短波上行辐照度进行校正,包括:辐射计观测余弦响应为对于太阳辐射计的任意观测天顶角,设观测天顶角为,设站点观测足迹辐射计余弦响应度为/>,辐射计观测像元/>中心位置相对垂直方向的角度为/>,观测角度为/>,像元/>在/>观测角度下对应的辐射计余弦响应误差为/>,辐射计观测足迹内高分辨率像元个数为/>,辐射计观测上行短波辐照度为/>,余弦响应校正后上行短波辐照度为/>,则:
;
。
山地复杂地形显著改变辐射传输过程,现有研究表明山地下行短波辐射空间分布受到区域地形的影响,忽略地形影响会造成下行短波辐射估算产生较大误差。因此,在进行站点观测下行辐射分解时需要考虑地形的影响,本发明引入典型地形辐射校正中的余弦校正模型,该模型可将山地观测短波辐射转换为水平地表辐射,实现山地辐射地形校正,因此该模型一定程度可表征山地下行辐射空间分布情况。综上,本发明将该模型校正系数作为下行短波辐射分解的权重因子,计算得到像元/>下行短波辐射占观测足迹短波下行辐射比例,基于辐射计观测下行辐照度计算观测范围下行辐射通量,与像元下行辐射比例相乘即可得到像元下行短波辐射通量。
可选的,基于辐射计观测下行短波辐照度,计算辐射计观测下行短波辐射通量,包括:设水平地面等效辐射通量为,倾斜地表观测的辐射通量为/>,坡面局地入射天顶角为/>,太阳天顶角为/>,坡度为/>,太阳方位角为/>,坡向角度为/>,像元/>下行辐射占整个观测足迹下行辐射的比例为/>,分解得到的像元/>下行辐射为/>,辐射计观测下行短波辐照度为/>,像元面积为/>;
水平地面等效辐射通量为:;
坡面局地入射天顶角为:;
像元下行辐射占整个观测足迹下行辐射的比例为:/>;
分解得到的像元下行辐射为:/>。
太阳天顶角、坡面局地入射天顶角/>与坡度/>的关系示意图如附图3所示。
辐射计观测到的地表上行辐射是观测足迹范围内所有地物反射太阳辐射的综合结果。为了获得观测足迹内对应高空间分辨率像元的地表反照率数据,将辐射计观测足迹内校正后的上行短波辐射通量分解到各像元尺度。考虑到卫星遥感观测获得的宽波段反射率数据与地表反照率具有较好的相关关系,本发明提出一种利用宽波段反射率数据作为像元短波波段内的反射辐射贡献度的分解模型。利用高分辨率Sentinel反射率数据,计算辐射计观测足迹内像元反射率在短波波段内的积分,用于表征该像元在短波波段内对太阳辐射的反射能力,然后,计算各像元反射率波段积分值占观测范围内所有像元反射率波段积分值总和的比例,表示该像元相对整个观测范围的反射辐射贡献度。像元反射辐射贡献度与经过校正后的辐射计观测上行短波辐照度计算得到的上行短波辐射通量相乘即可得到各像元上行短波辐射通量。
可选的,根据像元反射辐射贡献度与校正后的辐射计观测上行短波辐射通量,确定各像元上行短波辐射通量,包括:设像元反射率波段积分为/>,波长为/>,/>表示卫星第1个短波波段对应波长,卫星短波波段数量为/>,像元/>地表反射率为/>,像元/>反射率积分占整个观测足迹比例为/>,像元/>由辐射计观测分解所得上行短波辐射通量为/>,余弦响应校正后上行短波辐照度为/>,则:
像元反射率波段积分为:/>;
像元𝑖反射率积分占整个观测范围比例为:;
像元由辐射计观测分解所得短波上行辐射通量为:/>。
可选的,还包括根据三维激光雷达数据构建真实场景,并采用三维辐射传输模型对像元尺度上行辐射通量与像元尺度下行辐射通量进行验证,包括:
模拟10 m分辨率类Sentinel反射率数据,同时利用三维辐射传输模型模拟场景总的上行辐射、下行辐射及10 m分辨率像元尺度上行辐射、像元尺度下行辐射;分解场景总的上行辐射通量,计算得到10 m分辨率像元尺度地表反照率,并对模拟得到10 m分辨率像元尺度地表反照率进行比较,计算均方根误差和相关系数平方对高分辨率像元尺度地表反照率进行验证。
三维辐射传输模型是基于光线追踪三维真实结构辐射传输模型,它可以模拟入射光在场景中的传输过程(吸收、反射以及透射),并输出相应的模拟数据(例如反射率、反照率、FPAR等)。三维辐射传输模型充分利用了光线追踪的前向追踪模式模拟能量平衡问题以及后向追踪模式模拟大尺度(公里级)遥感影像,从而实现在同一三维辐射传输模型中实现多种遥感数据的模拟。
本发明获取激光雷达得到的真实场景的三维地形数据,构建一个相对真实的三维场景,场景元素包括树木、地形与草地。通过输入Sentinel-2波段设置,模拟了10 m分辨率类Sentinel反射率数据,同时利用三维辐射传输模型模拟场景总的上下行辐射及10 m分辨率像元尺度上下行辐射。分解场景总上行辐射通量,计算得到10 m分辨率像元尺度地表反照率,并与模拟得到10 m分辨率像元尺度地表反照率进行比较,以验证提出的方法精度、可靠性。验证指标选择均方根误差和相关系数平方/>。
可选的,计算均方根误差和相关系数平方,包括:
设模型模拟10 m分辨率像元尺度参考地表反照率为,降尺度后的目标地表反照率为/>,像元/>模型模拟参考地表反照率均值为/>,像元数量为/>,均方根误差为,相关系数平方为/>,则:
;
。
在实际应用过程中,通过模拟像元尺度的上下行短波辐射、地表反照率,三维辐射传输模型模拟地表反照率与降尺度地表反照率对比的散点密度图如附图4所示。可以看出,三维辐射传输模型模拟结果很好的反映了地表不同地物的反照率差异特征,本发明山地地表反照率站点观测遥感空间降尺度方法与三维辐射传输模型模拟像元尺度地表反照率目视结果具有非常好的一致性。
通过定量分析三维辐射传输模型模拟与分解的相关性,结果表明,显示两者的相关系数平方R2为0.8339,线性拟合斜率约接近1:1线,均方根误差RMSE为0.0169,两者的一致性较好,差异较小,表明本发明提出的山地地表反照率站点观测遥感空间降尺度方法效果良好,能够实现山地站点观测降尺度,得到较高精度高分辨率像元尺度地表反照率。
本发明通过视域分析方法提取山地辐射计站点观测范围,并基于辐射计观测特性和地表反射特性计算高分辨率像元反射辐射贡献度,分解站点观测短波上行辐射通量,基于局地地形特征分解站点观测下行短波辐射通量,得到高分辨率像元尺度地表反照率。
由于复杂地形的影响,山地站点辐射计观测范围显著区别于平坦地表,在利用辐射计站点观测进行山地能量平衡研究时考虑了地形及尺度效应的影响;
本发明提出的基于辐射计观测特性、地表反射特性及局地地形特征的山地地表反照率站点观测遥感空间降尺度方法能够有效实现山地站点观测降尺度,利用三维辐射传输模型模拟的验证结果表明降尺度地表反照率能够实现山地站点观测降尺度,得到较高精度、高分辨率像元尺度地表反照率。
以上仅为本发明的优选实施例而已,并不用于限制本发明,对于本领域的技术人员来说,本发明可以有各种更改和变化。凡在本发明的精神和原则之内,所作的任何修改、等同替换、改进等,均应包含在本发明的保护范围之内。
Claims (9)
1.山地地表反照率站点观测遥感空间降尺度方法,其特征在于,包括:
获取通量观测站点的辐射观测数据、卫星遥感反射率数据、数字高程模型数据;
根据辐射计观测特性与数字高程模型数据,利用视域分析方法提取山地辐射计站点观测足迹;
根据辐射计站点观测足迹各像元观测角度,计算辐射计观测余弦响应度,并对辐射计观测短波上行辐照度进行校正,得到校正后的辐射计观测上行短波辐射通量;
利用地形辐射校正中的余弦校正模型的校正系数作为下行短波辐射分解的权重因子;
基于辐射计观测下行辐照度,计算辐射计观测下行短波辐射通量;
根据辐射计观测下行短波辐射通量与下行短波辐射分解的权重因子,确定像元短波下行辐射通量;
利用卫星遥感反射率数据,计算辐射计观测范围内像元反射率在短波波段内的积分;
计算各像元反射率波段积分值占观测范围内所有像元反射率波段积分值总和的比例,确定像元反射辐射贡献度;
根据像元反射辐射贡献度与校正后的辐射计观测上行短波辐射通量,确定各像元上行短波辐射通量;
根据像元尺度上行辐射通量与像元尺度下行辐射通量的比值,得到降尺度后的目标像元地表反照率。
2.根据权利要求1所述山地地表反照率站点观测遥感空间降尺度方法,其特征在于,站点观测上行短波辐射通量与下行短波辐射通量的比值即为站点观测地表反照率,设为地表反照率,/>为辐射计观测下行太阳短波辐射,/>为辐射计观测上行太阳短波辐射,计算公式为:
。
3.根据权利要求1所述山地地表反照率站点观测遥感空间降尺度方法,其特征在于,卫星遥感反射率数据选择Sentinel-2多光谱高分辨率数据;数字高程模型数据由机载激光雷达数据生成。
4.根据权利要求1所述山地地表反照率站点观测遥感空间降尺度方法,其特征在于,根据辐射计观测特性与数字高程模型数据,利用视域分析方法提取山地辐射计站点观测足迹,包括:
设置视点参数、视场参数,通过求解视线与地形的交点得到视域范围;
设视域范围为,可视点为/>,视域范围为由所有以可视点/>为中心的像元组成的空间范围,观察点为/>,目标点分别为/>,观察点的高程为/>,目标点的高程为/>,/>与/>视线间的任意地形点为点/>,点/>的高程为/>,观察点/>与目标点/>之间的水平距离为/>,点/>与目标点/>之间的水平距离为,不遮挡观察点与目标点之间视线的点/>的最高高程为/>,点/>沿垂直于水平面方向与视线交点的高程为/>,则:
;
。
5.根据权利要求1所述山地地表反照率站点观测遥感空间降尺度方法,其特征在于,根据辐射计站点观测足迹各像元观测角度,计算辐射计观测余弦响应度,并对辐射计观测短波上行辐照度进行校正,包括:辐射计观测余弦响应为对于太阳辐射计的任意观测天顶角,设观测天顶角为,设站点观测足迹辐射计余弦响应度为/>,辐射计观测像元/>中心位置相对垂直方向的角度为/>,观测角度为/>,像元/>在/>观测角度下对应的辐射计余弦响应误差为/>,辐射计观测足迹内高分辨率像元个数为/>,辐射计观测上行短波辐照度为/>,余弦响应校正后上行短波辐照度为/>,则:
;
。
6.根据权利要求1所述山地地表反照率站点观测遥感空间降尺度方法,其特征在于,基于辐射计观测下行短波辐照度,计算辐射计观测下行短波辐射通量,包括:设水平地面等效辐射通量为,倾斜地表观测的辐射通量为/>,坡面局地入射天顶角为/>,太阳天顶角为,坡度为/>,太阳方位角为/>,坡向角度为/>,像元/>下行辐射占整个观测足迹下行辐射的比例为/>,分解得到的像元/>下行辐射为/>,辐射计观测下行短波辐照度为/>,像元面积为;
水平地面等效辐射通量为:;
坡面局地入射天顶角为:;
像元下行辐射占整个观测足迹下行辐射的比例为:/>;
分解得到的像元下行辐射为:/>。
7.根据权利要求1所述山地地表反照率站点观测遥感空间降尺度方法,其特征在于,根据像元反射辐射贡献度与校正后的辐射计观测上行短波辐射通量,确定各像元上行短波辐射通量,包括:设像元反射率波段积分为/>,波长为/>,/>表示卫星第1个短波波段对应波长,卫星短波波段数量为/>,像元/>地表反射率为/>,像元/>反射率积分占整个观测足迹比例为/>,像元/>由辐射计观测分解所得上行短波辐射通量为/>,余弦响应校正后上行短波辐照度为/>,则:
像元反射率波段积分为:/>;
像元𝑖反射率积分占整个观测范围比例为:;
像元由辐射计观测分解所得短波上行辐射通量为:/>。
8.根据权利要求1所述山地地表反照率站点观测遥感空间降尺度方法,其特征在于,还包括根据三维激光雷达数据构建真实场景,并采用三维辐射传输模型对像元尺度上行辐射通量与像元尺度下行辐射通量进行验证,包括:
模拟10 m分辨率类Sentinel反射率数据,同时利用三维辐射传输模型模拟场景总的上行辐射、下行辐射及10 m分辨率像元尺度上行辐射、像元尺度下行辐射;分解场景总的上行辐射通量,计算得到10 m分辨率像元尺度地表反照率,并对模拟得到10 m分辨率像元尺度地表反照率进行比较,计算均方根误差和相关系数平方对高分辨率像元尺度地表反照率进行验证。
9.根据权利要求1所述山地地表反照率站点观测遥感空间降尺度方法,其特征在于,计算均方根误差和相关系数平方,包括:
设模型模拟10 m分辨率像元尺度参考地表反照率为,降尺度后的目标地表反照率为/>,像元/>模型模拟参考地表反照率均值为/>,像元数量为/>,均方根误差为/>,相关系数平方为/>,则:
;
。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202410035742.1A CN117554300B (zh) | 2024-01-10 | 2024-01-10 | 山地地表反照率站点观测遥感空间降尺度方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202410035742.1A CN117554300B (zh) | 2024-01-10 | 2024-01-10 | 山地地表反照率站点观测遥感空间降尺度方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN117554300A true CN117554300A (zh) | 2024-02-13 |
CN117554300B CN117554300B (zh) | 2024-03-19 |
Family
ID=89821868
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202410035742.1A Active CN117554300B (zh) | 2024-01-10 | 2024-01-10 | 山地地表反照率站点观测遥感空间降尺度方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN117554300B (zh) |
Citations (9)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN104406686A (zh) * | 2014-12-10 | 2015-03-11 | 西北师范大学 | 复杂地形条件下太阳短波入射辐射估算方法 |
CN109885959A (zh) * | 2019-03-05 | 2019-06-14 | 中国科学院地理科学与资源研究所 | 一种地表温度鲁棒降尺度方法 |
CN110516816A (zh) * | 2019-08-30 | 2019-11-29 | 中国科学院、水利部成都山地灾害与环境研究所 | 基于机器学习的全天候地表温度生成方法及装置 |
CN111198162A (zh) * | 2020-01-09 | 2020-05-26 | 胡德勇 | 一种城区表面反射率遥感反演方法 |
CN111563228A (zh) * | 2020-05-07 | 2020-08-21 | 中国科学院、水利部成都山地灾害与环境研究所 | 基于地表入射短波辐射的山地地表反射率地形改正方法 |
GB2593592A (en) * | 2020-03-03 | 2021-09-29 | Nvidia Corp | Technique for Performing Bit-Linear Transformations |
WO2022166939A1 (zh) * | 2021-02-08 | 2022-08-11 | 南京农业大学 | 一种基于Sentinel-2卫星影像红边波段改进小麦生长早期叶面积指数估算的方法 |
CN115015147A (zh) * | 2022-04-21 | 2022-09-06 | 北京市遥感信息研究所 | 一种高空间分辨率高光谱热红外遥感影像模拟方法 |
CN117077437A (zh) * | 2023-10-12 | 2023-11-17 | 中科星图维天信科技股份有限公司 | 基于多源卫星的极区海表净辐射模型构建和确定方法 |
-
2024
- 2024-01-10 CN CN202410035742.1A patent/CN117554300B/zh active Active
Patent Citations (9)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN104406686A (zh) * | 2014-12-10 | 2015-03-11 | 西北师范大学 | 复杂地形条件下太阳短波入射辐射估算方法 |
CN109885959A (zh) * | 2019-03-05 | 2019-06-14 | 中国科学院地理科学与资源研究所 | 一种地表温度鲁棒降尺度方法 |
CN110516816A (zh) * | 2019-08-30 | 2019-11-29 | 中国科学院、水利部成都山地灾害与环境研究所 | 基于机器学习的全天候地表温度生成方法及装置 |
CN111198162A (zh) * | 2020-01-09 | 2020-05-26 | 胡德勇 | 一种城区表面反射率遥感反演方法 |
GB2593592A (en) * | 2020-03-03 | 2021-09-29 | Nvidia Corp | Technique for Performing Bit-Linear Transformations |
CN111563228A (zh) * | 2020-05-07 | 2020-08-21 | 中国科学院、水利部成都山地灾害与环境研究所 | 基于地表入射短波辐射的山地地表反射率地形改正方法 |
WO2022166939A1 (zh) * | 2021-02-08 | 2022-08-11 | 南京农业大学 | 一种基于Sentinel-2卫星影像红边波段改进小麦生长早期叶面积指数估算的方法 |
CN115015147A (zh) * | 2022-04-21 | 2022-09-06 | 北京市遥感信息研究所 | 一种高空间分辨率高光谱热红外遥感影像模拟方法 |
CN117077437A (zh) * | 2023-10-12 | 2023-11-17 | 中科星图维天信科技股份有限公司 | 基于多源卫星的极区海表净辐射模型构建和确定方法 |
Also Published As
Publication number | Publication date |
---|---|
CN117554300B (zh) | 2024-03-19 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
JP4425983B1 (ja) | 日射量の評価方法および評価装置 | |
Strahler et al. | MODIS BRDF/albedo product: algorithm theoretical basis document version 5.0 | |
Yin et al. | Optimization of L-band sea surface emissivity models deduced from SMOS data | |
CN108132220A (zh) | 林区机载推扫式高光谱影像的brdf归一化校正方法 | |
Gastellu-Etchegorry et al. | DART: A 3D model for remote sensing images and radiative budget of earth surfaces | |
Marzano et al. | Potential of high-resolution detection and retrieval of precipitation fields from X-band spaceborne synthetic aperture radar over land | |
Wang et al. | Modeling the angular effect of MODIS LST in urban areas: A case study of Toulouse, France | |
CN112798013B (zh) | 一种对光学载荷在轨绝对辐射定标结果进行验证的方法 | |
Liao et al. | Simplified vector-based model tailored for urban-scale prediction of solar irradiance | |
CN101876700B (zh) | 一种基于辐射度的复杂地形区域辐射传输模拟方法 | |
Reul et al. | Earth-viewing L-band radiometer sensing of sea surface scattered celestial sky radiation—Part II: Application to SMOS | |
CN104880701A (zh) | 一种星载传感器成像仿真方法及装置 | |
Ma et al. | Characterizing the three-dimensional spatiotemporal variation of forest photosynthetically active radiation using terrestrial laser scanning data | |
Corbett et al. | Accounting for the effects of sastrugi in the CERES clear-sky Antarctic shortwave angular distribution models | |
CN102073038B (zh) | 基于微小地形的遥感影像的地形校正的方法 | |
CN113091892B (zh) | 卫星遥感器在轨星上绝对辐射定标方法及系统 | |
Gastellu-Etchegorry et al. | Recent improvements in the dart model for atmosphere, topography, large landscape, chlorophyll fluorescence, satellite image inversion | |
CN117554300B (zh) | 山地地表反照率站点观测遥感空间降尺度方法 | |
CN113012276A (zh) | 基于辐射度的地表高分辨率光谱信息遥感反演方法 | |
CN116611352A (zh) | 极轨卫星地表温度时间可比性校正通用方法与系统 | |
Adavi et al. | Analyzing different parameterization methods in GNSS tomography using the COST benchmark dataset | |
Borfecchia et al. | Integrated GIS and remote sensing techniques to support pv potential assessment of roofs in urban areas | |
Miao et al. | Wet Tropospheric Correction Methods for Wide-Swath Altimeters | |
Schweitzer et al. | Validation of the background simulation model MATISSE: Comparing results with MODIS satellite images | |
Labarre et al. | MATISSE-v2. 0: new functionalities and comparison with MODIS satellite images |
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 |