一种夜光遥感影像相对辐射校正方法与流程

本发明涉及遥感影像校正技术领域，具体涉及一种夜光遥感影像相对辐射校正方法。

背景技术：

与多光谱的传统遥感影像相比，信息量较少的夜光遥感影像在反映人类社会经济活动方面更加突出，且由于其波段数较少，分辨率较低，影像数据量较小，大大减少了研究的计算量和存储容量要求。同时，由于时间跨度大，覆盖范围广等原因，美国军事气象卫星计划(defensemeteorologicalsatelliteprogram，dmsp)线性扫描业务系统(operationallinescansystem,ols)在夜光遥感领域得到广泛应用。

美国国防部发射的国防气象卫星dmsp提供从太阳同步卫星轨道拍摄的云层覆盖影像，主要用于监视气象、海洋和日地物理学。与此同时，它搭载了业务型线扫描系统(ols)，通过两个望远镜和光电倍增管(pmt)探测夜间影像，可见像素为从0到63。其独特的夜间光电放大能力，较长的时间覆盖范围及较大的空间覆盖范围使其在探索人类社会经济活动方面发挥着重大的作用。但由于卫星轨道的系统性偏差、传感器退化及缺乏星上校准，同一卫星获取的不同年份的同一地点的像元亮度值(dn)值存在波动，且由于大气层、地形起伏等的影响，不同卫星获取的同一年份的影像也存在差异。这些问题导致长时间序列的dmsp/ols影像一致性受到影响，因此，在利用长时间序列的dmsp/ols影像进行研究时，必须对其进行校正。

目前已经有很多学者研究了dmsp/ols夜光遥感影像的相对辐射校正方法，主要有不变目标区域的相互校正方法，稳健回归的全自动相对辐射校正，半自动的脊线采样回归方法。

不变目标区域的相互校正方法由elvidge等学者提出，该校正方法以灯光变化较少的西西里岛为校正的基准区域，以f121999影像作为校正的基准影像，通过二阶回归模型将其他年份的卫星影像校正以匹配f121999数据。在elvidge等人的基础上，liu等人通过应用年内和年际校正算法，利用土地利用和土地覆盖数据通过阈值分割的方法提取了城市信息。savory等人通过二阶回归模型校正dmsp/ols影像，使用高斯处理方法(gp)来平滑夜间灯光影像时间序列，描述了非洲城市增长的基本模式。吴建生等人首先对数据进行年际校正，再利用辅助数据ndvi来进行饱和校正。

为解决先验知识难以获得的问题，li等学者提出了稳健回归的全自动相对辐射校正算法，通过将两个不同日期影像的所有像素输入到线性回归模型中，迭代的丢弃异常值，最终获取稳定的像元并实施校正模型。在li等人基础上，bennie等人在方法细节方面做出了调整，使用通过中位数进行分位数回归估计的六阶多项式模型，减少了外围偏离值的影响。

本申请发明人在实施本发明的过程中，发现现有技术的方法，至少存在如下技术问题：

现有的方法有些需要较多的先验知识以确定校正的基准区域及基准影像，而这些先验知识是比较难获取的，部分自动校正方法，是在有限的地理范围内进行的，计算量较大，需要土地利用、土地覆盖和人口等多个参考数据集，并且没有充分验证。

由此可知，现有技术中的方法存在计算量大、校正效果不佳的技术问题。

技术实现要素：

有鉴于此，本发明提供了一种夜光遥感影像相对辐射校正方法，用以解决或者至少部分解决现有技术中的方法存在的计算量大、校正效果不佳的技术问题。

为了解决上述技术问题，本发明第一方面提供了一种夜光遥感影像相对辐射校正方法，包括：

步骤s1：获取研究区域所有dmsp/ols影像数据；

步骤s2：对获取的所有dmsp/ols影像数据进行重采样和图像配准，获得处理后的影像，处理后的影像包括校正基准影像和待校正影像；

步骤s3：建立相邻年份的待校正影像间密度图；

步骤s4：选取步骤s3中密度图的脊线点作为特征点，计算相邻年份的待校正影像的二阶多项式系数，其中，二阶多项式系数用于计算两相邻年份的待校正影像间dn值之间的关系；

步骤s5：基于二阶多项式系数和校正基准影像对后续待校正影像进行校正，获得校正后的第一影像。

在一种实施方式中，步骤s3中的密度图以校正基准影像的dn值范围作为x轴，以待校正影像的dn值范围作为y轴。

在一种实施方式中，步骤s4具体包括：

步骤s4.1：自动选取密度图中脊线上的点作为二阶回归模型的特征点，其中，特征点包括校正基准影像对应的第一特征点和待校正影像对应的第二特征点；步骤s4.2：利用公式(1)计算相邻年份的待校正影像的多项式系数，

y＝ax+bx²+c(1)

其中，y表示第一特征点的dn值，x表示第一特征点的dn值，a,b,c为利用最小二乘估计的二阶多项式系数。

在一种实施方式中，步骤s5具体包括：

步骤s5.1：获取校正基准影像所在年份和待校正影像所在年份；

步骤s5.2：当待校正影像所在年份小于校正基准影像所在年份时，计算待校正影像所在年份的原始影像与所在年份后一年的原始影像间的多项式系数，并将待校正影像所在年份的影像校正到所在年份后一年的辐射等级上，并将该经过校正的影像作为新的待校正影像；当待校正影像所在年份大于校正基准影像所在年份时，计算待校正影像所在年份的原始影像与所在年份前一年的原始影像间的多项式系数，并将待校正影像所在年份的影像校正到所在年份前一年的辐射等级上，并将该经过校正的影像作为新的待校正影像；

步骤s5.3：迭代执行步骤s5.2，直到将待校正影像所在年份的影像校正到校正基准影像所在年份的辐射等级上。

在一种实施方式中，所述方法还包括：

步骤s5.4：对所在年份小于校正基准影像年份的校正后影像，利用公式(2)进行连续性校正：

其中，n1＝1992,…,i年，dnn1为第n1年数据的dn值，dnn1-1为前一年数据的dn值；

步骤s5.5：对所在年份大于校正基准影像年份的校正后影像，利用公式(3)进行连续性校正：

其中，n2＝i,…,2013年，dnn2为第n2年数据的dn值，dnn2+1为后一年数据的dn值。

在一种实施方式中，所述方法还包括：

对于影像中对应同一地物的像元，对其附加限制条件，使得后一年的对应同一地物的像元，其dn值不小于前一年的dn值。

在一种实施方式中，在获得校正后的第一影像之后，所述方法还包括步骤s6：

对校正后的第一影像进行过饱和校正，获得最终校正影像。

在一种实施方式中，步骤s6具体包括：

将校正后的第一影像中像元的dn值大于dmsp/ols影像最大dn值63的影像校正为63，过饱和校正公式如下：

其中n3＝1992,…,2013年，dnn3为第n3年的dn值。

在一种实施方式中，所述方法还包括：

获取研究区域的区域矢量文件和gdp数据；

根据最终校正影像、研究区域的区域矢量文件以及gdp数据，计算最终校正影像研究区域内夜间灯光总量以及夜间灯光总量与gdp之间的相关性。

本申请实施例中的上述一个或多个技术方案，至少具有如下一种或多种技术效果：

本发明提供的一种夜光遥感影像相对辐射校正方法，首先获取研究区域所有dmsp/ols影像数据；并对获取的所有dmsp/ols影像数据进行重采样和图像配准，获得处理后的影像，处理后的影像包括校正基准影像和待校正影像；然后，建立相邻年份的待校正影像间密度图；接着选取密度图的脊线点作为特征点，计算相邻年份的待校正影像的二阶多项式系数；再基于二阶多项式系数和校正基准影像对后续待校正影像进行校正，获得校正后的第一影像。

与现有技术相比，本发明通过计算获得相邻两年待校正影像的伪不变特征点，并进行连续相对辐射校正，提高了相对辐射校正的精度。并且通过重采样及脊线采样回归的方法降低了运算量，实现了对全球数据的校正。

进一步地，对于影像中对应同一地物的像元，对其附加限制条件，使得后一年的对应同一地物的像元，其dn值不小于前一年的dn值，通过附加限制条件提高了相对辐射校正的连续性及其夜间灯光总量与社会经济活动之间的相关性。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明的夜光遥感影像相对辐射校正方法的流程示意图；

图2为本发明实施例中的密度图构建示意图；

图3为本发明实施例中校正关系示意图。

具体实施方式

本申请发明人通过大量的研究与实践发现：

zhang等学者提出的半自动的脊线采样回归相对辐射校正方法，以辐射亮度密度图脊线上的点作为特征点生成校正模型，从而在全球范围内创建一致的时间序列。该方法选取固定的基准影像，当待校正影像与基准影像相隔时间较长时，伪不变特征点将会减少，从而影响校正精度。

综合前人研究，目前dmsp/ols夜光遥感影像相对辐射校正方面算法都不够成熟，需要较多的先验知识，计算量大，地理范围有限，在距校正基准年份较远的情况下校正精度受到影响，校正结果的覆盖范围较小，夜间灯光总量的连续性及其与人类社会经济活动的相关性都不高。基于此，本发明提供了一种夜光遥感影像相对辐射校正方法，主要构思如下：

以位于原始影像时间序列中部的影像作为校正基准影像，以时间序列中的其他影像作为待校正影像，从校正基准影像开始，通过建立相邻两年待校正影像的密度图获得脊线，缩短两景影像间的时间差距，使地面上点的变化减少，进而可以获得更多的伪不变特征点，自动选取脊线上的点作为特征点构建二阶回归方程以使伪不变特征点数达到最多，利用构建的二阶回归方程对影像进行连续的相对辐射校正，同时通过对其施加限制条件保证影像响的连续性，最终提高影像相对辐射校正的精度和夜间灯光总量的连续性及其与社会经济活动的相关性。

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

实施例一

本实施例提供了一种夜光遥感影像相对辐射校正方法，该方法包括：

步骤s1：获取研究区域所有dmsp/ols影像数据。

具体实施时，dmsp/ols影像数据即为dmsp/ols稳定灯光数据。本步骤从ngdc网站下载1992年至2013年共34个卫星年(版本四)的数据。具体为：从f10,f12,f14,f15,f16,f18共六颗卫星获取，涵盖1992年至2013年，34卫星年(版本四)的数据，包括三种年度平均数据：无云覆盖，平均可见光和稳定光。选择覆盖城镇、乡镇和其他稳定发光光源并去除偶然噪声影响后的稳定灯光数据，数据的空间分辨率为30弧秒。

步骤s2：对获取的所有dmsp/ols影像数据进行重采样和图像配准，获得处理后的影像，处理后的影像包括校正基准影像和待校正影像。

具体实施时，采用最近邻域法将像元大小重采样为原来的0.5倍，并采用一次多项式进行图像配准。

步骤s3：建立相邻年份的待校正影像间密度图。

其中如图2所示，密度图的x轴为校正基准影像的dn值范围(0～63)，y轴为待校正影像的dn值范围(0～63)。密度图中密度最大的点构成的线称为脊线，y轴对应的一列中，密度最大的值为脊线上对应的点，在图二中表示为白色的点。

步骤s4：选取步骤s3中密度图的脊线点作为特征点，计算相邻年份的待校正影像的二阶多项式系数，其中，二阶多项式系数用于计算两相邻年份的待校正影像间dn值之间的关系。

具体来说，相邻年份的待校正影像的二阶多项式系数是指两个相邻年份的待校正影像(例如1992年的待校正影像和1993年的待校正影像)之间的二阶多项式系数，计算二阶多项式系数的目的是对两影像的dn值关系进行拟合，进而得到两影像dn值之间的关系，从而利用计算出的关系对后续影像进行校正。

步骤s5：基于二阶多项式系数和校正基准影像对后续待校正影像进行校正，获得校正后的第一影像。

具体来说，利用多项式系数进行的校正是将影像每一个像元对应的dn值进行多项式校正。根据计算出的二阶多项式系数，可以得到二阶回归模型，则可以对两相邻影像进行拟合，得到两影像dn值之间的关系，进而进行校正。

在一种实施方式中，步骤s4具体包括：

步骤s4.1：自动选取密度图中脊线上的点作为二阶回归模型的特征点，其中，特征点包括校正基准影像对应的第一特征点和待校正影像对应的第二特征点；

步骤s4.2：利用公式(1)计算相邻年份的待校正影像的多项式系数，

y＝ax+bx²+c(1)

其中，y表示第一特征点的dn值，x表示第一特征点的dn值，a,b,c为利用最小二乘估计的二阶多项式系数。

在一种实施方式中，步骤s5具体包括：

步骤s5.1：获取校正基准影像所在年份和待校正影像所在年份；

步骤s5.2：当待校正影像所在年份小于校正基准影像所在年份时，计算待校正影像所在年份的原始影像与所在年份后一年的原始影像间的多项式系数，并将待校正影像所在年份的影像校正到所在年份后一年的辐射等级上，并将该经过校正的影像作为新的待校正影像；当待校正影像所在年份大于校正基准影像所在年份时，计算待校正影像所在年份的原始影像与所在年份前一年的原始影像间的多项式系数，并将待校正影像所在年份的影像校正到所在年份前一年的辐射等级上，并将该经过校正的影像作为新的待校正影像；

步骤s5.3：迭代执行步骤s5.2，直到将待校正影像所在年份的影像校正到校正基准影像所在年份的辐射等级上。

在具体的实施过程中，

步骤5.1中，设i为校正的基准影像所在年份，x＜i，x为待校正影像所在年份小于i的年份，y＞i，且y表示待校正影像所在年份大于i的年份。对于小于基准影像所在年份的待校正年份，计算x与x+1之间的多项式系数，并将第x年的影像校正到x+1年的辐射等级上。对于大于基准影像所在年份的待校正年份，计算y与y-1之间的多项式系数，并将第y年的影像校正到y-1年的辐射等级上；

步骤5.2中，对于小于基准影像所在年份的待校正年份，将x+1视为步骤5.1的x，将x+2视为步骤5.1的x+1，计算得到x+2年辐射等级上的第x年的校正结果。对于大于基准影像所在年份的待校正年份，将y-1视为步骤5.1的y，将y-2视为步骤5.1的y-1，计算得到y-1年辐射等级上的第y年的校正结果；

步骤5.3，迭代进行步骤5.1，直到将第x年的影像校正到第i年的辐射等级上。

具体实施过程中，本步骤的具体操作如下：

如图3及公式5～8所示，选取f152002作为基准影像，计算f101992与f101993年的多项式系数为p1，计算f101993与f101994年的多项式系数为p2，……，f152001与f152002年的多项式系数为p10，通过p1将f101992的影像校正到f101993的辐射等级上，得到f101992⁽¹⁾，利用p2将f101992⁽¹⁾校正到f101994的辐射等级上，得到f101992⁽²⁾，……，利用p10将f101992⁽⁹⁾校正到f152002的辐射等级上，得到f101992⁽¹⁰⁾。

f101992⁽¹⁾＝fp1(f101992)(5)

f101992⁽²⁾＝fp2(fp1(f101992))(6)

f101992⁽⁹⁾＝fp9(fp8(fp7(fp6(fp5(fp4(fp3(fp2(fp1(f101992)))))))))(7)

f101992(¹⁰)＝fp10(fp9(fp8(fp7(fp6(fp5(fp4(fp3(fp2(fp1(f101992))))))))))(8)

在一种实施方式中，所述方法还包括：

步骤s5.4：对所在年份小于校正基准影像年份的校正后影像，利用公式(2)进行连续性校正：

其中，n1＝1992,…,i年，dnn1为第n1年数据的dn值，dnn1-1为前一年数据的dn值；

步骤s5.5：对所在年份大于校正基准影像年份的校正后影像，利用公式(3)进行连续性校正：

其中，n2＝i,…,2013年，dnn2为第n2年数据的dn值，dnn2+1为后一年数据的dn值。

在一种实施方式中，所述方法还包括：

对于影像中对应同一地物的像元，对其附加限制条件，使得后一年的对应同一地物的像元，其dn值不小于前一年的dn值。

在一种实施方式中，在获得校正后的第一影像之后，所述方法还包括步骤s6：

对校正后的第一影像进行过饱和校正，获得最终校正影像。

其中，步骤s6具体包括：

将校正后的第一影像中像元的dn值大于dmsp/ols影像最大dn值63的影像校正为63，过饱和校正公式如下：

其中n3＝1992,…,2013年，dnn3为第n3年的dn值。

在一种实施方式中，所述方法还包括：

获取研究区域的区域矢量文件和gdp数据；

根据最终校正影像、研究区域的区域矢量文件以及gdp数据，计算最终校正影像研究区域内夜间灯光总量以及夜间灯光总量与gdp之间的相关性。

具体来说，可以从网络下载中国地区的矢量文件，从中国科学院资源环境科学数据中心(http://www.resdc.cn/)下载2015年中国地市行政边界数据(矢量文件)。从世界银行组织(https://www.worldbank.org/)获取研究区域1992年至2013年gdp数据。

具体地，图1为夜光遥感影像相对辐射校正方法的流程示意图，回归方程即二阶回归模型，即公式(1)，相关性计算是指步骤s7。

夜间灯光总量指一景影像中所有像元dn值之和，利用arcgis统计中国行政区域内的夜间灯光总量，相关性由pearson相关系数表示。制作1992年至2013年夜间灯光总量折线图，计算夜间灯光总量与gdp间pearson相关系数。

下面将结合具体实例应用进一步说明本发明的技术方案及有益效果。

选择1992年至2013年dmsp/ols所有稳定灯光数据(即影像数据)，中国区域矢量数据及gdp数据。首先对获取的所有dmsp/ols稳定夜间灯数据进行重采样和图像配准，然后对经过重采样和图像配准的数据，建立相邻两幅影像间的密度图，选取密度图脊线点作为特征点，计算相邻两幅影像之间的二阶多项式系数，利用多项式系数对影像进行连续相对辐射校正并进行过饱和校正，最后，计算校正后中国区域内夜间灯光总量及其与gdp间相关性。结果显示其夜间灯光总量的连续性有了很大程度的提高，夜间灯光总量及其与gdp间相关性pearson系数由校正前的0.9444提高到0.9511，表明校正后夜光遥感影像与人类社会经济活动之间的相关性得到了提高。

尽管已描述了本发明的优选实施例，但本领域内的技术人员一旦得知了基本创造性概念，则可对这些实施例做出另外的变更和修改。所以，所附权利要求意欲解释为包括优选实施例以及落入本发明范围的所有变更和修改。

显然，本领域的技术人员可以对本发明实施例进行各种改动和变型而不脱离本发明实施例的精神和范围。这样，倘若本发明实施例的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也意图包含这些改动和变型在内。