本发明涉及城市不透水层的检测，更具体地说，涉及一种城市不透水层覆盖率遥感估算方法。
背景技术：
城市不透水层(urbanimpervioussurfaces，uis)是城市中的重要人造地表，它能够阻止地下水渗透到土壤中，割断了城市地表与地下水文的联系，主要由城市中的建筑物屋顶、广场、停车场和柏油道路等组成。城市不透水层不仅是衡量城市化程度的指示剂，也是衡量城市环境质量的重要指标。随着全球城市化进程的不断加快，城市不透水层所占的面积比例逐渐扩大，引发了城市热岛效应、环境污染、资源短缺等一系列生态、环境和社会问题。因此，提取城市不透水层，进而监测城市不透水层的变化对监测城市热环境变化和城市生态系统具有重要的现实意义，也是城市规划中要考虑的重要因素。
遥感技术由于具有快速获取地表信息的优势而成为进行uis提取和监测的重要手段。用于uis提取和监测的遥感数据源包括高分辨率遥感影像、中分辨率遥感影像、高光谱遥感影像、多源遥感影像的结合等。中分辨率遥感影像因其空间分辨率适中、覆盖范围大、成本低等优势而成为uis提取和监测的常用遥感数据源。由于城市景观复杂，地物具有强烈的空间异质性使得中分辨率遥感影像中存在了大量的混合像元。因此通常采用估算城市不透水层覆盖率(urbanimpervioussurfacespercent，uisp)—“单位遥感影像像元中不透水层所占的面积比例”作为衡量城市不透水层覆盖程度的指标。目前用于中分辨率遥感影像uisp估算的方法主要有混合像元分解、多元回归和机器学习等方法。混合像元分解方法假设地表各地物之间没有交叉反射，这不符合地表非朗伯性的现实情况，且像元分解的精度主要取决于端元选取的精度，因此该方法受人为因素影响大，稳定性差。普通的多元回归不能有效地描述遥感影像像元波谱特征与uisp之间复杂的非线性关系。目前人工神经网络是进行uisp估算的常用机器学习方法，但是该方法依赖训练样本太多、容易陷入局部最优、结果的泛化能力差。因此迫切需要应用新的方法和技术提高利用中分辨率遥感影像为主要数据源进行uisp遥感估算的精度。
支持向量机(svm)是建立在统计学习理论的vc维理论和结构风险最小原理基础上的一种机器学习方法，有严格的理论基础，其算法是一个凸二次优化问题，保证找到的解是全局最优解，能较好地解决小样本、非线性、高维数、局部极小点等实际问题，问题的复杂度不取决于特征的维数，且具有良好的推广能力,是复杂非线性科学和人工智能科学的研究前沿。
考虑到成本因素，本研究可用的高分辨率遥感影像数量较少，且uisp与影像像元的波谱信息及其衍生信息之间存在复杂的非线性关系，因此适用于利用支持向量回归机(supportvectorregression，svr)估算uisp。支持向量回归(svr)是建立在svm思想上的回归算法，具有坚实的理论基础，但是其理论优势得以实现的前提是要选取到合适的回归参数。因此迫切需要开发新的支持向量机应用于特定问题的参数优化算法。
支持向量回归机参数的优化对于其拟合精度的影响至关重要，目前常用的svr参数优化方法主要有实验法、留一法(loo)、网格搜索法(gridsearch)[14]、遗传算法(ga)[15]和粒子群算法(pso)[16]等。目前，ga和pso算法应用较多。但是ga算法存在实现比较复杂，对初始值依赖大、并行性有限等问题。pso算法虽然具有收敛快速、规则简单等特点，但是它的早熟收敛和易限于局部最优等缺点还是不容忽视的。鉴于此，本文将混沌算法融入到粒子群算法中，利用混沌算法的全局遍历特性，克服了粒子群算法的早熟收敛等缺点，使得混沌粒子群算法(cpso)兼具两者的优点。
因此本发明拟在：对支持向量回归机中的参数优化算法进行改进，提出一种应用于中分辨率遥感影像的改进型支持向量回归uisp遥感估算方法以提高支持向量回归在中分辨率遥感影像uisp估算中的精度。
技术实现要素：
本发明的方法支持向量回归机中的参数优化算法进行改进，提出一种应用于中分辨率遥感影像的改进型支持向量回归uisp遥感估算方法以提高支持向量回归在中分辨率遥感影像uisp估算中的精度。
为了达到上述目的，本发明提供一种城市不透水层覆盖率遥感估算方法，其特征在于，所述估算方法步骤如下：
s1、基于中分辨率遥感影像的多个波谱波段计算归一化植被指数(ndvi)、归一化建筑物指数(ndbi)、地表温度(lst)；
具体公式如下：
其中，nir是近红外波段的波谱值，r是红波段的波谱值，mir是中红外波段的波谱值；
s2、基于中分辨率遥感影像，利用灰度共生矩阵计算每个波段的纹理信息；
s3、选取与中分辨率遥感影像同时相，地理覆盖范围为中分辨率遥感影像覆盖范围之内且分散分布在东西南北四个方向的多块高分辨率遥感影像，利用面向对象的图像处理方法获取高分辨率遥感影像中每个像元的不透水性，其中，1表示不透水，0表示透水；
s4、统计每个中分辨率遥感影像像元对应的多个高分辨遥感影像像元中数值为1的像元数占总像元数的百分比，即为对应的中分辨率遥感影像像元的uisp；具体公式如公式(3)：
其中uispij表示中分辨率遥感影像上第i行第j列像元的uisp，uispij∈[0,1]，αpq表示高分辨率遥感影像上第p行第q列像元的不透水性，取值为0(透水)或1(不透水)；k表示每个中分辨率遥感影像像元中有k*k个高分辨率遥感影像像元；所述k取值如下：
其中，m为中分辨率遥感影像像元的大小，n为高分辨遥感影像像元的大小，如果m/n不等于整数，则需要把中分辨率遥感影像重采样为像元大小为r(r大于m，且是与m最靠近的整数)的图像，使得r/n等于整数；
s5、随机选取高分辨率遥感影像和与之空间对应的中分辨率遥感影像像元中的70％用作训练数据，30％用作验证数据；
s6、对s1和步s2中所得的数据进行归一化处理，具体公式如公式(5)：
其中xi为归一化前像元的值，xmin为整幅图像中像元的最小值，xmax为整幅图像中像元的最大值，xj为归一化后像元的值，介于0和1之间。
s7、对s2中经过s6处理后的纹理特征数据利用信息熵和相关系数进行初步约简，删除信息熵小和相关系数高的特征；
s8、采用改进的h离散化方法对步骤1中经过s6处理后的数据进行离散化处理；所述改进的h离散化方法为：
根据数据的统计规律分类数据，对于在取值空间分布较为均匀的数据，采用梯度函数对其进行离散化；对于空间分布不均匀，聚集效应很明显的属性数据采用指数函数对其进行离散化；
将离散化后的属性与uisp数据按顺序排列好，组成利用粗糙集进行遥感估算uisp的信息表，即为uisp估算的论域；
s9、采用基于差别矩阵的粗糙集约简算法对s8中经过离散化后的特征进行特征约简；
通过前面对样本数据的预处理，可以用样本数据组成一个信息系统，即，uisp＝(u,c,v,f)；
其中，u是信息系统的论域，是所有样本组成的集合。
c是信息系统中属性的集合，是从中分辨率遥感影像中提取出来的属性的集合，单个属性称为简单属性；
v是信息函数f的值域，vi为属性αi的值域，即各个属性的取值范围；
f是系统的信息函数，fi为属性αi的信息函数；
具体步骤如下：
step1：先求出其差别矩阵mn×n；
setp2：对差别矩阵中的所有取值不为空集合的元素cij(即：
step3：对step2得到的析取范式进行逻辑运算，得到内合取外析取的合取范式lv(∧)(is)如公式8所示：
step4：输出
s10、基于步骤9中约简后的数据建立用于uisp估算的支持向量回归机模型(svr-uisp)，svr-uisp模型的输入数据为中分辨率遥感影像像元的波谱信息及其衍生信息，输出数据为每个中分辨率像元的uisp。通过建立中分辨率影像像元的波谱信息及其衍生信息和对应uisp参考值之间的回归模型，从而可以根据训练数据的特征进而外推到中分辨率影像其它像元的uisp。
输入和输出数据的表达如下：
s11、基于s10，选用ε-支持向量回归机(ε-svr)和径向基核函数(brf)，采用混沌粒子群算法(cpso)优化不敏感系数ε、惩罚系数c和rbf核的宽度系数σ三个参数，建立cpso-ε-svr-uisp模型。cpso进行ε-svr-uisp参数优化的主要思想是将混沌算法嵌套进粒子群算法，利用混沌运动的随机性和遍历性以当前整个粒子群寻找到的最优位置为基点生成混沌序列，用生成的混沌序列中的适应度最高的粒子随机替换当前种群中的一个粒子。
svr拟合样本数据集{xi,yi}(i＝1,2,…,n；xi∈rd；yi∈r)得到的回归函数如公式9所示：
公式中，αi、
其中
s12、利用s10优化后的cpso-ε-svr-uisp模型，输入中分辨率遥感影像上的基本波谱信息及其衍生信息经过约简后的数据，模型的输出即为对应每个中分辨率像元中uisp的值，由此完成已知中分辨率遥感影像像元的波谱及其衍生信息估算uisp的任务。
所述地表温度(lst)的计算要求影像必须有热红外波段，采用单窗算法进行地表温度的计算；
所述纹理信息分别为：分别为中值(mean)、反差(contrast，con)、自相关(correlation,or)、dissimilarity(差异性)、entropy(熵，ent)、homogeneitymean(同质性)、angularsecondmoment(角二阶矩，asm)、variance(协方差)8种。
混沌粒子群优化支持向量回归机的具体步骤如下：
step1：初始化粒子群(ε，c，σ)，给定种群规模m，确定算法的最大、最小权重因子ωmax,ωmin值，预设该方法的最大迭代次数itermax；需要注意的是：由于是ε，c，σ三个参数同时优化，但它们的值不在相同量级上，因此在随机初始化粒子的速度v时，要为其添加调整系数；
step2：设置每个粒子的个体最优值pibest为当前最优值，使用拟合函数求出估计值，再利用适应度函数计算每个粒子的适应度，选择种群中适应度最好的个体所对应的最优值为初始全局极值gbest。
step3：按照公式(11)、(12)进行迭代计算，更新每个粒子的速度v和位置x；
vi+1＝w·vi+c1r1(pibest-vi)+c2r2(gbest-vi)(11)
xi+1＝xi+βvi+1(12)
step4：通过公式(9)、(10)求出粒子的个体适应度。
step5：将每个粒子的当前适应度值与它的pibest对应的适应度值比较，若当前适应度更好，令pibest等于当前适应度值，反之，保留原值。
step6：比较更新后的单个粒子的pibest与种群全局最优值gbest对应的适应度值，更新gbest。
step7：对gbest对应的最优位置xg＝(xg1,xg2,…,xgd)t进行混沌优化，通过公式
yq+1＝μyq(1-yq)(13)
其中yq为(0,1)的随机数，μ一般取4，此时logistic方程处于完全混沌的状态。然后通过下式：
进行迭代产生混沌变量序列
xqi＝ai+yqi(bi-ai)(14)
还原到原来的解空间，得，
在原来的解空间中计算混沌变量经历的每个可行解
step8：用p*替换从当前种群中随机选取的一个粒子。
step9：判断算法是否满足终止条件，如果迭代达到设定的最大次数或者得到的解停止变化，中止迭代；如果不满足终止条件，返回到step3。
本发明优点：提出了一种新的以性价比高的中分辨率遥感影像为主要数据源利用混沌粒子群优化支持向量回归估算城市不透水面覆盖率的新方法，对支持向量回归机中的参数优化算法进行改进，提出一种应用于中分辨率遥感影像的改进型支持向量回归uisp遥感估算方法以提高支持向量回归在中分辨率遥感影像uisp估算中的精度。进一步提高了以中分辨率遥感影像为主要数据源进行uisp估算的性能，为利用城市不透水面信息作为输入数据的相关研究提供更为精确的城市不透水面数据。
具体实施方式
针对以上技术问题，本发明提供了一种城市不透水层覆盖率遥感估算方法(uisp)估算方法，包括步骤如下：
s1、基于中分辨率遥感影像的多个波谱波段计算归一化植被指数(ndvi)、归一化建筑物指数(ndbi)、地表温度(lst)。
具体公式如下：
地表温度的计算要求影像必须有热红外波段，采用单窗算法进行地表温度的计算。具体算法参见以下文章。“谢远礼等.基于tm影像的兰州市地表温度反演及城市热岛效应研究”
其中，nir是近红外波段的波谱值，r是红波段的波谱值，mir是中红外波段的波谱值。
s2、基于中分辨率遥感影像，利用灰度共生矩阵计算每个波段的纹理信息。分别为中值(mean)、反差(contrast,con)、自相关(correlation,or)、dissimilarity(差异性)、entropy(熵,ent)、homogeneitymean(同质性)、angularsecondmoment(角二阶矩，asm)、variance(协方差)8种纹理信息。具体计算方法和操作步骤参见envi遥感图像处理软件的帮助。(需考虑不同滤波大小得到的纹理进行比较)
s3、选取与中分辨率遥感影像同时相，地理覆盖范围为中分辨率遥感影像覆盖范围之内且分散分布在东西南北四个方向的多块高分辨率遥感影像，利用面向对象的图像处理方法(参考文献为：)获取高分辨率遥感影像中每个像元的不透水性(1表示不透水，0表示透水)。
s4、统计每个中分辨率遥感影像像元对应的多个高分辨遥感影像像元中数值为1的像元数占总像元数的百分比，即为对应的中分辨率遥感影像像元的uisp。
具体公式如公式(3)：
其中uispij表示中分辨率遥感影像上第i行第j列像元的uisp，uispij∈[0,1]，αpq表示高分辨率遥感影像上第p行第q列像元的不透水性，取值为0(透水)或1(不透水)。k表示每个中分辨率遥感影像像元中有k*k个高分辨率遥感影像像元。k的取值如公式(5)。
其中，m为中分辨率遥感影像像元的大小，n为高分辨遥感影像像元的大小，如果m/n不等于整数，则需要把中分辨率遥感影像重采样为像元大小为r(r大于m，且是与m最靠近的整数)的图像，使得r/n等于整数。
s5、随机选取高分辨率遥感影像和与之空间对应的中分辨率遥感影像像元中的70％用作训练数据，30％用作验证数据。
s6、对s1和s2中所得的数据进行归一化处理，具体公式如公式(5)：
其中xi为归一化前像元的值，xmin为整幅图像中像元的最小值，xmax为整幅图像中像元的最大值，xj为归一化后像元的值，介于0和1之间。
s7、对s2中经过s6处理后的纹理特征数据利用信息熵和相关系数进行初步约简,删除信息熵小和相关系数高的特征。
s8、采用改进的h离散化方法对s1中经过s6处理后的数据进行离散化处理。
1)根据数据的统计规律分类数据，对于在取值空间分布较为均匀的数据，采用梯度函数对其进行离散化；对于空间分布不均匀，聚集效应很明显的属性数据采用指数函数对其进行离散化。
2)将离散化后的属性与uisp数据按顺序排列好，组成利用粗糙集进行遥感估算uisp的信息表，即为uisp估算的论域
s9、采用基于差别矩阵的粗糙集约简算法对步骤8中经过离散化后的特征进行特征约简。
通过前面对样本数据的预处理，可以用样本数据组成一个信息系统，即，
uisp＝(u,c,v,f)(6)
其中，u是信息系统的论域，是所有样本组成的集合。
c是信息系统中属性的集合，是从中分辨率遥感影像中提取出来的属性的集合，单个属性称为简单属性。
v是信息函数f的值域，vi为属性αi的值域，即各个属性的取值范围。
f是系统的信息函数，fi为属性αi的信息函数。
具体步骤如下：
step1：先求出其差别矩阵mn×n；
setp2：对差别矩阵中的所有取值不为空集合的元素cij(即：
step3：对step2得到的析取范式进行逻辑运算，得到内合取外析取的合取范式lv(∧)(is)如公式8所示：
step4：输出
s10、基于s9中约简后的数据建立用于uisp估算的支持向量回归机模型(svr-uisp)，svr-uisp模型的输入数据为中分辨率遥感影像像元的波谱信息及其衍生信息，输出数据为每个中分辨率像元的uisp。通过建立中分辨率影像像元的波谱信息及其衍生信息和对应uisp参考值之间的回归模型，从而可以根据训练数据的特征进而外推到中分辨率影像其它像元的uisp。
输入和输出数据的表达如下：
s11、基于s10，选用ε-支持向量回归机(ε-svr)和径向基核函数(brf)，采用混沌粒子群算法(cpso)优化不敏感系数ε、惩罚系数c和rbf核的宽度系数σ三个参数，建立cpso-ε-svr-uisp模型。cpso进行ε-svr-uisp参数优化的主要思想是将混沌算法嵌套进粒子群算法，利用混沌运动的随机性和遍历性以当前整个粒子群寻找到的最优位置为基点生成混沌序列，用生成的混沌序列中的适应度最高的粒子随机替换当前种群中的一个粒子。
svr拟合样本数据集{xi,yi}(i＝1,2,…,n；xi∈rd；yi∈r)得到的回归函数如公式9所示：
公式中，αi、
其中
混沌粒子群优化支持向量回归机的具体步骤如下：
step1：初始化粒子群(ε,c,σ)，给定种群规模m，确定算法的最大、最小权重因子ωmax,ωmin值，预设该方法的最大迭代次数itermax。需要注意的是：由于是ε,c,σ三个参数同时优化，但它们的值不在相同量级上，因此在随机初始化粒子的速度v时，要为其添加调整系数。
step2：设置每个粒子的个体最优值pibest为当前最优值，使用拟合函数求出估计值，再利用适应度函数计算每个粒子的适应度，选择种群中适应度最好的个体所对应的最优值为初始全局极值gbest。
step3：按照公式(11)、(12)进行迭代计算，更新每个粒子的速度v和位置x。
vi+1＝w·vi+c1r1(pibest-vi)+c2r2(gbest-vi)(11)
xi+1＝xi+βvi+1(12)
step4：通过公式(9)、(10)求出粒子的个体适应度。
step5：将每个粒子的当前适应度值与它的pibest对应的适应度值比较，若当前适应度更好，令pibest等于当前适应度值，反之，保留原值。
step6：比较更新后的单个粒子的pibest与种群全局最优值gbest对应的适应度值，更新gbest。
step7：对gbest对应的最优位置xg＝(xg1,xg2,…,xgd)t进行混沌优化，通过公式
yq+1＝μyq(1-yq)(13)
其中yq为(0,1)的随机数，μ一般取4，此时logistic方程处于完全混沌的状态。然后通过下式：
进行迭代产生混沌变量序列
xqi＝ai+yqi(bi-ai)(14)
还原到原来的解空间，得
在原来的解空间中计算混沌变量经历的每个可行解
step8：用p*替换从当前种群中随机选取的一个粒子。
step9：判断算法是否满足终止条件，如果迭代达到设定的最大次数或者得到的解停止变化，中止迭代；如果不满足终止条件，返回到step3。
s12、利用s10优化后的cpso-ε-svr-uisp模型，输入中分辨率遥感影像上的基本波谱信息及其衍生信息经过约简后的数据，模型的输出即为对应每个中分辨率像元中uisp的值，由此完成已知中分辨率遥感影像像元的波谱及其衍生信息估算uisp的任务。
以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明披露的技术范围内，根据本发明的技术方案及其发明构思加以等同替换或改变，都应涵盖在本发明的保护范围之内。