高斯过程回归的空间插值

Question

我有一个包含 140.000 点（行）的 csv 文件。它包括：

1. 经度值
2. 纬度值
3. 特定点的沉降值。我假设这些点在空间上是相关的。

我想对点的区域进行空间插值分析。意思是，我将使用例如克里金法（即高斯过程回归）进行地统计插值分析。

我正在 sci-kit learn 页面上阅读有关高斯回归的信息。但我不确定如何实现它。

哪些特性决定了我可以使用哪个内核？如何正确使用我的空间数据实现这一点？

Answer 1

首先，您应该将数据转换为投影坐标系。最好的取决于您的数据所在的位置；本质上，您希望在您的位置（例如赤道附近的墨卡托，或者如果您的数据都接近单个子午线，则为横轴墨卡托。您可以在geopandas 中实现这一点，例如：

import pandas as pd
import geopandas as gpd

data = {'latitude': [54, 56, 58], 'longitude': [-62, -63, -64], 'subsidence': [10, 20, 30]}
df = pd.DataFrame(data)

params ={
    'geometry': gpd.points_from_xy(df.longitude, df.latitude),
    'crs': 'epsg:4326',  # WGS84
}

gdf_ = gpd.GeoDataFrame(df, **params)
gdf = gdf_.to_crs('epsg:2961')   # UTM20N
gdf

此GeoDataFrame 现在位于投影坐标中。现在您可以进行一些空间预测：

import numpy as np
from sklearn.gaussian_process.kernels import RBF
from sklearn.gaussian_process import GaussianProcessRegressor

kernel = RBF(length_scale=100_000)
gpr = GaussianProcessRegressor(kernel=kernel)

X = np.array([gdf.geometry.x, gdf.geometry.y]).T
y = gdf.subsidence
gpr.fit(X, y)

现在您可以在某个位置进行预测，例如gpr.predict([(500_000, 5_900_000)]) 为我的玩具数据提供 array([22.86764555])。

要在网格上进行预测，您可以这样做：

x_min, x_max = np.min(gdf.geometry.x) - 10_000, np.max(gdf.geometry.x) + 10_000
y_min, y_max = np.min(gdf.geometry.y) - 10_000, np.max(gdf.geometry.y) + 10_000
grid_y, grid_x = np.mgrid[y_min:y_max:10_000, x_min:x_max:10_000]
X_grid = np.stack([grid_x.ravel(), grid_y.ravel()]).T
y_grid = gpr.predict(X_grid).reshape(grid_x.shape)

需要考虑的事情：

• 您应该阅读geopandas 和sklearn.gaussian_process 的文档
• 您应该使内核适合您的数据。
• 您可能想要使用各向异性内核。
• 估计器有几个超参数，您应该注意这些参数。
• 不要忘记对您的估计进行一些验证，检查残差的分布等。
• 您可能需要使用专业的地理统计软件包，例如 gstools，它会为您完成很多繁琐的工作。