基于随机森林与ARIMA模型的降水变化与灾害风险评估

Precipitation Variation and Disaster Risk Assessment Based on Random Forest and ARIMA Models

摘要: 随着对地观测技术的飞速发展，我们能够以前所未有的精度和频率获取地球表面的各种数据，通过进行更精细的空间分析和时间序列分析，可以揭示地理环境变化的深层次规律。本文旨在建立中国降水量变化趋势及其与海拔、坡度、土地利用之间的预测模型。通过对降水量、地形因素和五种主要土地覆盖类型的相关性分析，运用Logistic回归和随机森林模型探讨了这些因素对灾害发生的共同影响机制。此外，采用ARIMA时间序列模型预测了未来2025年到2035年间的降水量和土地利用格局，并结合随机森林模型评估了此期间各地区暴雨灾害风险的空间分布。研究结果揭示了在极端天气条件下最脆弱的地区，为灾害防范和土地规划提供了重要参考。

Abstract: With the rapid advancement of remote sensing technologies, we are now able to obtain various data on the Earth’s surface with unprecedented accuracy and frequency. Through more refined spatial and time series analyses, the underlying patterns of geographical environmental changes can be revealed. This study aims to establish predictive models for the trends in precipitation changes in China and their relationships with elevation, slope, and land use. By analyzing the correlations between precipitation, topographic factors, and five major land cover types, the study employs Logistic regression and Random Forest models to explore the joint impact of these factors on the occurrence of disasters. Additionally, the ARIMA time series model is utilized to forecast precipitation and land use patterns from 2025 to 2035, while the Random Forest model is applied to assess the spatial distribution of rainfall disaster risks during this period. The results of the study highlight the most vulnerable regions under extreme weather conditions, providing valuable insights for disaster prevention and land planning.

1. 前言

随着全球气候变化的加剧，极端天气事件如热浪、暴雨、干旱等频繁发生，对人类社会和自然环境造成了深远的影响。特别是中国，作为全球气候变化的敏感区域，正面临着极端天气事件频发、强度增强的严峻挑战。21世纪初以来，华北、华南和西北地区平均年降水量波动上升，而东北和华东地区降水量年际波动幅度增大，这些变化直接关联到极端天气事件的增加。全球变暖背景下，我国年均降水总量增加，降水结构发生了变化，极端降水事件的增加导致干旱风险和短时强降水发生的风险增加。这些变化不仅对农业、水资源管理和社会经济活动产生了重大影响，也对生态系统的稳定性和生物多样性构成了威胁。

大数据技术的引入为地理系统的研究提供了新的工具和方法[1]，使得研究人员能够进行更精细的空间分析和时间序列分析，从而揭示地理环境变化的深层次规律。例如，通过高分辨率的遥感影像[2]，可以详细追踪土地利用/土地覆被的变化；通过长时间序列的气象数据[3]，可以分析气候变化的趋势及其对生态系统的影响。因此，如何利用大数据手段对地理系统进行更深入的综合分析，探索全球气候变化下中国地理环境的演变，成为了当前地球科学研究的关键问题。

本文首先针对2009~2019年降雨、地形和土地利用指标分析以及暴雨频数等指标的统计，通过时间序列角度，绘制2009~2019年平均降雨量，清晰显示随时间变化我国平均降水量的升高或降低。这一时间段内的降水变化趋势对于分析气候变化影响具有重要意义，通过绘制时间序列图来展示这十年间中国平均降水量的变化趋势。然后，研究指标和暴雨频数的相关性，建立降水、地形、土地利用与成灾的相互作用模型，通过逻辑回归和随机森林算法分别建立相互作用模型，并分别进行残差分析和分析评价模型。通过设定输出概率结果大于零的情况下，出现成灾情况。并且利用相关性取最值取得八个指标的临界值。

其次，建立对降水和土地利用的时间预测模型：通过神经网络模型和ARIMA时间序列模型预测2025~2035年降水和土地利用数据的预测并且分析残差和对比模型优劣。最后，通过神经网络模型和ARIMA时间序列模型建立对降水和土地利用的时间预测模型对2025~2035年降水和土地利用数据的预测，形成2025~2035年成灾情况的热力分析图，找出应对暴雨灾害能力最为脆弱的地区。通过一系列建模研究，本文旨在为应对气候变化与极端天气带来的挑战提供数据驱动的科学分析，可以为各地区未来暴雨灾害的防控提供空间分布参考，帮助相关部门制定更具针对性的防灾减灾策略，具有重要的学术价值和现实意义。

2. 建立降水变化与灾害风险评估模型并求解

2.1. 问题概述

降雨、地形和土地利用对于暴雨等极端天气灾害的形成都具有不可忽视的影响。这其中，降雨的时空变异性和不可控性都最强；土地利用作为自然条件和人类活动的综合结果，虽然也随时空演化，但具有一定可控性；地形是最为稳定、不易改变的因素[4]。根据从“暴雨”到“灾害”中上述三方面因素的角色及其交互作用，确定暴雨成灾的临界条件；并结合降雨量和土地利用/土地覆被变化的历史时空演化特征，对2025~2035年间中国境内应对暴雨灾害能力最为脆弱的地区进行预测并可视化。

2.2. 降水、地形和土地利用指标分析

2.2.1. 降水指标分析

本文的高程数据[5]采用平均值重采样，因为降采样时通常要保持局部高度变化。气温数据采用双线性插值进行采样，适合连续的气温数据。将数据集中的正轴割圆锥等面积投影的图像进行栅格化。首先，使用ArcGIS软件，将图像逐像素地进行栅格化，即将每个点的信息转换成网格单元，每个网格单元代表一个固定大小的地理空间区域，通过内插算法完成。从栅格化的图像中可以提取出需要的特征数据，进行分类和统计分析。将不同坐标系的数据统一在同一张地图上展示，便于对比和分析。栅格化后的数据能够直接在地图上显示，方便后续指标的提取和模型的建立。

将NetCDF文件[6]中存储的降水量数据利用Python开源编程软件处理，提取和整合。首先计算全国平均降水量，设为$\overline{P}$ 。数据为每天固定经度纬度降水量，将一年中每天的全国的降水量求总和即为全国平均降水量。i表示第几天，N表示一年总共有多少天，正常一年为365天，闰年有366天。如下表1中计算出随年份变化全国降水量数值变化。

$\overline{P}={\displaystyle \sum _{i=1}^N{P}_i}$ (1)

这个公式用来表示每年全国范围内的降水量变化趋势，用折线图表示$\overline{P}$ 随着时间年份的变化情况，描述出降水量的整体变化趋势。

Table 1. National precipitation during 1990~2005

表1. 1990~2005年全国降水量

年份	全国降水量(毫米)	年份	全国降水量(毫米)
1990	9,834,608	2006	8,910,206
1991	9,287,716	2007	9,168,763
1992	8,859,833	2008	9,773,804
1993	9,733,307	2009	8,685,242
1994	9,499,462	2010	10,292,448
1995	9,236,443	2011	8,411,527
1996	9,493,789	2012	10,204,314
1997	9,022,143	2013	9,831,153
1998	10,510,959	2014	9,565,373
1999	9,170,193	2015	10,091,210
2000	9,195,043	2016	11,027,806
2001	9,009,950	2017	9,934,042
2002	9,868,190	2018	10,197,503
2003	9,581,776	2019	9,768,728
2004	9,007,960	2020	10,510,831
2005	9,442,227

下图中是平均降水量随时间变化而变化，图1中可以看到全国平均降水量在1990年至2020年间呈现出一定的波动，整体上没有明显的长期上升或下降趋势，但存在显著的年际间波动。在2015年左右，全国平均降水量出现了显著的高峰，尤其在2016年达到一个非常明显的峰值。2010年左右降水量也有一个显著的高峰，但在2011年迅速回落。全国平均降水量在过去30年间波动显著，没有明显的长期趋势，但存在一些年际间的极端事件，可能是由于极端降雨导致的。特别是2015年之后的高降水量，表明近年来的气候波动频繁，极端降水事件的频率和强度可能有所增加。

Figure 1. Line graph of national average precipitation over time

图1. 全国平均降水量随时间变化折线图

为了描述降水量长期变化趋势，可以计算降水量的年均变化率，见图2。设$P_{\text{t}}$ 是第t年的降水量变化，则年均变化率为：全国降水增长率的波动较为明显，增减幅度大，呈现出较大的年际变化。2010年降水增长率出现显著的峰值，而在1999年和2011年左右，降水增长率相对较低甚至为负值。2010年的全国降水增速达到了显著的峰值，这可能反映了当年极端气候事件的影响，例如强降雨、洪水等。这类事件可能导致降水量的显著增长。2011年之后降水增速出现了大幅度下降，甚至为负值，这表明全国降水量减少，可能伴随着干旱或降雨不足的气候现象。

$R_t=\frac{{\overline{P}}_t-{\overline{P}}_{t-1}}{{\overline{P}}_{t-1}}\times 100\%$ (2)

Figure 2. National precipitation growth rate during 1990~2020

图2. 1990~2020年全国降水增长率

通过年均变化率的时间序列分析，可以描述降水量的增加和减少趋势。通过折线图表示出来，更加清晰地分析出降水量变化率的增加或减少。

热力图(Heatmap)是一种二维图表，通过使用颜色编码来表示数据的相对值。它通常用于可视化矩阵或网格数据，其中每个单元格的数值通过颜色的深浅来表示。其中，颜色越深表示数值越大。热力图的主要作用是显示数据的模式、趋势和关联性。我们设定经度为横坐标，纬度为纵坐标，颜色数值表示降雨量。

通过对中国大陆0.25˚逐日降水数据集进行分析，可以快速发现数据集中的高值、低值、聚类模式、异常值等，很明显地看出中国南方的降水量普遍高于北方，这与中国的气候带分布有一致性，南方多为亚热带季风气候，北方则是温带大陆性气候。

图3是2010年至2018年间的降水年平均值变化情况，纵轴表示降水的年平均值，横轴表示年份。可以分析出，降水年平均值在2010年至2018年间呈现较为明显的波动，没有明显的长期上升或下降趋势。降水量在不同年份之间存在显著的变化。2010年和2016年是两个显著的高峰年份，降水年平均值接近或超过657。这表明在这两个年份内的降水较为集中，可能伴随着一些极端天气事件。2016年是近年来降水年平均值的最高点，这可能与极端气候事件有关，特别是在2016年全球范围内的气候现象，如厄尔尼诺可能对该年的降水量有重要影响。2011年出现了一个显著的低谷，降水年平均值较前后两年都明显减少，显示出该年的降水量相对较少，可能与干旱或其他气候异常有关。在2012年到2015年期间，降水年平均值在547.5到620.5之间小幅波动。整体上，降水量相对稳定，波动不大。在2016年达到高峰之后，降水年平均值在2017年和2018年出现了小幅度下降，逐步回到接近2014年的水平。

Figure 3. Line graph of annual average precipitation during 2009~2019

图3. 2009~2019降水量年均值折线图

2.2.2. 地形指标分析

1、海拔高度(E)

1) 定义：海拔高度是指地面某一点相对某一基准面的垂直距离，通常以米为单位。基准面一般为海平面，用h来表示。

2) 研究海拔高度意义：高度影响气候条件，如温度和降水。高度变化形成不同的生态系统，影响生物多样性。在地形图中，高度是表示地形起伏的重要因素。

3) 全国海拔高度的空间特征：对[5]中的高程图进行分析，可以看出中国不同地区的海拔高度分布，东部地区地势较为平坦，而西部和西南部地区为世界上最高的山区。此图有助于分析中国地形特征，并可用于进一步的地理、气候或生态环境分析。以地势相对较为平坦，主要为低海拔的平原、盆地和丘陵地带为主；中西部海拔逐渐升高，尤其是青藏高原、西藏和新疆部分地区，海拔显著升高，成为世界上最为高耸的高原和山脉区域。

2、坡度(S)

1) 定义：坡度是指地表某一斜面的倾斜程度，通常用角度或百分比表示，反映了地表的陡峭程度。

2) 研究坡度意义：坡度影响水流和土壤侵蚀，陡坡易导致水土流失。较大坡度地形不利于农业发展，而适中坡度更有利于农作。坡度影响植被分布，陡坡往往植被稀疏。

3) 全国坡度的空间特征：对[5]中的高程图进行分析，中国的东部和东北部区域的坡度相对较小，地势较为平缓。而在西部、西南部以及华北地区的部分区域的度较大，可能代表高原和山区，例如青藏高原、云南高原以及太行山脉等区域。值越大，表示地形越陡峭。

3、坡向(N)

1) 定义：坡向是指坡面的朝向，即坡面所面对的方向，通常以地理方位角度表示。

2) 研究坡向意义：坡向对于山地生态有着较大的作用。坡向影响阳光照射和气温，南坡通常阳光直射，温暖湿润，北坡相对阴凉。坡向影响降水的分布，某些方向可能因地形屏障而形成雨影区。

3) 全国坡向的空间特征：对[5]中的高程图进行分析，不同颜色代表的坡度方向能够帮助分析水流方向、土壤侵蚀、植被分布等与地形相关的自然现象。青藏高原(图中西南部地区)地形坡度变化剧烈，此区域的复杂地形，尤其是高山和峡谷，导致不同方向的坡度分布广泛。东北地区(图中的东北部)相对较平坦，坡度变化较少。华北平原和长江中下游平原(图中的东部和中部区域)相对平坦，地形坡度较小，平面区域占据主导地位。西南部和西北部(如云贵高原、四川盆地等地)也表现出坡度变化，显示了复杂的地形起伏。

坡度信息能够帮助研究水流在地表的汇聚与分散情况，尤其在坡度较大的区域(如青藏高原)，水土流失可能较为严重。了解坡度的方位能够帮助优化灌溉规划，减少水资源浪费，并确保农田能均匀分配水分。在山区修建公路或铁路时，坡度信息有助于优化路线选择，降低建设成本与风险。在滑坡易发区，坡度和坡向信息能够帮助预测和评估灾害风险。

2.2.3. 土地利用指标分析

土地利用/土地覆被：土地利用描述了“人类如何利用土地资源”，其结果可以是农业用地、工业用地、交通用地、居住用地等，提取耕地、林地、草地、灌木丛、湿地的占地面积[7]。为了观察近年来土地利用中五个类别的变化，绘制出2009~2019年各种类型的变化折线图。图4展示了2010年至2018年期间不同土地覆被类型(耕地、林地、草地、灌木丛、湿地)的年均值变化，纵轴表示每种土地覆被类型的比例，横轴表示年份。可以分析出：林地的年均值在图中表现为较高的比例，大约为0.45，说明林地是主要的土地覆被类型，且在2010年至2018年间保持了较为稳定的占比，略有下降趋势。耕地的比例在0.2附近，略低于林地，是第二大土地类型，且在这段时间内基本保持不变，略有波动，但变化不大。草地的比例大约为0.15，且从2010年到2018年保持了相对稳定的趋势。灌木丛的比例在0.05左右，且变化较为平缓。湿地的比例接近0，表明湿地占地面积非常小，且从2010年到2018年变化幅度极小，各类土地覆被类型的比例在2010年至2018年之间变化并不显著，显示出土地利用模式的稳定性。

Figure 4. Land cover type line graph during 2009~2019

图4. 2009~2019年土地覆被类型折线图

2.2.4. 暴雨频数分析

分析全国暴雨频次具有重要的实际意义，特别是在气候变化和极端天气事件频发的背景下，暴雨频次的变化能够为灾害进行预防。下面通过2010年至2018年期间暴雨总频次随年份的变化趋势，纵轴表示该年度的暴雨总频次，横轴表示年份，见图5。

暴雨总频次的波峰和波谷相差较为明显，显示了暴雨在不同年份间的强度差异。2010年和2016年是暴雨频发的高峰年份，其中2016年暴雨频次达到最高峰，接近16,000次。这表明2016年可能是一个暴雨频发的年份，可能伴随了极端天气现象，如强降雨、台风等。2010年也是暴雨频发的高峰年份，接近15,000次，显示当年暴雨较为集中。2011年和2015年出现了明显的低谷，暴雨频次显著减少，接近10,000次，这表明这两个年份暴雨发生的频次较少，可能是相对干旱或降雨强度较低的年份。尤其是2011年，暴雨频次从2010年的高点迅速下降，显示了暴雨事件的剧烈波动。从2016年高峰之后，暴雨频次在2017年和2018年出现下降趋势，但2018年有轻微回升，显示出暴雨频次在这些年间的波动性依然较强。

Figure 5. Storm frequency line graph during 2009~2019

图5. 2009~2019年暴雨频次折线图

2.3. 建立降水、地形、土地利用与成灾的相互作用模型

2.3.1. 降水、地形、土地利用与成灾中指标的皮尔逊相关性分析

进行降水、地形、土地利用与成灾中指标的皮尔逊相关性分析有助于理解这些因素之间的相互关系，并评估它们对灾害形成的影响。皮尔逊相关系数是一种衡量两个变量线性相关性强度的统计方法，其值在−1和1之间，表示从完全负相关到完全正相关的关系。值为0表示没有线性关系。

从图6中可以看出，降水、地形、土地利用与成灾指标之间存在显著的相关性。降水量与暴雨频次呈显著正相关(0.71)，表明降水是暴雨灾害的主要驱动因素。地形方面，海拔与降水量呈负相关(−0.24)，坡度与降水量也呈负相关(−0.3)，反映出地形对降水分布的影响。土地利用中，林地和草地占比与降水及暴雨频次呈正相关，尤其草地占比与暴雨频次的相关性较强(0.36)，而湿地占比与暴雨频次呈负相关(−0.14)，显示出湿地对成灾过程具有一定的调节作用。这表明降水、地形与土地利用类型共同影响着暴雨灾害的发生与分布。

Figure 6. Pearson correlation plot of precipitation, topography, land use and storm frequency

图6. 降水、地形、土地利用与暴雨频次皮尔逊相关性图

针对降水、地形、土地利用与暴雨成灾的相互作用模型，选择Logistic回归与随机森林两种方法进行对比分析，各自的优缺点使其在不同层面上发挥优势。Logistic回归作为一种线性概率模型，具有较强的可解释性，能够通过回归系数明确展示各个指标(如降水量、地形坡度、土地利用类型等)对暴雨成灾概率的影响方向及程度，帮助理解不同指标对成灾机制的贡献。此外，Logistic回归计算高效、适用于大规模数据，并能输出成灾概率，满足本研究中对概率判定和指标临界值的需求。然而，其局限性在于难以捕捉复杂的非线性关系，也无法自动识别各指标之间的交互效应，导致在多变量条件下模型性能可能受限。

相比之下，随机森林模型作为一种集成学习方法，能够有效处理非线性关系，自动捕捉降水、地形与土地利用之间复杂的高阶交互作用，提高模型的预测精度。此外，随机森林还能提供特征重要性分析，揭示各个指标对成灾概率的贡献度，有助于筛选关键指标并设定合理的临界值，从而增强模型在暴雨成灾预测中的科学性和实用性。随机森林在处理高维数据时表现出色，鲁棒性较强，对数据噪声和异常值不敏感，适合本研究中多指标、多区域的数据场景。然而，随机森林模型的不足之处在于可解释性较弱，其输出结果难以直观地揭示指标与成灾之间的具体关系，同时计算复杂度较高，在大数据场景下对计算资源要求较大。

为增强研究方法的合理性与说服力，本文通过对Logistic回归与随机森林的综合对比，利用精确率、召回率、ROC-AUC值等评价指标，量化两者的模型性能差异，弥补单一模型的局限性。具体而言，Logistic回归用于验证指标与成灾概率之间的线性关系，提供清晰的解释性结果，而随机森林则通过非线性建模捕捉更复杂的交互模式，提高预测准确度。此外，结合特征重要性分析与回归系数的方向性，筛选出影响成灾概率的关键指标及其临界值，为后续研究与应用提供科学依据。

2.3.2. 逻辑回归(Logistic Regression)算法介绍

逻辑回归也称作Logistic回归分析，是一种广义的线性回归分析模型，属于机器学习中的监督学习。其推导过程与计算方式类似于回归的过程，但实际上主要是用来解决二分类问题(也可以解决多分类问题)。通过给定的$n$ 组数据(训练集)来训练模型，并在训练结束后对给定的一组或多组数据(测试集)进行分类。

逻辑回归分类根据输入数据不同，可以分为线性可分和非线性可分两种分类。

线性分类与线性回归法不同，逻辑回归不会尝试在给定一组输入的情况下预测数值变量的值。相反，输出是给定输入点属于某个类的概率。为简单起见，假设我们只有两个类(对于多类问题，对应的是多项Logistic回归)，我们希望输出概率$P$ 大于某个值时，则判别为类别$A$ ，输出概率$P$ 小于某个值的时候判为类别$B$ 。因此，Logistic回归模型的输出值总是在[0, 1]中。

在实际数据中，有的数据并不能通过直线或者平面等方式进行直接区分。为了更好地拟合数据，这个时候就需要考虑利用现有的数据点创建更多的特征数据。在这个高维特征数据上训练的逻辑回归分类器将具有更复杂的决策边界，并且在我们的二维图中绘制时将呈现非线性，这就是逻辑回归中的线性不可分。

逻辑回归的预测函数：

$h_{\theta }\left(x\right)=g\left({\theta }^{T}x\right)=\frac{1}{1+e^{-{\theta }^{T}x}}$ (3)

其中：${\theta }_0+{\theta }_1{x}_1+\mathrm{...}+{\theta }_n{x}_n={\displaystyle {\sum }_{i=1}^n{\theta }_i{x}_i}={\theta }^{T}x$ ，对于二分类来说，

$P\left(y=1|x;\theta \right)=h_{\theta }\left(x\right)$ (4)

$P\left(y=0|x;\theta \right)=1-h_{\theta }\left(x\right)$ (5)

将上面两个式子整合可以得到：

$P\left(y|x;\theta \right)={\left(h_{\theta }\left(x\right)\right)}^y{\left(1-h_{\theta }\left(x\right)\right)}^{1-y}$ (6)

通过采用极大似然估计，寻找一组参数$\theta$ 使得每组数据发生概率最大：

$L\left(\theta \right)={\displaystyle \prod _{i=1}^{m}P\left(y_i|x_i;\theta \right)}={\displaystyle \prod _{i=1}^n{\left(h_{\theta }\left(x_i\right)\right)}^{y_i}}{\left(1-h_{\theta }\left(x_i\right)\right)}^{1-y_i}$ (7)

两边同时取对数并除以样本总数$n$ ：

$l\left(\theta \right)=-\frac{1}{n}\ln \left(L\left(\theta \right)\right)=-\frac{1}{n}{\displaystyle \sum _{i=1}^n\left(y_i\log \left(h_{\theta }\left(x_i\right)\right)\right)}+\left(1-y_i\right)\log \left(1-h_{\theta }\left(x_i\right)\right)$ (8)

通过对$\theta$ 求偏导：

$\frac{\partial l\left(\theta \right)}{\partial \theta }=-\frac{1}{n}{\displaystyle \sum _{i=1}^n\left(y_i-g\left({\theta }^T{x}_i\right)\right)}{x}_i$ (9)

正则化项主要是对各个参数的权重做控制的函数，也叫惩罚项，它可以让整个曲线更平滑，可以有效地抑制过拟合。比较常见的有$L1$ 正则化项、$L2$ 正则化项：$L1$ 正则有一个特殊的功能：它可以稀疏化$\theta$ 的取值，即会让某些${\theta }_i$ 值为0，特征权重趋向于稀疏化；$L2$ 正则假设参数的先验分布是$Gaussian$ 分布，可以保证模型的稳定性，也就是参数的值不会太大或太小。在实际使用中，如果特征是高维稀疏的，则使用$L1$ 正则；如果特征是低维稠密的，则使用$L2$ 正则。

因此，$Logistic$ 回归良好的稀疏泛化能力能够有效解决多特征变量影响的预测类问题，有利于信息建模，同时可以更好地建模具有时序关系和较多变量影响的暴雨成灾模型。

2.3.3. 基于Logistc回归算法的降水、地形、土地利用与成灾模型

使用逻辑回归算法[8]可以建立降水、地形、土地利用等特征与成灾情况之间的预测模型。逻辑回归是一个广泛应用的分类模型，适合用于预测未来是否发生暴雨成灾的情况，通过这个模型，可以评估各个特征对灾害发生的影响，并进行灾害风险预测。

1) 准备数据

建立逻辑回归模型变量的数据集：自变量(X)：预测成灾的特征，包括：降水量、地形特征(海拔、坡度)土地利用类型(耕地、林地、灌木丛、草地、湿地的占比)。因变量(y)：成灾标记(0表示未成灾，1表示成灾)。

2) 数据预处理

处理缺失值：清除或填补数据中的缺失值。特征标准化或归一化：由于不同特征的量纲可能差异较大(如降水量与坡度的单位不同)，标准化或归一化是一个常见的预处理步骤，以避免某些特征对模型产生过大的影响。划分训练集和测试集：将数据划分为训练集和测试集，以便于后续模型训练和评估。

3) 训练逻辑回归模型

使用Python的sklearn库，实现逻辑回归模型。得到函数：

$\begin{array}{l}logit\left(p\right)=-2.0914+2.8888X_1-1.4958X_2-0.1022X_3-0.3016X_4\\ -0.5071X_5-0.0953X_6-0.1831X_7-0.0642X_8\end{array}$ (10)

其中，$X_1$ 是降水量，$X_2$ 是海拔高度，$X_3$ 是坡度，$X_4$ ~$X_8$ 分别是耕地、林地、灌木丛、草地、湿地的占比面积。

4) 分析逻辑回归系数

逻辑回归模型通过系数来反映每个特征对成灾的影响。系数的符号和大小可以帮助我们了解哪些特征对成灾有更大的影响。逻辑回归模型通过对各个特征的系数进行分析，可以识别出哪些因素(降水量、地形特征、土地利用类型等)对成灾的影响最大。逻辑回归模型输出的预测值是成灾的概率，可以帮助政府和防灾部门预测特定地区在未来可能发生灾害的概率，进而采取相应的预防措施。

从图7中可以看出，通过混淆模型可以分析出模型表现较好，准确率达到了90%，模型整体上能较好地区分成灾和未成灾的情况。精确率高，召回率较高，模型在预测成灾的样本中，大多数预测是正确的(90.23%)，但在实际成灾的样本中，有19.63%被误分类为未成灾(召回率为90.36%)。其中，准确率(Accuracy)：模型在测试集上预测正确的比例。精确率(Precision)：预测为正例中，真正例的比例。召回率(Recall)：实际正例中，被正确预测为正例的比例。F1-score：精确率和召回率的调和平均，反映模型的整体表现。

Figure 7. Logistic regression confusion matrix plot

图7. Logistic回归混淆矩阵图

2.3.4. 随机森林(Random Forest, RF)算法介绍

随机森林是一种属于Bagging类型的集成算法，“随机”使它具有抗过拟合性，“森林”使模型具有更高的精确度。首先，RF使用CART决策树作为弱学习器，生成的每棵树仅随机选出少数特征，使用训练出来的分类器集合对新样本进行分类，使用多次投票统计所有分类器的分类结果，确定最高类别为最终标签。

CART树又称为分类回归树。当数据集因变量为连续性数值时，该算法可看成一个回归树，用叶节点观察的均值作为预测值；当数据集的因变量为离散型数值时，该算法就是一个分类树，可解决分类问题。此外，若某个非叶节点是连续变量时，决策树也将把它当作离散变量来处理。

特征选择常用方法是信息增益、增益率、基尼系数(GINI)和卡方检验。随机森林采用的CART决策树就是基于基尼系数选择特征的。对于一般的决策树，假如总共有$K$ 类，样本属于第$k$ 类的概率为$p_k$ ，则该概率分布的基尼系数为：

$Gini\left(p\right)={\displaystyle \sum _{k=1}^K{p}_k\left(1-p_k\right)}=1-{\displaystyle \sum _{k=1}^K{p}_k^2}$ (11)

基尼指数越大，说明不确定性就越大；基尼指数越小，不确定性越小，数据分割越彻底，对于CART树而言，可以通过下面的表示：

$Gini\left(p\right)=2p\left(1-p\right)$ (12)

在我们遍历每个特征的每个分割点时，当使用特征$A=a$ ，将$D$ 划分为两部分，即$D_1$ (满足$A=a$ 的样本集合)和$D_2$ (不满足$A=a$ 的样本集合)。则在特征$A=a$ 的条件下$D$ 的基尼指数为：

$Gini\left(D,A\right)=\frac{\left|D_1\right|}{\left|D\right|}Gini\left(D_1\right)+\frac{\left|D2\right|}{\left|D\right|}Gini\left(D_2\right)$ (13)

其中，$Gini\left(D\right)$ 表示集合$D$ 的不确定性；$Gini\left(A,D\right)$ 表示经过$A=a$ 分割后的集合$D$ 的不确定性。

随机森林中的每棵CART决策树都是通过不断遍历这棵树的特征子集的所有可能的分割点，寻找基尼系数最小的特征分割点，将数据集分成两个子集，直到满足停止条件为止，见图8。

Figure 8. Random forest algorithm flowchart

图8. 随机森林算法流程图

2.3.5. 随机森林算法建立降水、地形、土地利用与成灾模型

使用随机森林算法[9]可以建立降水、地形、土地利用与成灾的预测模型。随机森林是一种基于决策树的集成学习方法，能够有效应对多种特征，并具有较强的处理非线性关系的能力。它适合用于解决成灾与未成灾的二分类问题。

随机森林算法的主要步骤：

1) 准备数据

要建立逻辑回归模型烦人变量的数据集：自变量(X)：预测成灾的特征，包括：降水量、地形特征(海拔、坡度)土地利用类型(耕地、林地、灌木丛、草地、湿地的占比)。因变量(y)：成灾标记(0表示未成灾，1表示成灾)。

2) 数据预处理

检查并处理数据中的缺失值，可以使用均值填充或删除包含缺失值的样本。虽然随机森林不要求数据标准化，但根据实际情况，可以对特征进行标准化或归一化。划分训练集和测试集：将数据划分为训练集和测试集(例如80%的数据用于训练，20%的数据用于测试)。

3) 使用随机森林训练模型

随机森林通过构建多个决策树，并利用多数投票或平均结果来进行分类或回归。随机森林模型通过多个决策树的集成，能够准确预测降水、地形、土地利用等特征与成灾之间的关系。它能处理非线性特征之间的复杂关系，适用于灾害预测这种多因素综合影响的场景。通过分析随机森林模型中的特征重要性，可以识别出哪些因素(如降水量、坡度、耕地比例等)对成灾有更大的影响。这对灾害管理、土地利用规划以及政策制定提供了有力的数据支持。通过预测哪些地区在暴雨、特定地形条件下更容易发生灾害，相关部门可以提前制定防灾策略，建设防洪设施，优化土地利用规划，以减少暴雨等极端天气带来的灾害风险。随机森林能够处理多维数据，并且对线性关系和非线性关系都具有很好的处理能力。它不容易过拟合，适合大规模、复杂的数据集。

Table 2. Comparison table of logistic regression and random forest evaluation metrics

表2. 逻辑回归和随机森林评价指标对比表

评价指标	准确率	精确度	召回率	F1	ROC-AUC
Logistic回归	0.9036	0.9023	0.9036	0.9018	0.9593
随机森林	0.9055	0.9052	0.9055	0.9054	0.9627

从表2中可以发现，通过随机森林算法建立降水、地形、土地利用与成灾模型，可以实现较高的灾害预测准确率，并识别出影响成灾的关键特征。

Figure 9. Random forest confusion matrix plot

图9. 随机森林混淆矩阵图

通过随机森林混淆矩阵(见图9)可以推断，随机森林模型的准确率为90.55%，表明其大多数预测是正确的，整体表现较好。精确率与召回率平衡：精确率为90.52%，召回率为90.55%，两者比较接近，说明模型在成灾样本预测中的表现比较平衡，既能正确预测出多数成灾事件，同时保持较低的错误警报率。尽管模型表现不错，仍可以进一步提升召回率，以减少假负例(未能预测成灾)的数量，这对于防灾减灾来说尤为重要。如果需要进一步调整模型或改进性能，可以通过调优参数、处理数据不平衡等手段进行优化。

2.3.6. 预测临界值

确定降水量、海拔、坡度、耕地、林地、草地、灌木丛和湿地等特征的临界值是一个重要的步骤，这可以识别哪些特征的值在什么范围内会显著影响暴雨成灾的概率。

首先在Logistic模型中，特征的系数可以提供该特征与结果之间的关系强度。系数的正负表示特征与成灾风险的正负关系。

通过分析这些系数，可以初步判断特征的重要性和影响方向。之后我们利用分位数法来判定具分位数法是一种统计分析中常用的描述数据分布位置的数据划分方法，尤其适用于将数据集分为几个等距的部分。

分位数法的基本步骤：首先，选择分位数级别：常见的分位数有四分位数(Q1, Q2, Q3)，其中Q1是第一下四分位数，Q2是第二是中位数，Q3是第三上四分位数。对于每个分位数，找到对应数据集中特定比例的数值。因为本文的特征如降水量是无序的，因此需要先进行排序。通过Python编程得到临界，值如表3所示。

Table 3. Table of critical values by the quantile method

表3. 分位数法求临界值表

临界值	降水量	海拔	坡度	耕地	林地	草地	灌木丛	湿地
数值	1070	538.8	73.4911	0.005	0.015	0.015	0.015	0.01

2.4. 建立对降水和土地利用的时间预测模型

2.4.1. ARIMA模型

时间序列模型[10]是一种用于对时间序列数据进行建模和预测的统计模型。它基于过去观察到的时间序列数据的特征和模式，通过数学方法来描述和捕捉数据的趋势性、周期性、季节性等规律，从而实现对未来数值或趋势的预测。时间序列模型可以分为两大类：描述性模型和预测模型。

描述性模型：描述性模型旨在对时间序列数据进行总结和解释。常见的描述性模型包括时间序列的平均值、方差、自相关性、平稳性检验等。这些模型提供了对时间序列数据统计特性的度量和分析，帮助我们理解数据的基本性质。

预测模型：预测模型是在描述性模型的基础上，利用时间序列数据的历史信息，建立数学模型来预测未来的数值或趋势。预测模型可以用于短期预测和长期预测。常见的时间序列预测模型包括移动平均、自回归、移动平均自回归、差分自回归移动平均等。其中，ARIMA是一种常用的时间序列分析和预测方法。ARIMA模型[11]结合了自回归(AR)和移动平均(MA)的概念，并引入了差分(I)的操作来处理非平稳时间序列数据。

ARIMA模型的三个关键参数分别是p、d和q:p (AR阶数)：表示自回归模型中使用的历史观测值的数量。具体而言，ARIMA模型使用过去p个时刻的观测值作为预测当前时刻的值的基础。较高的p值可以提供更多的历史信息，但也可能增加模型的复杂度。

d (差分阶数)：表示对时间序列进行差分的次数。差分是为了使非平稳时间序列变得平稳。通过差分操作，我们可以消除趋势、季节性等影响，使时间序列更加稳定。

q (MA阶数)：表示移动平均模型中使用的滞后误差的数量。具体而言，ARIMA模型使用过去q个时刻的滞后误差作为预测当前时刻的值的基础。较高的q值可以捕捉到时间序列中的随机波动。

AR和MA模型的数学表达式如下：

$AR：Y_t=c+{\alpha }_1{Y}_{t-1}+{\alpha }_2{Y}_{t-2}+\cdots +{\alpha }_p{Y}_{t-p}+{\xi }_t$ (14)

$MR：Y_t=\mu +{\epsilon }_1{}_t+{\beta }_1{\epsilon }_{t-1}+\cdots +{\beta }_q{\epsilon }_{t-q}$ (15)

若暂时不考虑差分，那么ARIMA模型可以被看作是AR模型和MA模型的直接结合，形式上看，ARIMA模型的公式可以表示为：

$ARMA：Y_t=c+{\alpha }_1{Y}_{t-1}+{\alpha }_2{Y}_{t-2}+\cdots +{\alpha }_p{Y}_{t-p}+{\xi }_t+{\epsilon }_1{}_t+{\beta }_1{\epsilon }_{t-1}+\cdots +{\beta }_q{\epsilon }_{t-q}$ (16)

其中，$Y_t$ 使我们正在考虑的时间序列数据，${\alpha }_1$ 到${\alpha }_p$ 是AR模型的参数，这些参数用来描述当前值与过去p个时间点之间的值，${\beta }_1$ 和${\beta }_q$ 是MA模型的参数，这些参数用来描述当前值与过去q个时间点的误差之间的关系，${\xi }_t$ 是在t时间的误差项，c是一个常数项。

这个公式的基础是假设我们正在处理的时间序列是平稳的，这样我们可以直接应用AR和MA模型。如果时间序列是非平稳的，那么我们就需要考虑ARIMA模型中的I部分，也就是进行差分处理。

残差分析是统计学和回归分析中常用的一种方法，用于评估模型的拟合程度和检查模型的假设是否满足。在线性回归中，残差是指观测值与回归模型预测值之间的差异。

残差(Residual)是指观测值与回归模型预测值之间的差异，用于衡量模型的拟合程度和检测异常值。简单线性回归模型中的残差表示为：

$\text{R}=\text{Y}-\widehat{\text{Y}}$ (17)

式中，$\text{R}$ 为残差，$\text{Y}$ 为观测值，$\widehat{\text{Y}}$为回归模型对应的预测值。

标准化残差(Standardized Residual)是将残差除以其标准差，以消除量纲差异，并进行比较和分析。标准化残差通过除以残差的标准差来得到：

$\text{Z}=\text{R/}\sigma$ (18)

式中，Z为标准化残差，R为残差，$\sigma$ 的残差标准差。

残差分析的主要目的：检验模型的拟合优度：通过观察残差的分布，判断模型是否能够很好地解释数据的变异。检查模型的假设：如线性关系、误差项的独立性、方差齐性和正态性等。

进行残差分析的步骤通常包括以下几个方面：绘制残差图：将观测值的残差绘制在横坐标上，以检查其是否呈现特定的模式或结构。绘制预测值与残差的散点图：通过观察预测值与残差的关系，检查残差是否与预测值存在系统性的相关性。绘制残差的直方图或概率图：检查残差是否符合正态分布，可以通过绘制直方图或概率图来评估残差的分布情况。进行统计检验：如Ljung-Box检验、Durbin-Watson检验等，用于检验残差的自相关性和异方差性。通过残差分析，可以发现模型的不足之处，并进行修正或改进。同时，残差分析也有助于解释模型的预测能力和适用性，并帮助我们更好地理解数据和变量之间的关系。

ARIMA中的AR是自回归模型，利用自身变量历史出现的数据对自身未来时间进行预测，描述了目前情况下当前值和历史值之间的关系。自回归需要满足平稳性条件。时间序列的平稳性指的是在一段时间某个数据看起来没有大波动，很平坦。从数据上来说就是各阶统计特征不随时间的变化而变化。自回归模型也是有一定的限制，比如自回归模型需要用自身的数据进行预测并且需要数据具有平稳性和自相关性。ARIMA中的MA是移动平均模型，移动平均模型对于自回归模型中误差项的累加十分关注。并且通过移动平均法可以适当地消除在预测过程中随机的波动。因此，两部分组成的ARIMA为差分自回归移动平均模型，如图10所示。

Figure 10. Time series ARIMA flowchart

图10. 时间序列ARIMA流程图

2.4.2. 时间序列ARIMA预测降水和土地利用

在本文的研究中，选择ARIMA模型用于对未来2025~2035年降水量与土地利用数据的时间预测，具有清晰的理论依据和实际应用价值。ARIMA模型是一种经典的时间序列预测方法，能够有效捕捉数据的趋势性、周期性与随机扰动，尤其适用于处理具有线性趋势的时间序列数据。在降水量和土地利用变化中，时间维度的数据往往表现出长期趋势与年际波动，ARIMA通过差分平稳化处理和参数估计(自回归阶数p、差分阶数d、移动平均阶数q)，能提供数学上高度解释的预测结果，清晰揭示未来数据变化的整体规律。此外，ARIMA模型结构简单，计算效率高，非常适合长时间跨度的数据分析，如本文中对未来十年降水与土地利用变化的预测任务。

与Logistic回归和随机森林模型相比，ARIMA模型在时间序列预测任务中具有显著优势。Logistic回归主要用于解释降水、地形、土地利用等指标与与暴雨成灾概率之间的线性关系，具有结果可解释性强、计算高效等优点，但无法处理时间动态特征，且难以捕捉复杂的非线性关系。而随机森林模型作为一种集成学习方法，擅长处理非线性关系与多特征交互作用，能够揭示各指标对暴雨成灾的非线性贡献，并提供特征重要性分析。然而，随机森林主要用于静态数据建模，对时间序列数据的建模需额外调整，且计算复杂度较高，模型结果的可解释性相对较弱。因此，在本文的时间预测任务中，ARIMA模型弥补了Logistic回归和随机森林无法处理时间动态关系的不足，成为时间序列预测的优选工具。

在实际应用中，本文通过ARIMA模型预测2025~2035年间降水和土地利用的变化趋势，将结果与成灾模型相结合，分析未来降水变化对暴雨灾害的潜在影响，形成时间动态预测与空间分布分析相结合的研究框架。以下是详细的步骤和每个步骤的意义。

1) 数据准备

首先，获取并准备全国降水量的时间序列数据。时间序列按照年的时期对应每个时间点的降水量。

2) 可视化数据

在使用ARIMA之前，可视化数据，了解全国降水量随时间的变化趋势，见图11。可以帮助理解数据的基本趋势、季节性等特征。便于观察是否有明显的上升、下降或周期性波动。

Figure 11. Average annual precipitation during 2009~2019

图11. 2009~2019年年降水量平均值

3) 数据平稳性检验

ARIMA模型要求数据是平稳的，即数据的均值和方差随时间不发生显著变化。我们可以使用ADF (Augmented Dickey-Fuller)检验来检查数据的平稳性。表4以降水量为例。

Table 4. Smooth line test for precipitation data

表4. 降水数据平稳性检验

检验结果	ADF统计量	1%的显著性水平	5%的显著性水平	10%的显著性水平
降水量	−4.7116	−4.0689	−3.1271	−2.7017

4) 选择ARIMA模型参数(p, d, q)

ARIMA模型有三个参数：

p：自回归(AR)的阶数，表示使用多少个过去的值来预测当前值。

d：差分次数，表示进行了几次差分使数据平稳。

q：移动平均(MA)的阶数，表示模型使用过去的预测误差项的数量。

通过自相关图(ACF)和偏自相关图(PACF)来选择最佳的p和q值。

Figure 12. Autocorrelation (ACF) and partial autocorrelation (PACF) plots

图12. 自相关图(ACF)和偏自相关图(PACF)

根据图12可知，选取的模型为ARIMA(1, 0, 2)。

5) 构建ARIMA模型

基于前面选出的参数，构建ARIMA模型，并对降水量和土地利用进行建模。得到预测图可视化如图13、图14所示。

Figure 13. Time series projection diagram of precipitation during 2025~2035

图13. 2025~2035年降水量时间序列预测图

通过预测未来的降水趋势，可以帮助相关部门提前规划防洪、排涝等防灾措施，特别是在气候变化带来的极端天气事件增加的情况下。降水量预测对于农业生产非常重要，尤其是在农业活动依赖降水的地区。未来降水量的变化将直接影响作物的种植安排、灌溉计划等。水资源管理部门可以根据未来的降水预测，合理安排水库、河流的水量调度，确保在干旱或多雨季节实现水资源的合理分配。

Figure 14. Sequential projection diagram of cover types during 2025~2035

图14. 2025~2035年覆被类型序列预测图

ARIMA模型预测土地利用的意义在于理解和量化土地利用随时间的动态变化趋势，为政策制定、城市规划和自然资源管理等领域提供科学决策支持。土地利用的变化是一个涉及农业、林业、城市扩展等多因素驱动的复杂过程，通过时间序列分析能够捕捉这些动态变化，并提供未来的变化趋势预测。ARIMA模型凭借其高效性和准确性，可以为未来土地利用变化提供量化预测结果，揭示土地利用演变的长期趋势，有助于在气候变化背景下深入理解降水结构与土地利用演变的规律，为防灾减灾策略的制定提供科学依据和空间分布参考，从而为应对极端天气事件和资源管理挑战提供数据驱动的决策支持。

2.5. 预测成灾脆弱区

2.5.1. 预测结果

通过森林算法预测的在2025年至2035年期间的灾害成灾情况表示，尽管整体的成灾情况分布保持一定的一致性，但在某些年(2028年、2033年、2035年)可以看到某些地区的灾害成灾情况有所加剧，特别是在中国中部和东南沿海地区。并且统计出这11年期间的总体灾害成灾趋势，尤其是西南部和东南沿海地区在未来十多年里灾害成灾的风险较高。测结果表明，西南部和东南沿海地区在未来十多年里将面临较高的灾害风险，值得进一步关注和防范。

2.5.2. 预测成灾脆弱区与决策

1) 西南部地区：预测结果显示中国的西南部地区(如云南、贵州、四川部分地区)是未来成灾的高风险区。该区域多山，可能面临地质灾害(如泥石流、滑坡)或气象灾害(如强降雨、洪水)的高发风险。

2) 东南沿海地区：东南沿海区域(如广东、福建、浙江)在预测期内显示成灾风险较高，可能与强台风、暴雨和海岸线洪涝灾害相关。随着气候变化的加剧，沿海地区的成灾风险在未来可能进一步上升。

3) 中部地区：图中显示中国中部部分地区(如湖北、湖南)在某些年份的灾害风险显著，这可能与极端天气(如暴雨洪涝或干旱)相关，需进一步监测。

综上所述，在脆弱区(尤其是西南部和东南沿海地区)加强灾害监测与预警系统建设，确保在灾害发生前能够及时发布预警，减少人员伤亡和财产损失。西南部地区由于地质灾害风险高，需加强山区的基础设施建设，如修建抗灾能力强的道路、桥梁和住房，防止泥石流和滑坡带来的灾害。加强沿海地区的堤坝、防洪设施和台风防御措施，减轻台风和海潮灾害的影响，尤其是城市与工业集中区的防洪能力要进一步提升。在预测的高风险区，地方政府和相关机构应根据预测数据制定详细的应急预案，并定期组织应急演练，确保在灾害发生时能迅速响应。由于灾害的长期变化与气候变化密切相关，需制定适应气候变化的长期策略，例如加强气候变化监测、减少温室气体排放，防止极端天气频率加剧。在脆弱区推行生态修复项目，如植树造林、湿地保护等，利用生态系统缓冲灾害影响，减少洪水、干旱等灾害的破坏性。

3. 总结

本文聚焦全球气候变化背景下中国极端天气事件频发的现象，研究降水、地形、土地利用与暴雨成灾的相互关系，并通过模型分析和时间预测为防灾减灾提供科学依据。首先，基于2009~2019年的数据，本文通过时间序列分析展示了中国降水变化的趋势，揭示了不同地区降水波动与极端天气事件的关联性。随后，本文利用Logistic回归和随机森林模型分别建立了成灾风险评估模型。Logistic回归通过其线性解释性揭示了各指标对成灾概率的影响方向；随机森林则捕捉了非线性关系与指标间的交互效应，并通过特征重要性评估筛选出关键因素。通过残差分析和模型对比，明确了两者的性能差异和适用场景。此外，本文利用ARIMA时间序列模型和神经网络预测了2025~2035年间降水和土地利用的变化趋势，识别出高风险地区。本文通过静态指标建模与动态时间预测相结合，揭示了气候变化下降水与土地利用的演变规律，为防灾减灾策略的制定提供了数据支撑和科学参考。