CN105930867A - 一种基于时间序列夜光遥感影像的数据挖掘方法 - Google Patents

Abstract

本发明涉及一种基于时间序列夜光遥感影像的数据挖掘方法，其特征在于包括以下内容：1)读取时间序列夜光遥感影像，并从读取的时间序列夜光遥感影像中选择某几段时序的时间序列夜光遥感影像；2)对选取的某几段时间序列夜光遥感影像分别进行重投影和重采样，并获得重采样影像；3)将重采样影像合成具有多个波段的合成影像；4)对合成影像进行数据挖掘得到影像分类结果。本发明可以广泛应用于对夜光遥感影像的数据挖掘中。

Description

一种基于时间序列夜光遥感影像的数据挖掘方法

技术领域

本发明涉及一种基于时间序列夜光遥感影像的数据挖掘方法，属于遥感影像处理技术领域。

背景技术

夜光遥感影像是指在夜间探测到的地表微弱近红外辐射，从而可以获取不受日光影响的来自城市、乡镇等含有持久光源的夜光数据，夜光遥感影像将城镇表征为块状或片状亮斑区域(即灯光区)，使之明显区别于黑暗的乡村背景(即无灯光区)，有效避免了传统遥感影像产生的城镇区域与非城镇区域混淆的问题，为监测人类经济活动提供一种独特的数据获取手段。夜光遥感影像相对于白天遥感影像来说，具有不受日光、阴影、植被及其他地物影响的特征，并且可以探测低强度灯光，更加利于进行城市化强度及其时空分异分析与空间扩展分析的相关研究。相对于高分辨率数据，夜光遥感影像数据量小、易于获取且获取成本低，在长时间序列上能够客观、快速、准确的从整体把握区域的扩展过程，在研究上具有重要意义，同时也更适合在区域大尺度上的研究。随着技术不断发展，灯光数据产品精度不断提高，因而可以充分利用现有数据进行长时间序列制图，使得对城市发展研究的探测提供更多依据。目前，该数据已被广泛应用于大尺度城镇空间扩展研究、经济与人口估算、城市用电量和能源消耗分析、碳排放和光污染等环境问题评估。除此之外，夜光遥感影像还可以捕捉到夜间渔船、天然气燃烧、森林火灾发光等，因此广泛应用于社会经济参数估算、区域发展研究、重大事件评估、渔业监测等诸多研究领域。

美国国防气象卫星(DMSP)搭载的业务型线扫描传感器(OLS)以及美国新一代极轨运行环境卫星系统预备卫星Suomi国家极轨合作伙伴(Suomi-NPP)搭载的可见光红外成像辐射仪(VIIRS)获取的白天/夜间波段(DNB)夜间遥感影像由于提供了夜晚可见光波段内观测的定量测量方法，在特征物(例如烟尘、冰雪、火山、海洋表面粗糙度等)的探测以及光源监测(例如城市灯光、火情、船舶灯光、闪电、极光等)方面得到了广泛的应用。这些探测应用多数是在单时相影像产品上进行相应灯光检测，提取城镇等特定兴趣区，再在此基础上进行特征与其他相关经济社会参数的分析或提取特征间变化比较分析，例如选取一定的阈值对不同影像中的同类特征进行提取，然后利用提取的区域进行相关分析。由于夜光遥感影像为灰度影像，利用灰度分割方式取得的信息有限。另外，由于云及其他噪声因素的影响，使得一部分夜光遥感产品质量较低，不利于进行单幅影像分析，使得传统方法存在一定的局限性。更重要的是，上述方法往往忽略了不同特征在时间序列夜光遥感影像上灰度强度是变化的，而这些随时间变化的特性本身就具有很大的信息量，已有方法无法考虑这些信息，也不利于挖掘由于影像时间序列变化及相互影响而产生的影像特征。综上所述，针对夜光遥感影像特征提取与分析等深入的数据挖掘方法还有巨大的提升空间。

发明内容

针对上述问题，本发明的目的是提供一种能够提高特征提取类型数和特征分割精度的基于时间序列夜光遥感影像的数据挖掘方法。

为实现上述目的，本发明采取以下技术方案：一种基于时间序列夜光遥感影像的数据挖掘方法，其特征在于包括以下内容：1)读取时间序列夜光遥感影像，并从读取的时间序列夜光遥感影像中选择某几段时序的时间序列夜光遥感影像；2)对选取的某几段时间序列夜光遥感影像分别进行重投影和重采样，并获得重采样影像；3)将重采样影像进行合成获得具有多个波段的合成影像；4)对合成影像进行数据挖掘得到影像分类结果。

进一步，所述步骤2)对选取的某几段时间序列夜光遥感影像分别进行重投影和重采样，并获得重采样影像，具体过程为：2.1)根据选取的每一段时间序列夜光遥感影像的地理信息范围，确定重投影后影像范围，并获取影像覆盖区域；2.2)根据影像覆盖区域左上角点坐标及分辨率确定重投影后的影像大小；2.3)对影像进行重采样，即根据步骤2.1)及2.2)得到的重投影影像大小、分辨率和左上角点坐标，对选取的时间序列夜光遥感影像进行重采样，为构建的重投影影像各个像素点赋值。

进一步，所述步骤2.3)对影像进行重采样的具体过程为：2.31)从重投影后影像的左上角像素点坐标开始，逐像素遍历步骤2.2)中构建的重投影影像，将重投影影像上的左上角点坐标计算到重投影前时间序列夜光遥感影像上；2.32)在重投影前的时间序列夜光遥感影像上进行重采样，得到重投影影像中对应像元的DN值；2.33)对覆盖区域内的所有像素点重复步骤2.31)及2.32)，得到重采样影像。

进一步，所述步骤4)对合成影像采用非监督分类方法进行数据挖掘。

进一步，所述非监督分类方法为迭代自组织数据分析方法，对合成影像进行数据挖掘采用迭代自组织数据分析方法，得到影像分类结果的具体情况为：4.1)初始化，假设初始聚类数c以及期望得到的聚类数为K，确定初始中心m_i，i＝1,…,c：

$m_i=\frac{1}{N_i}\underset{{\mathrm{DN}}_k\∈{\mathrm{\Γ}}_i}{\Σ}{\mathrm{DN}}_k$

式中，N_i为第i类聚类Γ_i中的样本数目，DN_k为聚类Γ_i中的样本；4.2)将所有的样本分到距初始中心最近的类Γ_i中，若某个类Γ_j中的样本数量N_j小于一个聚类中的最小样本数θ_N，则去掉该类，令c＝c-1；4.3)重新计算聚类中心m_j为：

$\begin{array}{cc}{m}_j=\frac{1}{N_j}\underset{{\mathrm{DN}}_k\∈{\mathrm{\Γ}}_j}{\Σ}{\mathrm{DN}}_k,& j=1,\mathrm{...},c\end{array}$

式中，N_j为第j个聚类中的样本数目；4.4)计算第j类样本与其中心的平均距离和总平均距离

$\begin{array}{cc}{\stackrel{\‾}{\mathrm{\δ}}}_j=\frac{1}{N_j}\underset{{\mathrm{DN}}_k\∈{\mathrm{\Γ}}_j}{\Σ}||{\mathrm{DN}}_k-m_j||,& j=1,\mathrm{...},c\end{array}$

$\stackrel{\‾}{\mathrm{\δ}}=\frac{1}{N}\underset{j=1}{\overset{c}{\Σ}}{N}_j{\stackrel{\‾}{\mathrm{\δ}}}_j$

式中，N为总样本数目；4.5)设置迭代次数I并进行迭代，具体迭代过程为：i)若c≤K/2，或者2K≤c≤K/2且迭代次数为奇数，则进行分裂：对每个类求取各维的标准偏差σ_j＝[σ_j1,σ_j2,…,σ_jd]^T：

$\begin{array}{ccc}{\mathrm{\σ}}_{ji}=\sqrt{\frac{1}{N_j}\underset{{\mathrm{DN}}_k\∈{\mathrm{\Γ}}_j}{\Σ}{({\mathrm{DN}}_{ki}-m_{ji})}^2},& j=1,\mathrm{...},c,& i=1,\mathrm{...},d\end{array}$

式中，DN_ki为第k个样本的第i个分量，m_ji是当前第j个聚类中心的第i个分量，σ_ji是第j类第i个分量的标准差，d是样本维数；对每个类求出标准偏差最大的分量σ_jmax，j＝1,…,c；对各类的σ_jmax(j＝1,…,c)，若存在某个类的σ_jmax>标准偏差参数θ_s，且且N_j>2(θ_N+1)，或c≤K/2，则将Γ_j分裂为两类，中心分别为和置c＝c+1，分裂项γ_j＝[0,…,0,σ_jmax,0,…,0]^T；ii)若c≥2K，或者2K≤c≤K/2且迭代次数为偶数，则进行合并：计算各聚类中心两两之间距离δ_ij：δ_ij＝||m_i-m_j||，i,j＝1,…,c，i≠j，比较δ_ij与合并参数θ_C，l为小于θ_C的个数，对小于θ_C者排序：

${\mathrm{\&delta;}}_{i_1{j}_1}<{\mathrm{\&delta;}}_{i_2{j}_2}<...<{\mathrm{\&delta;}}_{i_l{j}_l}$

从最小的开始，把每个对应的和合并，组成新类，新类的中心为：

$m_l=\frac{1}{N_{i_l}+N_{j_l}}\&lsqb;N_{i_l}{m}_{i_l}+N_{j_l}{m}_{j_l}\&rsqb;$

并令c＝c-1，每次迭代中避免同一类被合并两次；4.6)若达到设置的最大迭代次数，则终止迭代；否则，迭代次数加1，进入步骤4.2)。

进一步，当提取特征类别已知的情况下，对合成影像采用基于机器学习方法进行数据挖掘，得出影像分类结果。

进一步，所述基于机器学习方法为支持向量机方法，对合成影像进行数据挖掘采用支持向量机方法，得到影像分类结果的具体过程为：4.1)选取特征样本数据集，样本数据集中的特征向量中记录了每类特征样本各个波段的影像亮度值及样本类别，记为：(dn₁,y₁),(dn₂,y₂),…,(dn_m,y_m)，m为样本个数；4.2)对合成影像采用支持向量机方法提取影像特征，并对合成影像进行分类，具体步骤如下：4.21)构建支持向量机预测函数，即：通过样本数据集构建分类支持向量机目标函数，求取a_i及b：

$\max Q\left(\mathrm{\&alpha;}\right)\underset{i=1}{\overset{n}{\&Sigma;}}{a}_i-\frac{1}{2}\underset{i,j=1}{\overset{n}{\&Sigma;}}{a}_i{a}_j{y}_j{y}_{i}K({\mathrm{dn}}_i\&CenterDot;{\mathrm{dn}}_j)$

$s.t.{\&Sigma;}_{i,j=1}^n{a}_i{y}_i=0$

核函数：K(dn_i·dn_j)＝((dn_i·dn_j)+1)²

式中，a_i，a_j为拉格朗日乘子，dn_i，dn_j为样本数据集特征向量，y_i，y_j＝±1为数据点类别，利用二次规划方法求解目标函数的最优解，得到最优拉格朗日乘子a_i，利用样本数据集中的支持向量dn，带入判别函数，求解偏差值b，其中，判别函数表达式为：

$f\left(x\right)=\mathrm{sgn}\left(\underset{i=1}{\overset{n}{\&Sigma;}}{a}_i{y}_{i}K\right({\mathrm{dn}}_i\&CenterDot;{\mathrm{dn}}_j)+b)$

式中，sgn为判别函数，若内部返回值大于0，则返回1；等于0，则返回0，小于0，则返回-1；根据上述方法分别构建n个支持向量机决策函数，其中，n为所划分的特征类别个数，即每个类别构建一个支持向量机决策函数；4.22)采用步骤4.21)构建出的支持向量机决策函数对特征向量进行分类：将特征向量带入第1类的支持向量机决策函数，进行分类，若f(x)为1，则该特征向量属于第1类，若f(x)为-1，则该特征向量不属于第1类，将不属于该类的特征向量带入第2类的支持向量机决策函数进行分类，以此类推，最终得出划分为n类特征的合成影像。

本发明由于采取以上技术方案，其具有以下优点：1、本发明以时间序列的夜光遥感影像进行数据挖掘，由于不同时间夜光遥感影像具有的信息不同，几何时间序列的夜光遥感影像可以形成信息互补，在集成了同一地区灯光数据的共性与不同时序所造成的灯光数据特性条件下，可以更加详尽地在夜光遥感影像上表征各种类型特征，具有更高的可挖掘性，从而有效提高特征提取类型数和特征的分割精度。2、本发明采用非监督分类方法或基于机器学习方法对时间序列夜光遥感影像合成后的产品进行数据挖掘，能够有效提取时间遥感影像中由于时序变化导致的灯光信息特征，由于时间序列图像的引入，使得整体信噪比增加，从而更加可以突出有用信息，消除噪声的影响，在对单幅夜光影像很难提取的弱特征(如道路、城区边界等)提取时具有较好的效果。综上所述，本发明可以广泛应用于不同时间序列的夜光遥感影像的数据挖掘中，例如美国国防气象卫星(Defense Meteorological SatelliteProgram,DMSP)搭载的业务型线扫描传感器(Operational Linescan System,OLS)年平均数据、美国新一代极轨运行环境卫星系统预备卫星Suomi国家极轨合作伙伴(Soumi NationalPolar-orbiting Partnership,Suomi-NPP)搭载的可见光红外成像辐射仪(VIIRS)的白天/夜间波段(DNB)获取的月平均数据等。

附图说明

图1是本发明的数据挖掘方法流程示意图。

具体实施方式

以下结合附图来对本发明进行详细的描绘。然而应当理解，附图的提供仅为了更好地理解本发明，它们不应该理解成对本发明的限制。

如图1所示，本发明提供的基于时间序列夜光遥感影像的数据挖掘方法，包括以下内容：

1、读取时间序列夜光遥感影像。

夜光遥感数据可以根据不同的获取时间制作成年平均夜光遥感影像、月平均夜光遥感影像产品等。年平均影像像元的DN(灰度)值代表了年平均灯光辐射亮度值(即面辐射源上某点在一定方向上辐射强弱的物理量)，是对每年中白天/夜晚波段数据中各像元无云观测次数及其对应进行过辐射校正的辐射亮度值进行平均化处理得到的。月平均影像像元的DN值代表了月平均灯光辐射亮度值，不同时间序列即为不同月份/年份的平均夜光遥感影像。

本发明以VIIRS/DNB月平均数据为具体实施例进行说明，但是不限于此，VIIRS/DNB月平均数据可以通过美国国家海洋和大气管理局网站下载获得，下载得到的数据为月平均栅格影像数据，参考系为WGS_84坐标系，夜光遥感影像的扫描带宽为3000km，水平空间分辨率为交轨<0.77km，沿轨<0.75km。下载得到的夜光遥感影像可以通过envi等软件或自行编程的方式读入多幅不同时间序列的同类型夜光遥感影像，获取夜光遥感图像上每个点的DN值，并获取每幅夜光遥感影像的地理坐标信息。

由于夜光遥感影像是以一定的地面分辨率进行排列的，所以影像的像平面坐标(x,y)与其投影坐标间(D_B，D_L)也具有线性的关系，具体为：

$\begin{array}{c}{D}_B=a_i+x*{\mathrm{dx}}_i\\ D_L=b_i+y*{\mathrm{dy}}_i\end{array}---\left(1\right)$

式中，a_i,b_i,dx_i,dy_i为影像存储的地理坐标信息，i为影像编号，a_i,b_i为影像左上角点的大地坐标，dx_i,dy_i为影像的列方向和行方向的分辨率。

2、根据读取的时间序列夜光遥感影像数量、时间范围信息并结合数据挖掘目的，从读取的时间序列夜光遥感影像中选择某几段时序合适的夜光遥感影像。

根据数据挖掘的具体目标，选择相应时间序列夜光遥感影像，例如需要提取某段时间城市边界的变化区域，则可以选择相应时间段内每个月的月平均数据进行处理，如果没有特定目的，则可以以一年为周期，选择一年内每个月的月平均数据进行处理，具体可以根据实际需要进行选取，在此不做限定。

3、对选取的某几段时间序列夜光遥感影像进行重投影和重采样，获得重采样影像，本发明中的影像重投影是指将地理坐标影像通过建立点与点之间函数关系的方法投影为平面坐标，即将曲面转化为平面的过程，具体过程为：

1)根据选取的时间序列夜光遥感影像的地理信息范围，确定重投影后影像范围，并获取影像覆盖区域。

根据公式(1)获取步骤2中选取的每幅夜光遥感影像的左上角像素点投影坐标及右下角点投影坐标i为影像编号，求取所有夜光遥感影像的最大值及的最大值作为夜光遥感影像公共覆盖区域的左上角点坐标同时求取所有夜光遥感影像的最小值及的最小值作为夜光遥感影像公共覆盖区域右下角点坐标

2)根据影像覆盖区域左上角点坐标及分辨率确定重投影后的影像大小。

将步骤1)得出的公共覆盖区域的左上角点及右下角点坐标带入公式(2)，求取重投影后的影像大小，假设重投影后影像宽度为w，高度为h：

$\begin{array}{c}w=(D_{L_i}^{rd}-D_{L_i}^{lt})/dx\\ h=(D_{B_i}^{rd}-D_{B_i}^{lt})/dy\end{array}---\left(2\right)$

式中，dx、dy为重投影后的影像分辨率，设定为原始输入影像分辨率的平均值。根据计算得出的重投影影像的高度、宽度及分辨率构建一幅重投影影像，此时影像中各像素点的值为空值。

3)影像重采样，即根据步骤1)、步骤2)得到的重投影影像大小、分辨率和左上角点坐标，对选取的时间序列夜光遥感影像进行重采样，为步骤2)构建的重投影影像各个像素点赋值，具体过程为：

3.1)从重投影后影像的左上角像素点坐标开始，逐像素遍历步骤2)中构建的重投影影像，将重投影影像上的左上角点坐标(m，n)计算到重投影前时间序列夜光遥感影像(x,y)上，计算公式如下：

$\begin{array}{c}{D}_B=D_{B_i}^{lt}+m*dx\\ D_L=D_{L_r}^{lt}+n*dy\\ x=(D_B-a_i)/{\mathrm{dx}}_i\\ y=(D_L-b_i)/{\mathrm{dy}}_i\end{array}---\left(3\right)$

3.2)在重投影前的时间序列夜间遥感影像上进行采样，得到重投影影像中对应像元的DN值。

由于(x,y)非整像素，本发明利用卷积方式进行采样，本发明实施例以双三次卷积方法进行说明，但是不限于此，具体为：根据步骤3.1)中重投影后影像左上角像素点坐标(m,n)计算得出的重投影前的时间序列夜光遥感影像中左上角像素点(x,y)，令i代表行数，j代表列数，α＝x-i,β＝y-j,i＝Integer(x),j＝Integer(y)，其中，Integer为取整。取与计算点(x,y)周围相邻的16个像素点参与计算，在x方向上，设m＝j-1,j,j+1,j+2，计算公式为：

$\begin{array}{c}f(x,m)={\mathrm{\&alpha;}}^2(\mathrm{\&alpha;}-1)f(i+2,m)+\mathrm{\&alpha;}(1+\mathrm{\&alpha;}-{\mathrm{\&alpha;}}^2)f(i+1,m)\\ +(1-2{\mathrm{\&alpha;}}^2+{\mathrm{\&alpha;}}^3)f(i,m)-\mathrm{\&alpha;}{(1-\mathrm{\&alpha;})}^{2}f(i-1,m)\end{array}---\left(4\right)$

采用同样的方法计算y方向的值f(x,y)，即该像素点的亮度值，计算公式为：

$\begin{array}{c}f(x,y)={\mathrm{\&beta;}}^2(\mathrm{\&beta;}-1)f(x,j+2)+\mathrm{\&beta;}(1+\mathrm{\&beta;}-{\mathrm{\&beta;}}^2)f(x,j+1)\\ +(1-2{\mathrm{\&beta;}}^2+{\mathrm{\&beta;}}^3)f(x,j)-\mathrm{\&beta;}{(1-\mathrm{\&beta;})}^{2}f(x,j-1)\end{array}---\left(5\right)$

3.3)对重采样后夜光遥感影像的每个像素点重复步骤3.1)及3.2)，对步骤2)中所有获得的覆盖区域内的所有像素点的亮度值分别进行计算，得到重采样影像。

4、通过影像波段合成方法合成得到具有多个波段的合成影像，波段数即为步骤2中选取的影像数，本发明以BSQ(Band sequential)数据格式存储合成影像，但是不限于此，具体过程为：按照波段顺序(步骤2中选取的影像顺序)依次排列的数据格式将所有影像存储为一幅合成影像，各波段影像内部还按照原始影像像素排列顺序进行存储。

5、对合成影像进行数据挖掘，对合成影像通过非监督分类方法进行影像数据挖掘，得到影像分类结果，具体情况为：

针对步骤4中产生的合成影像，逐像素遍历影像各像素点DN值，构建特征向量(DN₁，DN₂，……，DN_n)，然后对这些特征向量进行非监督分类。本发明的非监督分类方法以迭代自组织数据分析方法(ISODATA)为实施例进行说明，但是不限于此，具体过程为：

5.1)初始化，假设初始聚类数c以及期望得到的聚类数为K，本实施例中c取为6，但不限于此，采用下式确定初始中心m_i，i＝1,…,c。

$m_i=\frac{1}{N_i}\underset{{\mathrm{DN}}_k\&Element;{\mathrm{\&Gamma;}}_i}{\&Sigma;}{\mathrm{DN}}_k---\left(6\right)$

式中，N_i为第i类聚类Γ_i中的样本数目，DN_k为聚类Γ_i中的样本。

5.2)将所有的样本分到距距离中心最近的类Γ_i中，若某个类Γ_j中的样本数量N_j小于一个聚类中的最小样本数θ_N，则去掉该类(根据各样本到其它类中心的距离分别合入其他类)，令c＝c-1。

5.3)重新计算聚类中心m_j为：

$\begin{array}{cc}{m}_j=\frac{1}{N_j}\underset{{\mathrm{DN}}_k\&Element;{\mathrm{\&Gamma;}}_j}{\&Sigma;}{\mathrm{DN}}_k,& j=1,...,c\end{array}---\left(7\right)$

式中，N_j为第j个聚类中的样本数目。

5.4)计算第j类样本与其中心的平均距离和总平均距离

$\begin{array}{c}\begin{array}{cc}{\stackrel{\&OverBar;}{\mathrm{\&delta;}}}_j=\frac{1}{N_j}\underset{{\mathrm{DN}}_k\&Element;{\mathrm{\&Gamma;}}_j}{\&Sigma;}||{\mathrm{DN}}_k-m_j||,& j=1,\mathrm{...},c\end{array}\\ \stackrel{\&OverBar;}{\mathrm{\&delta;}}=\frac{1}{N}\underset{j=1}{\overset{c}{\&Sigma;}}{N}_j{\stackrel{\&OverBar;}{\mathrm{\&delta;}}}_j\end{array}---\left(8\right)$

式中，N为总样本数目。

5.5)设置迭代次数I进行迭代，本发明中I＝8，以此为例，不限于此，具体迭代过程为：

i)若c≤K/2，或者2K≤c≤K/2且迭代次数为奇数，则进行分裂：

对每个类求取各维的标准偏差σ_j＝[σ_j1,σ_j2,…,σ_jd]^T：

$\begin{array}{ccc}{\mathrm{\&sigma;}}_{ji}=\sqrt{\frac{1}{N_j}\underset{{\mathrm{DN}}_k\&Element;{\mathrm{\&Gamma;}}_j}{\&Sigma;}{({\mathrm{DN}}_{ki}-m_{ji})}^2},& j=1,\mathrm{...},c,& i=1,\mathrm{...},d\end{array}---\left(9\right)$

式中，DN_ki为第k个样本的第i个分量，m_ji是当前第j个聚类中心的第i个分量，σ_ji是第j类第i个分量的标准差，d是样本维数；

对每个类求出标准偏差最大的分量σ_jmax，j＝1,…,c；

对各类的σ_jmax(j＝1,…,c)，若存在某个类的σ_jmax>标准偏差参数θ_s，且且N_j>2(θ_N+1)，或c≤K/2，则将Γ_j分裂为两类，中心分别为和置c＝c+1，分裂项γ_j＝[0,…,0,σ_jmax,0,…,0]^T。

ii)若c≥2K，或者2K≤c≤K/2且迭代次数为偶数，则进行合并：

计算各聚类中心两两之间距离δ_ij：δ_ij＝||m_i-m_j||，i,j＝1,…,c，i≠j

比较δ_ij与合并参数θ_C，l为小于θ_C的个数，对小于θ_C者排序：

${\mathrm{\&delta;}}_{i_1{j}_1}<{\mathrm{\&delta;}}_{i_2{j}_2}<...<{\mathrm{\&delta;}}_{i_l{j}_l}$

从最小的开始，把每个对应的和合并，组成新类，新类的中心为：

$m_l=\frac{1}{N_{i_l}+N_{j_l}}\&lsqb;N_{i_l}{m}_{i_l}+N_{j_l}{m}_{j_l}\&rsqb;---\left(10\right)$

并令c＝c-1，每次迭代中避免同一类被合并两次。

5.6)若达到设置的最大迭代次数，则终止迭代；否则，迭代次数加1，进行步骤5.2)。

实际使用中，可以根据先验知识对合成影像的特征种类预先进行判断，当提取特征类别已知的情况下，对合成影像还可以采用基于机器学习方法进行模式识别，得出分类结果，具体为：针对步骤4中产生的合成影像，逐像素遍历影像各像素点DN值，构建特征向量(DN₁，DN₂，……，DN_n)，然后对这些特征向量利用机器学习的方法进行模式识别。本发明的机器学习方法以支持向量机方法为实施例进行说明，但是不限于此，具体过程为：

A)选取特征样本数据集。

从形成的特性向量集合中，勾选出一部分需要提取的已知特征(如道路、主城区、郊区)作为样本数据集特征向量，样本数据集特征向量中记录了每类特征样本各个波段的影像亮度值及样本类别，记为：(dn₁,y₁),(dn₂,y₂),…,(dn_m,y_m)，其中，m为样本个数。

B)采用支持向量机方法对合成影像进行特征提取，具体步骤如下：

B1)构建支持向量机决策函数，即：通过样本数据集构建分类支持向量机目标函数，求取a_i及b：

$\begin{array}{c}\max Q\left(\mathrm{\&alpha;}\right)=\underset{i=1}{\overset{n}{\&Sigma;}}{a}_i-\frac{1}{2}\underset{i,j=1}{\overset{n}{\&Sigma;}}{a}_i{a}_j{y}_j{y}_{i}K({\mathrm{dn}}_i\&CenterDot;{\mathrm{dn}}_j)\\ s.t.{\&Sigma;}_{i,j=1}^n{a}_i{y}_i=0\end{array}---\left(11\right)$

核函数：K(dn_i·dn_j)＝((dn_i·dn_j)+1)²(12)

式中，a_i，a_j为拉格朗日乘子，dn_i，dn_j为样本数据集特征向量，y_i，y_j＝±1为数据点类别。本发明利用二次规划方法求解目标函数的最优解，得到最优拉格朗日乘子a_i，利用样本数据集中的支持向量dn带入判别函数，求解偏差值b，其中，判别函数表达式为：

$f\left(x\right)=\mathrm{sgn}\left(\underset{i=1}{\overset{n}{\&Sigma;}}{a}_i{y}_{i}K\right({\mathrm{dn}}_i\&CenterDot;{\mathrm{dn}}_j)+b)---\left(13\right)$

式中，sgn为判别函数，若内部返回值大于0，则返回1；等于0，则返回0，小于0，则返回-1。

根据上述方法分别构建n个支持向量机决策函数，其中，n为所划分的特征类别个数，即每个类别构建一个支持向量机决策函数。

B2)采用步骤B1)构建出的支持向量机决策函数对特征向量进行分类。

将特征向量(DN₁，DN₂，……，DN_n)带入第1类的支持向量机决策函数，进行分类，若f(x)为1，则该特征向量属于第1类，若f(x)为-1，则该特征向量不属于第1类，将不属于该类的特征向量带入第2类的支持向量机决策函数进行分类，以此类推，最终得出划分为n类特征的合成影像。

6、获得分类完成的合成影像，根据分类结果，可以高精度的提取出感兴趣区域，同时也可以进行区域在时间序列上的扩展过程研究。另外本发明分类成果能区分出单幅夜光卫星影像上无法区分的信息，如商业区与居民区的区分，城市中心与边缘的区分等。这些结果更适合在区域大尺度上的精细研究，即通过进行多年度时间序列合成影像，应用于大尺度城镇空间扩展研究、道路、商业区、居民区提取分析等与人类活动相关的问题的研究。

上述各实施例仅用于说明本发明，其中方法的各实施步骤都是可以有所变化的，凡是在本发明技术方案的基础上进行的等同变换和改进，均不应排除在本发明的保护范围之外。

Claims (7)

1.一种基于时间序列夜光遥感影像的数据挖掘方法，其特征在于包括以下内容：

1)读取时间序列夜光遥感影像，并从读取的时间序列夜光遥感影像中选择某几段时序的时间序列夜光遥感影像；

2)对选取的某几段时间序列夜光遥感影像分别进行重投影和重采样，并获得重采样影像；

3)将重采样影像进行合成获得具有多个波段的合成影像；

4)对合成影像进行数据挖掘得到影像分类结果。

2.如权利要求1所述的一种基于时间序列夜光遥感影像的数据挖掘方法，其特征在于，所述步骤2)对选取的某几段时间序列夜光遥感影像分别进行重投影和重采样，并获得重采样影像，具体过程为：

2.1)根据选取的每一段时间序列夜光遥感影像的地理信息范围，确定重投影后影像范围，并获取影像覆盖区域；

2.2)根据影像覆盖区域左上角点坐标及分辨率确定重投影后的影像大小；

2.3)对影像进行重采样，即根据步骤2.1)及2.2)得到的重投影影像大小、分辨率和左上角点坐标，对选取的时间序列夜光遥感影像进行重采样，为构建的重投影影像各个像素点赋值。

3.如权利要求2所述的一种基于时间序列夜光遥感影像的数据挖掘方法，其特征在于，所述步骤2.3)对影像进行重采样的具体过程为：

2.31)从重投影后影像的左上角像素点坐标开始，逐像素遍历步骤2.2)中构建的重投影影像，将重投影影像上的左上角点坐标计算到重投影前时间序列夜光遥感影像上；

2.32)在重投影前的时间序列夜光遥感影像上进行重采样，得到重投影影像中对应像元的DN值；

2.33)对覆盖区域内的所有像素点重复步骤2.31)及2.32)，得到重采样影像。

4.如权利要求1或2或3所述的一种基于时间序列夜光遥感影像的数据挖掘方法，其特征在于，所述步骤4)对合成影像采用非监督分类方法进行数据挖掘。

5.如权利要求4所述的一种基于时间序列夜光遥感影像的数据挖掘方法，其特征在于，所述非监督分类方法为迭代自组织数据分析方法，对合成影像进行数据挖掘采用迭代自组织数据分析方法，得到影像分类结果的具体情况为：

4.1)初始化，假设初始聚类数c以及期望得到的聚类数为K，确定初始中心m_i，i＝1,…,c：

$m_i=\frac{1}{N_i}\underset{{\mathrm{DN}}_k\&Element;{\mathrm{\&Gamma;}}_i}{\&Sigma;}{\mathrm{DN}}_k$

式中，N_i为第i类聚类Γ_i中的样本数目，DN_k为聚类Γ_i中的样本；

4.2)将所有的样本分到距初始中心最近的类Γ_i中，若某个类Γ_j中的样本数量N_j小于一个聚类中的最小样本数θ_N，则去掉该类，令c＝c-1；

4.3)重新计算聚类中心m_j为：

$m_j=\frac{1}{N_j}\underset{{\mathrm{DN}}_k\&Element;{\mathrm{\&Gamma;}}_j}{\&Sigma;}{\mathrm{DN}}_k,j=1,...,c$

式中，N_j为第j个聚类中的样本数目；

4.4)计算第j类样本与其中心的平均距离和总平均距离

${\stackrel{\&OverBar;}{\mathrm{\&delta;}}}_j=\frac{1}{N_j}\underset{{\mathrm{DN}}_k\&Element;{\mathrm{\&Gamma;}}_j}{\&Sigma;}\left|\right|{\mathrm{DN}}_k-m_j\left|\right|,j=1,...,c$

$\stackrel{\&OverBar;}{\mathrm{\&delta;}}=\frac{1}{N}\underset{j=1}{\overset{c}{\&Sigma;}}{N}_j{\stackrel{\&OverBar;}{\mathrm{\&delta;}}}_j$

式中，N为总样本数目；

4.5)设置迭代次数I并进行迭代，具体迭代过程为：

i)若c≤K/2，或者2K≤c≤K/2且迭代次数为奇数，则进行分裂：

对每个类求取各维的标准偏差σ_j＝[σ_j1,σ_j2,…,σ_jd]^T：

${\mathrm{\&sigma;}}_{ji}=\sqrt{\frac{1}{N_j}\underset{{\mathrm{DN}}_k\&Element;{\mathrm{\&Gamma;}}_j}{\&Sigma;}{({\mathrm{DN}}_{ki}-m_{ji})}^2},j=1,...,c,i=1,...,d$

式中，DN_ki为第k个样本的第i个分量，m_ji是当前第j个聚类中心的第i个分量，σ_ji是第j类第i个分量的标准差，d是样本维数；

对每个类求出标准偏差最大的分量σ_jmax，j＝1,…,c；

对各类的σ_jmax(j＝1,…,c)，若存在某个类的σ_jmax>标准偏差参数θ_s，且且N_j>2(θ_N+1)，或c≤K/2，则将Γ_j分裂为两类，中心分别为和置c＝c+1，分裂项γ_j＝[0,…,0,σ_jmax,0,…,0]^T；

ii)若c≥2K，或者2K≤c≤K/2且迭代次数为偶数，则进行合并：

计算各聚类中心两两之间距离δ_ij：δ_ij＝||m_i-m_j||，i,j＝1,…,c，i≠j

比较δ_ij与合并参数θ_C，l为小于θ_C的个数，对小于θ_C者排序：

${\mathrm{\&delta;}}_{i_1{j}_1}<{\mathrm{\&delta;}}_{i_2{j}_2}<...<{\mathrm{\&delta;}}_{i_l{j}_l}$

从最小的开始，把每个对应的和合并，组成新类，新类的中心为：

$m_l=\frac{1}{N_{i_l}+N_{j_l}}\&lsqb;N_{i_l}{m}_{i_l}+N_{j_l}{m}_{j_l}\&rsqb;$

并令c＝c-1，每次迭代中避免同一类被合并两次；

4.6)若达到设置的最大迭代次数，则终止迭代；否则，迭代次数加1，进入步骤4.2)。

6.如权利要求1或2或3所述的一种基于时间序列夜光遥感影像的数据挖掘方法，其特征在于，当提取特征类别已知的情况下，对合成影像采用基于机器学习方法进行数据挖掘，得出影像分类结果。

7.如权利要求6所述的一种基于时间序列夜光遥感影像的数据挖掘方法，其特征在于，所述基于机器学习方法为支持向量机方法，对合成影像进行数据挖掘采用支持向量机方法，得到影像分类结果的具体过程为：

4.1)选取特征样本数据集，样本数据集中的特征向量中记录了每类特征样本各个波段的影像亮度值及样本类别，记为：(dn₁,y₁),(dn₂,y₂),…,(dn_m,y_m)，m为样本个数；

4.2)对合成影像采用支持向量机方法提取影像特征，并对合成影像进行分类，具体步骤如下：

4.21)构建支持向量机预测函数，即：通过样本数据集构建分类支持向量机目标函数，求取a_i及b：

$\max Q\left(\mathrm{\&alpha;}\right)=\underset{i=1}{\overset{n}{\&Sigma;}}{a}_i-\frac{1}{2}\underset{i,j=1}{\overset{n}{\&Sigma;}}{a}_i{a}_j{y}_j{y}_{i}K({\mathrm{dn}}_i\&CenterDot;{\mathrm{dn}}_j)$

$s.t.{\mathrm{\&Sigma;}}_{i,j=1}^n{a}_i{y}_i=0$

核函数：K(dn_i·dn_j)＝((dn_i·dn_j)+1)²

式中，a_i，a_j为拉格朗日乘子，dn_i，dn_j为样本数据集特征向量，y_i，y_j＝±1为数据点类别，利用二次规划方法求解目标函数的最优解，得到最优拉格朗日乘子a_i，利用样本数据集中的支持向量dn，带入判别函数，求解偏差值b，其中，判别函数表达式为：

$f\left(x\right)=\mathrm{sgn}\left(\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{a}_i{y}_{i}K\right({\mathrm{dn}}_i\&CenterDot;{\mathrm{dn}}_j)+b)$

式中，sgn为判别函数，若内部返回值大于0，则返回1；等于0，则返回0，小于0，则返回-1；

根据上述方法分别构建n个支持向量机决策函数，其中，n为所划分的特征类别个数，即每个类别构建一个支持向量机决策函数；

4.22)采用步骤4.21)构建出的支持向量机决策函数对特征向量进行分类：

将特征向量带入第1类的支持向量机决策函数，进行分类，若f(x)为1，则该特征向量属于第1类，若f(x)为-1，则该特征向量不属于第1类，将不属于该类的特征向量带入第2类的支持向量机决策函数进行分类，以此类推，最终得出划分为n类特征的合成影像。