灾害学研究引入空间信息科学的手段和方法是充分体现灾害事件空间行为特征, 进而深刻揭示灾害在一定空间地域范围内发生、发展、变化的必要要求, 也是向空间信息化方向发展的必然趋势。GIS技术已经大量应用于灾害信息管理和灾害损失评估当中^[1-4]。国外的一些相关研究更是充分利用空间信息科学的强大空间分析及空间模型演算能力, 结合多重、多源灾害信息数据及资料, 实现以综合减灾为目的的决策分析、灾害信息管理及空间模型演算等功能^[5-7]。对灾害事件的空间形态表现及特征的精确和完整定义是空间信息技术应用于灾害研究领域的关键, 也是充分发挥空间信息分析及挖掘能力的基础。记录灾害事件的空间数据需要描述事件发生的空间位置、范围、尺度等信息。按照空间信息科学对空间体的定义, 灾害事件的几何形态有点、线、面等3种形式, 以气象灾害为例, 雷击是以事件发生地点的“点形状”为空间特征表现; 大雾等恶劣天气条件对交通系统带来是线状影响, 或者是由多条道路所组成的线簇^[7-9]。气象灾害中最常见的空间影响形态是“面”, 高温、暴雨、大风等灾害性天气均是如此。灾害性天气在一定空间尺度范围内给城市不同地区造成的影响是不一样的, 而且对应下垫面的社会、经济、环境的特征与受灾害冲击的强度还不一一对应, 有较大的空间对应偏差。如何圈定出这种不对应情况下的灾害影响或危险程度的空间图谱化 (gridded atlas) 分布结果和形态特征, 是灾害评估的重要内容及防治关键, 也是灾害空间评估体系建设的基础和发展灾害学评估理论和方法的要点。

数学形态学方法为确定灾害影响及危险性的图谱化特征提供确实可行的方法和手段。数学形态学的实质是非线性滤波, 腐蚀 (erosion) 和膨胀 (dilation) 是滤波操作的两种数值化处理方法^[10-13]。气象要素值 (温度、湿度、降水量等) 的空间离散分布是度量气象灾害危险性及自然强度的定量化基础, 但是“承灾体”抗灾能力的差异是气象灾害给城市不同下垫面特征的地表造成不同程度影响的主要原因, 也是灾害评估不能简单地以某时刻观测到的气象要素数值大小作为评价依据的主要原因。

本文所采用的空间评估方法是从下垫面的社会、经济、环境构成情况及特征中提取数学形态学滤波中的目标函数, 函数换算值的高低确定“控制腐蚀和膨胀操作”的判断阈值和系数, 在观测要素值的空间离散化图谱上进行动态腐蚀或膨胀操作, 最终完成气象灾害给城市地区造成影响的空间图谱化特征图。讨论结果表明:该图比较贴切于灾害性天气的实际影响情况。

设E为二维离散的欧几里德空间, 图谱化图像A⊂E, 结构元素B⊂E, b ∈ B是欧氏空间中的一个点^[14]。

(1)

(2)

数学形态学的方法大致被分解为二值形态学 (Binary Morphology)、灰度形态学 (Grayscale Morphology)、软数学形态学 (Soft Mathematical Morphology) 等^[14-15]。按图像数据的格式及处理要求, 对灰度图像进行处理的灰度形态学方法是灾害应用研究的首选, 灰度形态学腐蚀运算的定义为:

(3)

灰度形态学膨胀运算的定义为:

(4)

以上定义的灰度形态学的腐蚀和膨胀运算是以结构元素g (i, j) 为模板的一个极值滤波, 从方法上看这种极值滤波显然不能作为灾害空间评估的基础, 因为下垫面的复杂情况决定结构元素g (i, j) 应该是一个动态的、自适应的可变模板, 模板系数的取值应受“下垫面”特征变量的控制。这里暂且将适应灾害空间评估要求的形态学腐蚀和膨胀运算定义如下:

腐蚀定义:

(5)

膨胀定义:

(6)

其中, O_neg和O_pos是根据下垫面的具体情况来决定Θ和⊕运算的选取。在灾害学意义上, 当对应于承灾体上的某一空间位置 (i, j) 处对灾害的承受能力较强时, 就要执行腐蚀运算, 即危险性减弱, 腐蚀程度由E_dyn (g (-i, -j)) 函数值决定; 反之, 则执行膨胀运算, 膨胀程度由E_dyn (g (i, j)) 函数值决定, 表明 (i, j) 处对灾害的承受能力较弱。实现E_dyn的函数表达是量化的根本及决定评估合理性的基础。

E_dyn的函数映射如式 (7)、(8) 所示, 其中t₁, t₂, …, t_n为 (i, j) 处的致灾因子参变量, 比如气象要素值 (观测值或其插值结果)、人口密度、土地覆盖类型、相对地势高差 (可从空间DEM数据中提取)、经济集中程度等, 变量选取由研究灾种的具体情况来定。g (i, j) 的设定也依据灾种的情况。E_dyn的换算函数的取值情况可由权重法或阈值判断的方法来建立, 即以 (i, j) 点处的换算值来确定结构元素g (i, j) 的构成, 建立腐蚀、膨胀的操作序列, 其表述如下:

①阈值法

(7)

②权重法

(8)

假设 (i, j) 处的某个致灾因子参变量t有N个不同权重的目标物类型 (比如下垫面的土地覆盖类型等) 参与计算, 每个目标物类型对应的腐蚀膨胀系数为R_k (N≤k≤1), 如果每个类型在 (i, j) 处所占的比重分别为D_k, 那么E_dyn函数的换算结果, 即结构元素g (i, j) 的取值为R_k×D_k, 其中D_i值可选用各自恰当的方法获取。

以1999年7月22日我国华北大部分地区所经历的高温天气为例, 该日北京市绝大部分地区的日最高气温超过40 ℃, 属于典型的灾害性高温天气, 彩图 1为该日的日最高气温分布示意图。一般来说, 高气温是高温天气对人体带来不舒适感的主要原因, 但不是唯一因素。城市地区不同下垫面类型所具备的潜热、显热交换能力和机制是不一样的, 因此, 其抵抗高温的能力和效果也是不同的, 一般而言, 植被的蒸腾对周边微环境有一定降温效果; 草地、湿地、水体的潜热交换强烈, 也能起较明显的抗高温作用。本案例采用GIS获取下垫面地表特征参数构建腐蚀、膨胀卷积模板, 完成对高温温度分布图基础上的数学形态方法处理, 生成“7·22”高温灾害性天气的“灾害严重度”图谱化分布图。

图 1

图 1. “7·22”北京最高气温分布色标示意图 Fig 1. Map of maximum temperature distribution in Beijing district

一般情况下, 数学形态学卷积有3×3的8邻域及4邻域模板 (图 2), 多数应用选用8邻域模板, 其中=1。

图 2

图 2. 数学形态学卷积模板 (a)8邻域, (b)4邻域 Fig 2. Template of morphological convolution (a)8 neighbors, (b)4 neighbors

假设待处理图像中网格化图谱为一定水平宽、一定垂直高的n×m的阵列, 即

(9)

式 (9) 中, 任一点V_ij处的周边取值与模板的卷积即为点V_ij的动态膨胀或腐蚀结果, 即如果i=1, j=1, i=n, j=m, 那么V_ij=0(边界上的点暂时不做卷积处理), 否则

卷积时, 模板中心从网格化图谱 (式 (9)) 的左上角开始从左到右, 从上到下移动到图的右下角, 每移动到一个新的位置, 完成一次动态膨胀或腐蚀操作, 模板的腐蚀、膨胀系数由所处位置的下垫面的土地覆盖类型动态确定。

步骤1:完成温度等气象要素值的全局网格化。即在北京市空间范围内, 将整个图幅划分为多行多列的网格点, 其中划分网格的空间范围是39.409352°~41.081614°N, 115.390701°~117.515098°E, 共划分了210行、268列 (彩图 3为划分网格的局部效果), 每个网格点大小为东西跨0.008经距, 南北跨0.008纬距, 并与“7·22”高温天气温度分布图 (面图) 进行“空间包含”计算 (rectangle in polygon), 对每个格点赋温度值, 最后网格化的温度图谱显示效果见彩图 4。

图 3

图 3. 模板在地理信息中的结构示意图 Fig 3. Map of template structure indicated in GIS

图 4

图 4. 高温网格化图谱（对应图 3） Fig 4. Grided altas of temperature (correspongding to Fig.3)

步骤2:利用GIS的空间叠置分析, 依据模板网格范围内的土地覆盖类型确定模板系数, 正常情况下, 取3×3大小卷积模板, 由于卷积处理的结果主要对应的是中心点处的值, 而且其取值还受周边数值的影响, 因此模板系数的初始值设定为

(10)

即中心点的权重系数为0.2, 周围8个点的权重系数分别为0.1, 总和为1。

但是, 不同下垫面类型对高温性天气所产生的作用不一样, 绿地、水体等对高温天气的图像形态起腐蚀作用, 而建筑物、裸地、道路水泥地面则起膨胀作用, 它们的膨胀和腐蚀系数取值如表 1所示 (取经验值, 其中系数大于1将执行膨胀, 等于1不执行任务操作, 小于1则执行腐蚀)。

表 1

表 1 不同土地覆盖类型下高温分布图腐蚀膨胀系数 Table 1 Erosion or dilation coefficient of temperature distributing on different types of land cover

表 1 不同土地覆盖类型下高温分布图腐蚀膨胀系数 Table 1 Erosion or dilation coefficient of temperature distributing on different types of land cover

假设在一个卷积模板网格内一共有N种不同的下垫面土地覆盖类型, 每种类型的腐蚀膨胀系数为R_k (N≤k≤1), 如图 3为3×3模板在GIS土地类型图上的空间叠置情况, 如果单个模板单元M_ij (i≤3, j≤3) 分别与N种不同的土地覆盖类型的“多边形图斑”相交, 相交面积为D₁, …, D_N, 那么模板单元M_ij的图像学腐蚀膨胀系数r′_ij=R_k×D_k, 其中D_i的值通过代表模板单元的网格状矩形几何体与单个土地类型图层作overlap空间叠置计算获得。

这时, 原来的图像操作卷积模板 (式 (10)) 取值即转变为下式

(11)

步骤3:利用经过空间下垫面类型换算了的3×3图像腐蚀膨胀模板 (式 (11)), 依次逐点逐行卷积原始高温图像值, 卷积公式也转变为

(12)

式 (12) 中, 当k=2, l=2时, R′_kl=0.2×r′_kl, 其他情况, R′_kl=0.1×r′_kl。

卷积后, 最终得到表示高温灾害影响程度大小的图谱化分布图, 即“7·22”北京高温天气影响状况图 (局部效果可见彩图 5~7)。

图 5

图 5. 根据下垫面类型数学形态动态膨胀腐蚀结果图（对应图 3） Fig 5. The result of dynamical erosion and dilation morphological method by surface underlying (corresponding to Fig.3)

图 6

图 6. 根据下垫面类型数学形态动态膨胀腐蚀结果图叠加上“绿地”（对应图 3） Fig 6. The result of dynamical erosion and dilation morphological method overlapped with the green by surface underlying (corresponding to Fig.3)

图 7

图 7. 根据下垫面类型数学形态动态膨胀腐蚀结果图叠加上“绿地”和“水体”（对应图 3） Fig 7. The result of dynamical erosion and dilation morphological method overlapped with the green and the water by surface underlying (corresponding to Fig.3)

实例针对1999年7月22日典型性高温天气的网格化温度分布图像进行数学形态学动态腐蚀和膨胀操作, 卷积模板的腐蚀和膨胀系数主要由下垫面类型来决定。彩图 3为在GIS上的局部网格划分情况, 各个网格矩形范围内的下垫面类型对灾害性高温天气的抵消和增强效果不一, 其中植被、水体、湿地、草地等通过植物蒸腾、潜热交换、植物光合作用等过程耗散部分热量, 阻止强高温情况出现, 而城市道路、水泥地面、高大建筑物等已被证明是导致城市热岛效应的主要成因^[16]。因此, 数学形态学卷积模板中的腐蚀、膨胀系数是依据网格中下垫面的土地覆盖类型的组成情况动态换算而成。以彩图 3的情况为例, 该图幅范围大致覆盖北京市海淀区颐和园、圆明园附近地区, 其空间地物组成既有水体较为宽广的昆明湖、史上的皇家园林区———颐和园、圆明园等, 又有人口集中、高楼密布、商业活动繁忙的中关村高科技园区, 这种两极分化的“下垫面”状况对研究高温天气及其造成的灾害性影响规律和空间形态分布特征非常具有代表性, 选用这一区域完成数学形态学操作, 可显示出较明显的图像动态腐蚀和膨胀效果。彩图 3处于“选中”状态的3×3网格点中, 西北角上下两个网格基本覆盖上了颐和园区内的“昆明湖”湖体和周边园林区, 其他的6~7个网格内的土地覆盖类型大多为水泥和建筑物用地, 这9个网格点内的土地覆盖类型差异很明显; 彩图 4为网格与“7·22”高温天气温度分布图 (面) 叠合后得到的温度网格图, 其3×3个网格 (选中状态) 的温度值均为41.5 ℃; 彩图 5为通过动态腐蚀和膨胀计算之后的“7·22”高温天气灾害严重度图谱分布, 处于“选中”状态的3×3个网格的处理结果值为

(13)

式 (13) 中, *号表示被执行了膨胀操作, 其他的表示被执行了腐蚀操作。即大致位于颐和园的2个网格的温度值在换算成“灾害严重度”(彩图 5~7中高温色标) 后, 由于下垫面的土地覆盖类型为水体和园林, 被执行了腐蚀操作 (Θ), 计算结果值分别为36.683476和35.605448, 而其他的6~7个网格, 由于下垫面以水泥路面和建筑物用地为主, 则被执行了膨胀操作 (⊕)。另外, 换算值大小还与在空间上距离园林区 (颐和园) 的远近有关, 距离近的网格膨胀轻微, 比如[行3, 列1]网格, 远的则膨胀较大。为了更加清晰地说明这个问题, 彩图 6和7分别为在处理结果图上叠合“绿地”和“水体”两个图层的情况, 两幅图显示结果表明:靠近或包含“水体”和“绿地”等的网格区域, 大多被执行了腐蚀操作 (Θ); 反之, 建筑用地、道路、水泥地面等分布较集中的区域, 则被执行了膨胀操作 (⊕), 腐蚀和膨胀是与网格区域下垫面类型相关的动态操作过程。这也间接证明了依据下垫面类型和特征等参数指标建立数学形态学方法的动态腐蚀和膨胀卷积操作的合理性、可行性及科学性。

数学形态学方法在诸多领域的数值处理和分析上得到较为广泛的应用。本文针对气象灾害的面状影响特征及灾害“影响面”的空间差异, 提出了利用气象灾害发生地的下垫面特征参数, 即能加剧或削弱气象灾害严重度的下垫面土地覆盖类型、人口集中程度、财产集中和分布等特征值, 来建立数学形态学方法的卷积模板中的系数值, 对引发气象灾害的气象要素值的空间离散化图谱进行动态腐蚀和膨胀操作, 获得气象灾害空间评估图谱化结果图。最后, 以1999年7月22日的高温天气为例, 依据下垫面的土地覆盖类型, 建立卷积模板系数, 对高温温度值的离散化结果图谱进行了图像动态腐蚀和膨胀操作。图谱化演算结果比较贴近城市下垫面特征下的高温天气影响情况, 其数字化处理结果有较好的可解释性, 可作为灾害性天气影响下定量化空间评估结果, 为针对气象灾害的减灾、防灾、治灾和相关预案制订提供重要参考。