卫星遥感夜光辐亮度的观测角度归一化方法

1.本发明属于统计回归模型在夜光遥感影像各向异性消除领域的应用，特别是涉及一种卫星遥感夜光辐亮度的观测角度归一化方法。


背景技术：

2.夜光遥感是获取夜间无云时地表发射的可见光信息，通过夜间灯光可以发现人类聚集区并记录地表灯光强度信息，进而可以通过夜光遥感影像间接反映居民在城镇中的行为活动与空间分布。深入挖掘夜光辐亮度时间序列的动态变化规律，对于准确监测城市化进程、电力供应以及评估自然灾害等重大突发事件有着重要的作用。
3.随着遥感技术的不断发展，夜光遥感产品逐渐丰富。特别是新一代夜光辐亮度数据npp/viirs(national polar
‑
orbiting partnership's visible infrared imaging radiometer suite)不仅具有检测微弱夜间灯光的能力，而且弥补了dmsp/ols数据低空间分辨率与低辐射分辨率的缺点，能更多地应用在建立夜光辐亮度与人类活动模式经验关系的研究中。由于夜间灯光信号会受到杂散光、大气成分、季节性等噪声信号干扰，从而造成夜光辐亮度时间序列的波动。为此，nasa团队生产出了一套全球可用的、包含丰富时间与空间信息的黑色大理石夜光影像产品(vnp46)。该产品去除了受云污染的像素，并对受大气、地形、植被、雪、月球以及杂散光影响的夜光辐亮度进行了校正，使得产品的夜光辐亮度更为可靠。但是li等学者发现，夜光辐亮度与卫星传感器观测天顶角存在非线性关系，特别是在高建筑物区域，受到城市三维形态的影响，夜光辐亮度会存在明显的冷点和热点效应。因此，在利用长时间序列的夜光影像进行研究时，为了更准确地反映城市的动态变化，消除夜光辐亮度的角度效应十分必要。
4.目前许多学者在光学各向异性的研究主要包括基于双向反射率分布函数(brdf)表征地球表面反射率各向异性的研究与模拟城市上空光分布的光污染研究。地球表面反射光会受到地物的物理结构、传感器视角以及太阳光照角度变化的影响，产生地球表面反射率各向异性特征。landsat卫星通过最低点一定范围的视角获取地球表面影像，会导致地表反射率产生较小的方向性影响。roy等人基于全球高质量的modis brdf产品，得出一组固定的brdf光谱模型参数，从而完成landsat数据地表反射率角度校正，但作为描述叶子空间分布格局的关键参数，聚集指数(ci)会受到brdf模型与太阳天顶角的影响。wei等人使用归一化冷热点差异(ndhd)评估brdf模型的不同参数配置与太阳天顶角对ci估计的影响，结果表明太阳天顶角为0
°
时会低估ci值，而60
°
时会高估ci值。此外，jiao等人在线性核驱动的brdf模型基础上，提出了一种表征纯雪反射率各向异性的雪核，从而更好地模拟雪的散特性。
5.在主动光源各向异性的研究领域，li等人利用卫星传感器提供的多角度观测信息，探究卫星传感器观测角度与夜光辐亮度之间的关系，结果表明，在城市区域二次模型可以较好地描述夜光辐亮度与观测天顶角之间的关系，并且城市内部空间结构会影响角度与辐亮度之间的关系。由于受到观测角度效应的影响，夜光辐亮度会造成错误检测，而无法反
映灯光的实际变化。elvidge等人针对夜光辐亮度随观测天顶角变化的四种趋势(平稳、凸起、凹陷、最低观测点峰值)，采用最接近观测天顶角0
°
簇内夜光辐亮度中位数与当前校正角度簇内夜光辐亮度中位数的比值作为系数，分别乘以校正角度簇内的夜光辐亮度观测值从而完成夜光辐亮度角度效应的归一化，进而有益于电力可靠性的时空模式挖掘。在光污染研究中，tong等人基于无暮光与月光的viirs夜光辐亮度数据绘制了北欧和北非向上人造光的角度分布，借此更好地了解光污染在夜间环境中的传播机制。kocifaj等人利用观测角度模拟夜间环境的人造光，并引入纠正功能改进城市发射函数，从而提高了夜晚城市上空光分布的估算结果。
6.综合前人研究，目前关于角度各向异性的研究主要针对地表反射率的遥感产品，很少考虑由主动光源生成的夜光影像角度各向异性，并且目前关于夜光遥感的研究多数聚集在城市化进程与社会经济参数估算等领域的应用，忽略了数据本身的质量存在不确定性与无效性，从而导致分析结果的偏差。


技术实现要素：

7.本发明针对现有技术的不足，提供一种卫星遥感夜光辐亮度的观测角度归一化方法。首先通过建立夜光辐亮度与观测天顶角之间的天顶
‑
辐亮度二次(zrq)模型，采用最小二乘多项式曲线拟合方式求解二次函数的系数，然后基于zrq模型，将观测到的夜光辐亮度分解为天顶角为0
°
时的夜光辐亮度与二次函数值的乘积，最后采用最小二乘法求解分解后的二次函数系数，并计算当前传感器获取的夜光辐亮度与二次函数值的比值，求解天顶角为0
°
时的夜光辐亮度，进而消除夜光辐亮度的角度效应。
8.为了达到上述目的，本发明提供的技术方案是一种卫星遥感夜光辐亮度的观测角度归一化方法，包括以下步骤：
9.步骤1，获取研究区域黑色大理石(black marble)夜光影像产品中vnp46a1与vnp46a2数据。
10.步骤2，对步骤1获取到的所有vnp46a1与vnp46a2数据进行数据格式转换，并得到需要的质量文件与影像文件数据。
11.步骤3，对步骤2所得到的质量标志数据，设置一定的筛选条件，得到符合需求的高质量夜光辐亮度数据。
12.步骤4，以像素为研究单位，从步骤3得到的高质量夜光辐亮度数据中提取像素的夜光辐亮度时间序列和传感器观测天顶角时间序列。
13.步骤5，对步骤4提取的夜光辐亮度与卫星传感器观测天顶角进行统计建模，确立两者之间的量化关系，并计算多项式系数，包括以下子步骤：
14.步骤5.1，以观测天顶角作为自变量，夜光辐亮度作为因变量，建立两者之间的天顶
‑
辐亮度二次(zrq)模型。
15.步骤5.2，计算zrq模型的二次函数系数。
16.步骤6，基于步骤5建立的天顶
‑
辐亮度二次(zrq)统计模型，消除夜光辐亮度角度效应，包括以下子步骤：
17.步骤6.1，基于步骤5建立zrq统计模型，将模型公式分解为天顶角为0
°
时的夜光辐亮度(常数项系数)与二次函数值的乘积形式。
18.步骤6.2，将步骤4提取的夜光辐亮度时间序列{r1,r2,...,r
n
}除以步骤5得到的二次模型的常数项系数c，获得新的夜光辐亮度时间序列并用新的夜光辐亮度时间序列与对应时刻的卫星传感器观测天顶角时间序列{z1,z2,...,z
n
}进行二次建模，使用如步骤5.2中的最小二乘法多项式曲线拟合求解二次函数系数a
′
、b
′
。
19.步骤6.3，依次计算时间序列每一时刻传感器观测到的夜光辐亮度与步骤6.1分解得到的二次函数的比值，从而求解天顶角为0
°
时的夜光辐亮度，最终得到的天顶角为0
°
时的夜光辐亮度时间序列为：{r1/(a
′
z
12
+b
′
z1+1)，r2/(a
′
z
22
+b
′
z2+1)，
…
，r
n
/(a
′
z
n2
+b
′
z
n
+1)}。
20.步骤7，计算夜光辐亮度时间序列角度归一化后的变异系数(cv)与归一化冷热点夜光辐亮度(ndhd
ntl
)，量化角度效应在夜光辐亮度时间序列中的影响程度。
21.而且，所述步骤5.1中以观测天顶角作为自变量，夜光辐亮度作为因变量，建立两者之间的天顶
‑
辐亮度二次(zrq)拟合关系如下：
22.r＝az2+bz+c
ꢀꢀ
(1)
23.式中，z为卫星观测天顶角，r为夜光辐亮度，a、b、c分别为二次模型系数。
24.而且，所述步骤5.2中假设夜光辐亮度时间序列为{r1,r2,...,r
n
}，传感器观测天顶角时间序列为{z1,z2,...,z
n
}，在夜光辐亮度与观测天顶角拟合二次曲线过程中定义拟合的目标函数使得拟合偏差最小，目标函数公式如下：
[0025][0026]
式中，z
i
与r
i
分别为时间序列中第i时刻的观测天顶角与夜光辐亮度，为拟合的夜光辐亮度，n为时间序列长度。
[0027]
建立误差方程，并求解二次曲线函数的系数矩阵，公式如下：
[0028]
v＝bx
‑
l
ꢀꢀ
(3)
[0029]
其中，
[0030]
x＝(b
t
b)
‑1b
t
l
ꢀꢀ
(4)
[0031]
式中，v表示拟合辐亮度与实测辐亮度之间的误差矩阵，b表示公式(1)中各系数对应的坐标组合项，x为公式(1)中待求的二次曲线函数的系数矩阵，l为实测辐亮度矩阵。
[0032]
而且，所述步骤6.1中将基于步骤5建立zrq统计模型分解为天顶角为0
°
时的夜光辐亮度(常数项系数)与二次函数值的乘积形式，公式如下：
[0033]
r＝c(t)
×
(a
′
z2+b
′
z+1)
ꢀꢀ
(5)
[0034]
式中，z代表传感器天顶角，c(t)代表分解成天顶角为0
°
的夜光辐亮度，(a
′
z2+b
′
z+1)代表夜光辐亮度分解后二次函数的值，a
′
、b
′
为二次函数的系数。
[0035]
而且，所述步骤7中cv与ndhd
ntl
计算公式如下：
[0036][0037][0038][0039][0040]
式中，ntl
hotspot
与ntl
darkspot
分别代表夜光辐亮度时间序列中最大、最小的夜光辐亮度，ntl
i
代表第i时刻的夜光辐亮度，n代表时间序列的个数。
[0041]
与现有技术相比，本发明具有如下优点：
[0042]
(1)本发明利用长时间序列的夜光辐亮度数据，建立夜光辐亮度与观测天顶角之间的统计模型，有助于更好地量化夜光辐亮度与观测天顶角之间的关系；
[0043]
(2)本发明基于天顶
‑
辐亮度二次模型构建夜光辐亮度角度归一化模型，考虑到了人类活动以及传感器测量误差等因素造成的时间序列波动，采用连续两次模型拟合并求解参数，提高了角度归一化的精度；
[0044]
(3)本发明通过最小二乘法多项式曲线拟合求解统计模型参数，运算效率较高，可以实现对长时间覆盖范围的夜光辐亮度角度归一化。
附图说明
[0045]
图1为本发明实施例的流程图。
[0046]
图2为本发明实施例中天顶
‑
辐亮度二次(zrq)模型的建立情况。
[0047]
图3为本发明实施例中洛杉矶研究点夜光辐亮度角度归一化前后时间序列对比示意图：其中图3(a)为洛杉矶1号研究点夜光辐亮度角度归一化前后时间序列对比示意图，图3(b)为洛杉矶2号研究点夜光辐亮度角度归一化前后时间序列对比示意图。
[0048]
图4为本发明实施例中拉斯维加斯研究点夜光辐亮度角度归一化前后时间序列对比示意图：其中图4(a)为拉斯维加斯1号研究点夜光辐亮度角度归一化前后时间序列对比示意图，图4(b)为拉斯维加斯2号研究点夜光辐亮度角度归一化前后时间序列对比示意图。
具体实施方式
[0049]
本发明提供一种卫星遥感夜光辐亮度的观测角度归一化方法，首先通过建立夜光辐亮度与观测天顶角之间的天顶
‑
辐亮度二次(zrq)模型，采用最小二乘多项式曲线拟合方式求解二次函数的系数，然后基于zrq模型，将观测到的夜光辐亮度分解为天顶角为0
°
时的夜光辐亮度与二次函数值的乘积，最后采用最小二乘法求解分解后的二次函数系数，并计算当前传感器获取的夜光辐亮度与二次函数值的比值，求解天顶角为0
°
时的夜光辐亮度，进而消除夜光辐亮度的角度效应。
[0050]
下面结合附图和实施例对本发明的技术方案作进一步说明。
[0051]
如图1所示，本发明实施例的流程包括以下步骤：
[0052]
步骤1，获取研究区域黑色大理石(black marble)夜光影像产品中vnp46a1与vnp46a2数据。
[0053]
步骤2，对步骤1获取到的所有vnp46a1与vnp46a2数据进行数据格式转换，并得到需要的质量文件与影像文件数据。
[0054]
步骤3，对步骤2所得到的质量标志数据，设置一定的筛选条件，得到符合需求的高质量夜光辐亮度数据。
[0055]
步骤4，以像素为研究单位，从步骤3得到的高质量夜光辐亮度数据中提取像素的夜光辐亮度时间序列和传感器观测天顶角时间序列。
[0056]
步骤5，对步骤4提取的夜光辐亮度与卫星传感器观测天顶角进行统计建模，确立两者之间的量化关系，并计算多项式系数，包括以下子步骤：
[0057]
步骤5.1，以观测天顶角作为自变量，夜光辐亮度作为因变量，建立两者之间的天顶
‑
辐亮度二次(zrq)拟合关系如下：
[0058]
r＝az2+bz+c
ꢀꢀ
(1)
[0059]
式中，z为卫星观测天顶角，r为夜光辐亮度，a、b、c分别为二次模型系数。
[0060]
步骤5.2，计算zrq模型的二次函数系数。
[0061]
假设夜光辐亮度时间序列为{r1,r2,...,r
n
}，传感器观测天顶角时间序列为{z1,z2,...,z
n
}，在夜光辐亮度与观测天顶角拟合二次曲线过程中定义拟合的目标函数使得拟合偏差最小，目标函数公式如下：
[0062][0063]
式中，z
i
与r
i
分别为时间序列中第i时刻的观测天顶角与夜光辐亮度，为拟合的夜光辐亮度，n为时间序列长度。
[0064]
建立误差方程，并求解二次曲线函数的系数矩阵，公式如下：
[0065]
v＝bx
‑
l
ꢀꢀ
(3)
[0066]
其中，
[0067]
x＝(b
t
b)
‑1b
t
l
ꢀꢀ
(4)
[0068]
式中，v表示拟合辐亮度与实测辐亮度之间的误差矩阵，b表示公式(1)中各系数对应的坐标组合项，x为公式(1)中待求的二次曲线函数的系数矩阵，l为实测辐亮度矩阵。
[0069]
步骤6，基于步骤5建立的天顶
‑
辐亮度二次(zrq)统计模型，消除夜光辐亮度角度效应，包括以下子步骤：
[0070]
步骤6.1，基于步骤5建立zrq统计模型，将模型公式分解为天顶角为0
°
时的夜光辐亮度(常数项系数)与二次函数值的乘积形式，公式如下：
[0071]
r＝c(t)
×
(a
′
z2+b
′
z+1)
ꢀꢀ
(5)
[0072]
式中，z代表传感器天顶角，c(t)代表分解成天顶角为0
°
的夜光辐亮度，(a
′
z2+b
′
z+1)代表夜光辐亮度分解后二次函数的值，a
′
、b
′
为二次函数的系数。
[0073]
步骤6.2，将步骤4提取的夜光辐亮度时间序列{r1,r2,...,r
n
}除以步骤5得到的二次模型的常数项系数c，获得新的夜光辐亮度时间序列并用新的夜光辐亮度时间序列与对应时刻的卫星传感器观测天顶角时间序列{z1,z2,...,z
n
}进行二次建模，使用如步骤5.2中的最小二乘法多项式曲线拟合求解二次函数系数a
′
、b
′
。
[0074]
步骤6.3，在求解得到系数之后，依次计算时间序列每一时刻传感器观测到的夜光辐亮度与步骤6.1分解得到的二次函数的比值，从而求解天顶角为0
°
时的夜光辐亮度。最终得到的天顶角为0
°
时的夜光辐亮度时间序列为：
[0075]
步骤7，计算夜光辐亮度时间序列角度归一化后的变异系数(cv)与归一化冷热点夜光辐亮度(ndhd
ntl
)，量化角度效应在夜光辐亮度时间序列中的影响程度。
[0076]
cv与ndhd
ntl
计算公式如下：
[0077][0078][0079][0080][0081]
式中，ntl
hotspot
与ntl
darkspot
分别代表夜光辐亮度时间序列中最大、最小的夜光辐亮度，ntl
i
代表第i时刻的夜光辐亮度，n代表时间序列的个数。
[0082]
下面以vnp46a1与vnp46a2实验数据为例，进一步阐述本发明的技术方案。
[0083]
步骤1，获取研究区域黑色大理石(black marble)夜光影像产品中vnp46a1与vnp46a2数据。
[0084]
从laads daac网站(https://ladsweb.modaps.eosdis.nasa.gov/)上下载viirs/npp卫星获取的每日夜光辐亮度产品vnp46a1与vnp46a2，数据的空间分辨率为15弧秒。影像涵盖2016年366天的数据，包括vnp46a1中的26个科学数据集以及vnp46a2中的7个科学数据集。本发明选择了包含洛杉矶与拉斯维加斯(h06v05)所在区域的夜光影像数据。
[0085]
步骤2，对下载获取到的所有vnp46a1与vnp46a2数据进行数据格式转换，并得到需要的质量文件与影像文件数据。
[0086]
将下载获取的hdf5格式科学数据集转换成geotiff格式，选取vnp46a1产品中传感
器观测天顶角、太阳天顶角、月光照度分数文件，以及vnp46a2产品中经过brdf校正的夜光影像、云层掩码、强制性质量指标文件。其中太阳天顶角、月光照度分数、云层掩码以及强制性质量指标文件作为后续像素筛选的质量标志数据，经过brdf校正的夜光影像与传感器观测天顶角将用于天顶角与夜光辐亮度的建模。
[0087]
步骤3，对步骤2中得到的质量标志数据，设置一定的筛选条件，得到符合需求的高质量夜光辐亮度数据。
[0088]
分别设置云层掩码指标、强制质量指标、太阳天顶角、月光照度分数作为逐像素夜光辐亮度质量的筛选条件，其中云层掩码是作为判别像素是否被云层遮挡的重要信息，其11位具体设置如表5.1所示。
[0089]
表5.1对云层掩码质量标志各个位数据设置规则
[0090][0091][0092]
但是在云层筛选的基础上，仍然存在部分像素值异常的情况，因此，将强制性质量指标设置为0，以保证所选像素均为高质量，从而提高数据质量。此外，剔除太阳天顶角大于108
°
或月光照度分数小于1％的像素，以保证获取像素的有效性和降低云层掩码标志数据质量错误的可能性。
[0093]
步骤4，以像素为研究单位，从步骤3得到的高质量夜光辐亮度数据中提取像素的夜光辐亮度时间序列和传感器观测天顶角时间序列。
[0094]
以当前像素3
×
3的窗口作为分析边界，首先保证3
×
3窗口内的9个像素经过条件筛选后均符合质量要求，然后再计算9个像素夜光辐亮度的均值作为当前像素的夜光辐亮度，传感器观测天顶角数据仅取当前像素位置对应的天顶角数据，最后分别读取2016年366天的夜光影像，提取出研究区域高质量的夜光辐亮度与观测角度时间序列。
[0095]
步骤5，对夜光辐亮度与卫星传感器观测天顶角进行统计建模，确立两者之间的量化关系，并计算多项式系数，包括以下子步骤：
[0096]
步骤5.1，以观测天顶角作为自变量，夜光辐亮度作为因变量，建立两者之间的天顶
‑
辐亮度二次(zrq)拟合关系如下：
[0097]
r＝az2+bz+c
ꢀꢀ
(1)
[0098]
式中，z为卫星观测天顶角，r为夜光辐亮度，a、b、c分别为二次模型系数。
[0099]
步骤5.2，计算zrq模型的二次函数系数。
[0100]
假设夜光辐亮度时间序列为{r1,r2,...,r
n
}，传感器观测天顶角时间序列为{z1,z2,...,z
n
}，在夜光辐亮度与观测天顶角拟合二次曲线过程中定义拟合的目标函数使得拟合偏差最小，目标函数公式如下：
[0101][0102]
式中，z
i
与r
i
分别为时间序列当中第i时刻的观测天顶角与夜光辐亮度，n为时间序列长度。
[0103]
建立误差方程，并求解二次曲线函数的系数矩阵，公式如下：
[0104]
v＝bx
‑
l
ꢀꢀ
(3)
[0105]
其中，
[0106]
x＝(b
t
b)
‑1b
t
l
ꢀꢀ
(4)
[0107]
式中，v表示拟合辐亮度与实测辐亮度之间的误差矩阵，b表示公式(1)中各系数对应的坐标组合项，x为公式(1)中待求的二次曲线函数的系数矩阵，l为实测辐亮度矩阵。
[0108]
步骤6，基于步骤5建立天顶
‑
辐亮度二次(zrq)统计模型，消除夜光辐亮度角度效应，包括以下子步骤：
[0109]
步骤6.1，基于步骤5建立zrq统计模型，将模型公式分解为天顶角为0
°
时的夜光辐亮度(常数项系数)与二次函数值的乘积形式，公式如下：
[0110]
r＝c(t)
×
(a
′
z2+b
′
z+1)
ꢀꢀ
(5)
[0111]
式中，z代表传感器天顶角，c(t)代表分解成天顶角为0
°
的夜光辐亮度，(a
′
z2+b
′
z+1)代表夜光辐亮度分解后二次函数的值，a
′
、b
′
为二次函数的系数。
[0112]
步骤6.2，将步骤4提取的夜光辐亮度时间序列{r1,r2,...,r
n
}除以步骤5得到的二次模型的常数项系数c，获得新的夜光辐亮度时间序列并用新的夜光辐亮度时间序列与对应时刻的卫星传感器观测天顶角时间序列{z1,z2,...,z
n
}进行二次建模，使用如步骤5.2中的最小二乘法多项式曲线拟合求解二次函数系数a
′
、b
′
。
[0113]
步骤6.3，在求解得到系数之后，依次计算时间序列每一时刻传感器观测到的夜光辐亮度与步骤6.1分解得到的二次函数的比值，从而求解天顶角为0
°
时的夜光辐亮度。最终得到的天顶角为0
°
时的夜光辐亮度时间序列为：
[0114]
步骤7，计算夜光辐亮度时间序列角度归一化后的变异系数(cv)与归一化冷热点夜光辐亮度(ndhd
ntl
)，量化角度效应在夜光辐亮度时间序列中的影响程度。
[0115]
cv与ndhd
ntl
计算公式如下：
[0116][0117][0118][0119][0120]
式中，ntl
hotspot
与ntl
darksport
分别代表夜光辐亮度时间序列中最大、最小的夜光辐亮度，ntl
i
代表第i时刻的夜光辐亮度，n代表时间序列的个数。
[0121]
角度校正前后不同区域研究点的归一化冷热点夜光辐亮度(ndhd
ntl
)与变异系数(cv)如表5.2结果所示，角度归一化使得时间序列更加平稳，总体cv下降超过30％，这表明归一化后的时间序列，由于消除了角度效应而降低了内部的差异性，从而提升了时间序列的稳定性。
[0122]
表5.2角度校正前后不同区域研究点的归一化冷热点夜光辐亮度(ndhd
ntl
)与变异系数(cv)
[0123][0124]
具体实施时，以上流程可采用计算机软件技术实现自动运行流程。
[0125]
本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。