一种融合地表温度和建筑指数的城市实体边界识别方法与流程

本发明属于国土空间规划和城乡规划领域，具体涉及一种融合地表温度和建筑指数的城市建成区实体边界识别方法。

背景技术：

伴随全球城市扩张难以避免、生态损毁困境不断加剧，我国城市建成区扩张呈现更加凸显的失衡三线困境(扩张弹性突破“合理界线”、城市环境逼近“生态底线”和发展态势逾越“指标上线”)，因此，如何科学识别城市扩张实体边界、调控凸显困境，成为我国城市可持续发展和社会管理的重大理论命题和技术难题。

目前，城市建成区实体边界提取方法通常采用多光谱遥感数据、夜间灯光数据和地价均衡模型等技术方法进行城市边界识别与界定，但以上方法存在较多不足之处：如多光谱遥感的光谱特征易将非城市不透水面误提取出来；夜间灯光遥感的分辨率较低，造成城市边界识别结果与实际存在一定的误差；地价均衡模型对地价空间异质性考虑不足，并且大都未能将城市实体边界扩张的生态环境效应纳入决策分析框架等。

技术实现要素：

发明目的：本发明公开的一种融合地表温度和建筑指数的城市实体边界识别方法，该方法通过结合遥感反演的地表温度和建筑指数可以更为精确地定位城市建成区实体边界，弥补了城市边缘工业区和城乡结合部识别难的缺陷，更为科学的识别城市建成区实体边界。

技术方案：一种融合地表温度和建筑指数的城市建成区边界识别方法，包括以下几个步骤：

步骤1：影像预处理

主要包括辐射定标、大气校正、图幅镶嵌、影像裁剪等，其中大气校正采用flaash大气校正。

步骤2：地表温度反演，如图2所示：

构建地表温度lst反演模型，利用单窗算法计算地表温度

lst＝{a(1-c-d)+[b(1-c-d)+c+d]tsenor-dta}/c

lst为地表温度；a、b为常数，a＝-67.3554，b＝0.45861；tsenor为亮温，ta为大气平均作用温度c＝ετ；，d＝(1-τ)[1+(1-ε)τ]，ε为地表比辐射率，τ为大气透射率。

步骤2.1：计算亮温tsenor

l(λ)＝gain·dn+bias

式中：gain、bias、k1和k2为卫星发射前预设常数，可在头条文件中获取。dn为遥感热红外影像的灰度值，lλ为遥感器所接收的辐射强度。

步骤2.2：计算地表比辐射率ε

地表比辐射率可以通过ndvi值估算求取，

当ndvi＞0.7时，为完全植被，地表比辐射较高，ε＝0.99

当ndvi＜0.05时，为完全裸土，地表比辐射率低，ε＝0.973

当0.05≤ndvi≤0.7时，为混合像元，ε＝0.004pv+0.986

其中植被指数ndvi能够很好的检测植被的生长状况和指数覆盖度，计算公式为：

pv＝((ndvi-ndvis))/((ndviv-ndvis)

ndvi＝((nir-red))/((nir-red))

式中：pv为植被覆盖度，取ndviv＝0.7，ndvis＝0.05，ndvi为归一化植被指数，red为红色波段，nir为近红外波段。

步骤2.3：大气透射率τ

当0.4≤w≤1.6时，τ＝0.974290-0.08007w

当1.6≤w≤3.0时，τ＝1.031412-0.11536w

式中，w为水分含量。

步骤2.4：大气平均作用温度ta

ta＝16.0110+0.92621t0

式中，t0是近地面气温(k)。

步骤3：归一化建筑指数与地表温度融合

首先利用波段计算器计算建筑指数ndbi，再与反演的地表温度数据结合，ndbi大于零且地表温度大于研究区内平均值的认为是建成区。

步骤3.1：对建筑指数ndbi进行阈值分割，得到二值分割后的二值化建筑指数i1；

ndbi＝(mir-nir)/(mir+nir)

其中：nir、mir分别代表近红外和中红外波段的像元值。

步骤3.2：对地表温度反演数据i2进行阈值分割，得到二值分割后的二值化地表温度指数i3；

其中，(i，j)表示图像中的像素坐标，th为阈值分割法选取的分割阈值；

其中，datamax为地表温度最大值；datamin为地表温度最小值；

步骤3.3：结合二值化建筑指数i1与二值化地表温度指数i3，获取综合建成区指数i4；

i4(i，j)＝i3(i，j)×i1(i，j)

其中，(i，j)表示图像中的像素坐标。

步骤4：城市建成区实体边界划定

步骤4.1：点转换。将i4值为1的栅格数据转换为点数据。

步骤4.2：建立不规则三角网。在tin模块中选择创建tin，然后选择描绘tin数据区，在这里面设置的最大边长为120m，此阈值选择主要是基于温度反演数据的分辨率为120m；

步骤4.3：确定tin范围。将创建的tin数据区生成tindomain文件；

步骤4.4：要素转换为筛选。将生成的tindomain文件转为面要素，筛选出城市建成区及其实体边界。

有益效果：本发明公开的一种融合地表温度和建筑指数的城市边界识别方法，通数据源只需要landsat数据，操作简单，正确率高，能够快速准确的得到城市建成区边界。大量实验结果表明，本发明边界识别正确率达到了94％以上，适用于国土空间规划、地理国情监测、城市蔓延治理分析等领域。

附图说明

图1本发明研究思路示意图；

图2本发明温度反演图；

图3本发明城市建成区实体边界识别结果。

具体实施方式

如图1所示，一种融合地表温度和建筑指数的城市建成区边界识别方法，包括以下几个步骤：

步骤1：影像预处理

主要包括辐射定标、大气校正、图幅镶嵌、影像裁剪等，其中大气校正采用flaash大气校正。

步骤2：地表温度反演

构建地表温度lst反演模型，利用单窗算法计算地表温度

lst＝{a(1-c-d)+[b(1-c-d)+c+d]tsenor-dta}/c

lst为地表温度；a、b为常数，a＝-67.3554，b＝0.45861；tsenor为亮温，ta为大气平均作用温度c＝ετ；，d＝(1-τ)[1+(1-ε)τ]，ε为地表比辐射率，τ为大气透射率。

步骤2.1：计算亮温tsenor

l(λ)＝gain·dn+bias

式中：gain、bias、k1和k2为卫星发射前预设常数，可在头条文件中获取。dn为遥感热红外影像的灰度值，lλ为遥感器所接收的辐射强度。

步骤2.2：计算地表比辐射率ε

地表比辐射率可以通过ndvi值估算求取，

当ndvi＞0.7时，为完全植被，地表比辐射较高，ε＝0.99

当ndvi＜0.05时，为完全裸土，地表比辐射率低，ε＝0.973

当0.05≤ndvi≤0.7时，为混合像元，ε＝0.004pv+0.986

其中植被指数ndvi能够很好的检测植被的生长状况和指数覆盖度，计算公式为：

pv＝((ndvi-ndvis))/((ndviv-ndvis)

ndvi＝((nir-red))/((nir-red))

式中：pv为植被覆盖度，取ndviv＝0.7，ndvis＝0.05，ndvi为归一化植被指数，red为红色波段，nir为近红外波段。

步骤2.3：大气透射率τ

当0.4≤w≤1.6时，τ＝0.974290-0.08007w

当1.6≤w≤3.0时，τ＝1.031412-0.11536w

式中，w为水分含量。

步骤2.4：大气平均作用温度ta

ta＝16.0110+0.92621t0

式中，t0是近地面气温(k)。

步骤3：归一化建筑指数与地表温度融合

首先利用波段计算器计算建筑指数ndbi，再与反演的地表温度数据结合，ndbi大于零且地表温度大于研究区内平均值的认为是建成区。

步骤3.1：对建筑指数ndbi进行阈值分割，得到二值分割后的二值化建筑指数i1；

ndbi＝(mir-nir)/(mir+nir)

其中：nir、mir分别代表近红外和中红外波段的像元值。

步骤3.2：对地表温度反演数据i2进行阈值分割，得到二值分割后的二值化地表温度指数i3；

其中，(i，j)表示图像中的像素坐标，th为阈值分割法选取的分割阈值；

其中，datamax为地表温度最大值；datamin为地表温度最小值；

步骤3.3：结合二值化建筑指数i1与二值化地表温度指数i3，获取综合建成区指数i4；

i4(i，j)＝i3(i，j)×i1(i，j)

其中，(i，j)表示图像中的像素坐标。

步骤4：城市建成区实体边界划定

步骤4.1：点转换。将i4值为1的栅格数据转换为点数据。

步骤4.2：建立不规则三角网。在tin模块中选择创建tin，然后选择描绘tin数据区，在这里面设置的最大边长为120m，此阈值选择主要是基于温度反演数据的分辨率为120m；

步骤4.3：确定tin范围。将创建的tin数据区生成tindomain文件；

步骤4.4：要素转换为筛选。将生成的tindomain文件转为面要素，筛选出城市建成区及其实体边界。如图3所示。