本文由SCY翻译，转载注明出处。

本章介绍了一系列与处理空间和时空数据相关的概念，并指向后续章节，其中这些概念将会被更详细地讨论。它还介绍了一系列构成所有空间数据科学语言实现基础的开源技术。

本文由SCY原创，转载注明出处。

时间序列趋势分析

时序趋势分析是一种重要的统计技术，用于分析随时间变化的数据模式。以下是一些常用的时序趋势分析方法：

• 线性和非线性趋势分析：

  通过线性回归或非线性回归模型，可以估计时间序列数据的趋势。例如，一元线性回归可以用来估计线性趋势，而多项式回归或指数回归可以用来估计非线性趋势。

• 季节分解（Seasonal Decomposition）：

  季节分解可以分离出时序数据中的趋势、季节性和随机组成部分。例如，可以使用STL（Seasonal and Trend decomposition using Loess）或季节分解的经典方法。

• 时间序列平滑（Time Series Smoothing）：

  时间序列平滑方法，如移动平均和指数平滑，可以用来估计时间序列数据的趋势。

• 自相关和偏自相关分析（Autocorrelation and Partial Autocorrelation Analysis）：

  通过计算时间序列数据的自相关和偏自相关函数，可以识别数据中的循环模式和趋势。

• 时间序列分解（Time Series Decomposition）：

  将时间序列分解为趋势、周期和随机噪声组件，通常使用模型如X-12-ARIMA或其变体。

• 时间序列模型（Time Series Modeling）：

  ARIMA（自回归整合移动平均模型）、ETS（错误、趋势、季节性模型）和状态空间模型等可以用于分析时间序列的趋势和周期性。

• 波段分析（Wavelet Analysis）：

  波段分析可以在不同的时间尺度上识别时序数据的趋势和周期性。

• Mann-Kendall趋势检验和Sen的斜率估计：

  这些非参数方法用于检测时间序列数据的趋势并估计变化的速度。

• 突变点分析（Changepoint Analysis）：

  用于识别时间序列数据中的结构变化，例如趋势的改变或方差的改变。

• 谱分析（Spectral Analysis）：

  通过分析时序数据的频率域特性来识别周期性和趋势。

一元线性回归和Mann-Kendall检验

本文重点介绍一元线性回归和Mann-Kendall检验的原理和R语言实现步骤，以1980——2020年的降雨数据为例。

一元线性回归和Mann-Kendall检验是分析时间序列趋势的两种不同方法，它们各自有其特点和适用场景。下面是对它们的比较和联系的说明：

1. 基本原理：

一元线性回归是基于参数的统计方法，它假设数据之间存在线性关系，并试图找到描述这种关系的线性方程。它提供了估计的斜率和截距，以及相关的统计测试，以评估这种关系的显著性和强度。

Mann-Kendall检验是一种非参数的统计方法，用于检测时间序列数据中的单调趋势，而不假设数据之间的特定关系。它基于数据的秩次，而不是数据的实际值。

2. 假设：

一元线性回归通常需要满足一些基本假设，例如误差的正态性和独立性，以及数据的线性关系。当数据不满足这些假设时，线性回归的结果可能会受到影响。

Mann-Kendall检验作为非参数检验，不需要数据满足正态分布或其他分布假设，因此它对异常值和非正态数据更具鲁棒性。

3. 输出：

一元线性回归提供了详细的模型参数（例如斜率和截距）和预测值，同时也提供了模型的显著性检验结果。

Mann-Kendall检验主要提供了趋势的显著性检验结果，但不提供具体的模型参数或预测值。

4. 适用场景：

当数据具有明确的线性关系，并且满足线性回归的基本假设时，一元线性回归是一个很好的选择。

当数据可能有单调趋势，但不一定是线性的，或者当数据可能包含异常值或不满足正态分布假设时，Mann-Kendall检验可能是一个更好的选择。

5. 联系：

两者都可以用于分析时间序列数据中的趋势，但方法和假设有很大的不同。

在某些情况下，它们可以互补使用。例如，可以首先使用Mann-Kendall检验来确定是否存在显著的趋势，然后使用一元线性回归来估计趋势的具体参数。

6. 扩展：

Mann-Kendall检验可以与Sen的斜率估计结合使用，以提供趋势的斜率估计，这在一定程度上类似于一元线性回归提供的斜率估计。

一元线性回归

一元线性回归用于研究一个变量（自变量）如何线性影响另一个变量（因变量）。以下是一元线性回归的基本原理和步骤：

数学模型：

一元线性回归假设两个变量之间存在线性关系，可以用下面的方程式表示： \[ Y = \beta_0 + \beta_1X + \varepsilon \]

其中： - ( Y ) 是因变量的值， - ( X ) 是自变量的值， - ( \(\beta_0\) ) 是截距项， - ( \(\beta_1\) ) 是斜率项， - ( \(\varepsilon\) ) 是误差项。

参数估计：

一元线性回归的目标是找到最佳的 ( \(\beta_0\) ) 和 ( \(\beta_1\) )，使得模型的预测值与实际值的差异（误差平方和，\(SSE\)）最小。这通常通过最小二乘法（OLS）实现，它的基本思想是最小化所有观测值的残差平方和，即： \[SSE = \sum_{i=1}^n (Y_i - (\beta_0 + \beta_1X_i))^2 \]

通过对 \(SSE\) 关于 \(\beta_0\) 和 \(\beta_1\) 的偏导数并令其为零，可以解得 \(\beta_0\) 和 \(\beta_1\) 的估计值。

显著性检验：

一旦得到了 \(\beta_0\) 和 \(\beta_1\) 的估计值，通常会进行假设检验来评估这些参数是否显著不为零。这涉及到计算\(t\)统计量和对应的\(p\)值。如果\(p\)值低于某个预定的显著性水平（通常为0.05），则认为参数是显著的。

模型评估：

模型评估通常包括计算 \(R^2\) （决定系数）和调整 \(R^2\)，以评估模型对数据的拟合程度。\(R^2\) 表示模型解释的数据变异的比例。

预测：

使用得到的 \(\beta_0\) 和 \(\beta_1\) 估计值，可以对新的 \(X\) 值做预测，并计算预测区间，以估计预测的不确定性。

模型诊断：

模型诊断是检查模型是否满足回归分析的基本假设，例如误差的正态性、独立性和方差齐性。这可以通过残差图、正态概率图和其他诊断图来完成。

案例分析

Mann-Kendall检验和Sen斜率估计

Mann-Kendall（MK）检验和Sen的斜率估计是两种常用的非参数方法，用于分析时间序列数据中的趋势。下面是它们的基本原理介绍：

Mann-Kendall检验：

Mann-Kendall检验是一种非参数检验，用于确定一个数据序列中是否存在单调的趋势。它不需要数据满足特定的分布假设，因此对于非正态分布的数据很有用。MK检验的基本步骤如下： 1. 对于序列中的每一对数据点(\(x_i\), \(x_j\))，计算符号统计量 (\(S\))： \[S = \sum_{i=1}^{n-1} \sum_{j=i+1}^n \text{sign}(x_j - x_i)\] 其中\(\text{sign}(x_j - x_i)\)是符号函数，如果 \(x_j\) > \(x_i\) 则值为 +1，如果 \(x_j\) < \(x_i\) 则值为 -1，如果 \(x_j\) = \(x_i\) 则值为 0。

2. 计算检验统计量 \(U\) 和对应的 \(p\) 值以判断趋势是否显著。

Sen的斜率估计：

Sen的斜率估计是一种非参数方法，用于估计数据序列中的趋势斜率。它通过计算所有可能的数据点对之间的斜率，然后取这些斜率的中位数作为趋势斜率的估计。Sen的斜率估计的基本步骤如下：

1. 对于序列中的每一对数据点\((x_i, y_i)\) 和 \((x_j, y_j)\)，计算斜率： \[d_k = \frac{(y_j - y_i)}{(x_j - x_i)}\] 其中 \(1 \leq i < j \leq n\)。

2. 从所有计算得到的斜率 \(d_k\) 中，取中位数作为Sen的斜率估计： \[b_{\text{Sen}} = \text{median}(d_k)\] Sen的斜率估计可以提供一个关于时间序列数据趋势的稳健（对异常值不敏感）的估计，而Mann-Kendall检验可以提供这种趋势是否显著的证据。通常，这两种方法可以结合使用，以提供对数据趋势的全面理解。在时间序列趋势分析中，MK检验通常用于检测趋势的显著性，而Sen的斜率估计用于量化趋势的大小。

案例分析

本文由SCY原创，转载注明出处。

本文主要讲解地图可视化

本文由SCY原创，转载注明出处。

本文主要讲解ggplot2包的绘图原理及案例。

本文由SCY原创，转载注明出处。

本文主要讲解R的terra包的使用技巧，包括数据导入，批量处理，绘图，导出等内容

本文由SCY翻译自《Geocomputation with R》第十六章

综合本书的内容，引用重复出现的主题/概念，并激发未来应用和发展的方向。本章不需要先修知识。然而，如果您已经阅读并尝试了第一部分（基础）中的练习，并考虑了地理计算如何帮助您解决工作、研究或其他问题，参考了第三部分（应用）中的章节，那么您可能会从中获得更多收益。

本文由SCY翻译自《Geocomputation with R》第十五章

本章通过模拟雾绿洲（也叫lomas）的植被分布，揭示了明显受到水源供应控制的不同植被区域。这个案例研究不仅整合了前几章中的核心观点，还拓展了这些观点，帮助你更全面地掌握使用R进行地理计算的相关技能。

本文由SCY翻译自《Geocomputation with R》第十四章

• T本章要求使用下列包 (tmaptools 必须安装）:

1
2
3
4
5
6

library(sf)
library(dplyr)
library(purrr)
library(terra)
library(osmdata)
library(spDataLarge)

• 必要的数据将会在适当的时候下载。

为了方便读者并确保易于复现，我们已将下载的数据放在 spDataLarge 包中供使用。

引言

本章展示了在第一部分和第二部分学到的技能如何应用于特定领域：地理营销（有时也称为位置分析或位置智能）。这是一个广泛的研究和商业应用领域。地理营销的典型例子是如何选择一个新店的位置。这里的目标是吸引最多的顾客，最终实现最大的利润。此外，还有许多非商业应用可以利用这种技术来造福公众，例如选择新的卫生服务设施的位置。

人类对于位置分析来说是基本要素，特别是人们可能会在哪里花费时间和其他资源。有趣的是，生态学的概念和模型与用于店铺选址分析的概念和模型非常相似。动植物可以在某些“最优”位置最好地满足其需求，这些位置是基于随空间变化的变量确定的。这是地理计算和地理信息科学的一个重要优势：概念和方法可以转移到其他领域。例如，北极熊更喜欢气温较低且食物（海豹和海狮）丰富的北纬地区。同样地，人类倾向于聚集在某些地方，创造出类似于北极地区的生态位的经济位。位置分析的主要任务是基于现有数据找出这些特定服务的“最优位置”在哪里。典型的研究问题包括：

• 目标群体居住在哪里，他们经常去哪些地区？
• 竞争的商店或服务位于哪里？
• 有多少人可以轻松到达特定的商店？
• 现有的服务是否过度或不足地利用了市场潜力？
• 公司在特定地区的市场份额是多少？

本章通过一个基于实际数据的假设案例研究，演示了地理计算如何回答这些问题。

案例研究：德国自行车商店

假设您正在德国开设一家自行车连锁店。这些店铺应该位于尽可能多的潜在顾客所在的城市地区。此外，一个虚构的调查（仅用于本章，非商业用途！）表明，单身年轻男性（年龄在20到40岁之间）最有可能购买您的产品：这就是目标受众。幸运的是，您有足够的资本来开设多家店铺。但是，它们应该放在哪里呢？咨询公司（雇佣地理营销分析师的公司）通常会收取高额费用来回答此类问题。幸运的是，借助开放数据和开源软件，我们可以自己解决这些问题。以下几节将演示如何将本书前几章学到的技术应用于执行服务位置分析中的常见步骤：

• 整理来自德国人口普查的输入数据
• 将表格化的人口普查数据转换为栅格对象
• 识别人口密度较高的大都市地区
• 为这些地区下载详细的地理数据（使用osmdata从OpenStreetMap下载）
• 使用地图代数创建评分栅格，以评估不同位置的相对吸引力

尽管我们将这些步骤应用于一个特定的案例研究，但它们可以推广到许多店铺选址或公共服务提供的情景。

整理数据

德国政府提供了分辨率为1公里或100米的栅格化人口普查数据。下面的代码块用于下载、解压缩和读取1公里分辨率的数据。

1
2
3
4

download.file("https://tinyurl.com/ybtpkwxz", 
              destfile = "census.zip", mode = "wb")
unzip("census.zip") # unzip the files
census_de = readr::read_csv2(list.files(pattern = "Gitter.csv"))

请注意，census_de 数据也可以从 spDataLarge 包中获取：

1

data("census_de", package = "spDataLarge")

census_de 对象是一个包含13个变量的数据框，涵盖德国境内超过360,000个栅格单元。对于我们的工作，我们只需要其中的一个子集：东向坐标 (x)、北向坐标 (y)、居民数量（人口；pop）、平均年龄（mean_age）、女性比例（women）和平均家庭规模（hh_size）。下面的代码块从这些变量中选择并将其从德文重命名为英文，并在表格 @ref(tab:census-desc) 中进行了总结。此外，mutate() 函数被用于将值-1和-9（表示“未知”）转换为 NA。

1
2
3
4
5
6
7

# pop = population, hh_size = household size
input = select(census_de, x = x_mp_1km, y = y_mp_1km, pop = Einwohner,
                      women = Frauen_A, mean_age = Alter_D, hh_size = HHGroesse_D)
# set -1 and -9 to NA
input_tidy = dplyr::mutate(input,
  dplyr::across(.cols = c(pop, women, mean_age, hh_size), 
                .fns =  ~ifelse(.x %in% c(-1, -9), NA, .x)))

创建人口普查栅格

在预处理之后，可以使用rast()函数将数据转换为SpatRaster对象。当将其 type 参数设置为 xyz 时，输入数据框的 x 和 y 列应该对应于正规栅格上的坐标。所有其余列（在这里是 pop、women、mean_age、hh_size）将用作栅格图层的值（请参阅图 @ref(fig:census-stack)；还可以在我们的 GitHub 存储库中的code/14-location-figures.R 文件中找到相关代码）。

1

input_ras = terra::rast(input_tidy, type = "xyz", crs = "EPSG:3035")

1
2
3
4
5
6
7
8
9
10

input_ras
#> class       : SpatRaster 
#> dimensions  : 868, 642, 4  (nrow, ncol, nlyr)
#> resolution  : 1000, 1000  (x, y)
#> extent      : 4031000, 4673000, 2684000, 3552000  (xmin, xmax, ymin, ymax)
#> coord. ref. : ETRS89-extended / LAEA Europe (EPSG:3035) 
#> source(s)   : memory
#> names       : pop, women, mean_age, hh_size 
#> min values  :   1,     1,        1,       1 
#> max values  :   6,     5,        5,       5

📌请注意，我们使用的是等面积投影（EPSG:3035；欧洲兰伯特等面积投影），即每个栅格单元的面积相同，这里是1000 x 1000平方米。由于我们主要使用诸如每个栅格单元的居民数量或女性比例等密度数据，因此每个栅格单元的面积相同非常重要，以避免“苹果与橙子比较”。请注意，在地理坐标系（CRS）中，栅格单元的面积在向极地方向不断减小，因此需要谨慎处理。

Gridded German census data of 2011 (see Table @ref(tab:census-desc) for a description of the classes).

下一步是根据前面章节提到的调查，使用 terra 函数 classify() 对存储在 input_ras 中的栅格图层的值进行重新分类。对于人口数据，我们使用类均值将类别转换为数值数据类型。假设栅格单元的值为1（'class 1' 中的单元包含3到250名居民），则栅格单元的人口被假定为127；如果值为2（包含250到500名居民的单元），则人口被假定为375，依此类推（请参阅上表）。对于'类别6'，栅格单元的值被选为8000名居民，因为这些单元包含的人口超过8000人。当然，这些都是对真实人口的近似值，而不是精确值。¹然而，这个详细级别足以划定大都市区域（请参阅下一节）。

与变量 pop 不同，该变量表示总人口的绝对估计，其他变量被重新分类为与调查中使用的权重相对应的权重。例如，变量 women 中的'类别1'代表的是人口中女性占0到40%的地区；这些地区被重新分类为较高的权重3，因为目标人群主要是男性。同样地，包含最年轻人口和最高比例的单身家庭的类别被重新分类为高权重。

1
2
3
4
5
6
7
8
9

rcl_pop = matrix(c(1, 1, 127, 2, 2, 375, 3, 3, 1250, 
                   4, 4, 3000, 5, 5, 6000, 6, 6, 8000), 
                 ncol = 3, byrow = TRUE)
rcl_women = matrix(c(1, 1, 3, 2, 2, 2, 3, 3, 1, 4, 5, 0), 
                   ncol = 3, byrow = TRUE)
rcl_age = matrix(c(1, 1, 3, 2, 2, 0, 3, 5, 0),
                 ncol = 3, byrow = TRUE)
rcl_hh = rcl_women
rcl = list(rcl_pop, rcl_women, rcl_age, rcl_hh)

请注意，我们确保了列表中重新分类矩阵的顺序与 input_ras 的元素顺序相同。例如，第一个元素在两种情况下都对应于人口。随后，for 循环 将重新分类矩阵应用于相应的栅格图层。最后，下面的代码块确保 reclass 图层与 input_ras 的图层具有相同的名称。

1
2
3
4
5

reclass = input_ras
for (i in seq_len(terra::nlyr(reclass))) {
  reclass[[i]] = terra::classify(x = reclass[[i]], rcl = rcl[[i]], right = NA)
}
names(reclass) = names(input_ras)

1
2
3
4
5

reclass
#> ... (full output not shown)
#> names       :  pop, women, mean_age, hh_size 
#> min values  :  127,     0,        0,       0 
#> max values  : 8000,     3,        3,       3

定义都市区

我们特意将大都市区域定义为面积为20平方公里且居住人口超过500,000人的像素。在这种粗分辨率下，可以通过使用 aggregate()快速创建像素。下面的命令使用参数 fact = 20 将结果的分辨率降低了20倍（回想一下，原始栅格的分辨率是1平方公里）：

1
2
3
4
5
6
7
8
9
10

pop_agg = terra::aggregate(reclass$pop, fact = 20, fun = sum, na.rm = TRUE)
summary(pop_agg)
#>       pop         
#>  Min.   :    127  
#>  1st Qu.:  39886  
#>  Median :  66008  
#>  Mean   :  99503  
#>  3rd Qu.: 105696  
#>  Max.   :1204870  
#>  NA's   :447

接下来的步骤是仅保留居住人口超过500,000人的单元格。您可以使用以下代码来实现：

1

pop_agg = pop_agg[pop_agg > 500000, drop = FALSE] 

绘制这些数据将显示出八个大都市区域。每个区域由一个或多个栅格单元组成。如果我们能够将属于同一区域的所有单元格连接起来，那将会很好。terra 的 patches() 命令正是这样的功能。随后，as.polygons() 将栅格对象转换为空间多边形，而 st_as_sf() 将其转换为 sf 对象。

1
2
3
4

metros = pop_agg |> 
  terra::patches(directions = 8) |>
  terra::as.polygons() |>
  sf::st_as_sf()

The aggregated population raster (resolution: 20 km) with the identified metropolitan areas (golden polygons) and the corresponding names.

生成的八个适合自行车店的大都市区域；请参阅 code/14-location-figures.R 以创建图表）仍然缺少名称。逆地理编码方法可以解决这个问题：根据坐标找到对应的地址。因此，提取每个大都市区域的中心坐标可以作为逆地理编码API的输入。这正是 tmaptools 包中的 rev_geocode_OSM() 函数所期望的。此外，将 as.data.frame 设置为 TRUE 将返回一个包含多列关于位置的 data.frame，包括街道名称、门牌号和城市。然而，在这里，我们只关注城市的名称。

1
2
3
4
5
6

metro_names = sf::st_centroid(metros, of_largest_polygon = TRUE) |>
  tmaptools::rev_geocode_OSM(as.data.frame = TRUE) |>
  select(city, town, state)
# smaller cities are returned in column town. To have all names in one column,
# we move the town name to the city column in case it is NA
metro_names = dplyr::mutate(metro_names, city = ifelse(is.na(city), town, city))

为了确保读者使用完全相同的结果，我们已将它们存储在 spDataLarge 中，对象名称为 metro_names。

总体而言，我们对 city 列作为大都市名称（上表）感到满意，除了一个例外，即属于杜塞尔多夫大区的费尔伯特。因此，我们将 Velbert 替换为 Düsseldorf（上图）。像 ü 这样的特殊字符可能会在后续的操作中引起问题，例如在使用 opq() 确定大都市区域的边界框时（请参见下面的内容），因此我们避免使用它们。

1
2
3
4

metro_names = metro_names$city |> 
  as.character() |>
  {\(x) ifelse(x == "Velbert", "Düsseldorf", x)}() |>
  {\(x) gsub("ü", "ue", x)}()

兴趣点

osmdata包提供了易于使用的访问OSM数据的方法。我们不是下载整个德国的商店数据，而是将查询限制在了定义好的大都市区域，以减少计算负担并且只获取感兴趣区域内的商店位置。下面的代码块使用了许多函数，包括：

• map()（tidyverse 中的 lapply()等效函数），它遍历了八个大都市名称，随后在OSM查询函数 opq() 中定义了边界框。
• add_osm_feature() 用于指定带有键值为 shop 的OSM元素（请参阅 wiki.openstreetmap.org 以获取常见的键值对列表）。
• osmdata_sf() 将OSM数据转换为空间对象（sf 类）。
• while()，如果第一次下载失败，将重复尝试（在本例中为三次）。²

在运行此代码之前，请考虑它将下载近2GB的数据。为了节省时间和资源，我们已将名为 shops 的输出放入 spDataLarge 包中。要在您的环境中使用它，请运行 data("shops", package = "spDataLarge")。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

shops = purrr::map(metro_names, function(x) {
  message("Downloading shops of: ", x, "\n")
  # give the server a bit time
  Sys.sleep(sample(seq(5, 10, 0.1), 1))
  query = osmdata::opq(x) |>
    osmdata::add_osm_feature(key = "shop")
  points = osmdata::osmdata_sf(query)
  # request the same data again if nothing has been downloaded
  iter = 2
  while (nrow(points$osm_points) == 0 && iter > 0) {
    points = osmdata_sf(query)
    iter = iter - 1
  }
  # return only the point features
  points$osm_points
})

在我们定义的任何大都市区域中几乎不可能没有商店。以下的 if 条件仅仅检查每个区域是否至少有一家商店。如果没有的话，我们建议尝试再次为该/这些特定区域下载商店数据。

1
2
3
4
5
6

# checking if we have downloaded shops for each metropolitan area
ind = purrr::map_dbl(shops, nrow) == 0
if (any(ind)) {
  message("There are/is still (a) metropolitan area/s without any features:\n",
          paste(metro_names[ind], collapse = ", "), "\nPlease fix it!")
}

为确保每个列表元素（一个 sf 数据框）具有相同的列³，我们只保留 osm_id 和 shop 列，使用 map_dfr 循环将所有商店合并为一个大的 sf 对象。

1
2

# select only specific columns
shops = purrr::map_dfr(shops, select, osm_id, shop)

注意：shops 已经在 spDataLarge 中提供，并且可以通过以下方式访问：

1

data("shops", package = "spDataLarge")

唯一剩下的任务是将空间点对象转换为栅格。sf 对象 shops 将被转换为一个栅格，其参数（维度、分辨率、CRS）与 reclass 对象相同。重要的是，在此处使用 length() 函数来计算每个单元格中的商店数量。

因此，下面的代码块的结果是商店密度的估计（商店/平方公里）。在使用 rasterize() 之前，使用 st_transform() 来确保两个输入的CRS匹配。

1
2
3

shops = sf::st_transform(shops, st_crs(reclass))
# create poi raster
poi = terra::rasterize(x = shops, y = reclass, field = "osm_id", fun = "length")

与其他栅格图层（人口、女性、平均年龄、户籍人口）一样，poi 栅格也被重新分类为四个类别。在一定程度上，定义类别间隔是一个主观的任务。可以使用等距断点、分位数断点、固定值或其他方法。在这里，我们选择了费舍尔-詹金斯自然断点法，该方法最小化了类内方差，其结果为重新分类矩阵提供了一个输入。

1
2
3
4
5
6
7
8
9

# construct reclassification matrix
int = classInt::classIntervals(terra::values(poi), n = 4, style = "fisher")
int = round(int$brks)
rcl_poi = matrix(c(int[1], rep(int[-c(1, length(int))], each = 2), 
                   int[length(int)] + 1), ncol = 2, byrow = TRUE)
rcl_poi = cbind(rcl_poi, 0:3)  
# reclassify
poi = terra::classify(poi, rcl = rcl_poi, right = NA) 
names(poi) = "poi"

明确合适的位置

在将所有图层组合在一起之前，只剩下几个步骤：将 poi 添加到 reclass 栅格堆叠中，并将人口图层从中移除。后者的原因有两点。首先，我们已经勾勒出了大都市区域，即人口密度高于德国其他地区平均水平的地区。其次，虽然在特定的服务区域内有许多潜在的顾客可能是有优势的，但仅仅数量本身可能并不真正代表所需的目标群体。例如，高层住宅区是人口密度较高的地区，但不一定具有购买昂贵自行车配件的高购买力。

1
2
3

# remove population raster and add poi raster
reclass = reclass[[names(reclass) != "pop"]] |>
  c(poi)

与其他数据科学项目一样，数据检索和“整理”在整个工作量中占据了很大一部分。有了干净的数据，最后一步——通过将所有栅格图层相加来计算最终得分——可以在一行代码中完成。

1
2

# calculate the total score
result = sum(reclass)

例如，得分大于9的分数可能是一个适当的阈值，表示可以放置自行车店的栅格单元格；请参阅 code/14-location-figures.R）。

1
2
3
4
5
6
7
8

#> The legacy packages maptools, rgdal, and rgeos, underpinning the sp package,
#> which was just loaded, will retire in October 2023.
#> Please refer to R-spatial evolution reports for details, especially
#> https://r-spatial.org/r/2023/05/15/evolution4.html.
#> It may be desirable to make the sf package available;
#> package maintainers should consider adding sf to Suggests:.
#> The sp package is now running under evolution status 2
#>      (status 2 uses the sf package in place of rgdal)

1
2

<div class="leaflet html-widget html-fill-item-overflow-hidden html-fill-item" id="htmlwidget-841de324d41155df19a0" style="width:100%;height:415.296px;"></div>
<script type="application/json" data-for="htmlwidget-841de324d41155df19a0">{"x":{"options":{"crs":{"crsClass":"L.CRS.EPSG3857","code":null,"proj4def":null,"projectedBounds":null,"options":{}}},"calls":[{"method":"addTiles","args":["https://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png",null,null,{"minZoom":0,"maxZoom":18,"tileSize":256,"subdomains":"abc","errorTileUrl":"","tms":false,"noWrap":false,"zoomOffset":0,"zoomReverse":false,"opacity":1,"zIndex":1,"detectRetina":false,"attribution":"&copy; <a href=\"https://openstreetmap.org/copyright/\">OpenStreetMap<\/a>,  <a href=\"https://opendatacommons.org/licenses/odbl/\">ODbL<\/a>"}]},{"method":"addRasterImage","args":["data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAACgAAAAoCAYAAACM/rhtAAAAM0lEQVRYhe3SMQ0AMAwEsQcb/hTSuQSaVLIR3HAJAADfqfR0wm1dEC9VevcCq+MAAIBhB20vBvBK3JZrAAAAAElFTkSuQmCC",[[52.69660085729196,13.08200261479863],[52.32107408861835,13.69815913809763]],0.8,null,null,null]},{"method":"addLegend","args":[{"colors":["darkgreen"],"labels":["potential locations"],"na_color":null,"na_label":"NA","opacity":0.5,"position":"bottomright","type":"unknown","title":"Legend","extra":null,"layerId":null,"className":"info legend","group":null}]}],"limits":{"lat":[52.32107408861835,52.69660085729196],"lng":[13.08200261479863,13.69815913809763]}},"evals":[],"jsHooks":[]}</script>

讨论和下一步

所呈现的方法是GIS的典型应用示例。我们将调查数据与基于专家知识和假设的方法相结合（定义大都市区域，定义类别间隔，定义最终得分阈值）。这种方法不如应用分析适用于科学研究，因为它可以提供基于证据的适合自行车店的区域，应该与其他信息来源进行比较。对方法的一些变更可以改进分析：

• 在计算最终得分时，我们使用了相等的权重，但其他因素，如户籍人口，可能与女性比例或平均年龄一样重要。
• 我们使用了所有兴趣点，但只有与自行车店相关的兴趣点，如自行车店、五金店、自行车、钓鱼、狩猎、摩托车、户外和运动用品店（请参阅 OSM Wiki 上可用的店铺值范围）可能会产生更精细的结果。
• 更高分辨率的数据可能会改善输出（请参阅练习）。
• 我们仅使用了有限的变量集和来自其他来源的数据，如 INSPIRE geoportal 或来自 OpenStreetMap 的自行车路径数据，可以丰富分析。
• 未考虑交互作用，如男性比例和单身户之间可能存在的关系。

简而言之，分析可以在多个方向上进行扩展。然而，它应该给您对如何在geomarketing背景下在R中获取和处理空间数据的第一印象和理解。

最后，我们必须指出，所呈现的分析仅仅是找到合适位置的第一步。到目前为止，我们已经确定了大小为1x1公里的区域，代表根据我们的调查可能适合自行车店的位置。分析的后续步骤可以是：

• 基于特定服务区域内的居民数量找到最佳位置。例如，在骑自行车的15分钟行程范围内，店铺应该为尽可能多的人提供可达性（服务区域路由）。在此过程中，我们应该考虑到远离店铺的人越远，他们实际访问店铺的可能性就越小（距离衰减函数）。
• 此外，考虑竞争对手也是一个好主意。也就是说，如果已经有一家自行车店在所选位置附近，可能的顾客（或销售潜力）应该在竞争对手之间分配。
• 我们需要找到适合并且价格合理的房地产，例如在可访问性、停车位的可用性、过路人的期望频率、有大窗户等方面。

练习

E1. 首先，您需要从提供的链接下载包含居民信息的CSV文件（以100米单元格分辨率）。请注意，解压缩后的文件大小为1.23 GB。您可以使用readr::read_csv将其读入R中。在具有16 GB RAM的机器上，这需要30秒钟。data.table::fread()可能会更快，它返回一个data.table()类的对象。使用dplyr::as_tibble()将其转换为tibble。构建一个居民栅格，将其聚合到1 km的单元格分辨率，并将其与我们使用类平均值创建的居民栅格（inh）进行比较。

E2. 假设我们的自行车店主要向老年人销售电动自行车。相应地更改年龄栅格，重复剩余的分析，并将更改与我们的原始结果进行比较。

本文由SCY翻译自《Geocomputation with R》第十三章

前提条件

• 本章使用下列包：¹

²：nabor 必须安装。

1
2
3
4
5
6
7

library(sf)
library(dplyr)
library(spDataLarge)
library(stplanr)      # for processing geographic transport data
library(tmap)         # map making (see Chapter 9)
library(ggplot2)      # data visualization package
library(sfnetworks)   # spatial network classes and functions 

引言

很少有其他部门的地理空间比交通更有形。移动(克服距离)的努力是地理学“第一定律”的核心，Waldo Tobler 在1970年定义如下(Tobler 1970) :

一切事物都与其他事物相关，但是近的事物比远的事物更相关。

这个“定律”是空间自相关和其他关键地理概念的基础。它适用于各种各样的现象，如友谊网络和生态多样性，可以用交通成本来解释——在时间、精力和金钱方面——这些构成了“距离摩擦”。从这个角度来看，运输技术具有颠覆性，改变了包括移动人员和货物在内的地理实体之间的空间关系: “运输的目的是克服空间”(Rodrigue，Comtois，and Slack 2013)。

运输是一个固有的空间活动，包括从一个起点‘ A’到目的地‘ B’，通过无限的位置之间。因此，交通运输研究人员长期以来一直采用地理和计算方法来理解运动模式，以及干预措施如何能够提高他们的表现，这并不令人惊讶(Lovelace 2021)。

本章介绍了不同地理层次运输系统的地理分析:

• Areal units: 交通模式可以通过参考区域总量来理解，比如主要的出行方式(例如，汽车、自行车或步行) ，以及生活在特定区域的人们的平均出行距离(见第13.3节)
• Desire lines: 表示“起点-目的地”数据的直线，记录有多少人在地理空间的不同地点(点或区域)之间旅行(或可能旅行) ，这是第13.4节的主题
• Nodes: 这些是交通系统中可以代表共同起点和目的地的点，以及公共交通站点，如公共汽车站和火车站，这是第13.5节的主题
• Routes: 这些线路表示沿着路由网络沿着所需的线路和节点之间的路径。路由(可以表示为单行字符串或多个短段)和生成它们的路由引擎，在第13.6节中介绍
• Route networks: 这些代表了一个区域内的道路、路径和其他线性特征的系统，在第13.6节中有介绍。 它们可以表示为地理特征（通常是组成完整网络的短路段）或结构化为一个相互连接的图，不同段上的交通流量被运输模型师称为“流量”(Hollander 2016)。

This highlights an important feature of transport systems: they are closely linked to broader phenomena and land-use patterns. 另一个关键层次是代理，如你我以及使我们能够移动的交通工具，比如自行车和公交车。这些可以在像MATSim和A/B Street这样的软件中以计算方式表示，它们使用基于代理的建模（ABM）框架，通常具有很高的空间和时间分辨率，来表示交通系统的动态性(Horni，Nagel 和 Axhausen 2016)。尽管ABM是一种具有很大整合潜力的强大的交通研究方法，尤其是与R的空间类(Thiele 2014; Lovelace and Dumont 2016) ，但它不在本章的范围内。除了地理层次和代理之外，许多交通模型中的基本分析单位是行程，即从一个起点'A'到一个终点'B'的单一目的地之旅(Hollander 2016)。行程连接了交通系统的不同层次，并可以简单地表示为连接区域质心（节点）的地理需求线，或作为遵循交通路线网络的路线。在这种情况下，代理人通常是在交通网络内移动的点实体。

交通系统是动态的(Xie and Levinson 2011)。尽管本章的重点是交通系统的地理分析，但它也提供了如何使用这种方法来模拟变化情景的见解，在第13.8节中有介绍。地理交通建模的目的可以理解为以捕捉其本质的方式简化这些时空系统的复杂性。选择适当的地理分析层次可以在不失去其最重要的特征和变量的情况下，简化这种复杂性，从而实现更好的决策和更有效的干预(Hollander 2016)。

通常，模型被设计用来解决特定问题，比如如何提高交通系统的安全性或环境性能。因此，本章围绕一个政策情境展开，将在下一节中介绍，即：如何增加布里斯托尔市的自行车使用？第14章演示了地理计算的一个相关应用：优先考虑新自行车店的位置。章节之间有一个联系：新的和有效定位的自行车基础设施可以使人们开始骑自行车，从而增加对自行车店和当地经济活动的需求。这突出了交通系统的一个重要特点：它们与更广泛的现象和土地利用模式密切相关。

布里斯托尔案例研究

本章使用的案例研究位于英格兰西部的布里斯托尔市，距离威尔士首都加的夫约30公里以东。该地区的交通网络概览在下图中展示，图中展示了多样的交通基础设施，包括自行车、公共交通和私家车。

布里斯托尔是一个充满活力和多样性的城市，拥有各种各样的交通选项。从自行车道到公交线路，再到繁忙的高速公路，布里斯托尔的交通网络反映了其居民多样化的出行需求。这为进行地理和交通研究提供了一个极为丰富的背景。

Bristol's transport network represented by colored lines for active (green), public (railways, black) and private motor (red) modes of travel. Black border lines represent the inner city boundary (highlighted in yellow) and the larger Travel To Work Area (TTWA).

布里斯托尔是英格兰第十大的市议会，拥有50万人口，尽管其交通服务覆盖区域更广泛（参见13.3节）。该市拥有一个充满活力的经济体，包括航空航天、媒体、金融服务和旅游公司，以及两所主要大学。尽管布里斯托尔的人均收入较高，但也有严重贫困的地区(布里斯托尔市议会2015年)。

从交通角度看，布里斯托尔的铁路和公路连接非常便利，活跃出行的水平相对较高。根据Active People Survey，19%的市民每月至少骑一次自行车，88%的人每月至少走一次路（全国平均分别为15%和81%）。2011年的人口普查数据显示，8%的人口自行车上班，而全国平均仅为3%。

与许多城市一样，布里斯托尔面临严重的交通拥堵、空气质量和体力不活跃的问题。自行车出行能有效地解决这些问题：与步行相比，自行车更有可能替代汽车出行，因为其典型速度为15-20 km/h，而步行为4-6 km/h。因此，布里斯托尔的交通战略对自行车出行有着雄心勃勃的计划。

为了强调政策考虑在交通研究中的重要性，本章旨在为那些负责让人们从汽车转向更可持续出行方式（特别是步行和自行车）的人（如交通规划师、政治家和其他利益相关者）提供证据。更广泛的目标是演示地理计算如何支持基于证据的交通规划。

在本章中，你将学习如何：

• 描述城市交通行为的地理模式
• 识别支持多模式出行的关键公共交通节点
• 分析旅行"愿望线"，找出哪里有许多人驾车短途出行
• 识别鼓励少开车、多骑自行车的自行车路线位置

为了在实践方面展开本章的内容，下一节将开始加载关于旅行模式的区域数据。这些区域级别的数据集虽小，但对于获取对城市整体交通系统的基本了解至关重要。

交通区域

尽管交通系统主要基于线性特征和节点——包括路径和车站——从面数据开始通常更有意义，以将连续空间划分为可触及的单位 (Hollander 2016)。除了定义研究区域的边界（在这种情况下是布里斯托尔）之外，对交通研究人员尤为感兴趣的两种区域类型是：起源区域和目的区域。通常，相同的地理单位用于起点和目的地。然而，不同的分区系统，例如 '工作区'，可能更适合代表在有许多'旅行吸引点'的区域内旅行目的地的增加密度 (国家统计局2014年)。

最简单定义研究区域的方式通常是第一个由OpenStreetMap返回的匹配边界。这可以通过诸如 osmdata::getbb("Bristol", format_out = "sf_polygon", limit = 1) 的命令来完成。这返回一个 sf 对象（或者如果没有指定 limit = 1，则返回一组 sf 对象），代表最大匹配城市区域的范围，要么是边界框的矩形多边形，要么是一个详细的多边形边界。³对于布里斯托尔，返回了一个详细的多边形，由 spDataLarge 包中的 bristol_region 对象表示。请参见上图中的内部蓝色边界：这种方法有几个问题：

• 第一个由OSM返回的OSM边界可能不是当地政府使用的官方边界
• 即使OSM返回官方边界，这也可能不适合交通研究，因为它们与人们的旅行路线关系不大

工作通勤区域（TTWAs）通过创建一个类似于水文流域的分区系统来解决这些问题。TTWAs最初被定义为连续区域，在其中75%的人口前往工作(Coombes, Green, and Openshaw 1986)，这是本章使用的定义。由于布里斯托尔是一个主要的雇主，吸引了周边城镇的旅行，因此其TTWA比城市界限要大得多（见上图）。代表这一以交通为导向的边界的多边形存储在由章节开始时加载的 spDataLarge 包提供的 bristol_ttwa 对象中。

本章使用的起点和终点区域是相同的：官方定义的中等地理分辨率区域（它们的官方名称是中层超级输出区域或MSOAs）。 每个区域大约有8,000人。这种行政区域可以为交通分析提供重要的背景信息，例如可能最受特定干预措施影响的人群类型(例如, Moreno-Monroy, Lovelace, and Ramos 2017)。

这些区域的地理分辨率很重要：小区域具有高地理分辨率通常更可取，但它们在大区域中的高数量可能对处理（尤其是对于起点-终点分析，在其中可能性作为区域数量的非线性函数增加）有影响 (Hollander 2016)。

📌与小区域相关的另一个问题与匿名规则有关。为了使在区域内无法推断个人身份，详细社会人口统计变量通常仅在低地理分辨率下可用。例如，旅行方式的年龄和性别细分在英国的地方当局层面上可用，但不在更高的输出区域层面上，其中每个区域大约包含100户家庭。有关更多详细信息，请参阅 www.ons.gov.uk/methodology/geography。

本章中使用的102个区域存储在 bristol_zones 中，如下图所示。注意，在人口密集地区，区域变得更小：每个区域都有类似数量的人。bristol_zones 不包含关于交通的属性数据，只有每个区域的名称和代码：

1
2

names(bristol_zones)
#> [1] "geo_code" "name"     "geometry"

为了添加旅行数据，我们将执行属性连接，这是在矢量属性连接节中描述的常见任务。我们将使用来自英国2011年人口普查关于工作通勤的问题的旅行数据，该数据存储在 bristol_od 中，由 ons.gov.uk 数据门户提供。bristol_od 是一个关于英国2011年人口普查中各区域间通勤旅行的起点-终点（OD）数据集。第一列是起始区域的ID，第二列是目的区域。bristol_od 的行数比 bristol_zones 多，代表的是区域之间的旅行，而不是区域本身：

1
2
3
4

nrow(bristol_od)
#> [1] 2910
nrow(bristol_zones)
#> [1] 102

前一个代码块的结果显示，每个区域有超过10个起点-终点（OD）对，这意味着在将其与 bristol_zones 进行连接之前，我们需要对起点-终点数据进行聚合，如下图所示 ：

1
2
3
4

zones_attr = bristol_od |> 
  group_by(o) |> 
  summarize(across(where(is.numeric), sum)) |> 
  dplyr::rename(geo_code = o)

前面的代码块执行了以下操作：

• 按照起点区域（包含在o列中）对数据进行了分组
• 如果在bristol_od 数据集中的变量是数值型的，就对它们进行了聚合，以找出每个区域中按交通方式分布的人口总数⁴
• 将分组变量o重命名，使其与bristol_zones对象中的ID列geo_code匹配

由此产生的对象zones_attr是一个数据框，其行代表各个区域和一个ID变量。我们可以使用%in%操作符验证这些ID是否与zones数据集中的ID匹配，操作如下：

1
2
3

summary(zones_attr$geo_code %in% bristol_zones$geo_code)
#>    Mode    TRUE 
#> logical     102

结果显示，新对象中包含所有102个区域，而zone_attr具有可以与区域连接的形式⁵。这是通过使用连接函数 left_join() 来完成的（注意，在这里inner_join()也会产生相同的结果）。

1
2
3
4
5
6

zones_joined = left_join(bristol_zones, zones_attr, by = "geo_code")
sum(zones_joined$all)
#> [1] 238805
names(zones_joined)
#> [1] "geo_code"   "name"       "all"        "bicycle"    "foot"      
#> [6] "car_driver" "train"      "geometry"

结果是zones_joined，它包含新的列，代表研究区域内每个区域起始的总出行次数（接近四分之一百万次）以及它们的出行方式（自行车、步行、汽车和火车）。出行起点的地理分布在下图的左侧地图中有所体现。这显示了大多数区域在研究区内有0至4000次的出行起始。生活在布里斯托尔市中心附近的人进行的出行次数更多，而在城市边缘的则较少。这是为什么呢？记住我们只在研究区域内处理出行：边缘区域低出行次数可以归因于这些周边区域的许多人会前往研究区域之外的其他地区。通过一个特殊的目的地ID，研究区域之外的出行可以被包括在区域模型中，该ID涵盖了任何前往模型中未表示的区域的出行(Hollander 2016)。然而，bristol_od中的数据简单地忽略了这样的出行：它是一个“区域内”的模型。

与OD（出发-目的地）数据集可以聚合到出发区域的方式相同，它们也可以被聚合以提供关于目的地区域的信息。人们倾向于聚集在中心地区。这解释了为什么下图右侧面板中表示的空间分布相对不均匀，最常见的目的地区域集中在布里斯托尔市中心。结果是zones_od，其中包含一个新列，报告了任何出行方式的目的地数量，如下所创建：

1
2
3
4
5

zones_destinations = bristol_od |> 
  group_by(d) |> 
  summarize(across(where(is.numeric), sum)) |> 
  select(geo_code = d, all_dest = all)
zones_od = inner_join(zones_joined, zones_destinations, by = "geo_code")

下面的代码创建了一个简化版本的下图（参见书籍GitHub仓库的 code 文件夹中的 12-zones.R 文件以重现该图，以及有关使用tmap创建分面地图的详细信息）：

1
2

qtm(zones_od, c("all", "all_dest")) +
  tm_layout(panel.labels = c("Origin", "Destination"))

Number of trips (commuters) living and working in the region. The left map shows zone of origin of commute trips; the right map shows zone of destination (generated by the script 13-zones.R).

欲望线

欲望线连接起点和终点，表示人们想要去的地方，通常是在各个区域之间。它们代表了从A点到B点最快的“直线”或“鸟飞”路线，如果没有建筑物和弯曲道路的阻碍的话（我们将在下一节看到如何将欲望线转化为实际路线）。通常，欲望线在地理上表示为每个区域的地理（或人口加权）中心的起点和终点。这是我们将在本节中创建和使用的欲望线类型，尽管值得注意的是，“抖动”技术能够增加多个起点和终点，以提高基于OD数据构建的分析的空间覆盖和准确性 (Lovelace，Félix 和 Carlino 2022)。

我们已经在数据集bristol_od中加载了代表欲望线的数据。这个起点-终点（OD）数据框对象表示在o和d所代表的区域之间旅行的人数，如下表所示。要按所有旅行对OD数据进行排序，然后只过滤出前5名，请输入：

1
2

od_top5 = bristol_od |> 
  slice_max(all, n = 5)

生成的表格提供了关于布里斯托尔通勤（上班往返）出行模式的快照。它表明，在前5大起点-终点（OD）组合中，步行是最受欢迎的出行方式，而E02003043 区域是一个受欢迎的目的地（布里斯托尔市中心，所有前5大OD组合的目的地），并且区内 出行，从E02003043 区域的一个部分到另一个部分（表的第一行），是数据集中出行最多的OD组合。但从政策角度看，表中呈现的原始数据用处有限：除了它仅包含了2,910个OD组合中的一小部分之外，它几乎没有告诉我们在哪里需要政策措施，或者有多少比例的出行是通过步行和骑自行车完成的。 以下命令计算了由这些积极模式完成的每条欲望线的百分比：

1
2

bristol_od$Active = (bristol_od$bicycle + bristol_od$foot) /
  bristol_od$all * 100

有两种主要类型的OD（起点-终点）组合：区间（Interzonal） 和 区内（Intrazonal）。区间的OD组合代表了目的地与起点不同区域之间的出行。区内的OD组合代表了在同一区域内的出行（参见表的顶行）。以下代码块将 od_bristol 分为这两种类型：

1
2

od_intra = filter(bristol_od, o == d)
od_inter = filter(bristol_od, o != d)

下一步是将区间（interzonal）的OD组合转换为一个表示“欲望线（desire lines）”的sf对象，这可以通过使用 stplanr 函数 od2line() 来在地图上进行绘制。⁶

1
2

desire_lines = od2line(od_inter, zones_od)
#> Creating centroids representing desire line start and end points.

下图展示了结果的示例，其简化版本可以通过以下命令创建（要完全复制该图形，请查看 13-desire.R 中的代码，有关使用 tmap 进行可视化的详细信息，请参见章节高级制图）：

1

qtm(desire_lines, lines.lwd = "all")

Desire lines representing trip patterns in Bristol, with width representing number of trips and color representing the percentage of trips made by active modes (walking and cycling). The four black lines represent the interzonal OD pairs in Table 7.1.

该地图显示，市中心主导了该地区的交通模式，这暗示政策应当优先考虑这里，尽管也可以看到一些边缘的次中心。欲望线是交通系统中重要的概括性组成部分。更具体的组成部分包括节点，这些节点有特定的目的地（而不是欲望线中所表示的假设直线）。下一节将涵盖节点。

节点

在地理交通数据集中，节点是组成交通网络的主要线性特征之间的点。大致上，有两种主要类型的交通节点：

1. 不直接位于网络上的节点，例如区域质心或个体的起点和终点（如住宅和工作场所）。
2. 是交通网络一部分的节点。 从技术上讲，一个节点可以位于交通网络上的任何点，但实际上它们通常是特殊类型的顶点，例如路径（交叉口）之间的交点和进入或退出交通网络的点，如公交站和火车站[^13-transport-6]。

交通网络可以表示为图形，其中每个段通过地理线（代表边）与网络中的一个或多个其他边相连。可以通过“质心连接器”添加网络之外的节点，这些新的路线段连接到网络上附近的节点(Hollander 2016)[^13-transport-7]。网络中的每个节点然后通过一个或多个代表网络上个别段的“边”进行连接。我们将在路线网络节中看到如何将交通网络表示为图。

公共交通站点是特别重要的节点，可以表示为两种类型的节点：是道路的一部分的公交站，或由距离铁路轨道数百米的行人入口点代表的大型铁路站。我们将使用火车站来说明与布里斯托尔市增加自行车使用相关的公共交通节点。这些站由 spDataLarge 在 bristol_stations 中提供。

一个常见的障碍是，从家到工作的距离太远，无法步行或骑自行车。公共交通可以通过为进入城市的常见路线提供快速和高容量的选择来减少这一障碍。从积极出行的角度来看，长途旅行的公共交通“腿”将旅行分为三部分：

• 起点腿，通常从住宅区到公共交通站
• 公共交通腿，通常从离旅行起点最近的站到离目的地最近的站
• 目的地腿，从下车站到目的地

在欲望线节进行的分析的基础上，公共交通节点可用于为可乘坐公共汽车和（在本例中使用的）火车的旅行构建三部分欲望线。第一阶段是确定具有最多公共交通出行的欲望线，在我们的情况下这很容易，因为我们之前创建的数据集 desire_lines 已经包含了描述火车旅行数量的变量（也可以使用如 OpenTripPlanner 等公共交通路线服务来估计公共交通潜力）。为了使该方法更容易遵循，我们将仅选择关于铁路使用方面的前三条欲望线。

1

desire_rail = top_n(desire_lines, n = 3, wt = train)

现在的挑战是将每一条这样的线分解为三部分，代表通过公共交通节点的出行。这可以通过将一条欲望线转换为一个由三个线几何体组成的多线字符串对象来实现，这三个线几何体分别代表行程的起点、公共交通和终点阶段。这个操作可以分为三个阶段：矩阵创建（代表起点、终点和表示铁路站的“途经”点）、最近邻标识和转换为多线字符串。这些都由line_via()函数来完成。这个stplanr函数接受输入线和点，并返回一份欲望线的副本——具体细节请参见?line_via()。 输出与输入线相同，只不过它有新的几何列，代表通过公共交通节点的行程，如下所示：

1
2
3
4
5

ncol(desire_rail)
#> [1] 9
desire_rail = line_via(desire_rail, bristol_stations)
ncol(desire_rail)
#> [1] 12

如下图所示，初始的 desire_rail 线现在有三个额外的几何列表列，分别代表从家到起点站、从那里到目的地，以及最后从目的地站到目的地的旅行。在这种情况下，目的地段非常短（步行距离），但起点段可能足够远，以证明需要在骑行基础设施上进行投资，以鼓励人们在上班途中到居住区周围的三个起点站骑自行车，如下图所示。

Station nodes (red dots) used as intermediary points that convert straight desire lines with high rail usage (thin green lines) into three legs: to the origin station (orange) via public transport (blue) and to the destination (pink, not visible because it is so short).

路线

从地理角度看，路线是不再直线的需求线：起点和终点与旅行的需求线表示相同，但从A到B的路径更为复杂。路线的几何形状通常（但不总是）由交通网络决定。

虽然需求线仅包含两个顶点（其起点和终点），但路线可以包含任意数量的顶点，代表由直线连接的A和B之间的点：这就是linestring几何形状的定义。覆盖大距离或遵循复杂网络的路线可能有数百个顶点；基于网格或简化道路网络的路线通常较少。

路线是由需求线生成的，或更常见地，由包含代表需求线的坐标对的矩阵生成。这个路由过程是由一系列广义定义的路由引擎完成的：软件和Web服务，它们返回描述如何从起点到终点的几何形状和属性。根据它们相对于R运行的位置，路由引擎可以被分类为：

• 内存路由，使用使路线计算成为可能的R包
• 本地托管的、外部于R的路由引擎，可以从R中调用
• 由外部实体远程托管的路由引擎，提供一个可以从R中调用的Web API

在描述每一个之前，值得概述分类路由引擎的其他方式。路由引擎可以是多模式的，意味着它们可以计算由多种交通方式组成的行程，或者不是。多模式路由引擎可以返回由不同交通方式组成的多个段的结果。从住宅区到商业区的最佳路线可能涉及1）步行到最近的公交车站，2）乘公交车到离目的地最近的节点，以及3）步行到目的地，给定一组输入参数。这三段之间的转换点通常被称为“进入”和“出口”，意味着上下公共交通工具。像R5这样的多模式路由引擎比如OSRM这样的'单一模式'路由引擎更为复杂，并且有更大的输入数据要求。

多模式引擎的一个主要优势是它们能够表示由火车、公共汽车等组成的“公共交通”行程。多模型路由引擎需要代表公共交通网络的输入数据集，通常以General Transit Feed Specification（GTFS）文件形式，这些可以用tidytransit和gtfstools包（也有其他用于处理GTFS文件的包和工具）中的函数进行处理。对于专注于特定（非公共）交通方式的项目，单一模式路由引擎可能就足够了。另一种分类路由引擎（或设置）的方式是通过输出的地理级别：路线、段和环节。

路线、段和环节

路由引擎可以在三个地理级别生成输出：路线、段和环节：

• 路线级别的输出每个起点-终点对包含一个单一特征（通常是数据框表示中的多线几何和关联行），意味着每次行程的数据有一行。
• 段级别的输出在每个起点-终点对内的每个模式都包含一个单一特征和相关属性。对于仅涉及一种模式的行程（例如，从家到工作的驾驶，忽略走到汽车的短距离），段与路线相同：即汽车行程。对于涉及公共交通的行程，段提供关键信息。r5r函数detailed_itineraries()返回的段有时令人困惑地被称为“环节”。
• 环节级别的输出提供了有关路线的最详细信息，具有交通网络的每个小段的记录。通常，环节在长度上类似于，或与，OpenStreetMap中的道路相同。cyclestreets函数journey()返回在环节级别的数据，这些数据可以通过对stplanr中route()函数返回的起点和终点级别数据进行分组而聚合。

大多数路由引擎默认返回路线级别的输出，尽管多模式引擎通常在段级别提供输出（每个连续运动的单一交通模式一个特征）。 环节级别的输出有提供更多细节的优点。 cyclestreets包返回每条路线的多个“安静度”级别，使得可以识别自行车网络中的“最薄弱环节”。 环节级别输出的劣势包括增加的文件大小和与额外细节相关的复杂性。

使用函数stplanr::overline()[@morgan_travel_2020]，路线级别的结果可以转换为环节级别的结果。在使用环节或段级别的数据时，通过对表示行程起点和终点的列进行分组，并汇总/聚合包含环节级别数据的列，可以返回路线级别的统计数据。

使用R进行内存内路由

R中的路由引擎允许将存储为R对象的内存内路由网络用作路由计算的基础。选项包括sfnetworks、dodgr 和 cppRouting 包，每个包都提供了它们自己的类系统来表示路由网络，这是下一节的主题。

尽管快速和灵活，使用原生R的路由选项通常比用于现实路由计算的专用路由引擎更难设置。路由是一个复杂问题，开源路由引擎已经投入了很多时间，这些引擎可以下载并在本地主机上运行。另一方面，基于R的路由引擎可能非常适用于模型实验和网络上变化影响的统计分析。在单一语言中更改路由网络特性（或与不同路由段类型相关联的权重）、重新计算路由和分析多个场景下的结果对研究应用有益。

本地托管的专用路由引擎

本地托管的路由引擎包括OpenTripPlanner、Valhalla 和R5（它们是多模式的），以及OpenStreetMap路由机器（OSRM）（它是“单模式”的）。这些可以通过opentripplanner、valhalla、r5r 和 osrm 包从R中访问(Morgan 等，2019; Pereira 等，2021)。本地托管的路由引擎在用户的计算机上运行，但在一个与R分开的进程中。它们的优点包括执行速度和对不同运输模式的权重配置文件的控制。劣势包括在本地表示复杂网络的困难；时间动态（主要是由于交通）；以及需要专用外部软件。

远程托管的专用路由引擎

远程托管的路由引擎使用web API发送有关起点和终点的查询并返回结果。基于开源路由引擎的路由服务，例如OSRM的公开可用服务，在从R调用时与本地托管实例的工作方式相同，只需更新指定“基础URL”的参数即可。然而，由于外部路由服务托管在专用机器上（通常由有激励生成准确路由的商业公司资助），这可能给它们带来优势，包括：

• 提供全球（或通常至少在大区域内）的路由服务
• 已建立的路由服务通常会定期更新，并且经常能够响应交通水平
• 路由服务通常在专用硬件和软件上运行，包括像负载均衡器这样的系统以确保性能稳定

远程路由服务的劣势包括批量作业不可能时的速度（它们通常依赖于一条条路由的基础上通过互联网进行数据传输）、价格（例如，谷歌路由API限制了免费查询的数量）和许可问题。

googleway 和 mapbox 包通过提供对Google和Mapbox的路由服务的访问来展示这种方法。免费（但有速率限制）的路由服务包括 OSRM 和 openrouteservice.org，它们可以通过 osrm 和 openrouteservice 包从R中访问，后者不在CRAN上。还有更具体的路由服务，例如由 CycleStreets.net 提供的，这是一个自行车旅行规划器和非营利性交通技术公司“为骑自行车者，由骑自行车者”。虽然R用户可以通过 cyclestreets 包访问CycleStreets路由，但许多路由服务缺乏R接口，这代表了一个巨大的软件包开发机会：构建一个R包以提供一个接口到web API可能是一种有益的经历。

R包用于计算和导入代表交通网络上路由的数据的广泛范围是一种优势，这意味着该语言近年来越来越多地用于交通研究。然而，这种包和方法的激增的一个小缺点是有很多包和函数名称需要记住。 stplanr 包通过提供一个用于生成路由的统一接口 route() 函数来解决这个问题。该函数接受广泛的输入，包括地理欲望线（使用 l = 参数）、坐标甚至表示独特地址的文本字符串，并返回作为一致 sf 对象的路由数据。

路径规划：一个实际示例

与其对部分生成的所有期望线进行路径规划，我们专注于高度政策相关的一个子集。在尝试处理整个数据集之前，先对一个子集进行计算密集型操作通常是明智的，这同样适用于路径规划。路径规划可能需要消耗大量时间和内存，并生成大型对象，这是因为路线对象的详细几何形状和额外属性。因此，我们将在本节中过滤期望线，然后计算路径。

当骑行替代开车出行时，效益最大。短距离（大约5公里，以20公里/小时的速度可以在15分钟内骑完）有相对较高的被骑行的概率，当使用电动自行车进行出行时，最大距离会增加(Lovelace 等，2017)。以下代码片段考虑了这些因素，过滤了期望线，并返回了表示OD对（起点到终点对）的对象 desire_lines_short，在这些OD对之间有很多（100+）短距离（2.5至5公里的欧几里得距离）的车程：

1
2
3

desire_lines$distance_km = as.numeric(st_length(desire_lines)) / 1000
desire_lines_short = desire_lines |> 
  filter(car_driver >= 100, distance_km <= 5, distance_km >= 2.5)

在上面的代码中，st_length() 计算了每个期望线的长度。dplyr 的 filter() 函数基于上述标准过滤了 desire_lines 数据集。下一步是将这些期望线转换为路线。这是通过下面的代码块中使用公开可用的OSRM服务和 stplanr 的 route() 和 route_osrm() 函数来完成的：

1
2

routes_short = route(l = desire_lines_short, route_fun = route_osrm,
                     osrm.profile = "bike")

输出是 routes_short，一个代表适用于骑自行车的交通网络上的路线的 sf 对象（至少根据 OSRM 路由引擎是这样），每个期望线对应一个。注意：像上面的命令中那样调用外部路由引擎只有在有互联网连接的情况下才能工作（有时还需要存储在环境变量中的 API 密钥，尽管在这种情况下不需要）。除了 desire_lines 对象中包含的列之外，新的路线数据集还包含 distance（这次是指路线距离）和 duration（以秒为单位）列，这些提供了有关每条路线性质的潜在有用额外信息。我们将绘制沿着这些线路进行许多短途汽车行驶的期望线和骑行路线。通过使路线的宽度与可能被替换的汽车行驶数量成比例，提供了一种有效的方式来优先考虑对道路网络进行干预(Lovelace 等 2017)。 下面的代码块绘制了期望线和路线，该图显示了人们驾驶短距离的路线：⁷

Routes along which many (100+) short (<5km Euclidean distance) car journeys are made (red) overlaying desire lines representing the same trips (black) and zone centroids (dots).

通过在交互式地图上绘制结果，例如使用 mapview::mapview(st_geometry(routes_short))，可以看到许多短途汽车行程发生在布拉德利斯托克（Bradley Stoke）及其周围，距离布里斯托尔（Bristol）市中心大约10公里以北。找出该地区高度依赖汽车的原因很容易：根据维基百科，布拉德利斯托克是“欧洲最大的由私人投资建设的新城镇”，这暗示了有限的公共交通供应。此外，该城镇被包括M4和M5高速公路在内的大型（对自行车不友好）的道路结构所环绕（tallon 2007)。

将旅行需求线转换为路线有很多好处。重要的是要记住，我们不能确定有多少（如果有的话）行程会沿着路由引擎计算出的确切路线进行。然而，路线和街道/道路/段级别的结果可能具有很高的政策相关性。路段结果可以根据可用数据(Lovelace 等 2017)，实现在最需要的地方优先进行投资。

路线网络

虽然路线通常包含有关旅行行为的数据，在与需求线和OD（起点-终点）对相同的地理层级上，路线网络数据集通常代表物理交通网络。路线网络中的每个段大致对应于交叉口之间的连续街道段，并且只出现一次，尽管段的平均长度取决于数据源（本节中使用的由OSM（开放街道地图）派生的bristol_ways数据集中的段平均长度刚好超过200米，标准差接近500米）。段长度的可变性可以由这样一个事实来解释：在一些农村地区，交叉口相距很远，而在密集的城市区域，每隔几米就有交叉口和其他段的断点。

路线网络可以是交通数据分析项目的输入或输出，或者两者都是。任何涉及路线计算的交通研究都需要内部或外部路由引擎中的路线网络数据集（在后一种情况下，路线网络数据不一定需要导入到R中）。然而，路线网络在许多交通研究项目中也是重要的输出：总结诸如特定段上可能进行的潜在行程数量的数据，并以路线网络的形式表示，可以帮助在最需要的地方优先进行投资。

为了演示如何从路线级别的数据中创建路线网络作为输出，想象一个简单的模式转移场景。假设0到3公里的路线距离之间的50%的汽车行程被自行车取代，这一百分比每增加1公里的路线距离就减少10个百分点，以至于6公里的汽车行程有20%被自行车取代，而8公里或更长的汽车行程没有被自行车取代。当然，这是一个不现实的情景(Lovelace 等 2017)，但它是一个有用的起点。 在这种情况下，我们可以按照以下方式模拟从汽车到自行车的模式转移：

1
2
3
4
5
6
7
8
9
10
11
12
13

uptake = function(x) {
  case_when(
    x <= 3 ~ 0.5,
    x >= 8 ~ 0,
    TRUE ~ (8 - x) / (8 - 3) * 0.5
  )
}
routes_short_scenario = routes_short |> 
  mutate(uptake = uptake(distance / 1000)) |> 
  mutate(bicycle = bicycle + car_driver * uptake,
         car_driver = car_driver * (1 - uptake))
sum(routes_short_scenario$bicycle) - sum(routes_short$bicycle)
#> [1] 3733

在创建了一个大约有4000次行程从驾驶转向骑自行车的场景后，我们现在可以模拟这个更新后的模拟自行车活动将在哪里发生。为此，我们将使用stplanr包中的overline()函数。该函数在交叉口（两条或更多的折线几何相交的地方）处断开折线，并计算每个唯一路段的聚合统计信息(Morgan 和 Lovelace 2020)，它将包含路线的对象和要汇总的属性的名称作为第一和第二个参数：

1

route_network_scenario = overline(routes_short_scenario, attrib = "bicycle")

前面两个代码块的输出总结如下面的图所示。

Illustration of the percentage of car trips switching to cycling as a function of distance (left) and route network level results of this function (right).

具有在路段级别的记录的交通网络，通常带有诸如道路类型和宽度之类的属性，构成了一种常见的路线网络。这种路线网络数据集在全球范围内都可以从OpenStreetMap获得，并可以使用诸如osmdata和osmextract这样的软件包进行下载。为了节省下载和准备OSM的时间，我们将使用来自spDataLarge包的bristol_ways对象，这是一个带有LINESTRING几何图形和属性的sf对象，代表了案例研究区域内交通网络的一个样本（详见?bristol_ways以获取详细信息），如下面的输出所示：

1
2
3
4
5

summary(bristol_ways)
#>      highway       maxspeed             ref                     geometry   
#>  cycleway:1721   Length:6160        Length:6160        LINESTRING   :6160  
#>  rail    :1017   Class :character   Class :character   epsg:4326    :   0  
#>  road    :3422   Mode  :character   Mode  :character   +proj=long...:   0

输出显示bristol_ways代表了交通网络上刚刚超过6千个路段。这种和其他地理网络可以表示为数学图，网络上的节点由边连接。一些R软件包已经为处理这样的图而开发，尤其是igraph。您可以手动将路线网络转换为igraph对象，但地理属性将会丢失。为了克服igraph的这一局限性，开发了sfnetworks软件包，该软件包能够同时将路线网络表示为图和地理线。我们将在bristol_ways对象上演示sfnetworks的功能。

1
2
3
4

bristol_ways$lengths = st_length(bristol_ways)
ways_sfn = as_sfnetwork(bristol_ways)
class(ways_sfn)
#> [1] "sfnetwork" "tbl_graph" "igraph"

1
2
3
4
5
6
7
8
9

ways_sfn
#> # A sfnetwork with 5728 nodes and 4915 edges
#> # A directed multigraph with 1013 components with spatially explicit edges
#> # Node Data:     5,728 × 1 (active)
#> # Edge Data:     4,915 × 7
#>    from    to highway maxspeed ref                              geometry lengths
#>   <int> <int> <fct>   <fct>    <fct>                    <LINESTRING [°]>     [m]
#> 1     1     2 road    <NA>     B3130 (-2.61 51.4, -2.61 51.4, -2.61 51.…    218.
#> # … 

由于空间考虑，前一个代码块的输出（最后的输出被缩短，仅包含最重要的8行）显示ways_sfn是一个复合对象，以图和空间形式包含节点和边。ways_sfn是sfnetwork类，该类基于igraph包中的igraph类。在下面的示例中，计算了'边介数'，也就是经过每条边的最短路径的数量（有关更多详细信息，请参见?igraph::betweenness）。边介数计算的输出显示在下图，该图以用overline()函数计算的自行车路线网络数据集作为叠加层进行比较。结果表明，每个图的边代表一个路段：靠近道路网络中心的路段具有最高的介数值，而靠近布里斯托尔中心的路段基于这些简单数据集具有更高的自行车潜力。

1
2
3

ways_centrality = ways_sfn |> 
  activate("edges") |>  
  mutate(betweenness = tidygraph::centrality_edge_betweenness(lengths)) 

Illustration of route network datasets. The grey lines represent a simplified road network, with segment thickness proportional to betweenness. The green lines represent potential cycling flows (one way) calculated with the code above

人们还可以使用sfnetworks包中的这种路由网络的图表示法来找到起点和终点之间的最短路线。与可能的情况相比，本节介绍的方法相对简单。sfnetworks开放的双重图/空间功能使得许多新的强大技术得以实现，这些在本节中无法完全涵盖。然而，本节确实为进一步探索和研究该领域提供了一个坚实的起点。最后一点是，我们上面使用的示例数据集相对较小。也可能值得考虑如何将工作适应到更大的网络：在数据的子集上测试方法，并确保您有足够的RAM将有所帮助，尽管也值得探索其他优化用于大型网络的交通网络分析工具，比如R5 (Alessandretti 等 2022)。

优先考虑新基础设施

本节演示了地理计算如何在交通规划领域创造与政策相关的成果。我们将使用一个简单的方法（出于教育目的）来识别可持续交通基础设施投资的有希望的地点。

本章概述的数据驱动方法的一个优点是其模块性：每个方面本身就很有用，并可以融入更广泛的分析中。使我们达到这一阶段的步骤包括在路线节中识别短途但依赖汽车的通勤路线（由欲望线生成），以及使用sfnetworks包在路线网络中分析路由网络特性。本章的最后一个代码块将这些分析纵向结合，通过在代表距离自行车基础设施很近的区域的新数据集上覆盖上一节中的自行车潜力估计。

1
2
3
4

existing_cycleways_buffer = bristol_ways |> 
  filter(highway == "cycleway") |>    # 1) filter out cycleways
  st_union() |>                       # 2) unite geometries
  st_buffer(dist = 100)               # 3) create buffer

下一步是创建一个数据集，该数据集代表网络上存在高自行车潜力但几乎没有为自行车提供便利设施的点：

1
2
3
4
5

route_network_no_infra = st_difference(
  route_network_scenario,
  route_network_scenario |> st_set_crs(st_crs(existing_cycleways_buffer)),
  existing_cycleways_buffer
)

1
2

#> Warning: attribute variables are assumed to be spatially constant throughout
#> all geometries

前面代码块的结果显示在下图中，该图显示了具有高度依赖汽车、具有高自行车潜力但没有自行车道的路线。

1
2
3

tmap_mode("view")
qtm(route_network_no_infra, basemaps = leaflet::providers$Esri.WorldTopoMap,
    lines.lwd = 5)

1
2

#> Warning: Some legend items or map compoments do not fit well (e.g. due to the
#> specified font size).

Potential routes along which to prioritise cycle infrastructure in Bristol to reduce car dependency. The static map provides an overview of the overlay between existing infrastructure and routes with high car-bike switching potential (left). The screenshot the interactive map generated from the qtm() function highlights Whiteladies Road as somewhere that would benefit from a new cycleway (right).

该方法存在一些局限性：实际上，人们并不总是前往区域中心或始终使用特定模式的最短路径算法。然而，结果展示了如何使用地理数据分析来突出可能特别有益于新投资自行车道的地方，尽管这种方法很简单。为了在实践中指导交通规划设计，这种分析需要大幅扩展——包括使用更大的输入数据集。

旅行的未来方向

这一章提供了使用地理计算进行交通研究的可能性的一个概览，并使用开放数据和可复制代码探讨了构成城市交通系统的一些关键地理元素。这些结果有助于规划需要投资的地方。

交通系统在多个相互作用的层面上运行，这意味着地理计算方法有很大的潜力，可以生成有关它们如何运作以及不同干预措施可能产生的影响的见解。在这一领域还有更多可以做的事情：基于本章提供的基础，有很多方向可以进行拓展。

交通是许多国家温室气体排放增长最快的来源，并将成为“尤其是在发达国家中最大的温室气体排放部门”。交通相关的排放在社会上分布不均，但（与食物和供暖不同）对幸福来说并非必需品。通过需求减少、车队电气化以及积极采取像步行和骑自行车这样的活动出行方式，该部门有很大的快速减碳潜力。新技术可以通过更多的汽车共享来减少对汽车的依赖。诸如无桩自行车和电动滑板车方案这样的“微移动”系统也正在出现，以General Bikeshare Feed Specification（GBFS）格式创建有价值的数据集，可以用gbfs包进行导入和处理。

从方法论的角度看，本章提供的基础可以通过在分析中包括更多变量来进行扩展。例如，通过使用R语言的统计建模能力，这可以用来预测当前和未来的自行车使用水平。

这种类型的分析是Propensity to Cycle Tool（PCT）的基础，这是一个用R语言开发的，用于优先考虑英格兰各地自行车投资的公开可访问的映射工具（lovelace 2017）。类似的工具可以用于鼓励与空气污染和公共交通准入等其他主题相关的基于证据的交通政策。

练习

E1. 在本章的大部分分析中，我们关注的是主动模式，但是驾驶行程呢？

• 在desire_lines对象中，有多大比例的行程是通过驾驶完成的？
• 有多大比例的desire_lines的直线长度为5公里或更远？
• 在长度超过5公里的desire_lines中，有多大比例的行程是通过驾驶完成的？
• 绘制那些长度小于5公里且其中超过50%的行程是通过汽车完成的desire_lines。
• 你注意到这些依赖汽车但desire_lines较短的地方有什么特点？

E2. 如果在离现有自行车道超过100米的区段上建造了上图中展示的所有路线，会增加多少长度的自行车道？

E3. desire_lines中表示的行程有多大比例在routes_short_scenario对象中得到了体现？

• 奖励：所有行程中有多大比例发生在穿越routes_short_scenario的desire_lines上？

E4. 本章展示的分析旨在教授如何将地理计算方法应用于交通研究。 如果你是在政府或交通咨询公司中真正进行这项工作，你会有哪三个不同的做法？

E5. 显然，上图中识别的路线只是部分情况。 你会如何扩展这个分析？

E6. 假设你想通过创建投资基于场所的自行车政策的关键区域（而不是路线）来扩展场景，比如无车区、自行车停车点和减少汽车停车策略。 栅格数据集如何协助完成这项工作？

• 奖励：开发一个光栅层，将布里斯托地区划分为100个单元格（10x10），并从bristol_ways数据集（参见章节 @ref(location)）估算每个单元格中道路的平均速度限制。

本文由SCY翻译自《Geocomputation with R》第十二章

前提条件

本章假设您已经具备地理数据分析的熟练技能，例如通过学习并完成章节空间数据至地理数据重投影中的练习来获得。强烈推荐熟悉广义线性模型（GLM）和机器学习。

本章使用以下的R包：¹

1
2
3
4
5
6
7
8
9
10
11
12

library(sf)
library(terra)
library(dplyr)
library(future)
library(lgr)
library(mlr3)
library(mlr3learners)
library(mlr3extralearners)
library(mlr3spatiotempcv)
library(mlr3tuning)
library(mlr3viz)
library(progressr)

当然，所需数据将在适当的时候附上。

引言

统计学习主要关注使用统计和计算模型来识别数据中的模式，并根据这些模式进行预测。由于其起源，统计学习是 R 语言的一大优势（见章节地理计算软件）。²。

本章重点介绍有训练数据集的监督技术，而非像聚类这样的无监督技术。响应变量可以是二元的（如山体滑坡发生与否）、分类的（土地利用）、整数（物种丰富度计数）或数值的（以 pH 衡量的土壤酸度）。监督技术用于建模——响应与一个或多个预测变量之间的关系——这些响应是基于一组观测样本而已知的。

大多数机器学习研究的主要目的是进行良好的预测。机器学习在“大数据”时代蓬勃发展，因为其方法对输入变量几乎没有假设，并且能处理巨大的数据集。机器学习适用于预测未来客户行为、推荐服务（音乐、电影、接下来购买什么）、面部识别、自动驾驶、文本分类和预测性维护（基础设施、产业）等任务。

本章基于一个案例研究：（空间）预测山体滑坡。 该应用与定义在引言章中的地理计算的应用性有关，并说明了机器学习在唯一目的是预测时如何借鉴统计学领域。因此，本章首先使用广义线性模型介绍建模和交叉验证的概念。在此基础上，该章实施了一个更典型的机器学习算法，即支持向量机（SVM）。模型的预测性能将使用空间交叉验证（CV）进行评估，该方法考虑到地理数据是特殊的。

CV确定模型泛化到新数据的能力，通过将数据集（反复）分割为训练集和测试集。它使用训练数据来拟合模型，并检查其在对测试数据进行预测时的性能。CV 有助于检测过拟合，因为过于紧密地预测训练数据（噪声）的模型通常会在测试数据上表现不佳。

随机分割空间数据可能导致训练点与测试点在空间上是邻近的。由于空间自相关，在这种情况下，测试数据集和训练数据集将不是独立的，因此交叉验证无法检测出可能的过拟合。空间交叉验证缓解了这个问题，并且是本章的核心主题。

案例: 滑坡敏感性

本案例研究基于南厄瓜多尔滑坡地点的一个数据集，如下图所示，并在@muenchow_geomorphic_2012中详细描述。该论文中使用的数据集的一个子集提供在spDataLarge包中，可以如下加载：

1
2

data("lsl", "study_mask", package = "spDataLarge")
ta = terra::rast(system.file("raster/ta.tif", package = "spDataLarge"))

T上面的代码加载了三个对象：一个名为lsl的data.frame，一个名为study_mask的sf对象，以及一个名为ta的SpatRaster，其中包含地形属性栅格。lsl包含一个因子列lslpts，其中TRUE对应于观察到的滑坡'发起点'，坐标存储在列x和y中。³ 有175个滑坡点和175个非滑坡点，如summary(lsl$lslpts)所示。这175个非滑坡点是从研究区域随机抽样的，但必须落在滑坡多边形周围的小缓冲区之外。

Landslide initiation points (red) and points unaffected by landsliding (blue) in Southern Ecuador.

lsl的前三行，四舍五入为两位有效数字，可以在下表中找到。

要对滑坡易发性进行建模，我们需要一些预测因子。由于地形属性经常与滑坡有关[@muenchow_geomorphic_2012]，我们已经从ta提取到lsl中以下的地形属性：

• slope：坡度角（°）
• cplan：平面曲率（rad m⁻¹）表示坡度的汇聚或发散，从而表示水流
• cprof：剖面曲率（rad m^-1）作为流速加速的衡量标准，也称为坡度角的下坡变化
• elev：海拔高度（m a.s.l.）作为研究区不同海拔带植被和降水的表示
• log10_carea：集水区面积（log10 m²）的常用对数，表示流向某个位置的水量

使用R-GIS桥计算地形属性并将其提取到滑坡点（参见本章末尾的练习部分）可能是一个值得尝试的练习。

在R中的传统建模方法

在介绍mlr3包之前，这是一个提供统一接口以访问数十种学习算法的 umbrella 包，值得先看一下R传统的建模接口。这一介绍到有监督的统计学习为进行空间CV提供了基础，并有助于更好地理解随后呈现的mlr3方法。

监督学习涉及将响应变量作为预测因子的函数进行预测。在R中，建模函数通常使用公式来指定（有关R公式的更多详细信息，请参阅?formula）。以下命令指定并运行一个广义线性模型：

1
2
3

fit = glm(lslpts ~ slope + cplan + cprof + elev + log10_carea,
          family = binomial(),
          data = lsl)

理解这三个输入参数各自的意义是非常有价值的：

• 一个公式，用于指定由预测因子决定的滑坡发生（lslpts）
• 一个家族（family），用于指定模型的类型，在这个例子中是binomial，因为响应是二元的（见?family）
• 包含响应和预测因子（作为列）的数据框

可以如下打印这个模型的结果（summary(fit)提供了关于结果的更详细的描述）：

1
2
3
4
5
6
7
8
9
10
11
12
13
14

class(fit)
#> [1] "glm" "lm"
fit
#> 
#> Call:  glm(formula = lslpts ~ slope + cplan + cprof + elev + log10_carea, 
#>     family = binomial(), data = lsl)
#> 
#> Coefficients:
#> (Intercept)        slope        cplan        cprof         elev  log10_carea  
#>    2.511364     0.079011   -28.941961   -17.563601     0.000179    -2.274877  
#> 
#> Degrees of Freedom: 349 Total (i.e. Null);  344 Residual
#> Null Deviance:	    485 
#> Residual Deviance: 373 	AIC: 385

模型对象fit，其类别为glm，包含了定义响应和预测因子之间拟合关系的系数。这个对象也可以用于预测。这是通过通用的predict()方法来完成的，在这个案例中，它调用了predict.glm()函数。将type设置为response会返回lsl中每个观测的预测概率（即滑坡发生的概率），如下图所示（见?predict.glm）：

1
2
3
4

pred_glm = predict(object = fit, type = "response")
head(pred_glm)
#>      1      2      3      4      5      6 
#> 0.1901 0.1172 0.0952 0.2503 0.3382 0.1575

通过将系数应用于预测因子的栅格数据，可以创建空间分布图。这可以手动完成，也可以使用terra::predict()来完成。除了一个模型对象（fit）外，这个函数还需要一个具有预测因子（栅格层）的SpatRaster，这些预测因子的名称应与模型输入数据框中的名称相同（见下图）。

1
2

# making the prediction
pred = terra::predict(ta, model = fit, type = "response")

Spatial distribution mapping of landslide susceptibility using a GLM.

在进行预测时，我们忽略了空间自相关，因为我们假设平均预测准确性在有或没有空间自相关结构的情况下保持不变。然而，也可以将空间自相关结构纳入模型和预测中。 尽管这超出了本书的范围，但我们为感兴趣的读者提供了一些查阅的方向：

1. 回归克里金（Regression Kriging）的预测结合了回归的预测和回归残差的克里金。
2. 人们还可以向广义最小二乘模型[nlme::gls()]添加空间相关（依赖）结构。
3. 人们也可以使用混合效应建模方法。

基本上，随机效应对响应变量施加了依赖结构，从而使某一类的观察结果与另一类的观察结果更相似。类别可以是蜜蜂巢、猫头鹰巢、植被样带或高度分层。这种混合建模方法假设随机截距是正态和独立分布的。这甚至可以通过使用正态和空间依赖的随机截距来扩展。然而，为此，你可能需要采用贝叶斯建模方法，因为频率论软件工具在这方面尤其有限，尤其是对于更复杂的模型。

空间分布图是模型（上图）的一个非常重要的结果。更重要的是，基础模型在制作这些图时有多好，因为如果模型的预测性能差，预测图就毫无用处。评估二项模型预测性能的最流行指标是接收器操作特性曲线下的面积（AUROC）。这是一个介于0.5和1.0之间的值，其中0.5表示模型不比随机模型好，1.0表示两个类的完美预测。因此，AUROC越高，模型的预测能力越好。以下代码块使用roc()计算模型的AUROC值，该函数将响应和预测值作为输入。auc()返回曲线下的面积。

1
2

pROC::auc(pROC::roc(lsl$lslpts, fitted(fit)))
#> Area under the curve: 0.8216

AUROC 值为 0.82 代表了一个良好的拟合。然而，这是一个过于乐观的估计，因为我们是在完整数据集上进行的计算。为了得到一个偏差较小的评估，我们需要使用交叉验证，并且在处理空间数据的情况下，应使用空间CV。

介绍（空间）交叉验证

交叉验证属于重采样方法的家族 [@james_introduction_2013]。基本思想是将数据集（反复）分为训练集和测试集，其中训练数据用于拟合模型，然后将该模型应用于测试集。通过使用如二项式案例中的AUROC这样的性能度量，将预测值与测试集中已知的响应值进行比较，从而得到模型将学习到的关系泛化到独立数据的能力的偏差较小的评估。例如，100次重复的5折交叉验证意味着将数据随机分为五个部分（折叠），每个折叠都被用作一次测试集（见下图的上行）。这保证了每个观测值都在测试集中使用一次，并需要拟合五个模型。随后，这一过程重复100次。当然，每次重复中的数据分割都会有所不同。总体而言，这合计为500个模型，而所有模型的平均性能度量（AUROC）即为模型的整体预测能力。

然而，地理数据是特殊的。正如我们将在运输章节中看到的，地理学的“第一定律”指出，彼此靠近的点通常比远离的点更相似 [@miller_tobler_2004]。这意味着这些点在统计上并不独立，因为在传统的CV中，训练点和测试点通常过于靠近（见下图的第一行）。靠近“测试”观测值的“训练”观测值可能提供一种“偷窥预览”：这些信息本应对训练数据集不可用。为了缓解这一问题，“空间划分”被用于将观测值分为空间上不相交的子集（使用观测值的坐标进行k-均值聚类；@brenning_spatial_2012；见下图的第二行）。这种划分策略是空间CV和传统CV之间的唯一差异。因此，空间CV导致了对模型预测性能的偏差较小的评估，从而有助于避免过拟合。

Spatial visualization of selected test and training observations for cross-validation of one repetition. Random (upper row) and spatial partitioning (lower row).

使用mlr3进行空间交叉验证

如CRAN 机器学习任务视图所描述的，有数十个用于统计学的包。熟悉其中每一个包，包括如何进行交叉验证和超参数调整，可能是一个耗时的过程。比较来自不同包的模型结果可能更为繁琐。mlr3 包和生态系统是为了解决这些问题而开发的。它充当一个“元包”，为包括分类、回归、生存分析和聚类在内的流行的有监督和无监督统计学习技术提供统一的接口 [@lang_mlr3_2019; @becker_mlr3_2022]。标准化的 mlr3 接口基于八个“构建块”。 如下图所示，这些构建块有明确的顺序。

(ref:building-blocks) mlr3包的基础构建块。来源：@becker_mlr3_2022。（得到了友好地重新使用此图的许可。）

(ref:building-blocks)

mlr3 建模过程包含三个主要阶段。 首先，任务（task）指定了数据（包括响应变量和预测变量）和模型类型（如回归或分类）。 其次，学习器（learner） 定义了应用于创建的任务的特定学习算法。 第三，重采样（resampling）方法评估了模型的预测性能，即其泛化到新数据的能力。

广义线性模型

要在 mlr3中实现GLM，我们必须创建一个包含滑坡数据的任务（task）。由于响应是二元（两类变量）且具有空间维度，我们使用 mlr3spatiotempcv 包的TaskClassifST$new()创建一个分类任务[@schratz_mlr3spatiotempcv_2021，对于非空间任务，使用 `mlr3::TaskClassif$new()` 或对于回归任务使用`mlr3::TaskRegr$new()`，其他任务类型请参见 `?Task`]。⁴ 这些 Task*$new() 函数的第一个重要参数是 backend。backend 期望输入数据包括响应变量和预测变量。target 参数指示响应变量的名称（在我们的例子中是 lslpts），而 positive 决定响应变量的两个因子水平中哪一个表示滑坡起点（在我们的例子中是 TRUE）。lsl 数据集的所有其他变量将作为预测因子。对于空间交叉验证，我们需要提供一些额外的参数。coordinate_names 参数期望坐标列的名称（参见空间交叉验证章节和图3）。另外，我们应该指示使用的坐标参考系统（crs）并决定是否希望在建模中使用坐标作为预测因子（coords_as_features）。

1
2
3
4
5
6
7
8
9
10

# create task
task = mlr3spatiotempcv::TaskClassifST$new(
  id = "ecuador_lsl",
  backend = mlr3::as_data_backend(lsl), 
  target = "lslpts", 
  positive = "TRUE",
  coordinate_names = c("x", "y"),
  coords_as_features = FALSE,
  crs = "EPSG:32717"
  )

注意，mlr3spatiotempcv::as_task_classif_st() 也接受一个 sf-对象作为 backend 参数的输入。在这种情况下，你可能只想额外指定 coords_as_features 参数。我们没有将 lsl 转换为一个 sf-对象，因为 TaskClassifST$new() 在背景中会将其转换回一个非空间的 data.table 对象。对于简短的数据探索，mlr3viz 包的 autoplot() 函数可能会很方便，因为它会绘制响应与所有预测因子的关系，以及所有预测因子之间的关系（未显示）。

1
2
3
4

# plot response against each predictor
mlr3viz::autoplot(task, type = "duo")
# plot all variables against each other
mlr3viz::autoplot(task, type = "pairs")

创建了任务之后，我们需要选择一个学习器（learner），以确定要使用的统计学方法。所有分类的学习器以 classif. 开头，所有回归的学习器以 regr. 开头（详见 ?Learner）。mlr3extralearners::list_mlr3learners() 列出了所有可用的学习器以及 mlr3 从哪个软件包中导入了它们（表 @ref(tab:lrns)）。要了解能够对二元响应变量进行建模的学习器，我们可以运行：

1
2
3
4

mlr3extralearners::list_mlr3learners(
  filter = list(class = "classif", properties = "twoclass"), 
  select = c("id", "mlr3_package", "required_packages")) |>
  head()

这将产生所有能够处理两类问题（有或没有滑坡）的学习器。我们选择在传统建模节中使用的二项分类方法，该方法在 mlr3learners 中实现为 classif.log_reg。另外，我们需要指定 predict.type，它决定了预测的类型，prob 会得出滑坡发生概率在0到1之间的预测值（这对应于 predict.glm 中的 type = response）。

1

learner = mlr3::lrn("classif.log_reg", predict_type = "prob")

要访问学习器的帮助页面并找出它来自哪个包，我们可以运行：

1

learner$help()

1
2

# performance estimation level
perf_level = mlr3::rsmp("repeated_spcv_coords", folds = 5, repeats = 100)

请注意，这与节glm中用于GLM的重采样所用的代码完全相同；我们只是在这里重复它作为提醒。

到目前为止，这个过程与节glm中描述的过程相同。然而，下一步是新的：调整超参数。使用相同的数据进行性能评估和调优可能会导致过于乐观的结果[@cawley_overfitting_2010]。这可以通过使用嵌套的空间交叉验证（CV）来避免。

Schematic of hyperparameter tuning and performance estimation levels in CV. [Figure was taken from Schratz et al. (2019). Permission to reuse it was kindly granted.]

这意味着我们将每个折叠再次分成五个空间上不重叠的子折叠，这些子折叠用于确定最佳的超参数（下面代码块中的tune_level对象；见上图以获取视觉表示）。为了找到最优的超参数组合，我们在这些子折叠中分别拟合50个模型（下面代码块中的terminator对象），并为超参数C和Sigma随机选择值。C和Sigma的随机选择值还受到预定义调优空间（search_space对象）的额外限制。调优空间的范围是根据文献中推荐的值来选择的[@schratz_hyperparameter_2019]。

1
2
3
4
5
6
7
8
9
10

# five spatially disjoint partitions
tune_level = mlr3::rsmp("spcv_coords", folds = 5)
# use 50 randomly selected hyperparameters
terminator = mlr3tuning::trm("evals", n_evals = 50)
tuner = mlr3tuning::tnr("random_search")
# define the outer limits of the randomly selected hyperparameters
search_space = paradox::ps(
  C = paradox::p_dbl(lower = -12, upper = 15, trafo = function(x) 2^x),
  sigma = paradox::p_dbl(lower = -15, upper = 6, trafo = function(x) 2^x)
)

下一步是根据所有定义超参数调优的特性，使用AutoTuner$new()来修改学习器lrn_ksvm。

1
2
3
4
5
6
7
8

at_ksvm = mlr3tuning::AutoTuner$new(
  learner = lrn_ksvm,
  resampling = tune_level,
  measure = mlr3::msr("classif.auc"),
  search_space = search_space,
  terminator = terminator,
  tuner = tuner
)

现在调优设置已完成，用于确定一个折叠的最佳超参数将拟合250个模型。重复这个过程针对每个折叠，我们最终得到125,000（250 * 5）个模型用于每次重复。重复100次意味着拟合总共125,000个模型以确定最佳超参数（见图）。这些将用于性能评估，这需要拟合另外500个模型（5个折叠 * 100次重复；参见图。 为了更清晰地描述性能评估处理链，让我们记录下我们给计算机的命令：

1. 性能层次（图的左上部分）- 将数据集分成五个空间上不交叠的（外部）子折叠。
2. 调优层次（图的左下部分）- 使用性能层次的第一个折叠，并再次将其空间上分成五个（内部）子折叠进行超参数调优。在这些内部子折叠中使用50个随机选定的超参数，即拟合250个模型。
3. 性能估计 - 使用上一步（调优层次）中的最佳超参数组合，并将其应用于性能层次中的第一个外部折叠以估计性能（AUROC）。
4. 重复步骤2和3，针对其余四个外部折叠。
5. 重复步骤2至4，共100次。

超参数调优和性能估计的过程在计算上是非常密集的。为了减少模型运行时间，mlr3提供了使用future包进行并行化的可能性。由于我们即将进行嵌套交叉验证，我们可以决定是否希望并行化内部循环或外部循环（见图的左下部分）。由于前者将运行125,000个模型，而后者仅运行500个，很明显我们应该并行化内部循环。要设置内部循环的并行化，我们执行：

1
2
3
4

library(future)
# execute the outer loop sequentially and parallelize the inner loop
future::plan(list("sequential", "multisession"), 
             workers = floor(availableCores() / 2))

确实，我们指示future包仅使用可用核心的一半而非全部（默认设置），这样做是为了在使用高性能计算集群的情况下允许其他可能的用户也能在同一个集群上工作。

现在，我们已经为执行嵌套空间交叉验证做好了准备。指定resample()参数的过程与使用GLM时完全相同，唯一的区别是store_models和encapsulate参数。将前者设置为TRUE将允许我们提取超参数调优结果，这在我们计划对调优进行后续分析时非常重要。后者确保即使其中一个模型出错，处理过程也会继续。这避免了由于一个失败的模型而导致整个过程停止，这在大规模模型运行中是非常需要的。一旦处理完成，可以查看失败的模型。处理完成后，最好使用future:::ClusterRegistry("stop")明确地停止并行化。最后，我们将输出对象（result）保存到磁盘，以防我们希望在另一个R会话中使用它。

在运行以下代码之前，请注意，由于它将执行包括125,500个模型在内的空间交叉验证，所以这是一个耗时的过程。在现代笔记本电脑上轻松运行半天也是可能的。需要注意的是，运行时间取决于多个方面：CPU速度、所选算法、所选核心数以及数据集。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

progressr::with_progress(expr = {
  rr_spcv_svm = mlr3::resample(task = task,
                               learner = at_ksvm, 
                               # outer resampling (performance level)
                               resampling = perf_level,
                               store_models = FALSE,
                               encapsulate = "evaluate")
})

# stop parallelization
future:::ClusterRegistry("stop")
# compute the AUROC values
score_spcv_svm = rr_spcv_svm$score(measure = mlr3::msr("classif.auc")) 
# keep only the columns you need
score_spcv_svm = score_spcv_svm[, .(task_id, learner_id, resampling_id, classif.auc)]

如果你不想在本地运行代码，我们已经将 score_svm 保存在本书的GitHub仓库中。你可以按照以下方式加载它们：

1
2
3

score = readRDS("extdata/12-bmr_score.rds")
score_spcv_svm = score[learner_id == "classif.ksvm.tuned" & 
                         resampling_id == "repeated_spcv_coords"]

让我们看一下最终的AUROC：模型区分两个类别的能力。

1
2
3

# final mean AUROC
round(mean(score_spcv_svm$classif.auc), 2)
#> [1] 0.74

在这个特定情况下，似乎GLM（汇总的AUROC是 r score[resampling_id == "repeated_spcv_coords" & learner_id == "classif.log_reg", round(mean(classif.auc), 2)]）比SVM略好一些。为了保证一个绝对公平的比较，还应确保两个模型使用完全相同的划分——这是我们在这里没有展示，但在背后默默使用的（更多信息请参见本书GitHub仓库中的code/12_cv.R）。为此，mlr3提供了benchmark_grid()和benchmark()函数[参见https://mlr3book.mlr-org.com/perf-eval-cmp.html#benchmarking，@becker_mlr3_2022]。我们将在练习中更详细地探讨这些函数。还请注意，在SVM的随机搜索中使用超过50次迭代可能会产生导致模型具有更好AUROC的超参数 [@schratz_hyperparameter_2019]。另一方面，增加随机搜索迭代次数也会增加总模型数量，从而增加运行时间。

到目前为止，空间CV已被用于评估学习算法泛化到未见数据的能力。对于预测性制图目的，人们会在完整的数据集上调整超参数。这将在生态章中讨论。

结论

重采样方法是数据科学家工具箱中的一个重要组成部分[@james_introduction_2013]。本章使用交叉验证来评估各种模型的预测性能。具有空间坐标的观测值可能由于空间自相关而不是统计上独立的，这违反了交叉验证的一个基本假设。空间CV通过减少空间自相引入的偏见来解决这个问题。

mlr3软件包便于使用（空间）重采样技术，结合最受欢迎的统计学习技术，包括线性回归，如广义加性模型之类的半参数模型，以及如随机森林，SVMs和增强回归树[@bischl_mlr:_2016;@schratz_hyperparameter_2019]等机器学习技术。机器学习算法通常需要超参数输入，其中最优的“调整”可能需要数千次模型运行，这需要大量的计算资源，消耗大量时间、RAM和/或核心。mlr3通过启用并行化来解决这个问题。

总体而言，机器学习及其用于理解空间数据是一个大领域，本章提供了基础知识，但还有更多需要学习。我们推荐以下方向的资源：

• mlr3 book [@becker_mlr3_2022; https://mlr-org.github.io/mlr-tutorial/release/html/]，特别是关于处理时空数据的章节
• 关于超参数调整的学术论文[@schratz_hyperparameter_2019]
• 关于如何使用mlr3spatiotempcv的学术论文[@schratz_mlr3spatiotempcv_2021]
• 在处理时空数据的情况下，应当在进行CV时考虑空间和时间自相关[@meyer_improving_2018]

练习

E1. 使用R-GIS桥梁（参见GIS软件桥梁章节），从通过terra::rast(system.file("raster/ta.tif", package = "spDataLarge"))$elev加载的elev数据集中计算以下地形属性：

• 坡度
• 平面曲率
• 轮廓曲率
• 流域面积

E2. 从相应的输出栅格中提取值到lsl数据框（通过data("lsl", package = "spDataLarge")），并添加名为slope、cplan、cprof、elev和log_carea的新变量。

E3. 结合GLM使用派生的地形属性栅格，制作一个与图12.2中显示的相似的空间预测地图。 运行data("study_mask", package = "spDataLarge")将学习区域的掩码附加上。

E4. 基于GLM学习器计算100次重复的5折非空间交叉验证和空间CV，并利用箱线图比较两种重采样策略的AUROC值。

提示：你需要指定一个非空间重采样策略。

另一个提示：你可能想借助mlr3::benchmark()和mlr3::benchmark_grid()（更多信息，请参考https://mlr3book.mlr-org.com/performance.html#benchmarking）一次性解决练习4至6。 在做此操作时，请记住计算可能需要很长时间，可能需要几天。 这当然取决于你的系统。 你拥有的RAM和核心越多，计算时间就会越短。

E5. 使用二次判别分析（QDA）对滑坡敏感性进行建模。 评估QDA的预测性能。 QDA和GLM的空间交叉验证平均AUROC值之间有什么区别？

E6. 不调整超参数运行SVM。 使用\(\sigma\) = 1 和 C = 1的rbfdot核。 否则，在kernlab的ksvm()中不指定超参数将初始化自动的非空间超参数调整。