地质灾害（滑坡、崩塌等均可）易发性评价

一、前言

仅代表个人见解，如有错误，纯属我个人问题！！！

如有不对，还请不吝赐教！

方法

💡 Tips：目前由于个人技术受限+见识浅薄，个人操作下来机器学习方法性能效果好一点（因为样本<灾害点>数量少导致）

• 常见的机器学习方法-->针对灾害点点（即样本大小为1*1）上的数据进行训练预测
• 深度学习方法-->针对灾害点+周边(即比如说7*7的样本大小)环境进行预测

数据

• 灾害点数据
• 各类因子数据（如DEM、坡度、坡向等等） 不排除因子越多越好，当然有的有因子加了也可能会导致准确度变差

流程

---

下面进行详细讲解

二、具体步骤

1.数据准备

1. 收集各种数据 无需多言
2. 将以上收集到的各种数据，统一分辨率，如 30 m * 30 m ， 90 m * 90 m ，统一坐标系 注：矢量数据也要转成栅格数据
3. 灾害点数据，暂时是不用转成栅格

2.数据预处理

此块是将 前面前面准备好的各种因子数据，由于 它们的 单位 什么的 还有就是属性值什么的都不一样，因此需要将数据进行转换

2.1 影响因子数据处理

！！！切记最后的一定影响因子是只有里面的数据不一样，其它均一样--地理参考、分辨率，tif的大小等

转换的常见方法有 直接重分类、频率比、信息量等

频率比 和 信息量 都是在 重分类 的基础之上进行的

重分类的目的是将拥有 某一条件 的地物分成一类 如：将高程0~200（注：我随意选的规则，打胡乱说的）的划分为平原

以下进行重分类展示，频率比和信息量 建议百度 B站上也有视频 可直接搜 易发性评价

将前期准备的因子数据按照一定规则进行重分类

如以下，都是我自己找的别人论文的重分类的规则

重分类完一般是将分类的结果值命名为数字 如 1、2、3、... 如0~200为1、200~400为2

（但是出图的时候记得修改标注 如将1改回为0~200）

基于以上步骤你就可以得到 各个因子数据 因子的数据均是由 1 2 3 4 5 .... 构成

信息量其实就是和在不同分区中，和灾害点建立关系。

2.2 灾害点数据处理

由于机器学习/深度学习均需要正、负样本进行训练 而这里我们只有灾害点的数据（即为正样本） 无负样

现在绝大多数论文所采用的方法 就是 随机采样负样本

具体方法为：

在灾害点的附近生成一点范围大小（如1km或者500m 视情况而定）的缓冲区，然后在缓冲区外&研究区内 生成相同灾害点数量的随机点 将这些随机点 命为 负样本点 相同灾害点数量是因为要照顾 正负样本均衡

而后将灾害点的属性赋值为 1 ，非灾害点（负样本点）赋值为 0

2.3 多值提取至点

利用arcgis 多值提取至点工具 ，将因子数据提取到灾害点和非灾害点的属性表上

后续如果是机器学习的话，用属性表里面的excel进行训练就行了

shp的属性表转excel 可以用表转EXCEL工具实现

3 影响因子筛选

因此有的因子会有共线性会影响结果，用共线性分析以下更为严谨

相关性分析

共线性分析

将前后保存的excel 数据表的表头如下：

以下代码用的皮尔逊相关性分析、和VIF进行评价，一般来说VIF>10的值就要剔除，本文这里给定的阈值是5，请注意修改，这里实现了VIF值大于5则打印应删除的因子

import numpy as np
import seaborn as sns
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.stats.outliers_influence import variance_inflation_factor
from statsmodels.tools.tools import add_constant
import re


# 计算VIF的函数
def calculate_vifs(dataframe):
    X = add_constant(dataframe)
    return pd.Series([variance_inflation_factor(X.values, i)
                      for i in range(X.shape[1])], index=X.columns)


# 移除高VIF值的变量的函数
def remove_high_vif_factors(dataframe, vif_threshold=5.0):
    removed_factors = []  # 用于记录被移除的变量
    while True:
        VIFs = calculate_vifs(dataframe).drop('const')
        max_vif = VIFs.max()
        if max_vif <= vif_threshold:
            break
        factor_to_drop = VIFs.idxmax()
        removed_factors.append(factor_to_drop)  # 添加被移除的变量到列表
        dataframe = dataframe.drop(columns=[factor_to_drop])
        print(f"Removed {factor_to_drop} with VIF={max_vif:.2f}")  # 打印被移除的变量
    return dataframe, removed_factors


# 读取数据
df = pd.read_csv(r"./Landslide-susceptibility-mapping/数据集16个因子.csv")

factors = df[list(df.columns)[:-1]]  # 选择您的自变量列
print(factors.columns)

# 移除高VIF值的变量，并获取被移除变量的列表
factors_with_low_vif, removed_factors_list = remove_high_vif_factors(factors, vif_threshold=7.0)

# 打印被移除的变量的列表
print("\n打印被移除的变量列表：")
print(removed_factors_list)

# 计算最终的VIF和TOL值
final_vifs = calculate_vifs(factors_with_low_vif).drop('const')
final_tols = 1 / final_vifs

# 创建结果表
results = pd.DataFrame({
    'Factor': final_vifs.index,
    'VIF': final_vifs.values,
    'TOL': final_tols.values
})

# 打印结果表
print("VIF和TOL结果表：")
print(results)
results.to_excel('.\结果分析存储\VIF和TOL结果表.xlsx')


# 计算皮尔逊相关性矩阵
pearson_correlation_matrix = factors_with_low_vif.corr()
print("\n皮尔逊相关性矩阵：")
print(pearson_correlation_matrix)

# 如果需要将皮尔逊相关性矩阵保存到Excel文件
pearson_correlation_matrix.to_excel('.\结果分析存储\pearson_correlation_matrix.xlsx')

4 机器学习

就是将前期的excel，进行训练，预测

def parse_args():
    parser = argparse.ArgumentParser(description="Train ML models Processes on data")
    parser.add_argument("--feature_path", default='origin_data2/feature/', type=str)
    parser.add_argument("--label_path", default='origin_data2/label/label1.tif', type=str)
    parser.add_argument("--model", default='ModelRF', type=str)
    args = parser.parse_args()
    return args

def main():
    data = pd.read_csv("数据集14个因子.csv")
    X, y = data.iloc[:, :-1], data.iloc[:, -1]
    y = np.array(y).ravel()
    y = pd.DataFrame(y, columns=['label'])

    x_train, x_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=42)  # 70%训练，30%测试
    print("x_train，y_train，x_test，y_test：", x_train.shape, y_train.shape, x_test.shape, y_test.shape)

    #创建模型
    models = {
        'LR': LogisticRegression(),
        "RF": RandomForestClassifier(),
        "ET": ExtraTreesClassifier(),
        "KNN": KNeighborsClassifier(n_neighbors=2),
        "SVM": SVC(probability=True),
        "DC": DecisionTreeClassifier(),
        "CatBoost": CatBoostClassifier(verbose=0),
        "LGBM": lgb.LGBMClassifier(verbose=-1),
    }

    # 格式化表头
    print(
        "{:<20} | {:<9} | {:<9} | {:<9} | {:<9} | {:<9}".format("Model", "Accuracy", "Precision", "Recall", "F1 Score",
                                                                "AUC"))
    results_df = pd.DataFrame(columns=["Model", "Accuracy", "Precision", "Recall", "F1 Score", "AUC"])

    for name, model in models.items():
        model.fit(x_train, y_train) # 模型训练
        y_pred = model.predict(x_test)
        y_pred_prob = model.predict_proba(x_test)[:, 1]  # 一列是类别为 0 的概率，另一列是类别为 1 的概率
        # print('Model Name: ', name)
        # tp, fn, fp, tn = confusion_matrix(y_test,y_pred,labels=[1,0]).reshape(-1)
        # # accuracy = (tp + tn) / (tp + tn + fp + fn)
        # # precision = tp / (tp + fp)
        # # recall = tp / (tp + fn)
        # # f1 = 2 * (precision * recall) / (precision + recall)
        #
        # accuracy = round((tp+tn)/(tp+fp+tn+fn), 3)
        # print('Accuracy :', round(accuracy*100, 2),'%')

        accuracy = accuracy_score(y_test, y_pred)
        precision, recall, f1_score, _ = precision_recall_fscore_support(y_test, y_pred, average='binary')
        auc1 = roc_auc_score(y_test, y_pred)
        fpr, tpr, thresholds = roc_curve(y_test, y_pred_prob)
        auc2 = auc(fpr, tpr)
        print("{:<20} | {:<9.2%} | {:<9.2%} | {:<9.2%} | {:<9.2%} | {:<9.2%}".format(name, accuracy, precision, recall,
                                                                                     f1_score, auc2))
        # # 将结果添加到DataFrame
        # results_df = results_df.append(
        #     {"Model": name, "Accuracy": accuracy, "Precision": precision, "Recall": recall, "F1 Score": f1_score},
        #     ignore_index=True)

    drawAUC_TwoClass2(models, x_test, y_test)
    # predict_tif3(models, X.columns) # 预测并保存tif文件

# 用窗口预测
def predict_tif3(models, feature_names):
    args = parse_args()
    # os.listdir(args.feature_path)[0]
    # 假设config['data_path']包含所有tif文件的路径
    # 预先定义一个变量来存储nodata值
    nodata_value = None

    # 以第一个tif文件为基准，获取元数据和窗口大小
    with rasterio.open(os.path.join(args.feature_path, os.listdir(args.feature_path)[0])) as src:
        meta = src.meta
        if nodata_value is None:
            nodata_value = src.nodatavals[0]
        # 窗口大小，这里可以根据您的内存大小进行调整
        width, height = src.width // 2, src.height // 2

    # 更新元数据
    meta.update(count=1, dtype=rasterio.float32, nodata=nodata_value)

    # 循环执行预测并保存tif文件
    for name, model in models.items():
        start_time = time.time()  # Start timing
        # 定义结果数组，用nodata值初始化
        prediction_image = np.full((meta['height'], meta['width']), nodata_value, dtype=rasterio.float32)

        # 按窗口读取和预测
        for row_off in range(0, meta['height'], height):
            for col_off in range(0, meta['width'], width):
                window = Window(col_off, row_off, width, height)
                arrays = []
                nodata_mask = None
                for path in os.listdir(args.feature_path):
                    # num_class = float(re.search(r'_(\d+)\.', path).group(1))   # 断层_5.tif --> 5
                    path = os.path.join(args.feature_path, path)
                    with rasterio.open(path) as src:
                        band = src.read(1, window=window).astype('float32')  # / num_class

                        # 初始化nodata掩膜
                        if nodata_mask is None:
                            nodata_mask = band == nodata_value
                        # 掩膜无效数据
                        band[nodata_mask] = np.nan
                        arrays.append(band)

                # 将窗口内的所有数组合并成一个3D数组
                stacked_arrays = np.dstack(arrays)
                # 提取有效像素
                valid_mask = ~np.isnan(stacked_arrays).all(axis=2)
                valid_pixels_features = stacked_arrays[valid_mask]

                # 如果窗口内没有有效像素，则跳过
                if valid_pixels_features.size == 0:
                    continue

                # 使用特征名称创建DataFrame
                valid_pixels_features_df = pd.DataFrame(valid_pixels_features, columns=feature_names)

                # print(valid_pixels_features_df.min())

                # 对有效像素特征进行预测
                predictions = model.predict_proba(valid_pixels_features_df)[:, 1]

                # 将预测填充到相应位置
                prediction_image[row_off:row_off + height, col_off:col_off + width][valid_mask] = predictions

        # 保存预测图像为新的tif文件
        save_path = f'./易发性图/{name}_predict.tif'
        with rasterio.open(save_path, 'w', **meta) as dst:
            dst.write(prediction_image, 1)
        print(name, "保存成功！ ", save_path)

        # Calculate and print the total time spent
        end_time = time.time()
        print(f"time for {name}: {end_time - start_time:.2f} seconds")

        # 清理变量
        prediction_image = None


if __name__ == '__main__':
    main()

5 深度学习

代码的话参考如下：

https://github.com/scar-on/Landslide-susceptibility-mapping

此项目已经把深度学习写得比较完善了，有能力的人，还可在模型上修改。

另：Landslide-susceptibility-mapping/train_CNN.py文件中 第103行 scheduler.step()，应该再往里缩进一下

---

三、结语

到这里大概就结束了

后续分享一些小工具，感谢他们无私分享，向他们致敬！！！

SVM-LSM-Toolbox 这个确实方便了我很多，可以做了很多的预处理工作，DEM生成其它相关因子（坡度、坡向、地形起伏度等，一键生成），还有就是非灾害点生成那里其实蛮麻烦的，先要弄缓冲区，然后再裁剪缓冲区，最后生成随机点，这个工具一下就使用好了，当然你也可以使用他这个工具去进行SVM预测。

注：我记得这个工具好像是要用UTM投影，记得处理一下。

中文简介：

2022（一等奖）B360基于Arcpy与SVM的滑坡易发性评价_svm-lsm工具箱-CSDN博客

详细操作视频：

2022（一等奖）B360基于Arcpy与SVM的滑坡易发性评价_哔哩哔哩_bilibili

工具样式如下：

另外就是B站上由很多关于易发性评价的视频，有用信息量、频率比的，等等等等

如

ArcGIS滑坡易发性评价_哔哩哔哩_bilibili

建议阅读论文：

[1]孔嘉旭,庄建琦,彭建兵,等.基于信息量和卷积神经网络的黄土高原滑坡易发性评价[J].地球科学,2023,48(05):1711-1729.

中国知网-登录