使用随机森林回归预测森林生物量（python）

随机森林算法的基本思想是基于多颗决策树的集成学习过程，使用场景广泛，一般的分类回归问题都可以使用。我们以光学影像为例，来估测森林生物量。

建立回归关系需要满足的条件

1、线性关系：回归关系应该是线性的，即自变量和因变量之间的关系应该是线性的。
2、独立性：自变量之间应该是独立的，即自变量之间的相关性应该尽可能小，自变量之间不应该存在多重共线性。
3、正态性：残差应该是正态分布的，即残差应该符合正态分布的假设。
4、同方差性：残差的方差应该是恒定的，即残差的方差应该在自变量的不同取值下是相同的。
5、随机性：误差项应该是随机的，即误差项应该是不可预测的，不能被自变量解释。

使用传统机器学习方法估测森林生物量的过程

基本思想：由于遥感影像一般可以覆盖整个研究区，而地面调查工作一般不会对整个研究区进行展开，所以就可以利用遥感影像和地面调查数据建立回归关系，然后再利用所建立的这个关系对未进行地面调查的地方进行预测。

指标筛选：建立回归关系首先要满足指标变量和预测变量之间存在相关性。先从光学影像中提取指标变量，单波段或者波段组合运算啥的，可以参考相关论文来选。

指标提取：在同一张影像中提取的自变量之间难免存在共线性问题，但也不是不能减少其对模型的影响，例如：使用PCA对指标进行降维；对指标进行逐步回归或者相关性分析提取对生物量解释程度很高的几个变量作为自变量；在损失函数中添加惩罚项来降低模型复杂度，从而缓解指标共线性问题；使用树模型，例如典型的随机森林，由于随机森林在建立决策树过程是随机选取指标以及特征，这种随机性可以使得模型对于共线性较弱的特征更加敏感，从而避免了共线性的影响；随机森林模型在建立过程中是通过计算特征调用率，建立分裂树信息熵减少量等来判断特征对于预测指标的重要性，可以剔除一些重要性较低的指标从而提高模型的泛化能力。

建立随机森林回归

提取样地点对应指标变量的像元值作为自变量，样地点生物量作为因变量后就可以建立随机森林模型并对整个研究区生物量就行预测，具体实现过程如下：

数据准备

import pandas as pd
import numpy as np
from matplotlib import pyplot as plt
from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error, r2_score
from pyswarm import pso
import rasterio

# 读取Excel表格数据
data = pd.read_excel(r'D:\code\3s\redata.xlsx', sheet_name='Sheet3')
y = data.iloc[:, 0].values  # 生物量
X = data.iloc[:, 1:].values  # 指标变量
# 分割数据集为训练集和测试集
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=42)

定义随机森林模型

def random_forest_regression(X, y, n_estimators, max_depth, max_features, min_samples_split, min_samples_leaf):
    model = RandomForestRegressor(n_estimators=n_estimators, max_depth=max_depth, max_features=max_features,
                                  min_samples_split=min_samples_split, min_samples_leaf=min_samples_leaf)
    model.fit(X, y)
    return model

定义粒子群算法优化函数

def objective_function(params):
    n_estimators = int(params[0])
    max_depth = int(params[1])
    max_features = int(params[2])
    min_samples_split = int(params[3])
    min_samples_leaf = int(params[4])
    model = random_forest_regression(X_train, y_train, n_estimators, max_depth, max_features, min_samples_split,
                                     min_samples_leaf)
    y_pred = model.predict(X_test)
    mse = mean_squared_error(y_test, y_pred)
    return mse
   

使用粒子群算法优化随机森林模型

lb = [50, 2, 1, 2, 1]  # 参数下限
ub = [200, 10, len(X[0]), 10, 10]  # 参数上限
bounds = [(low, high) for low, high in zip(lb, ub)]
xopt, fopt = pso(objective_function, lb, ub, swarmsize=20, maxiter=50)

训练最优模型

当然，没有一个模型或者优化算法是万能的，建立的模型要针对自己的数据进行调整，包括0.3这种验证集比例都会影响模型结果。而我们这里只是用了粒子群算法对超参数进行了简单优化，也没有做交叉验证来测试模型性能。需要对自己的数据结构以及模型原理有更深入的理解，再加上不断地调试模型，才能得到一个比较理想的结果。

n_estimators = int(xopt[0])
max_depth = int(xopt[1])
max_features = int(xopt[2])
min_samples_split = int(xopt[3])
min_samples_leaf = int(xopt[4])
rf_model = random_forest_regression(X_train, y_train, n_estimators, max_depth, max_features, min_samples_split,
                                    min_samples_leaf)

对模型进行精度验证

y_pred = rf_model.predict(X_test)
mse = mean_squared_error(y_test, y_pred)
rmse = np.sqrt(mse)
r2 = r2_score(y_test, y_pred)
print('最优模型参数: n_estimators={}, max_depth={}, max_features={}, min_samples_split={}, min_samples_leaf={}'.format(
    n_estimators, max_depth, max_features, min_samples_split, min_samples_leaf))
print('RMSE={:.2f}'.format(rmse))
print('R2分数={:.2f}'.format(r2))

预测结果可视化

#定义折线图绘制函数，只取前一百条可视化，因为数据太多了
def linechart_plot(TureValues,PredictValues):
    plt.plot(np.arange(100), TureValues[:100], "go-", label="True value")
    plt.plot(np.arange(100), PredictValues[:100], "ro-", label="Predict value")
    plt.title(f"RandomForest--score:{r2}")
    plt.legend(loc="best")
    plt.show()
#定义散点图函数
def scatter_plot(TureValues,PredictValues):
    fig, ax = plt.subplots()
    ax.scatter(TureValues, PredictValues,s=20 , c='r' , edgecolors='k' , marker='o' , alpha=0.8)
    ax.plot([y_test.min(), y_test.max()], [y_test.min(), y_test.max()], c='0' , linewidth=1 , linestyle=':' , marker='.' , alpha=0.3)
    plt.title('RandomForestRegressionScatterPlot')
    plt.show()
    
#调用函数进行绘图
linechart_plot(y_test,y_pred)
scatter_plot(y_test,y_pred)

至此，随机森林模型已经建立好了，接着就是对没有进行实地调查的位置进行生物量预测了，如果建立的是参数回归模型，则可以根据回归函数直接进行栅格运算得到预测的生物量，而随机森林模型数据非参数模型，建立的是自变量与因变量的回归关系，没有具体回归函数，但是可以使用python完成这个任务。

单波段栅格数据的数值实质上就是一个二维矩阵，使用rasterio和numpy库读取栅格数据存为二维数组，这样就可以拿到所有指标变量的数据了，再调用训练好的模型就可以计算出来预测生物量的数值，最后在使用rasterio将计算出来的预测数值写入栅格中，得到最终结果！！！
需要注意的是：所有的栅格数据必须保持形状一致，即每个栅格指标的行列像元数要相等，可以使用研究对象的shp文件掩膜掉不需要计算的部分，以减少计算冗余。最重要的是：预测过程中指标变量的排列顺序需要和训练模型过程中自变量的顺序保持一致，否则可能会导致预测结果出现偏差。

预测生物量

# 读取五个栅格指标变量数据并整合为一个矩阵
with rasterio.open(r'F:\raster\dem_h1.tif') as src:
    data1 = src.read(1)
    meta = src.meta
with rasterio.open(r'F:\raster\te_inter1.tif') as src:
    data2 = src.read(1)
with rasterio.open(r'F:\raster\te_best1.tif') as src:
    data3 = src.read(1)
with rasterio.open(r'F:\raster\n_ca_photo.tif') as src:
    data4 = src.read(1)
with rasterio.open(r'D:\code\3s\30m\n_ca_photo_tif.tif') as src:
    data5 = src.read(1)

X = np.stack((data1, data2, data3, data4,data5), axis=-1)
rows, cols, bands = X.shape

# 使用已经训练好的随机森林模型对整合后的栅格指标变量数据进行生物量的预测
X_2d = X.reshape(rows * cols, bands)
y_pred = rf_model.predict(X_2d)

# 将预测结果保存为一个新的栅格数据
with rasterio.open(r'F:\raster\result\biomass.tif', 'w', **meta) as dst:
    dst.write(y_pred.reshape(rows, cols), 1)
    
print("预测结束")

结束语

在使用模型进行生物量或者蓄积量预测等作业中，其实最影响预测结果精度的是样本数据的精度，分类问题也同理。对于模型的精度其实影响不是那么高，只有底子好了，再加上一点点适当的方法，做的成果才会好。但是如果你有用不完的数据可以用来进行高强度深度学习训练，那当我没说。