刚好我之前处理过类似规模的空间数据集，给你梳理一套实操性强的解决方案，完美适配你18万多边形的场景：

第一步：定位孤立多边形 #

首先得从已生成的皇后邻接列表里揪出那些没任何邻接对象的单元：

# 加载必备包
library(spdep)
library(sp)

# 找出邻接列表长度为0的多边形索引
isolated_idx <- which(sapply(sample_queen_nb, length) == 0)
# 拿到对应的多边形ID
isolated_ids <- row.names(sp.sample)[isolated_idx]

第二步：为孤立多边形计算最近邻 #

划重点！18万多边形的规模，绝对别用两两距离计算，会慢到怀疑人生。推荐用spdep里的k近邻搜索，专门针对大数据集优化过：

# 获取所有多边形的质心坐标（用质心算距离足够精准，还快）
coords <- coordinates(sp.sample)

# 为所有多边形计算1近邻的邻接列表（后续只取孤立单元的结果）
knn_nb <- knn2nb(knearneigh(coords, k=1), row.names = row.names(sp.sample))

# 提取每个孤立多边形的最近邻ID
isolated_neighbors <- sapply(isolated_idx, function(i) knn_nb[[i]])

如果你的数据是地理坐标系（经纬度），记得给knearneigh加个参数longlat=TRUE，用球面距离计算更准确。

第三步：补全皇后邻接列表 #

把找到的最近邻塞进原来的皇后邻接列表里就行：

# 遍历每个孤立多边形，把最近邻加入它的邻接列表
for (i in isolated_idx) {
  sample_queen_nb[[i]] <- isolated_neighbors[which(isolated_idx == i)]
}

第四步：重新生成权重矩阵并验证 #

现在可以重新生成权重矩阵，再绘图确认连通性：

# 重新创建权重列表
sample_W.list.queen <- nb2listw(sample_queen_nb, style="W", zero.policy = TRUE)
# 转换为矩阵格式
sample_W.queen <- listw2mat(sample_W.list.queen)

# 再次可视化检查，孤立多边形现在应该有绿色连线了
plot(sp.sample)
plot(sample_queen_nb, coords, add=TRUE, col="green", cex=0.5)

如果你想用sf包的替代方案 #

要是你已经转用sf对象，流程差不多，只是坐标提取和邻接计算的方式稍有不同：

library(sf)
library(spdep)

# 假设你的sf对象是sf.sample
# 提取多边形质心坐标
coords_sf <- st_coordinates(st_centroid(sf.sample))

# 生成皇后邻接列表（先转成Spatial对象给poly2nb用）
sf_queen_nb <- poly2nb(as(sf.sample, "Spatial"), queen=TRUE)

# 定位孤立多边形
isolated_idx_sf <- which(sapply(sf_queen_nb, length) == 0)

# 计算1近邻
knn_nb_sf <- knn2nb(knearneigh(coords_sf, k=1), row.names = sf.sample$ID)

# 补全邻接列表
for (i in isolated_idx_sf) {
  sf_queen_nb[[i]] <- knn_nb_sf[[i]]
}

# 生成权重矩阵
sf_W.list.queen <- nb2listw(sf_queen_nb, style="W", zero.policy = TRUE)
sf_W.queen <- listw2mat(sf_W.list.queen)

小提醒 #

1. 补全后可以检查权重矩阵的行和，因为用了style="W"，每行的和应该都是1，孤立单元补全后也会符合这个规则。
2. 如果实在担心质心距离不够精准，也可以用多边形边界距离，但计算量会暴增，18万数据真心不推荐。

内容的提问来源于stack exchange，提问作者Tom