2024年全国研究生数学建模竞赛华为杯D题大数据驱动的地理综合问题求解全过程文档及程序
2024年全国研究生数学建模竞赛华为杯
D题 大数据驱动的地理综合问题
原题再现：
地理系统是自然、人文多要素综合作用的复杂巨系统[1-2]，地理学家常用地理综合的方式对地理系统进行主导特征的表达[3]。如以三大阶梯概括中国的地形特征，以秦岭—淮河一线和其它地理区划的方式揭示中国气温、降水、植被、土壤及生态环境在水平和垂直方向上的地带性与非地带性规律，利用胡焕庸线[4]、T型开发结构等描绘我国人口、社会和经济发展的总体格局。这些方法早期以宏观结构和定性分析为主体，对我国生态保护、社会经济发展和国家安全保障起到了巨大的支撑作用。伴随着对地观测体系的快速发展，当前已经积累了巨量的对地观测数据。如何利用大数据的手段对地理系统进行综合[5-6]，探索全球气候变化下中国地理环境的演化，是当前地球科学研究的关键问题。
请利用“附件数据”一栏中列出的描述全球和中国地理不同方面特征的数据集，回答以下问题：
1.在众多描述地理环境的变量中，一些简单的指标背后蕴藏了深厚的内涵，对人类的生存发展具有重大深远的影响，如大气中二氧化碳的浓度、全球年平均气温等。降水量是一个连续变化的变量，而土地利用/土地覆被类型则是一个存在突变和离散分布的变量。同时，它们都具有时空分布不均匀的特征。请从附件数据中选取相关数据集，为这两个变量分别构建一套描述性统计方法，用13个较为简洁的统计指标或统计图表，对这两个变量在19902020年间中国范围内的时空演化特征进行描述和总结。
2.近年来，以暴雨为代表的极端天气事件对人类的生产生活造成了越来越难以忽视的影响。请结合附件中所给的数据，建立数学模型，说明地形-气候相互作用在极端天气形成过程中的作用。
3.降雨、地形和土地利用对于暴雨等极端天气灾害的形成都具有不可忽视的影响。这其中，降雨的时空变异性和不可控性都最强；土地利用作为自然条件和人类活动的综合结果，虽然也随时空演化，但具有一定可控性；地形是最为稳定、不易改变的因素。请考虑第2问所反映的从“暴雨”到“灾害”中上述三方面因素的角色及其交互作用，确定暴雨成灾的临界条件；并结合第1问中降雨量和土地利用/土地覆被变化的历史时空演化特征，对2025~2035年间中国境内应对暴雨灾害能力最为脆弱的地区进行预测。请以地图的形式呈现你们的预测结果。
在中国级别的尺度上，描述自然地理特征的地形可以概括为“三级阶梯”，而降水中具有标志性意义的“800mm等降水量线”则与区分我国南北方的“秦岭—淮河”一线大体重合；描述人文地理特征的人口分布及其社会经济活动总量等指标，则被由连接黑龙江黑河与云南腾冲的“胡焕庸线”清晰地划分成东密西疏的两部分。那么，对于自然地理和人文地理交汇点的土地利用/土地覆被情况，结合其在前三问中描述、估计和预测任务中的“特性”，利用地理大数据，建立相应的数学模型，对数据进行简化和综合，描述中国土地利用变化的特征与结构。从准确性和有用性两个方面解释验证你们的总结。
整体求解过程概述(摘要)
地理系统是多要素作用的复杂系统。根据中央气象台和国家气象信息中心国家级气象观测站数据统计，以暴雨为代表的极端天气事件给人类的生产生活带来了巨大的影响，2012—2023 年（截至 2023 年 9 月 15 日），我国发生暴雨天气事件共 503 次。本文利用地理大数据，使用统计测度、机器学习和统计模型等方法，系统地研究了极端天气成因及其影响因素、地形—气候间的交互作用和我国地理变化的特征与结构。
针对问题一，降水量和土地利用（土地覆被类型）时空演化特征建模。首先，从数据集 3（日降水量数据）和数据集 4（年土地利用类型数据）中提取 1990-2019 年的经纬度、降水量和土地利用类型的数据。其次，对降水数据定义不同年份的年均降水量、月均降水量、降水距平百分率和区域方差指标，对土地利用数据定义每年的动态度和成分度。再次，对月均降水量进行时间序列分解，提取其长期趋势、循环因素和随机误差。最后，对上述原始数据以及各种时空指标使用箱线图、折线图、地图、堆砌条形图进行可视化，得出如下结论：
针对降水量：（1）西北地区和青藏地区是全国年均降水量最少的区域。华中地区与东北地区降水量适中。（2）华南地区、华东沿地区海和西藏南部地区降水量变化大，容易发生极端降水事件。（3）自 2014 年，华中地区容易发生极端降水事件。针对土地利用/土地覆被：（1）草地利用分布范围最广，湿地最少。耕地利用率最高在东北地区、华北地区和华中地区；黑龙江地区和华东沿海地区的林地利用率最高；内蒙古中部地区的灌木丛利用率最高；西北中部地区的湿地利用率最高。（2）湿地和灌木丛的利用率最低，且呈现逐年递减的趋势。土地类型利用率最高的是草地。（3）各土地利用类型的动态度均呈现不同程度的递减趋势，其中灌木丛和湿地利用率持续呈现递减趋势。更详细的时空结果见正文。
针对问题二，地形—气候在极端天气形成中的交互作用研究。首先，提取数据 1（山影、坡度、坡向、海拔和高程）和数据 4（土地利用）中的地形数据以及数据 2（日气温）和数据 3（日降水）的气候数据。通过重抽样将不同数据源的像素统一成 0.5。，根据经纬度、时间将这些数据连接到一起，共得到 6169803 条数据，并将降水量通过等级划分对数据进行打标签。其次，检验整合后数据库的缺失值并进行删除和插值得到 5990100 条数据库（具体数据变量见正文），通过对降水数据划分等级来探究数据的不平衡特征，其中特大暴雨仅有 92 条，占全数据的 0.0025%，可见该数据库中的数据呈现出极其不平衡状态。再次，采用欠采样技术抽取“无雨”数据 6000 条，“小雨”数据 5000 条，“中雨”数据 4000条，“大雨”数据 3000 条，“暴雨”数据 1500 条，“大暴雨”数据 1208 条，“特大暴雨”数据 92 条建立抽样数据库。最后，基于抽样数据对 92 条特大暴雨的样本点做统计分析，对全部抽样数据做相关性检验，建立联合均值与方差模型、决策树模型和 MCP 变量选择模型，得出具体结论如下：
地形方面：（1）海拔和坡度与降水等级呈负相关，山影和坡向与降水等级呈正相关。（2）草地和灌木丛与降水等级呈负相关，耕地和林地与降水等级呈正相关。（3）针对特大暴雨数据分析，在暴雨发生区域，林地利用率是 0.797，然而耕地、草地、湿地和灌木丛利用率低，表明不是暴雨发生的主要地形区域。气候方面：（1）气温和降水量对降水等级呈正相关。暴雨天气的形成与特定的环境条件紧密相关，暴雨天气的形成与气温（超过 8°）和海拔（超过 1084 米）紧密相关，这些条件有利于降水形成。通过决策树模型可知经度大于等于 97.2°是划分暴雨和大暴雨的第一个分界点，其次是气温<20.6°、经度<109°和坡度大于等于0.112°。综合均值和方差模型的结果结论一致，即纬度、气温和土地利用率对极端天气形成有较大的影响。更详细的地形—气候相互作用结果见正文。
针对问题三，降雨、地形和土地利用对极端天气的影响因素和预测研究。首先，根据问题二的抽样数据将数据变量离散化后建立 XGBoost 模型，探究雨量等级的影响因素和不同变量间的交互效应。分别对降水量和气温使用高斯混合模型进行空间聚类，将空间数据划分成9份不同的区域，并根据第二问的结论对不同的空间数据进行打标签，得出2025—2035 年应对暴雨灾害最脆弱的地区。其次，通过经纬度把不同源数据的数据聚合到空间上，识别删除聚合后的缺失值，建立支持向量机和随机森林模型。最后，通过 ROC 曲线、混淆矩阵、精度和 AUC 值衡量不同模型的好坏，其中支持向量机的精度和 AUC 值为 0.993 和0.671；随机森林的精度和 AUC 值为 0.999 和 0.999。得出主要结论如下：（1）影响极端天气的前三个因素为气温、海拔和草地的利用率。（2）聚类模型表明东南沿海地区、长江中下游地区、西南地区附近在 2025—2035 年可能成为暴雨重灾区或暴雨高风险地区。（3）预测模型表明大珠江三角洲地区、西藏自治区和云南交界线、长江中下游（长沙、武汉和南昌）地区附近在 2025-2035 年可能会成为暴雨重灾区。针对问题四，中国土地利用变化结构特征建模分析。首先，读取数据 5（人口空间分布数据）和数据 6（GDP 空间分布数据）的年空间数据，关于人口数据选取 5000，3000，1000，300 进行截断，对 GDP 数据选取 100 进行截断。其次，对上述数据使用地图进行可视化，分别从三级阶梯、南北地区、东西地区进行分析。最后，对地形数据、土地利用类型数据和气候数据分别根据南北部和东西部的分界线进行聚合统计，使用假设检验（独立
t 检验）分析不同变量是否存在显著差异，得出如下结论：（1）东西部地区存在显著差异，降水量差幅为 1.604，气温差幅为 1.353。除此之外，草地、林地、灌木丛和湿地利用率均存在显著差异。耕地的利用率差幅中等。（2）南北方地区在降水量上差幅为 1.452，气温上差幅为 0.776。在耕地、林地、草地、灌木丛和湿地的利用率上均存在显著差异。
模型假设：
（1） 假设优化过程中各观测数据样本是独立的，不会相互影响。
（2） 假设人口密度越高，城市化越严重，会增加暴雨成灾的风险。
（3） 假设我们所使用的数据可靠，具有合理的分辨率和精度。
（4） 假设实验过程中不会有极端天气影响。
（5） 假设数据变化趋势不受政策影响。
（6） 假设我们所使用的数据具有代表性，模型经过训练后具有可推广性。
问题分析：
针对问题一，对降水量和土地利用/土地覆被两个变量分别构建描述性统计分析。首先根据给定的中国大陆 0.25°逐日降水数据集和土地利用和覆盖变化数据集，提取出 1990 年-2020 年间的每日降水量和土地利用/土地覆被数据依次按照月和年进行数据聚合，定义出13 个指标，分时间和空间两个维度进行可视化。
针对问题二，极端天气中地形-气候影响作用的建模分析。本题主要以暴雨为代表的极端天气事件，根据题目给出的数据集，选择前 4 个数据集（即中国数字高程图、中国 1°近地表气温数据集、中国大陆 0.25°逐日降水数据集和中国 0.5°土地利用和覆盖变化数据集）。在提取数据之前需对所有图片文件的像素精度进行统一。统一像素精度之后，以空间维度整合地形数据，以时间维度整合气候数据，将时空数据合并在一起。结合国家发布的降水量划分标准，给每条数据打上雨量等级的标签，若出现不平衡数据，则需进行重采样。最后，分别对抽样的地形—气候数据建立联合均值和方差模型以及决策树模型，研究地形和气候的变量对极端天气产生的影响大小。
针对问题三，地形-气候对极端天气的影响和交互作用建模分析与预测。直接使用问题二的数据库，建立混合高斯模型对各地区的降水量进行聚类，对全国地区进行划分，尤其是划分出重灾区，对这部分数据进行重点分析。随后，将各连续变量数据进行划分，将地形—气候数据转换为分类数据，再建立 XGBoost 模型，使用 Shap 方法细致化研究地形—气候变量对灾难的影响以及各变量间的交互作用。最后，建立聚类模型，根据以上的分析结果建立模型预测各地区发生灾难的概率，并对结果进行可视化。
针对问题四，中国土地利用变化的特征与结构分析。本题使用的数据在问题二和问题三的基础上，加入中国大陆 1km 逐年历史人口空间分布公里网格数据集和中国大陆 1km逐年历史 GDP 空间分布公里网格数据集。提取人口空间分布数据和历史 GDP 空间分布数据之后，需对两个原始数据集进行描述性统计分析以及可视化。随后，结合前三个问题得到的关于地形、气候和土地利用类型对极端天气的影响结论，根据南北结构，东西差异，三级阶梯进行数据整合，并进行假设检验。
模型的建立与求解整体论文缩略图
全部论文请见下方“ 只会建模 QQ名片” 点击QQ名片即可
程序代码：
部分程序代码如下:
install.packages("lattice")
install.packages("versions")
install.packages("rasterVis")
# 1. library packages
library(raster)
library(ncdf4)
library(lubridate)
library(data.table)
library(ggplot2)
library(patchwork)
library(rasterVis)
# 1.1 Definition of rainfall metrics
path1 = ".\\data\\"
Names1 <- list.files(path = path1, "Data3")
Data3 <- stack(x = paste0(path1, Names1))
x <- names(Data3[[1]])
NameToDate <- function(x){ x <- strsplit(x, split = "\\.")[[1]] x <- c(unlist(strsplit(x[1], "X"))[[2]], x[2:3]) x <- as.Date(paste0(x, collapse = "/")) return(x)
}
Names <- names(Data3)
Namesi <- c()
for(i in 1:22645){ Namesi[i] <- NameToDate(Names[i]) cat("i = ,", i, "\n")
}
Namesi <- as.Date("1970/01/01") + Namesi
Start <- which(Namesi == as.Date("1990/1/1"))
End <- which(Namesi == as.Date("2020/12/31"))
Namesi <- Namesi[Start:End]
Data3 <- Data3[[Start:End]]
Ye <- 1990:2020
Year <- 1990
Day3 <- c()
Dam3 <- c()
DancM <- Data3[[1:length(Ye)]]
DancS <- Data3[[1:length(Ye)]]
DancSD <- Data3[[1:length(Ye)]]
names(DancM) <- paste0("Year", Ye)
names(DancSD) <- paste0("Year", Ye)
n <- length(values(Data3[[3]]))
for(Year in Ye){ da3 <- c() Pos <- which(year(Namesi) == Year) for(pos in Pos){ Da <- Data3[[pos]] da3 <- rbind(da3, data.table(Pos = 1:n, Time = Namesi[pos], rainfall = values(Da))) cat("Year = ", Year, "pos = ,", pos, "\n") } DancD <- da3[, .(M = mean(rainfall), S = sum(rainfall), SD = sd(rainfall)), .(Pos)] da3 <- da3[is.na(rainfall) != T,] Day3 <- rbind(Day3, da3[,.(RYm = mean(rainfall), RYs = sum(rainfall), RYsd = sd(rainfall)), by = .(year(Time))]) Dam3 <- rbind(Dam3, da3[,.(Rmm = mean(rainfall, na.rm = T), Rms = sum(rainfall, na.rm = T), Rmss = sd(rainfall, na.rm = T)), by = .(month(Time))]) values(DancM[[which(Year == Ye)]]) <- DancD$M values(DancS[[which(Year == Ye)]]) <- DancD$S values(DancSD[[which(Year == Ye)]]) <- DancD$SD
}
write.csv(Day3, "Day3.csv")
write.csv(Dam3, "Dam3.csv")
p <- ggplot(data = Day3, aes(x = year , y = RYm, col = "red"))+
guides(col = FALSE)+ theme_bw()
P1 <- p + geom_point(aes(x = year , y = RYm, col = "red"),size=5)+ geom_errorbar(aes(x = year,ymin = RYm - RYsd, ymax = RYm + RYsd), width = 0.6, size = 2, col = "red" )
P2 <- p + geom_point(size=5) + geom_line()
P3 <- p + geom_point(aes(y = RYsd), size=5) + geom_line(aes(y = RYsd)) + ylab("RYsd")
P1 / P2 + P3
Dam3 <- cbind(year = rep(Ye, each = 12), Dam3)
p <- ggplot(data = Dam3, aes(x = year , y = Rmm, col = factor(month)))+
guides(col = FALSE)+ theme_bw()
P1 <- p + geom_point(size=5)+ geom_errorbar(aes(x = year,ymin = Rmm - Rmss, ymax = Rmm + Rmss), width = 0.6, size = 2, col = "red" ) +
facet_wrap(month~., 3, 4, scales="fixed")
P2 <- p + geom_point(size=5)+
facet_wrap(month ~ ., 3, 4, scales="free") + geom_smooth()
P2 <- p + geom_point(size=5) + geom_line()
Ts <- ts(data = Dam3$Rmm, start = c(1990, 1), frequency = 12)
Res <- decompose(Ts, type = "additive")
plot(Res)
# 没有划分等级
i = 1
DancM1 <- DancM
for(i in 1:31){ Tem <- values(DancM1[[i]]) values(DancM1[[i]]) <- (Tem > 1)+ (Tem > 3) + (Tem > 5) + (Tem > 7) + (Tem > 9)
}
plot(DancM1[[1:16]])
plot(DancM1[[17:31]])
DancSD1 <- DancSD
for(i in 1:31){ Tem <- values(DancSD1[[i]]) values(DancSD1[[i]]) <- (Tem > 5)+ (Tem > 10) + (Tem > 15) + (Tem > 20)
}
plot(DancSD1[[1:16]])
plot(DancSD1[[17:31]])
YeSelec <- seq(1, 31, by = 6)
plot(DancM1[[YeSelec]])
plot(DancSD1[[YeSelec]])
DancM2 <- DancM
Vals <- c()
for(i in 1:31){ Tem <- values(DancM2[[i]]) Vals <- cbind(Vals, Tem)
}
DancM2 <- DancM2[[1]]
values(DancM2) <- apply(Vals, 1 ,sd)
plot(DancM2)
test <- DancM[[1]]
test_spdf <- as(test, "SpatialPixelsDataFrame")
test_df <- as.data.table(test_spdf)
# 赋值列名
colnames(test_df) <- c("value", "x", "y")
#开始绘图
land_use <- ggplot() + geom_raster(data = test_df , aes(x = x, y = y,fill = factor(value))) + scale_fill_manual(name = "Category",values = classcolor,labels = landclass) + labs(x ='Longitude(E)',y="Latitude(N)", title = "Raster data Charts Exercise", subtitle = "land use data", caption = 'Visualization by DataCharm')+theme_bw()+ theme(text = element_text(family = "Times_New_Roman",face='bold'), title = element_text(family = 'Times_New_Roman',size = 16,face = 'bold'), axis.text = element_text(family = 'Times_New_Roman',size = 12,face = 'bold'), #修改刻度线内axis.ticks.length=unit(-0.22, "cm"), #设置刻度 label 的边距axis.text.x = element_text(margin=unit(c(0.5,0.5,0.5,0.5), "cm")), axis.text.y = element_text(margin=unit(c(0.5,0.5,0.5,0.5), "cm")) )
land_use