基于Cesium和Kriging插值的三维地形生成：原理、实现与性能优化

1次阅读

共计 3365 个字符，预计需要花费 9 分钟才能阅读完成。

传统地形生成方法（如规则网格插值）存在两个致命缺陷：

精度不足：在采样点稀疏区域会产生明显的阶梯状 artifacts
效率低下 ：全局插值计算导致 O(n³) 时间复杂度，无法处理大规模数据

Kriging 作为地质统计学经典算法，其核心优势在于：

考虑空间自相关性（通过半变异函数建模）
提供最优无偏估计（BLUE 性质）
可输出预测误差面（用于置信度评估）

算法类型	计算复杂度	适用场景	平滑效果	参数敏感性
IDW	O(n)	均匀采样	中等	低
RBF	O(n³)	不规则采样	强	高
Kriging	O(n³)	空间相关性明显数据	可调节	中等

关键结论：当处理具有明显空间连续性的地理数据时，Kriging 在精度 / 效率平衡上表现最优。

Kriging 的核心是求解以下方程组：

\begin{cases}
\sum_{j=1}^n w_j γ(h_{ij}) + μ = γ(h_{i0}) & i=1,...,n \\
\sum_{j=1}^n w_j = 1
\end{cases}

其中：
– γ(h) 为半变异函数
– w_j 是权重系数
– μ 为拉格朗日乘子

常用半变异函数模型：

// 球状模型
function spherical(h, range, sill) {
  return h < range 
    ? sill * (1.5*(h/range) - 0.5*(h/range)**3) 
    : sill;
}

DEM 数据转换：

gdal_translate -of XYZ input.tif output.xyz

归一化处理：

function normalize(points) {const xs = points.map(p => p[0]);
  const ys = points.map(p => p[1]);
  const zs = points.map(p => p[2]);

  return {
    points: points.map(p => [(p[0]-Math.min(...xs))/(Math.max(...xs)-Math.min(...xs)),
      (p[1]-Math.min(...ys))/(Math.max(...ys)-Math.min(...ys)),
      (p[2]-Math.min(...zs))/(Math.max(...zs)-Math.min(...zs))
    ]),
    scales: [xs, ys, zs].map(arr => ({min: Math.min(...arr),
      max: Math.max(...arr)
    }))
  };
}

const viewer = new Cesium.Viewer('cesiumContainer');

// 创建自定义地形 provider
class KrigingTerrainProvider extends Cesium.TerrainProvider {constructor(options) {
    super({
      // 必须配置的属性
      tilingScheme: new Cesium.GeographicTilingScheme(),
      availability: new Cesium.TerrainAvailability(),
      hasVertexNormals: true
    });

    this._krigingResult = options.krigingData;
  }

  requestTileGeometry(x, y, level) {
    return Cesium.TerrainProvider.
      createGriddedTerrainProvider({
        // 将 kriging 结果转换为网格
        grid: this._convertToGrid(x, y, level)
      });
  }
}

// 核心 kriging 计算类
class OrdinaryKriging {constructor(points) {
    this.points = points;
    this.variogram = this._computeVariogram();}

  predict(target) {const X = this._buildKrigingMatrix(target);
    const weights = numeric.solve(X.matrix, X.rhs);

    return {value: this._calculatePrediction(weights),
      variance: this._calculateVariance(weights)
    };
  }

  _computeVariogram() {
    // 实现半变异函数计算
    const lags = this._createLags(20); // 分 20 个滞后距
    return {
      model: 'spherical',
      parameters: this._fitModel(lags)
    };
  }
}

// Cesium 地形生成入口
async function generateTerrain() {
  // 1. 加载原始数据
  const demData = await loadDEM('data.xyz');

  // 2. 执行 kriging 插值
  const kriging = new OrdinaryKriging(demData);
  const grid = kriging.generateGrid(256); // 256x256 网格

  // 3. 创建地形 provider
  const provider = new KrigingTerrainProvider({krigingData: grid});

  // 4. 应用到场景
  viewer.terrainProvider = provider;
}

// 将大区域划分为 16x16 的子块
const CHUNK_SIZE = 16;

function processChunked(data) {const results = [];

  for(let x=0; x<data.width; x+=CHUNK_SIZE) {for(let y=0; y<data.height; y+=CHUNK_SIZE) {
      const chunk = data.getSubset(
        x, y, 
        Math.min(CHUNK_SIZE, data.width-x),
        Math.min(CHUNK_SIZE, data.height-y)
      );

      results.push(processChunk(chunk));
    }
  }

  return mergeResults(results);
}

// 主线程
const workerPool = new WorkerPool(4); // 4 个 worker

workerPool.dispatch({
  task: 'kriging', 
  data: chunkData
}, result => {// 处理计算结果});

// worker.js
self.onmessage = (e) => {if(e.data.task === 'kriging') {const result = performKriging(e.data.data);
    self.postMessage(result);
  }
};

使用 Float32Array 替代常规数组
及时释放中间计算矩阵
利用 Cesium 的 TerrainProvider 缓存机制

参数配置黄金法则：

半变异函数 range 参数应设为采样点平均间距的 1.5 倍
nugget 值建议设置为测量误差方差
对高程数据使用对数变换提升稳定性

数值稳定性保障：

// 添加 jitter 避免矩阵奇异
function addJitter(matrix, epsilon=1e-6) {return matrix.map((row, i) => 
    row.map((val, j) => 
      i === j ? val + epsilon : val
    )
  );
}