CN113569386B - 卫星遥感夜光辐亮度的观测角度归一化方法 - Google Patents
卫星遥感夜光辐亮度的观测角度归一化方法 Download PDFInfo
- Publication number
- CN113569386B CN113569386B CN202110743139.5A CN202110743139A CN113569386B CN 113569386 B CN113569386 B CN 113569386B CN 202110743139 A CN202110743139 A CN 202110743139A CN 113569386 B CN113569386 B CN 113569386B
- Authority
- CN
- China
- Prior art keywords
- radiance
- luminous
- angle
- zenith
- noctilucent
- 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
- 238000010606 normalization Methods 0.000 title claims abstract description 26
- 238000012887 quadratic function Methods 0.000 claims abstract description 38
- 238000000034 method Methods 0.000 claims abstract description 27
- 230000000694 effects Effects 0.000 claims abstract description 23
- 238000011160 research Methods 0.000 claims description 19
- 239000011159 matrix material Substances 0.000 claims description 16
- 238000012216 screening Methods 0.000 claims description 11
- 230000008569 process Effects 0.000 claims description 6
- PEDCQBHIVMGVHV-UHFFFAOYSA-N Glycerine Chemical compound OCC(O)CO PEDCQBHIVMGVHV-UHFFFAOYSA-N 0.000 claims description 5
- 230000005855 radiation Effects 0.000 claims description 5
- 238000004364 calculation method Methods 0.000 claims description 4
- 238000006243 chemical reaction Methods 0.000 claims description 4
- 238000013139 quantization Methods 0.000 claims description 4
- 238000000354 decomposition reaction Methods 0.000 claims description 3
- 230000002159 abnormal effect Effects 0.000 claims description 2
- 230000001419 dependent effect Effects 0.000 claims description 2
- 238000005259 measurement Methods 0.000 abstract description 2
- 238000013179 statistical model Methods 0.000 description 9
- 238000012937 correction Methods 0.000 description 6
- 238000010586 diagram Methods 0.000 description 6
- 238000002310 reflectometry Methods 0.000 description 6
- 230000008859 change Effects 0.000 description 4
- 238000004458 analytical method Methods 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 2
- 230000007547 defect Effects 0.000 description 2
- 230000008030 elimination Effects 0.000 description 2
- 238000003379 elimination reaction Methods 0.000 description 2
- DFPOZTRSOAQFIK-UHFFFAOYSA-N S,S-dimethyl-beta-propiothetin Chemical compound C[S+](C)CCC([O-])=O DFPOZTRSOAQFIK-UHFFFAOYSA-N 0.000 description 1
- 238000007792 addition Methods 0.000 description 1
- 238000012512 characterization method Methods 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000005315 distribution function Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000005286 illumination Methods 0.000 description 1
- 238000003384 imaging method Methods 0.000 description 1
- 238000003331 infrared imaging Methods 0.000 description 1
- 239000004579 marble Substances 0.000 description 1
- 230000007246 mechanism Effects 0.000 description 1
- 238000005065 mining Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000012544 monitoring process Methods 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 230000001932 seasonal effect Effects 0.000 description 1
- 230000003595 spectral effect Effects 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 230000002123 temporal effect Effects 0.000 description 1
- 238000012876 topography Methods 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
- G01J1/00—Photometry, e.g. photographic exposure meter
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/02—Reliability analysis or reliability optimisation; Failure analysis, e.g. worst case scenario performance, failure mode and effects analysis [FMEA]
-
- 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
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Mathematical Physics (AREA)
- Mathematical Analysis (AREA)
- Data Mining & Analysis (AREA)
- Pure & Applied Mathematics (AREA)
- Mathematical Optimization (AREA)
- General Engineering & Computer Science (AREA)
- Computational Mathematics (AREA)
- Computing Systems (AREA)
- Geometry (AREA)
- Evolutionary Computation (AREA)
- Computer Hardware Design (AREA)
- Algebra (AREA)
- Spectroscopy & Molecular Physics (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- Photometry And Measurement Of Optical Pulse Characteristics (AREA)
Abstract
本发明涉及一种卫星遥感夜光辐亮度的观测角度归一化方法。首先建立夜光辐亮度与观测天顶角之间的天顶‑辐亮度二次(ZRQ)模型，采用最小二乘法求解二次函数的系数，然后基于ZRQ模型将观测到的夜光辐亮度分解为天顶角为0°时的夜光辐亮度与二次函数值的乘积，采用最小二乘法求解分解后的二次函数系数，并计算当前传感器获取的夜光辐亮度与二次函数值的比值，求解天顶角为0°时的夜光辐亮度，进而消除夜光辐亮度的角度效应。本发明基于ZRQ模型构建夜光辐亮度角度归一化模型，顾及人类活动和传感器测量误差等因素造成的时间序列波动，连续两次采用模型拟合求解参数，提高了角度归一化的精度，可实现对长时间覆盖范围的夜光辐亮度角度归一化。
Description
技术领域
本发明属于统计回归模型在夜光遥感影像各向异性消除领域的应用，特别是涉及一种卫星遥感夜光辐亮度的观测角度归一化方法。
背景技术
夜光遥感是获取夜间无云时地表发射的可见光信息，通过夜间灯光可以发现人类聚集区并记录地表灯光强度信息，进而可以通过夜光遥感影像间接反映居民在城镇中的行为活动与空间分布。深入挖掘夜光辐亮度时间序列的动态变化规律，对于准确监测城市化进程、电力供应以及评估自然灾害等重大突发事件有着重要的作用。
随着遥感技术的不断发展，夜光遥感产品逐渐丰富。特别是新一代夜光辐亮度数据NPP/VIIRS(National Polar-Orbiting Partnership's Visible Infrared ImagingRadiometer Suite)不仅具有检测微弱夜间灯光的能力，而且弥补了DMSP/OLS数据低空间分辨率与低辐射分辨率的缺点，能更多地应用在建立夜光辐亮度与人类活动模式经验关系的研究中。由于夜间灯光信号会受到杂散光、大气成分、季节性等噪声信号干扰，从而造成夜光辐亮度时间序列的波动。为此，NASA团队生产出了一套全球可用的、包含丰富时间与空间信息的黑色大理石夜光影像产品(VNP46)。该产品去除了受云污染的像素，并对受大气、地形、植被、雪、月球以及杂散光影响的夜光辐亮度进行了校正，使得产品的夜光辐亮度更为可靠。但是Li等学者发现，夜光辐亮度与卫星传感器观测天顶角存在非线性关系，特别是在高建筑物区域，受到城市三维形态的影响，夜光辐亮度会存在明显的冷点和热点效应。因此，在利用长时间序列的夜光影像进行研究时，为了更准确地反映城市的动态变化，消除夜光辐亮度的角度效应十分必要。
目前许多学者在光学各向异性的研究主要包括基于双向反射率分布函数(BRDF)表征地球表面反射率各向异性的研究与模拟城市上空光分布的光污染研究。地球表面反射光会受到地物的物理结构、传感器视角以及太阳光照角度变化的影响，产生地球表面反射率各向异性特征。Landsat卫星通过最低点一定范围的视角获取地球表面影像，会导致地表反射率产生较小的方向性影响。Roy等人基于全球高质量的MODIS BRDF产品，得出一组固定的BRDF光谱模型参数，从而完成Landsat数据地表反射率角度校正，但作为描述叶子空间分布格局的关键参数，聚集指数(CI)会受到BRDF模型与太阳天顶角的影响。Wei等人使用归一化冷热点差异(NDHD)评估BRDF模型的不同参数配置与太阳天顶角对CI估计的影响，结果表明太阳天顶角为0°时会低估CI值，而60°时会高估CI值。此外，Jiao等人在线性核驱动的BRDF模型基础上，提出了一种表征纯雪反射率各向异性的雪核，从而更好地模拟雪的散特性。
在主动光源各向异性的研究领域，Li等人利用卫星传感器提供的多角度观测信息，探究卫星传感器观测角度与夜光辐亮度之间的关系，结果表明，在城市区域二次模型可以较好地描述夜光辐亮度与观测天顶角之间的关系，并且城市内部空间结构会影响角度与辐亮度之间的关系。由于受到观测角度效应的影响，夜光辐亮度会造成错误检测，而无法反映灯光的实际变化。Elvidge等人针对夜光辐亮度随观测天顶角变化的四种趋势(平稳、凸起、凹陷、最低观测点峰值)，采用最接近观测天顶角0°簇内夜光辐亮度中位数与当前校正角度簇内夜光辐亮度中位数的比值作为系数，分别乘以校正角度簇内的夜光辐亮度观测值从而完成夜光辐亮度角度效应的归一化，进而有益于电力可靠性的时空模式挖掘。在光污染研究中，Tong等人基于无暮光与月光的VIIRS夜光辐亮度数据绘制了北欧和北非向上人造光的角度分布，借此更好地了解光污染在夜间环境中的传播机制。Kocifaj等人利用观测角度模拟夜间环境的人造光，并引入纠正功能改进城市发射函数，从而提高了夜晚城市上空光分布的估算结果。
综合前人研究，目前关于角度各向异性的研究主要针对地表反射率的遥感产品，很少考虑由主动光源生成的夜光影像角度各向异性，并且目前关于夜光遥感的研究多数聚集在城市化进程与社会经济参数估算等领域的应用，忽略了数据本身的质量存在不确定性与无效性，从而导致分析结果的偏差。
发明内容
本发明针对现有技术的不足，提供一种卫星遥感夜光辐亮度的观测角度归一化方法。首先通过建立夜光辐亮度与观测天顶角之间的天顶-辐亮度二次(ZRQ)模型，采用最小二乘多项式曲线拟合方式求解二次函数的系数，然后基于ZRQ模型，将观测到的夜光辐亮度分解为天顶角为0°时的夜光辐亮度与二次函数值的乘积，最后采用最小二乘法求解分解后的二次函数系数，并计算当前传感器获取的夜光辐亮度与二次函数值的比值，求解天顶角为0°时的夜光辐亮度，进而消除夜光辐亮度的角度效应。
为了达到上述目的，本发明提供的技术方案是一种卫星遥感夜光辐亮度的观测角度归一化方法，包括以下步骤：
步骤1，获取研究区域黑色大理石(Black Marble)夜光影像产品中VNP46A1与VNP46A2数据。
步骤2，对步骤1获取到的所有VNP46A1与VNP46A2数据进行数据格式转换，并得到需要的质量文件与影像文件数据。
步骤3，对步骤2所得到的质量标志数据，设置一定的筛选条件，得到符合需求的高质量夜光辐亮度数据。
步骤4，以像素为研究单位，从步骤3得到的高质量夜光辐亮度数据中提取像素的夜光辐亮度时间序列和传感器观测天顶角时间序列。
步骤5，对步骤4提取的夜光辐亮度与卫星传感器观测天顶角进行统计建模，确立两者之间的量化关系，并计算多项式系数，包括以下子步骤：
步骤5.1，以观测天顶角作为自变量，夜光辐亮度作为因变量，建立两者之间的天顶-辐亮度二次(ZRQ)模型。
步骤5.2，计算ZRQ模型的二次函数系数。
步骤6，基于步骤5建立的天顶-辐亮度二次(ZRQ)统计模型，消除夜光辐亮度角度效应，包括以下子步骤：
步骤6.1，基于步骤5建立ZRQ统计模型，将模型公式分解为天顶角为0°时的夜光辐亮度(常数项系数)与二次函数值的乘积形式。
步骤6.2，将步骤4提取的夜光辐亮度时间序列{R1,R2,...,Rn}除以步骤5得到的二次模型的常数项系数c，获得新的夜光辐亮度时间序列并用新的夜光辐亮度时间序列与对应时刻的卫星传感器观测天顶角时间序列{z1,z2,...,zn}进行二次建模，使用如步骤5.2中的最小二乘法多项式曲线拟合求解二次函数系数a′、b′。
步骤6.3，依次计算时间序列每一时刻传感器观测到的夜光辐亮度与步骤6.1分解得到的二次函数的比值，从而求解天顶角为0°时的夜光辐亮度，最终得到的天顶角为0°时的夜光辐亮度时间序列为：{R1/(a′z1 2+b′z1+1)，R2/(a′z2 2+b′z2+1)，…，Rn/(a′zn 2+b′zn+1)}。
步骤7，计算夜光辐亮度时间序列角度归一化后的变异系数(CV)与归一化冷热点夜光辐亮度(NDHDNTL)，量化角度效应在夜光辐亮度时间序列中的影响程度。
而且，所述步骤5.1中以观测天顶角作为自变量，夜光辐亮度作为因变量，建立两者之间的天顶-辐亮度二次(ZRQ)拟合关系如下：
R＝az2+bz+c (1)
式中，z为卫星观测天顶角，R为夜光辐亮度，a、b、c分别为二次模型系数。
而且，所述步骤5.2中假设夜光辐亮度时间序列为{R1,R2,...,Rn}，传感器观测天顶角时间序列为{z1,z2,...,zn}，在夜光辐亮度与观测天顶角拟合二次曲线过程中定义拟合的目标函数使得拟合偏差最小，目标函数公式如下：
式中，zi与Ri分别为时间序列中第i时刻的观测天顶角与夜光辐亮度，为拟合的夜光辐亮度，n为时间序列长度。
建立误差方程，并求解二次曲线函数的系数矩阵，公式如下：
V＝BX-L (3)
其中，
X＝(BTB)-1BTL (4)
式中，V表示拟合辐亮度与实测辐亮度之间的误差矩阵，B表示公式(1)中各系数对应的坐标组合项，X为公式(1)中待求的二次曲线函数的系数矩阵，L为实测辐亮度矩阵。
而且，所述步骤6.1中将基于步骤5建立ZRQ统计模型分解为天顶角为0°时的夜光辐亮度(常数项系数)与二次函数值的乘积形式，公式如下：
R＝c(t)×(a′z2+b′z+1) (5)
式中，z代表传感器天顶角，c(t)代表分解成天顶角为0°的夜光辐亮度，(a′z2+b′z+1)代表夜光辐亮度分解后二次函数的值，a′、b′为二次函数的系数。
而且，所述步骤7中CV与NDHDNTL计算公式如下：
式中，NTLhOtspOt与NTLdarkSpOt分别代表夜光辐亮度时间序列中最大、最小的夜光辐亮度，NTLi代表第i时刻的夜光辐亮度，n代表时间序列的个数。
与现有技术相比，本发明具有如下优点：
(1)本发明利用长时间序列的夜光辐亮度数据，建立夜光辐亮度与观测天顶角之间的统计模型，有助于更好地量化夜光辐亮度与观测天顶角之间的关系；
(2)本发明基于天顶-辐亮度二次模型构建夜光辐亮度角度归一化模型，考虑到了人类活动以及传感器测量误差等因素造成的时间序列波动，采用连续两次模型拟合并求解参数，提高了角度归一化的精度；
(3)本发明通过最小二乘法多项式曲线拟合求解统计模型参数，运算效率较高，可以实现对长时间覆盖范围的夜光辐亮度角度归一化。
附图说明
图1为本发明实施例的流程图。
图2为本发明实施例中天顶-辐亮度二次(ZRQ)模型的建立情况。
图3为本发明实施例中洛杉矶研究点夜光辐亮度角度归一化前后时间序列对比示意图：其中图3(a)为洛杉矶1号研究点夜光辐亮度角度归一化前后时间序列对比示意图，图3(b)为洛杉矶2号研究点夜光辐亮度角度归一化前后时间序列对比示意图。
图4为本发明实施例中拉斯维加斯研究点夜光辐亮度角度归一化前后时间序列对比示意图：其中图4(a)为拉斯维加斯1号研究点夜光辐亮度角度归一化前后时间序列对比示意图，图4(b)为拉斯维加斯2号研究点夜光辐亮度角度归一化前后时间序列对比示意图。
具体实施方式
本发明提供一种卫星遥感夜光辐亮度的观测角度归一化方法，首先通过建立夜光辐亮度与观测天顶角之间的天顶-辐亮度二次(ZRQ)模型，采用最小二乘多项式曲线拟合方式求解二次函数的系数，然后基于ZRQ模型，将观测到的夜光辐亮度分解为天顶角为0°时的夜光辐亮度与二次函数值的乘积，最后采用最小二乘法求解分解后的二次函数系数，并计算当前传感器获取的夜光辐亮度与二次函数值的比值，求解天顶角为0°时的夜光辐亮度，进而消除夜光辐亮度的角度效应。
下面结合附图和实施例对本发明的技术方案作进一步说明。
如图1所示，本发明实施例的流程包括以下步骤：
步骤1，获取研究区域黑色大理石(Black Marble)夜光影像产品中VNP46A1与VNP46A2数据。
步骤2，对步骤1获取到的所有VNP46A1与VNP46A2数据进行数据格式转换，并得到需要的质量文件与影像文件数据。
步骤3，对步骤2所得到的质量标志数据，设置一定的筛选条件，得到符合需求的高质量夜光辐亮度数据。
步骤4，以像素为研究单位，从步骤3得到的高质量夜光辐亮度数据中提取像素的夜光辐亮度时间序列和传感器观测天顶角时间序列。
步骤5，对步骤4提取的夜光辐亮度与卫星传感器观测天顶角进行统计建模，确立两者之间的量化关系，并计算多项式系数，包括以下子步骤：
步骤5.1，以观测天顶角作为自变量，夜光辐亮度作为因变量，建立两者之间的天顶-辐亮度二次(ZRQ)拟合关系如下：
R＝az2+bz+c (1)
式中，z为卫星观测天顶角，R为夜光辐亮度，a、b、c分别为二次模型系数。
步骤5.2，计算ZRQ模型的二次函数系数。
假设夜光辐亮度时间序列为{R1,R2,...,Rn}，传感器观测天顶角时间序列为{z1,z2,...,zn}，在夜光辐亮度与观测天顶角拟合二次曲线过程中定义拟合的目标函数使得拟合偏差最小，目标函数公式如下：
式中，zi与Ri分别为时间序列中第i时刻的观测天顶角与夜光辐亮度，为拟合的夜光辐亮度，n为时间序列长度。
建立误差方程，并求解二次曲线函数的系数矩阵，公式如下：
V＝BX-L (3)
其中，
X＝(BTB)-1BTL (4)
式中，V表示拟合辐亮度与实测辐亮度之间的误差矩阵，B表示公式(1)中各系数对应的坐标组合项，X为公式(1)中待求的二次曲线函数的系数矩阵，L为实测辐亮度矩阵。
步骤6，基于步骤5建立的天顶-辐亮度二次(ZRQ)统计模型，消除夜光辐亮度角度效应，包括以下子步骤：
步骤6.1，基于步骤5建立ZRQ统计模型，将模型公式分解为天顶角为0°时的夜光辐亮度(常数项系数)与二次函数值的乘积形式，公式如下：
R＝c(t)×(a′z2+b′z+1) (5)
式中，z代表传感器天顶角，c(t)代表分解成天顶角为0°的夜光辐亮度，(a′z2+b′z+1)代表夜光辐亮度分解后二次函数的值，a′、b′为二次函数的系数。
步骤6.2，将步骤4提取的夜光辐亮度时间序列{R1,R2,...,Rn}除以步骤5得到的二次模型的常数项系数c，获得新的夜光辐亮度时间序列并用新的夜光辐亮度时间序列与对应时刻的卫星传感器观测天顶角时间序列{z1,z2,...,zn}进行二次建模，使用如步骤5.2中的最小二乘法多项式曲线拟合求解二次函数系数a′、b′。
步骤6.3，在求解得到系数之后，依次计算时间序列每一时刻传感器观测到的夜光辐亮度与步骤6.1分解得到的二次函数的比值，从而求解天顶角为0°时的夜光辐亮度。最终得到的天顶角为0°时的夜光辐亮度时间序列为：
步骤7，计算夜光辐亮度时间序列角度归一化后的变异系数(CV)与归一化冷热点夜光辐亮度(NDHDNTL)，量化角度效应在夜光辐亮度时间序列中的影响程度。
CV与NDHDNTL计算公式如下：
式中，NTLhOtSpOt与NTLdarkSpOt分别代表夜光辐亮度时间序列中最大、最小的夜光辐亮度，NTLi代表第i时刻的夜光辐亮度，n代表时间序列的个数。
下面以VNP46A1与VNP46A2实验数据为例，进一步阐述本发明的技术方案。
步骤1，获取研究区域黑色大理石(Black Marble)夜光影像产品中VNP46A1与VNP46A2数据。
从LAADS DAAC网站(https://ladsweb.modaps.eosdis.nasa.gov/)上下载VIIRS/NPP卫星获取的每日夜光辐亮度产品VNP46A1与VNP46A2，数据的空间分辨率为15弧秒。影像涵盖2016年366天的数据，包括VNP46A1中的26个科学数据集以及VNP46A2中的7个科学数据集。本发明选择了包含洛杉矶与拉斯维加斯(h06v05)所在区域的夜光影像数据。
步骤2，对下载获取到的所有VNP46A1与VNP46A2数据进行数据格式转换，并得到需要的质量文件与影像文件数据。
将下载获取的HDF5格式科学数据集转换成GeoTiff格式，选取VNP46A1产品中传感器观测天顶角、太阳天顶角、月光照度分数文件，以及VNP46A2产品中经过BRDF校正的夜光影像、云层掩码、强制性质量指标文件。其中太阳天顶角、月光照度分数、云层掩码以及强制性质量指标文件作为后续像素筛选的质量标志数据，经过BRDF校正的夜光影像与传感器观测天顶角将用于天顶角与夜光辐亮度的建模。
步骤3，对步骤2中得到的质量标志数据，设置一定的筛选条件，得到符合需求的高质量夜光辐亮度数据。
分别设置云层掩码指标、强制质量指标、太阳天顶角、月光照度分数作为逐像素夜光辐亮度质量的筛选条件，其中云层掩码是作为判别像素是否被云层遮挡的重要信息，其11位具体设置如表5.1所示。
表5.1对云层掩码质量标志各个位数据设置规则
但是在云层筛选的基础上，仍然存在部分像素值异常的情况，因此，将强制性质量指标设置为0，以保证所选像素均为高质量，从而提高数据质量。此外，剔除太阳天顶角大于108°或月光照度分数小于1％的像素，以保证获取像素的有效性和降低云层掩码标志数据质量错误的可能性。
步骤4，以像素为研究单位，从步骤3得到的高质量夜光辐亮度数据中提取像素的夜光辐亮度时间序列和传感器观测天顶角时间序列。
以当前像素3×3的窗口作为分析边界，首先保证3×3窗口内的9个像素经过条件筛选后均符合质量要求，然后再计算9个像素夜光辐亮度的均值作为当前像素的夜光辐亮度，传感器观测天顶角数据仅取当前像素位置对应的天顶角数据，最后分别读取2016年366天的夜光影像，提取出研究区域高质量的夜光辐亮度与观测角度时间序列。
步骤5，对夜光辐亮度与卫星传感器观测天顶角进行统计建模，确立两者之间的量化关系，并计算多项式系数，包括以下子步骤：
步骤5.1，以观测天顶角作为自变量，夜光辐亮度作为因变量，建立两者之间的天顶-辐亮度二次(ZRQ)拟合关系如下：
R＝az2+bz+c (1)
式中，z为卫星观测天顶角，R为夜光辐亮度，a、b、c分别为二次模型系数。
步骤5.2，计算ZRQ模型的二次函数系数。
假设夜光辐亮度时间序列为{R1,R2,...,Rn}，传感器观测天顶角时间序列为{z1,z2,...,zn}，在夜光辐亮度与观测天顶角拟合二次曲线过程中定义拟合的目标函数使得拟合偏差最小，目标函数公式如下：
式中，zi与Ri分别为时间序列当中第i时刻的观测天顶角与夜光辐亮度，n为时间序列长度。
建立误差方程，并求解二次曲线函数的系数矩阵，公式如下：
V＝BX-L (3)
其中，
X＝(BTB)-1BTL (4)
式中，V表示拟合辐亮度与实测辐亮度之间的误差矩阵，B表示公式(1)中各系数对应的坐标组合项，X为公式(1)中待求的二次曲线函数的系数矩阵，L为实测辐亮度矩阵。
步骤6，基于步骤5建立天顶-辐亮度二次(ZRQ)统计模型，消除夜光辐亮度角度效应，包括以下子步骤：
步骤6.1，基于步骤5建立ZRQ统计模型，将模型公式分解为天顶角为0°时的夜光辐亮度(常数项系数)与二次函数值的乘积形式，公式如下：
R＝c(t)×(a′z2+b′z+1) (5)
式中，z代表传感器天顶角，c(t)代表分解成天顶角为0°的夜光辐亮度，(a′z2+b′z+1)代表夜光辐亮度分解后二次函数的值，a′、b′为二次函数的系数。
步骤6.2，将步骤4提取的夜光辐亮度时间序列{R1,R2,...,Rn}除以步骤5得到的二次模型的常数项系数c，获得新的夜光辐亮度时间序列并用新的夜光辐亮度时间序列与对应时刻的卫星传感器观测天顶角时间序列{z1,z2,...,zn}进行二次建模，使用如步骤5.2中的最小二乘法多项式曲线拟合求解二次函数系数a′、b′。
步骤6.3，在求解得到系数之后，依次计算时间序列每一时刻传感器观测到的夜光辐亮度与步骤6.1分解得到的二次函数的比值，从而求解天顶角为0°时的夜光辐亮度。最终得到的天顶角为0°时的夜光辐亮度时间序列为：
步骤7，计算夜光辐亮度时间序列角度归一化后的变异系数(CV)与归一化冷热点夜光辐亮度(NDHDNTL)，量化角度效应在夜光辐亮度时间序列中的影响程度。
CV与NDHDNTL计算公式如下：
式中，NTLhotspot与NTLdarksport分别代表夜光辐亮度时间序列中最大、最小的夜光辐亮度，NTLi代表第i时刻的夜光辐亮度，n代表时间序列的个数。
角度校正前后不同区域研究点的归一化冷热点夜光辐亮度(NDHDNTL)与变异系数(CV)如表5.2结果所示，角度归一化使得时间序列更加平稳，总体CV下降超过30％，这表明归一化后的时间序列，由于消除了角度效应而降低了内部的差异性，从而提升了时间序列的稳定性。
表5.2角度校正前后不同区域研究点的归一化冷热点夜光辐亮度(NDHDNTL)与变异系数(CV)
具体实施时，以上流程可采用计算机软件技术实现自动运行流程。
本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。
Claims (8)
1.一种卫星遥感夜光辐亮度的观测角度归一化方法，其特征在于，包括如下步骤：
步骤1，获取研究区域Black Marble夜光影像产品中VNP46A1与VNP46A2数据；
步骤2，对步骤1获取到的所有VNP46A1与VNP46A2数据进行数据格式转换，并得到需要的质量文件与影像文件数据；
将下载获取的HDF5格式科学数据集转换成GeoTiff格式，选取VNP46A1产品中传感器观测天顶角、太阳天顶角、月光照度分数文件，以及VNP46A2产品中经过BRDF校正的夜间灯光影像、云层掩码、强制性质量指标文件，其中太阳天顶角、月光照度分数、云层掩码以及强制性质量指标文件作为后续像素筛选的质量标志数据，经过BRDF校正的夜间灯光影像与传感器观测天顶角将用于天顶角与夜光辐亮度的建模；
步骤3，对步骤2所得到的质量标志数据，设置一定的筛选条件，得到符合需求的高质量夜间灯光数据；
步骤4，以像素为研究单位，从步骤3得到的高质量夜间灯光数据中提取像素的夜光辐亮度时间序列和传感器观测天顶角时间序列；
步骤5，对步骤4提取的夜光辐亮度与卫星传感器观测天顶角进行统计建模，确立两者之间的量化关系，并计算多项式系数；
步骤6，基于步骤5建立的天顶-辐亮度二次模型，消除夜光辐亮度角度效应；
步骤7，计算夜光辐亮度时间序列角度归一化后的变异系数与归一化冷热点夜光辐亮度，量化角度效应在夜间灯光时间序列中的影响程度。
2.如权利要求1所述的一种卫星遥感夜光辐亮度的观测角度归一化方法，其特征在于：所述步骤3中设置云层掩码指标、强制质量指标、太阳天顶角和月光照度分数作为逐像素夜间灯光质量的筛选条件，但是在云层筛选的基础上，仍然存在部分像素值异常的情况，因此，将强制性质量指标设置为0，以保证所选像素均为高质量，从而提高数据质量，剔除太阳天顶角大于108°或月光照度分数小于1％的像素，以保证获取像素的有效性和降低云层掩码标志数据质量错误的可能性。
3.如权利要求1所述的一种卫星遥感夜光辐亮度的观测角度归一化方法，其特征在于：所述步骤5包括以下子步骤：
步骤5.1，以观测天顶角作为自变量，夜光辐亮度作为因变量，建立两者之间的天顶-辐亮度二次模型；
步骤5.2，计算天顶-辐亮度二次模型的二次函数系数。
4.如权利要求3所述的一种卫星遥感夜光辐亮度的观测角度归一化方法，其特征在于：所述步骤5.1是以观测天顶角作为自变量，夜光辐亮度作为因变量，建立两者之间的天顶-辐亮度二次模型拟合关系如下：
R＝az2+bz+c (1)
式中，z为卫星观测天顶角，R为夜光辐亮度，a、b、c分别为二次模型系数。
5.如权利要求3所述的一种卫星遥感夜光辐亮度的观测角度归一化方法，其特征在于：所述步骤5.2中中假设夜光辐亮度时间序列为{R1,R2,...,Rn}，传感器观测天顶角时间序列为{z1,z2,...,zn}，在夜光辐亮度与观测天顶角拟合二次曲线过程中定义拟合的目标函数使得拟合偏差最小，目标函数公式如下：
式中，zi与Ri分别为时间序列中第i时刻的观测天顶角与夜光辐亮度，为拟合的夜光辐亮度，n为时间序列长度；
建立误差方程，并求解二次曲线函数的系数矩阵，公式如下：
V＝BX-L (3)
其中，
X＝(BTB)-1BTL (4)
式中，V表示拟合辐亮度与实测辐亮度之间的误差矩阵，B表示公式(1)中各系数对应的坐标组合项，X为公式(1)中待求的二次曲线函数的系数矩阵，L为实测辐亮度矩阵。
6.如权利要求5所述的一种卫星遥感夜光辐亮度的观测角度归一化方法，其特征在于：所述步骤6包括以下子步骤：
步骤6.1，基于步骤5建立天顶-辐亮度二次模型，将模型公式分解为天顶角为0°时的夜光辐亮度与二次函数值的乘积形式；
步骤6.2，将步骤4提取的夜光辐亮度时间序列{R1,R2,...,Rn}除以步骤5得到的二次模型的常数项系数c，获得新的夜光辐亮度时间序列并用新的夜光辐亮度时间序列与对应时刻的卫星传感器观测天顶角时间序列{z1,z2,...,zn}进行二次建模，用步骤5.2中的最小二乘法多项式曲线拟合求解二次函数系数a′、b′；
步骤6.3，依次计算时间序列每一时刻传感器观测到的夜光辐亮度与步骤6.1分解得到的二次函数的比值，从而求解天顶角为0°时的夜光辐亮度，最终得到的天顶角为0°时的夜光辐亮度时间序列为：{R1/(a′z1 2+b′z1+1)，R2/(a′z2 2+b′z2+1)，…，Rn/(a′zn 2+b′zn+1)}。
7.如权利要求6所述的一种卫星遥感夜光辐亮度的观测角度归一化方法，其特征在于：所述步骤6.1中将基于步骤5建立天顶-辐亮度二次模型分解为天顶角为0°时的夜光辐亮度与二次函数值的乘积形式，公式如下：
R＝c(t)×(a′z2+b′z+1) (5)
式中，z代表传感器天顶角，c(t)代表分解成天顶角为0°的夜光辐亮度，(a′z2+b′z+1)代表夜光辐亮度分解后二次函数的值，a′、b′为二次函数的系数。
8.如权利要求1所述的一种卫星遥感夜光辐亮度的观测角度归一化方法，其特征在于：所述步骤7中变异系数CV与归一化冷热点夜光辐亮度NDHDNTL计算公式如下：
式中，NTLhotspot与NTLdarkspot分别代表夜间灯光时间序列中最大、最小的夜光辐亮度，NTLi代表第i时刻的夜光辐亮度，n代表时间序列的个数。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110743139.5A CN113569386B (zh) | 2021-07-01 | 2021-07-01 | 卫星遥感夜光辐亮度的观测角度归一化方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110743139.5A CN113569386B (zh) | 2021-07-01 | 2021-07-01 | 卫星遥感夜光辐亮度的观测角度归一化方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN113569386A CN113569386A (zh) | 2021-10-29 |
| CN113569386B true CN113569386B (zh) | 2023-08-22 |
Family
ID=78163309
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202110743139.5A Active CN113569386B (zh) | 2021-07-01 | 2021-07-01 | 卫星遥感夜光辐亮度的观测角度归一化方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN113569386B (zh) |
Families Citing this family (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115439642B (zh) * | 2022-08-26 | 2024-04-19 | 武汉大学 | 一种夜光遥感能量影像合成方法、装置及设备 |
| CN115880583B (zh) * | 2022-08-26 | 2024-03-19 | 武汉大学 | 一种夜光遥感影像的农田火识别与去除方法 |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110120018A (zh) * | 2019-04-10 | 2019-08-13 | 武汉大学 | 一种面阵高动态范围夜光成像卫星在轨相对辐射定标方法 |
| CN111561936A (zh) * | 2020-05-19 | 2020-08-21 | 中国科学院微小卫星创新研究院 | 旋转大幅宽光学卫星精准处理方法及系统 |
| CN111563962A (zh) * | 2020-04-09 | 2020-08-21 | 中国科学院空天信息创新研究院 | 一种基于几何辐射一体化采样的遥感图像仿真方法 |
| CN111680659A (zh) * | 2020-06-17 | 2020-09-18 | 中国科学院空天信息创新研究院 | 国际空间站rgb夜间灯光图像的相对辐射归一化方法 |
| CN112924391A (zh) * | 2021-02-03 | 2021-06-08 | 重庆邮电大学 | 基于遥感大数据的fy-4a/agri交叉辐射定标方法 |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US8165813B2 (en) * | 2011-07-25 | 2012-04-24 | Clean Power Research, L.L.C. | Computer-implemented system and method for efficiently performing area-to-point conversion of satellite imagery for photovoltaic power generation fleet output estimation |
-
2021
- 2021-07-01 CN CN202110743139.5A patent/CN113569386B/zh active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110120018A (zh) * | 2019-04-10 | 2019-08-13 | 武汉大学 | 一种面阵高动态范围夜光成像卫星在轨相对辐射定标方法 |
| CN111563962A (zh) * | 2020-04-09 | 2020-08-21 | 中国科学院空天信息创新研究院 | 一种基于几何辐射一体化采样的遥感图像仿真方法 |
| CN111561936A (zh) * | 2020-05-19 | 2020-08-21 | 中国科学院微小卫星创新研究院 | 旋转大幅宽光学卫星精准处理方法及系统 |
| CN111680659A (zh) * | 2020-06-17 | 2020-09-18 | 中国科学院空天信息创新研究院 | 国际空间站rgb夜间灯光图像的相对辐射归一化方法 |
| CN112924391A (zh) * | 2021-02-03 | 2021-06-08 | 重庆邮电大学 | 基于遥感大数据的fy-4a/agri交叉辐射定标方法 |
Non-Patent Citations (1)
| Title |
|---|
| 杨磊 等.高分一号卫星相机的辐射交叉定标研究.《红外与激光工程》.2015,第44卷(第08期),2456-2460. * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN113569386A (zh) | 2021-10-29 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Zheng et al. | Developing a new cross-sensor calibration model for DMSP-OLS and Suomi-NPP VIIRS night-light imageries | |
| Zhao et al. | Building a series of consistent night-time light data (1992–2018) in Southeast Asia by integrating DMSP-OLS and NPP-VIIRS | |
| Xu et al. | An integrated method for validating long-term leaf area index products using global networks of site-based measurements | |
| CN109919875B (zh) | 一种高时频遥感图像特征辅助的居民地提取与分类方法 | |
| CN113569386B (zh) | 卫星遥感夜光辐亮度的观测角度归一化方法 | |
| Ma et al. | Quantifying spatiotemporal patterns of urban impervious surfaces in China: An improved assessment using nighttime light data | |
| Liu et al. | Evaluation of consistency among three NDVI products applied to High Mountain Asia in 2000–2015 | |
| Liu et al. | Application of a new leaf area index algorithm to China's landmass using MODIS data for carbon cycle research | |
| CN103824077A (zh) | 一种基于多源遥感数据的城市不透水层率息提取方法 | |
| Zhao et al. | Spatiotemporal characteristics of urban surface temperature and its relationship with landscape metrics and vegetation cover in rapid urbanization region | |
| Shen et al. | Spatial–temporal land-use/land-cover dynamics and their impacts on surface temperature in Chongming Island of Shanghai, China | |
| Xiang et al. | Integration of tillage indices and textural features of Sentinel-2A multispectral images for maize residue cover estimation | |
| Ma et al. | Mapping fine-scale building heights in urban agglomeration with spaceborne lidar | |
| Zhang et al. | Development of the direct-estimation albedo algorithm for snow-free Landsat TM albedo retrievals using field flux measurements | |
| Zhang et al. | Methods for automatic identification and extraction of terraces from high spatial resolution satellite data (China-GF-1) | |
| CN115204691A (zh) | 基于机器学习和遥感技术的城市人为热排放量估算方法 | |
| Hu et al. | Modeling the spatiotemporal dynamics of global electric power consumption (1992–2019) by utilizing consistent nighttime light data from DMSP-OLS and NPP-VIIRS | |
| Zhao et al. | Evaluating the potential of H8/AHI geostationary observations for monitoring vegetation phenology over different ecosystem types in northern China | |
| Xin et al. | Seasonal differences in the dominant factors of surface urban heat islands along the urban-rural gradient | |
| Rakuasa et al. | Analysis of Vegetation Index in Ambon City Using Sentinel-2 Satellite Image Data with Normalized Difference Vegetation Index (NDVI) Method based on Google Earth Engine | |
| Geng et al. | Vegetation coverage of desert ecosystems in the Qinghai-Tibet Plateau is underestimated | |
| Nisar et al. | Assessment and Monitoring of VIIRS-DNB and SQML-L light Pollution in Lahore-Pakistan | |
| Yang et al. | Influence of urban extent discrepancy on the estimation of surface urban heat island intensity: A global-scale assessment in 892 cities | |
| CN114969425A (zh) | 一种基于地理加权回归的夜光遥感数据空间降尺度方法 | |
| CN115439642B (zh) | 一种夜光遥感能量影像合成方法、装置及设备 |
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 |