SPSS筛选因子之后如何采用python建模和反演整个研究区?（以反演生物量为例）

引言

在遥感数据分析中，因子筛选和建模是关键步骤。筛选出与目标变量（如生物量）显著相关的因子，不仅可以提高模型的预测精度，还能简化模型结构。

通常我们使用SPSS软件进行因子的初步筛选，通过相关性分析选出与目标变量相关性较高的因子。在SPSS中，我们可以使用皮尔逊相关系数来评估因子的显著性，通常选择带有两颗星号（表示相关性显著）的变量用于后续建模。

本教程将以反演森林生物量为例，详细介绍如何使用Python对筛选后的因子进行建模，并应用于整个研究区的反演任务。这一流程具有较强的通用性，适用于其他类似的遥感反演任务。

任务描述

本文以建模和反演森林生物量为例。我们将使用随机森林算法构建预测模型，并将其应用于整个研究区的生物量反演。

本文所用的示例数据和最终的反演结果已上传至百度网盘，整理不易，如有需要公众号（Python与遥感）后台回复 49 获取。

通过本文，读者可以完整地掌握从因子筛选到建模、再到研究区生物量反演的整个流程。

文件描述

1. 随机森林反演.ipynb：Jupyter Notebook代码文件，用于实现建模和反演。
2. jianmo.csv：用于建模的数据文件，包含反演变量（生物量）和5个最重要的影响因子。
3. 5个.tif文件：建模和反演所用的栅格变量数据。
4. 生物量.tif：反演得到的研究区生物量结果，并已导出为png方便查看。

完成反演后，我们将得到整个研究区的生物量分布图。下图展示了根据森林区域掩膜文件提取后的研究区生物量分布结果。

代码文件截图

如何打开代码运行，请看这篇文章

代码详解

1. 导入所需的库

首先，导入进行数据处理、建模和可视化所需的库：

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
import scipy.stats as stats
from sklearn import metrics
from sklearn.model_selection import GridSearchCV, RandomizedSearchCV
from joblib import dump

• pandas和numpy：数据处理和数值计算。
• matplotlib：绘图和可视化。
• sklearn：机器学习建模、模型评估和超参数调整。
• scipy：统计计算。

2. 读取数据并进行预处理

我们从CSV文件中读取数据，并将其分为自变量（影响因子）和因变量（生物量）：

data = pd.read_csv(r'.\jianmo.csv', na_values=[' '])
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=44)

• pd.read_csv读取CSV文件，na_values=[' ']将空白值视为缺失值。
• y为因变量（生物量），X为自变量（影响因子）。
• train_test_split将数据分为训练集和测试集，测试集占比30%。

3. 构建随机森林模型并进行训练

使用随机森林算法进行建模：

regressor = RandomForestRegressor(n_estimators=500, max_features=5)
regressor.fit(X_train, y_train)

• RandomForestRegressor构建随机森林回归模型，设置n_estimators=500构建500棵树，max_features=5表示每次分裂时最多考虑5个特征。
• fit函数训练模型。

4. 误差分析

接下来，对模型的预测误差进行分析，并计算Pearson相关系数、R²和均方根误差（RMSE）：

random_forest_predict = regressor.predict(X_test)
random_forest_error = random_forest_predict - y_testplt.figure(1)
plt.clf()
plt.hist(random_forest_error)
plt.xlabel('Prediction Error')
plt.ylabel('Count')
plt.grid(False)
plt.show()random_forest_pearson_r = stats.pearsonr(y_test, random_forest_predict)
random_forest_R2 = metrics.r2_score(y_test, random_forest_predict)
random_forest_RMSE = metrics.mean_squared_error(y_test, random_forest_predict) ** 0.5
print('Pearson correlation coefficient = {0}、R2 = {1} and RMSE = {2}.'.format(random_forest_pearson_r[0], random_forest_R2, random_forest_RMSE))

• regressor.predict对测试集进行预测，并计算误差。
• 绘制误差分布图，帮助理解模型的误差情况。
• 计算Pearson相关系数、R²和RMSE，评估模型的预测性能。

5. 栅格数据的读取与处理

读取研究区的栅格变量数据，并整合为一个矩阵：

import rasterio
from tqdm import tqdmwith rasterio.open(r'.\podu.tif') as src:data1 = src.read(1)meta = src.meta
with rasterio.open(r'.\dem.tif') as src:data2 = src.read(1)
with rasterio.open(r'.\lai.tif') as src:data3 = src.read(1)
with rasterio.open(r'.\ireci.tif') as src:data4 = src.read(1)
with rasterio.open(r'.\mcari.tif') as src:data5 = src.read(1)X = np.stack((data1, data2, data3, data4, data5), axis=-1)

• 使用rasterio库读取栅格文件，将每个变量的数据存储在data1至data5中。
• meta保存栅格文件的元数据，用于保存预测结果。
• np.stack将所有栅格数据整合为一个三维矩阵。

6. 数据清洗

对栅格数据进行清洗，去除缺失值（NaN）：

X_2d = X.reshape(-1, X.shape[-1])
print(np.isnan(X_2d).any())  # 检查是否存在 NaN 值X_2d[np.isnan(X_2d)] = 0  # 将 NaN 值替换为0

• 将三维数据转换为二维数据，便于处理。
• 检查并处理缺失值，确保数据完整性。

7. 使用模型进行预测

利用训练好的随机森林模型对栅格数据进行生物量预测：

y_pred = []
for i in tqdm(range(0, X_2d.shape[0], 10000)):y_pred_chunk = regressor.predict(X_2d[i:i+10000])y_pred.append(y_pred_chunk)
y_pred = np.concatenate(y_pred)

• 将栅格数据分块输入模型预测，减少内存占用。
• 使用tqdm显示进度条，实时监控预测进度。

8. 保存预测结果

将预测得到的生物量数据保存为新的栅格文件：

with rasterio.open(r'.\生物量.tif', 'w', **meta) as dst:dst.write(y_pred.reshape(X.shape[:-1]), 1)
print("预测结束")

• 创建新的栅格文件，并将预测结果写入其中。
• 确保预测结果与原始栅格的空间信息一致。

9. 结果展示

最后，我们得到了整个研究区的生物量反演图。这一过程不仅适用于生物量的反演，也适用于其他多光谱遥感数据的反演任务。

总结

通过这篇文章，我们展示了如何使用Python对筛选后的影响因子进行随机森林建模，并应用于遥感数据反演森林生物量的任务。该流程不仅可以提高预测精度，还能有效处理大规模的遥感数据，为生态监测和森林管理提供重要支持。

如有更多问题或需要进一步指导，欢迎公众号（Python与遥感）后台联系作者获取更多信息。