CN116151610B - 一种非均质城市下垫面承灾体风险暴露空间模拟方法 - Google Patents

Abstract

本发明公开了一种非均质城市下垫面承灾体风险暴露空间模拟方法，首先采集城市的多源异构数据并进行预处理，得到预处理数据；再根据预处理数据计算城市下垫面承灾体风险暴露参数，构建城市下垫面承灾体风险暴露特征集；然后采集城市下垫面承灾体的网格样本数据，构建承灾体重置成本暴露网格样本；最后根据城市下垫面承灾体风险暴露特征集和承灾体重置成本暴露网格样本，采用随机森林算法进行非均质城市下垫面承灾体风险暴露空间模拟。本发明突破了以行政边界为统计单元的承灾体风险暴露估算误差较大的瓶颈问题，有助于提升城市下垫面风险暴露模拟的准确性，适用于城市防灾减灾和应急管理等领域。

Description

一种非均质城市下垫面承灾体风险暴露空间模拟方法

技术领域

本发明属于灾害风险模拟技术领域，具体涉及一种非均质城市下垫面承灾体风险暴露空间模拟方法的设计。

背景技术

科学准确的城市下垫面承灾体风险暴露空间分布是自然灾害灾前防损、灾中救援和灾后评估的基础数据，对于水利、建筑和保险等领域风险预警和灾损评估具有重要价值。

目前城市下垫面承灾体风险暴露信息获取多以行政边界为单元，采用定期风险普查方式调查统计，对于国家和地区宏观了解自然灾害风险暴露空间分布具有一定参考价值，然而，这种调查统计方式存在以下问题：(1)因统计口径差异和人为主观因素，导致城市承灾体风险暴露估算信息误差较大；(2)难以反映行政边界内部，尤其是城市区域的非均质承灾体风险暴露信息空间分布差异；(3)统计调查的周期较长，难以实现城市承灾体风险暴露信息快速动态更新。因此，传统的统计调查方式一定程度上制约着城市区域承灾体风险暴露科学准确性和时效性，难以满足城市大范围自然灾害风险动态预警和精细评估的需求。

由于城市下垫面承灾体类型复杂、财产密度高及其脆弱性存在空间异质性，其遭受自然灾害的承灾体风险暴露存在显著的非均质特点。当前多采用统计年鉴和土地利用数据进行城市承灾体风险暴露空间模拟，难以展现土地利用单位内部城市下垫面非均质承灾体风险暴露空间差异。为了突破城市承灾体风险暴露空间高精度模拟的技术瓶颈，需要融合夜光灯光卫星遥感、Landsat8/9卫星遥感、土地利用、网络兴趣点、城市建筑轮廓等多源异构数据，开展非均质城市下垫面承灾体风险暴露空间网格化模拟，以期为自然灾害风险预警和灾损理赔提供支撑。

发明内容

本发明的目的是为了解决现有采用统计年鉴和土地利用数据进行城市承灾体风险暴露空间模拟的方法难以展现土地利用单位内部城市下垫面非均质承灾体风险暴露空间差异的问题，提出了一种非均质城市下垫面承灾体风险暴露空间模拟方法。

本发明的技术方案为：一种非均质城市下垫面承灾体风险暴露空间模拟方法，包括以下步骤：

S1、采集城市的多源异构数据并进行预处理，得到预处理数据。

S2、根据预处理数据计算城市下垫面承灾体风险暴露参数，构建城市下垫面承灾体风险暴露特征集。

S3、采集城市下垫面承灾体的网格样本数据，构建承灾体重置成本暴露网格样本。

S4、根据城市下垫面承灾体风险暴露特征集和承灾体重置成本暴露网格样本，采用随机森林算法进行非均质城市下垫面承灾体风险暴露空间模拟。

进一步地，步骤S1中的多源异构数据包括夜间灯光遥感数据、Landsat8/9遥感数据和网络兴趣点数据。

进一步地，夜间灯光遥感数据的采集和预处理方法为：

A1、选取珞珈一号和可持续发展卫星的原始高分辨率夜间灯光遥感数据，选取要求如下：

(1)影像为无云数据；

(2)经过几何和辐射精校正；

(3)无太阳和月亮炫光；

(4)空间分辨率优于130m。

A2、采用多阈值法，识别原始高分辨率夜间灯光遥感数据中背景噪声和短暂性光源的临界域阈值，将对应影像像素值修改为0，消除原始高分辨率夜间灯光遥感数据中的背景噪声和短暂性光源，得到过滤后的夜间灯光遥感数据。

A3、采集多期过滤后的夜间灯光遥感数据，采用多期影像中值合成，获取稳定夜间灯光遥感数据，并采用最近邻像元法重采样为30m分辨率，实现对夜间灯光遥感数据的预处理。

进一步地，Landsat8/9遥感数据的采集和预处理方法为：

B1、根据夜间灯光遥感数据的监测时间，选取时相相近的Landsat8/9遥感数据，选取要求如下：

(1)影像为无云和无冰雪覆盖的数据；

(2)经过几何和辐射精校正；

(3)成像时间为4-10月。

B2、对Landsat8/9遥感数据进行地表反射率计算、空间镶嵌和裁剪，实现对Landsat8/9遥感数据的预处理。

进一步地，网络兴趣点数据的采集和预处理方法为：

C1、通过资源环境数据云平台和高德地图获取相同时间的网络兴趣点数据。

C2、对网络兴趣点数据进行数据去重和坐标转换，获取网络兴趣点位置和属性信息。

C3、根据网络兴趣点位置和属性信息，将网络兴趣点数据分类成交通设施、居民住宅和企业厂矿3个一级类别。

进一步地，步骤S2中的城市下垫面承灾体风险暴露参数包括夜光强度、交通设施的核密度、居民住宅的核密度、企业厂矿的核密度、归一化植被指数、归一化建筑指数和归一化不透水面指数。

夜光强度的计算公式为：

其中P_i，j为(i，j)坐标位置的夜间灯光像素值，A_i，j为(i，j)坐标位置的夜间灯光像素面积，n₁为像素个数。

交通设施的核密度、居民住宅的核密度和企业厂矿的核密度的计算公式为：

其中为归一化后m类别的核密度，取值范围是(0，1)，K^m为m类别的核密度，为m类别核密度的最小值，/>为m类别核密度的最大值，m＝1，2，3，m＝1时表示交通设施，m＝2时表示居民住宅，m＝3时表示企业厂矿，每个类别的核密度计算公式为：

其中f(·)为核密度函数，h为距离阈值，n₂为距离阈值范围内包含的对应交通设施、居民住宅或企业厂矿的网络兴趣点数量，a为任一网络兴趣点，a_i为n₂范围内的网络兴趣点，K(·)为核函数并且有：

其中d_i表示两个网络兴趣点之间的欧式距离。

归一化植被指数NDVI的计算公式为：

其中DN_NIR为Landsat8/9影像的近红外波段的像元值，对应影像中的波段5，DN_RED为Landsat8/9影像的红光波段的像元值，对应影像中的波段4。

归一化建筑指数NDBI的计算公式为：

其中DN_SWIR为Landsat8/9影像的短波红外的像元值，对应影像中的波段6。

归一化不透水面指数NDISI的计算公式为：

其中MNDWI为归一化水体指数，DN_TIR为Landsat8/9影像的热红外波段的像元值，对应影像中的波段10，DN_G为Landsat8/9影像的绿光波段的像元值，对应影像中的波段3。

进一步地，步骤S2中构建城市下垫面承灾体风险暴露特征集的方法为：

通过最近邻重采样方法，将夜光强度、交通设施的核密度、居民住宅的核密度、企业厂矿的核密度、归一化植被指数、归一化建筑指数和归一化不透水面指数重采样为30m×30m，构建城市下垫面承灾体风险暴露特征集Δ：

Δ＝S(V₁，V₂，V₃，V₄，V₅，V₆，V₇)

其中S为集合特征，V₁为重采样后的夜光强度，V₂为重采样后的交通设施的核密度，V₃为重采样后的居民住宅的核密度，V₄为重采样后的企业厂矿的核密度，V₅为重采样后的归一化不透水面指数，V₆为重采样后的归一化植被指数，V₇为重采样后的归一化建筑指数。

进一步地，步骤S3中的城市下垫面承灾体的网格样本数据包括城市建筑物轮廓数据、土地利用数据和承灾体造价数据。

城市建筑物轮廓数据为矢量数据，其属性值包括建筑物名称、位置和楼层高度。

土地利用数据选用公开发布的10m分辨率土地利用二级精细分类数据。

承灾体造价数据根据统计年鉴数据计算得到。

进一步地，步骤S3中构建承灾体重置成本暴露网格样本的方法为：

D1、将城市研究区划分成单元网格尺度为30m×30m的网格，从研究区中随机选取3000个样本网格，满足覆盖12类土地利用类型。

D2、针对每个网格，采用空间聚合方法计算网格内的承灾体重置成本，构建承灾体重置成本暴露网格样本：

T_i＝∑_{s＝1，2，3}V_s×A_s×H_s

其中T_i表示第i个网格内的承灾体重置成本，s表示土地类型，s＝1时表示交通设施，s＝2时表示居民住宅，s＝3时表示企业厂矿，V_s为s类型土地的造价数据，A_s为s类型土地的面积，H_s为s类型土地的建筑物轮廓高度。

进一步地，步骤S4包括以下分步骤：

S41、将3000个样本网格对应位置的城市下垫面承灾体风险暴露特征集作为训练样本集X＝{x_j1，x_j2，...，x_jn}，将承灾体重置成本暴露网格样本作为目标集Y＝{y₁，y₂，...，y_n}，其中j为特征数，n为样本数。

S42、采用二分法对训练样本集X中的连续特征值进行离散化，对特征x₁的连续特征值，即当j＝1时的所有训练样本的特征值由小到大升序排列，记作{F₁，F₂，...，F_n}，得到特征x₁的n-1个切分点：

其中表示第i个训练样本的特征值，/>表示特征x₁的n-1个切分点的集合。

S43、根据特征x₁进行切分，得到n-1个切分点的集合S_X＝{s_i|1≤i≤n-1}。

S44、采用具有K个叶子节点的回归树将训练样本集X划分为K个子集R_k，1≤k≤K。

S45、选取第j个特征x_j和第i个切分点s_i，得到训练样本集X中被切分的两个特征空间子集：

R₁(j，s)＝{x|x_j≤s_i}，R₂(j，s)＝{x|x_j>s_i}

其中x表示样本点。

S46、遍历训练样本集X中的所有特征，寻找最优切分特征x_i和最优切分点s：

其中，x⁽ⁱ⁾为第i个样本点的第j个特征值，y⁽ⁱ⁾为第i个样本点第j个特征对应的目标值；c₁为选择切分特征x_j和切分点s时，特征空间子集R₁对应的所有目标的平均值；c₂为选择切分特征x_j和切分点s时，特征空间子集R₂对应的所有目标的平均值。

S47、将每个特征空间子集所有目标的平均值作为该特征空间子集的输出值c_k：

c_k＝mean(y|x∈R_k)

其中mean(·)表示平均值函数。

S48、重复步骤S46～S47，直到步骤S46算出的平方误差最小，得到随机森林回归模型F(x)：

其中表示第k个样本空间子集的估计值，K表示回归树叶子节点数量，I(x∈R_k)是一个指数，当x∈R_k时I(x∈R_k)＝1，否则I(x∈R_k)＝0。

S49、将城市下垫面承灾体风险暴露特征集输入随机森林回归模型F(x)，输出得到承灾体的风险暴露数值，估算研究区范围内的单元土地利用网格的承灾体重置成本信息，实现非均质城市下垫面承灾体风险暴露的空间模拟。

本发明的有益效果是：

(1)本发明可实现对非均质城市下垫面承灾体的风险暴露空间网格化模拟，突破了以行政边界为统计单元的承灾体风险暴露估算误差较大的瓶颈问题，有助于提升城市下垫面风险暴露模拟的准确性，适用于城市防灾减灾和应急管理等领域。

(2)本发明采用卫星遥感等多源异构数据，构建了非均质城市下垫面风险暴露特征参量，可以客观反映城市下垫面承灾体风险暴露的空间异质性。

(3)本发明构建了网格尺度的承灾体重置成本暴露样本，结合随机森林的智能模型，实现了现城市下垫面风险暴露高精度空间网格化模拟。

附图说明

图1所示为本发明实施例提供的一种非均质城市下垫面承灾体风险暴露空间模拟方法流程图。

图2所示为本发明实施例提供的非均质城市下垫面承灾体风险暴露的空间模拟效果示意图。

具体实施方式

现在将参考附图来详细描述本发明的示例性实施方式。应当理解，附图中示出和描述的实施方式仅仅是示例性的，意在阐释本发明的原理和精神，而并非限制本发明的范围。

本发明实施例提供了一种非均质城市下垫面承灾体风险暴露空间模拟方法，如图1所示，包括以下步骤S1～S4：

S1、采集城市的多源异构数据并进行预处理，得到预处理数据。

本发明实施例中，多源异构数据包括夜间灯光遥感数据、Landsat8/9遥感数据和网络兴趣点(POI)数据。

其中，夜间灯光遥感数据的采集和预处理方法为：

A1、选取珞珈一号和可持续发展卫星的原始高分辨率夜间灯光遥感数据，选取要求如下：

(1)影像为无云数据；

(2)经过几何和辐射精校正；

(3)无太阳和月亮炫光；

(4)空间分辨率优于130m。

A2、采用多阈值法，识别原始高分辨率夜间灯光遥感数据中背景噪声和短暂性光源的临界域阈值，将对应影像像素值修改为0，消除原始高分辨率夜间灯光遥感数据中的背景噪声和短暂性光源，得到过滤后的夜间灯光遥感数据。本发明实施例中，短暂性光源包括火灾、气体火炬、火山等。

A3、由于夜间灯光遥感辐射亮度受季节和成像角度影响，因此采集多期过滤后的夜间灯光遥感数据，采用多期影像中值合成，获取稳定夜间灯光遥感数据，并采用最近邻像元法重采样为30m分辨率，实现对夜间灯光遥感数据的预处理。

Landsat8/9遥感数据的采集和预处理方法为：

B1、根据夜间灯光遥感数据的监测时间，选取时相相近的Landsat8/9遥感数据，选取要求如下：

(1)影像为无云和无冰雪覆盖的数据；

(2)经过几何和辐射精校正；

(3)成像时间为4-10月。

B2、对Landsat8/9遥感数据进行地表反射率计算、空间镶嵌和裁剪，实现对Landsat8/9遥感数据的预处理。

网络兴趣点数据的采集和预处理方法为：

C1、通过资源环境数据云平台和高德地图获取相同时间的网络兴趣点数据。

网络兴趣点(POI)数据可从资源环境数据云平台和高德地图下载和更新，当前，资源环境数据云平台的POI数据时效性较低，可根据高德地图提供的API对下载的数据进行更新。

C2、对网络兴趣点数据进行数据去重和坐标转换，获取网络兴趣点位置和属性信息。

C3、根据网络兴趣点位置和属性信息，将网络兴趣点数据分类成交通设施、居民住宅和企业厂矿3个一级类别。

S2、根据预处理数据计算城市下垫面承灾体风险暴露参数，构建城市下垫面承灾体风险暴露特征集。

本发明实施例中，城市下垫面承灾体风险暴露参数包括夜光强度、交通设施的核密度、居民住宅的核密度、企业厂矿的核密度、归一化植被指数、归一化建筑指数和归一化不透水面指数。

其中，夜光强度的计算公式为：

其中P_i，j为(i，j)坐标位置的夜间灯光像素值，A_i，j为(i，j)坐标位置的夜间灯光像素面积，n₁为像素个数。

交通设施的核密度、居民住宅的核密度和企业厂矿的核密度的计算公式为：

其中为归一化后m类别的核密度，取值范围是(0，1)，K^m为m类别的核密度，为m类别核密度的最小值，/>为m类别核密度的最大值，m＝1，2，3，m＝1时表示交通设施，m＝2时表示居民住宅，m＝3时表示企业厂矿，每个类别的核密度计算公式为：

其中f(·)为核密度函数，h为距离阈值，n₂为距离阈值范围内包含的对应交通设施、居民住宅或企业厂矿的网络兴趣点数量，a为任一网络兴趣点，a_i为n₂范围内的网络兴趣点，K(·)为核函数并且有：

其中d_i表示两个网络兴趣点之间的欧式距离。

归一化植被指数NDVI的计算公式为：

其中DN_NIR为Landsat8/9影像的近红外波段的像元值，对应影像中的波段5，DN_RED为Landsat8/9影像的红光波段的像元值，对应影像中的波段4。

归一化建筑指数NDBI的计算公式为：

其中DN_SWIR为Landsat8/9影像的短波红外的像元值，对应影像中的波段6。

归一化不透水面指数NDISI的计算公式为：

其中MNDWI为归一化水体指数，DN_TIR为Landsat8/9影像的热红外波段的像元值，对应影像中的波段10，DN_G为Landsat8/9影像的绿光波段的像元值，对应影像中的波段3。

通过最近邻重采样方法，将夜光强度、交通设施的核密度、居民住宅的核密度、企业厂矿的核密度、归一化植被指数、归一化建筑指数和归一化不透水面指数重采样为30m×30m，构建城市下垫面承灾体风险暴露特征集Δ：

Δ＝S(V₁，V₂，V₃，V₄，V₅，V₆，V₇)

其中S为集合特征，V₁为重采样后的夜光强度，V₂为重采样后的交通设施的核密度，V₃为重采样后的居民住宅的核密度，V₄为重采样后的企业厂矿的核密度，V₅为重采样后的归一化不透水面指数，V₆为重采样后的归一化植被指数，V₇为重采样后的归一化建筑指数。

S3、采集城市下垫面承灾体的网格样本数据，构建承灾体重置成本暴露网格样本。

本发明实施例中，城市下垫面承灾体的网格样本数据包括城市建筑物轮廓数据、土地利用数据和承灾体造价数据。

其中，城市建筑物轮廓数据为矢量数据，可通过高德或百度平台抓取，其属性值包括建筑物名称、位置和楼层高度。

土地利用数据选用公开发布的10m分辨率土地利用二级精细分类数据。

承灾体造价数据根据统计年鉴数据计算得到，包括交通设施、居民住宅、企业厂矿等每平米的造价数据。

采集得到城市下垫面承灾体的网格样本数据后，构建承灾体重置成本暴露网格样本的方法为：

D1、将城市研究区划分成单元网格尺度为30m×30m的网格，从研究区中随机选取3000个样本网格，满足覆盖12类土地利用类型。

D2、针对每个网格，采用空间聚合方法计算网格内的承灾体重置成本，构建承灾体重置成本暴露网格样本：

T_i＝∑_{s＝1，2，3}V_s×A_s×H_s

其中T_i表示第i个网格内的承灾体重置成本，s表示土地类型，s＝1时表示交通设施，s＝2时表示居民住宅，s＝3时表示企业厂矿，V_s为s类型土地的造价数据，A_s为s类型土地的面积，H_s为s类型土地的建筑物轮廓高度，特别地当土地类型为草地时，H_s＝1。

S4、根据城市下垫面承灾体风险暴露特征集和承灾体重置成本暴露网格样本，采用随机森林算法进行非均质城市下垫面承灾体风险暴露空间模拟。

步骤S4包括以下分步骤S41～S49：

S41、将3000个样本网格对应位置的城市下垫面承灾体风险暴露特征集作为训练样本集X＝{x_j1，x_j2，...，x_jn}，将承灾体重置成本暴露网格样本作为目标集Y＝{y₁，y₂，...，y_n}，其中j为特征数，n为样本数。

S42、采用二分法对训练样本集X中的连续特征值进行离散化，对特征x₁的连续特征值，即当j＝1时的所有训练样本的特征值由小到大升序排列，记作{F₁，F₂，...，F_n}，得到特征x₁的n-1个切分点：

其中表示第i个训练样本的特征值，/>表示特征x₁的n-1个切分点的集合。

S43、根据特征x₁进行切分，得到n-1个切分点的集合S_X＝{s_i|1≤i≤n-1}。

S44、采用具有K个叶子节点的回归树将训练样本集X划分为K个子集R_k，1≤k≤K。

S45、选取第j个特征x_j和第i个切分点s_i，得到训练样本集X中被切分的两个特征空间子集：

R₁(j，s)＝{x|x_j≤s_i}，R₂(j，s)＝{x|x_j>s_i}

其中x表示样本点。

S46、遍历训练样本集X中的所有特征，寻找最优切分特征x_j和最优切分点s：

其中，x⁽ⁱ⁾为第i个样本点的第j个特征值，y⁽ⁱ⁾为第i个样本点第j个特征对应的目标值；c₁为选择切分特征x_j和切分点s时，特征空间子集R₁对应的所有目标的平均值；c₂为选择切分特征x_j和切分点s时，特征空间子集R₂对应的所有目标的平均值。

S47、将每个特征空间子集所有目标的平均值作为该特征空间子集的输出值c_k：

c_k＝mean(y|x∈R_k)

其中mean(·)表示平均值函数。

S48、重复步骤S46～S47，直到步骤S46算出的平方误差最小，得到随机森林回归模型F(x)：

其中表示第k个样本空间子集的估计值，K表示回归树叶子节点数量，I(x∈R_k)是一个指数，当x∈R_k时I(x∈R_k)＝1，否则I(x∈R_k)＝0。

S49、将城市下垫面承灾体风险暴露特征集输入随机森林回归模型F(x)，输出得到承灾体的风险暴露数值，估算研究区范围内的单元土地利用网格的承灾体重置成本信息，实现非均质城市下垫面承灾体风险暴露的空间模拟，模拟效果如图2所示。

本发明实施例中，对非均质城市下垫面承灾体风险暴露进行空间模拟后，采用与承灾体重置成本暴露网格样本相同的构建方法，另选取500个下垫面承灾体风险暴露验证样本点，对随机森林回归模型F(x)估算的值和实际承灾体的造价值进行验证分析，采用均方根误差计算模型精度。

本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。本领域的普通技术人员可以根据本发明公开的这些技术启示做出各种不脱离本发明实质的其它各种具体变形和组合，这些变形和组合仍然在本发明的保护范围内。

Claims (1)

1.一种非均质城市下垫面承灾体风险暴露空间模拟方法，其特征在于，包括以下步骤：

S1、采集城市的多源异构数据并进行预处理，得到预处理数据；

S2、根据预处理数据计算城市下垫面承灾体风险暴露参数，构建城市下垫面承灾体风险暴露特征集；

S3、采集城市下垫面承灾体的网格样本数据，构建承灾体重置成本暴露网格样本；

S4、根据城市下垫面承灾体风险暴露特征集和承灾体重置成本暴露网格样本，采用随机森林算法进行非均质城市下垫面承灾体风险暴露空间模拟；

所述步骤S1中的多源异构数据包括夜间灯光遥感数据、Landsat8/9遥感数据和网络兴趣点数据；

所述夜间灯光遥感数据的采集和预处理方法为：

A1、选取珞珈一号和可持续发展卫星的原始高分辨率夜间灯光遥感数据，选取要求如下：

(1)影像为无云数据；

(2)经过几何和辐射精校正；

(3)无太阳和月亮炫光；

(4)空间分辨率优于130m；

A2、采用多阈值法，识别原始高分辨率夜间灯光遥感数据中背景噪声和短暂性光源的临界域阈值，将对应影像像素值修改为0，消除原始高分辨率夜间灯光遥感数据中的背景噪声和短暂性光源，得到过滤后的夜间灯光遥感数据；

A3、采集多期过滤后的夜间灯光遥感数据，采用多期影像中值合成，获取稳定夜间灯光遥感数据，并采用最近邻像元法重采样为30m分辨率，实现对夜间灯光遥感数据的预处理；

所述Landsat8/9遥感数据的采集和预处理方法为：

B1、根据夜间灯光遥感数据的监测时间，选取时相相近的Landsat8/9遥感数据，选取要求如下：

(1)影像为无云和无冰雪覆盖的数据；

(2)经过几何和辐射精校正；

(3)成像时间为4-10月；

B2、对Landsat8/9遥感数据进行地表反射率计算、空间镶嵌和裁剪，实现对Landsat8/9遥感数据的预处理；

所述网络兴趣点数据的采集和预处理方法为：

C1、通过资源环境数据云平台和高德地图获取相同时间的网络兴趣点数据；

C2、对网络兴趣点数据进行数据去重和坐标转换，获取网络兴趣点位置和属性信息；

C3、根据网络兴趣点位置和属性信息，将网络兴趣点数据分类成交通设施、居民住宅和企业厂矿3个一级类别；

所述步骤S2中的城市下垫面承灾体风险暴露参数包括夜光强度、交通设施的核密度、居民住宅的核密度、企业厂矿的核密度、归一化植被指数、归一化建筑指数和归一化不透水面指数；

所述夜光强度的计算公式为：

其中P_i,j为(i,j)坐标位置的夜间灯光像素值，A_i,j为(i,j)坐标位置的夜间灯光像素面积，n₁为像素个数；

所述交通设施的核密度、居民住宅的核密度和企业厂矿的核密度的计算公式为：

其中为归一化后m类别的核密度，取值范围是(0,1)，K^m为m类别的核密度，/>为m类别核密度的最小值，/>为m类别核密度的最大值，m＝1,2,3，m＝1时表示交通设施，m＝2时表示居民住宅，m＝3时表示企业厂矿，每个类别的核密度计算公式为：

其中f(·)为核密度函数，h为距离阈值，n₂为距离阈值范围内包含的对应交通设施、居民住宅或企业厂矿的网络兴趣点数量，a为任一网络兴趣点，a_i为n₂范围内的网络兴趣点，K(·)为核函数并且有：

其中d_i表示两个网络兴趣点之间的欧式距离；

所述归一化植被指数NDVI的计算公式为：

其中DN_NIR为Landsat8/9影像的近红外波段的像元值，对应影像中的波段5，DN_RED为Landsat8/9影像的红光波段的像元值，对应影像中的波段4；

所述归一化建筑指数NDBI的计算公式为：

其中DN_SWIR为Landsat8/9影像的短波红外的像元值，对应影像中的波段6；

所述归一化不透水面指数NDISI的计算公式为：

其中MNDWI为归一化水体指数，DN_TIR为Landsat8/9影像的热红外波段的像元值，对应影像中的波段10，DN_G为Landsat8/9影像的绿光波段的像元值，对应影像中的波段3；

所述步骤S2中构建城市下垫面承灾体风险暴露特征集的方法为：

通过最近邻重采样方法，将夜光强度、交通设施的核密度、居民住宅的核密度、企业厂矿的核密度、归一化植被指数、归一化建筑指数和归一化不透水面指数重采样为30m×30m，构建城市下垫面承灾体风险暴露特征集Δ：

Δ＝S(V₁,V₂,V₃,V₄,V₅,V₆,V₇)

其中S为集合特征，V₁为重采样后的夜光强度，V₂为重采样后的交通设施的核密度，V₃为重采样后的居民住宅的核密度，V₄为重采样后的企业厂矿的核密度，V₅为重采样后的归一化不透水面指数，V₆为重采样后的归一化植被指数，V₇为重采样后的归一化建筑指数；

所述步骤S3中的城市下垫面承灾体的网格样本数据包括城市建筑物轮廓数据、土地利用数据和承灾体造价数据；

所述城市建筑物轮廓数据为矢量数据，其属性值包括建筑物名称、位置和楼层高度；

所述土地利用数据选用公开发布的10m分辨率土地利用二级精细分类数据；

所述承灾体造价数据根据统计年鉴数据计算得到；

所述步骤S3中构建承灾体重置成本暴露网格样本的方法为：

D1、将城市研究区划分成单元网格尺度为30m×30m的网格，从研究区中随机选取3000个样本网格，满足覆盖12类土地利用类型；

D2、针对每个网格，采用空间聚合方法计算网格内的承灾体重置成本，构建承灾体重置成本暴露网格样本：

T_i＝∑_s＝1,2,3V_s×A_s×H_s

其中T_i表示第i个网格内的承灾体重置成本，s表示土地类型，s＝1时表示交通设施，s＝2时表示居民住宅，s＝3时表示企业厂矿，V_s为s类型土地的造价数据，A_s为s类型土地的面积，H_s为s类型土地的建筑物轮廓高度；

所述步骤S4包括以下分步骤：

S41、将3000个样本网格对应位置的城市下垫面承灾体风险暴露特征集作为训练样本集X＝{x_j1,x_j2,…,x_jn}，将承灾体重置成本暴露网格样本作为目标集Y＝{y₁,y₂,…,y_n}，其中j为特征数，n为样本数；

S42、采用二分法对训练样本集X中的连续特征值进行离散化，对特征x₁的连续特征值，即当j＝1时的所有训练样本的特征值由小到大升序排列，记作{F₁,F₂,…,F_n}，得到特征x₁的n-1个切分点：

其中表示第i个训练样本的特征值，/>表示特征x₁的n-1个切分点的集合；

S43、根据特征x₁进行切分，得到n-1个切分点的集合S_X＝{s_i|1≤i≤n-1}；

S44、采用具有K个叶子节点的回归树将训练样本集X划分为K个子集R_k，1≤k≤K；

S45、选取第j个特征x_j和第i个切分点s_i，得到训练样本集X中被切分的两个特征空间子集：

R₁(j,s)＝{x|x_j≤s_i},R₂(j,a)＝{x|x_j>s_i}

其中x表示样本点；

S46、遍历训练样本集X中的所有特征，寻找最优切分特征x_j和最优切分点s：

其中，x⁽ⁱ⁾为第i个样本点的第j个特征值，y⁽ⁱ⁾为第i个样本点第j个特征对应的目标值；c₁为选择切分特征x_j和切分点s时，特征空间子集R₁对应的所有目标的平均值；c₂为选择切分特征x_j和切分点s时，特征空间子集R₂对应的所有目标的平均值；

S47、将每个特征空间子集所有目标的平均值作为该特征空间子集的输出值c_k：

c_k＝mean(y|x∈R_k)

其中mean(·)表示平均值函数；

S48、重复步骤S46～S47，直到步骤S46算出的平方误差最小，得到随机森林回归模型F(x)：

其中表示第k个样本空间子集的估计值，K表示回归树叶子节点数量，I(x∈R_k)是一个指数，当x∈R_k时I(x∈R_k)＝1，否则I(x∈R_k)＝0；

S49、将城市下垫面承灾体风险暴露特征集输入随机森林回归模型F(x)，输出得到承灾体的风险暴露数值，估算研究区范围内的单元土地利用网格的承灾体重置成本信息，实现非均质城市下垫面承灾体风险暴露的空间模拟。