如何在实际地图上构建灾难模拟场景：利用 Mapbox GL JS 和 Three.js 技术实现动画效果的画布层以及可坍塌的 3D 建筑模型

步骤8：将建筑物模型从地图自带的矢量瓦片中分离出来

你其实并不需要专门的建筑物数据集——因为地图本身已经在下载这样的数据了。Mapbox Streets的矢量瓦片包含了building这个源层，其中记录了建筑物的多边形轮廓、高度、最低高度以及是否需要进行三维建模等信息。querySourceFeatures函数可以读取当前已加载的所有这些数据。

function buildingSource(map) {
  const sources = map.style().sources;
  if (sources.composite) return "composite";         // Mapbox样式
  if (sources.openmaptiles) return "openmaptiles";   // MapLibre + OpenMapTiles
  for (const id in sources) if (sources[id].type === "vector") return id;
  return null;
}

function collectBuildings(map, center, kx, ky) {
  const src = buildingSource(map);
  if (!src) return [];
  const feats = map.querySourceFeatures(src, { sourceLayer: "building" });
  const out = [], seen = new Set();

  for (const f of feats) {
    const p = f.properties || {}, g = f.geometry;
    if (!g || p.extrude === "false" || p.underground === "true") continue;

    const polys = g.type === "Polygon" ? [g_coordinates]
                : g.type === "MultiPolygon" ? gcoordinates : null;
    if (!polys || !polys[0] || polys[0][0].length < 4) continue;

    // 外环的质心和边界坐标
    const ring = polys[0][0], n = ring.length - 1;
    let sx = 0, sy = 0, x0 = Infinity, x1 = -Infinity, y0 = Infinity, y1 = -Infinity;
    for (let j = 0; j < n; j++) {
      const [x, y] = ring[j];
      sx += x; sy += y;
      if (x < x0) x0 = x; if (x > x1) x1 = x;
      if (y < y0) y0 = y; if (y > y1) y1 = y;
    }
    sx /= n; sy /= n;

    // 同一建筑物在它所覆盖的每个图块中只会出现一次
    const key = sx.toFixed(6) + "," + sy.toFixed(6) + "," + n;
    if (seen.has(key)) continue;
    seen.add(key);

    let h = +(p.height ?? p.render_height);
    if (!(h > 0)) h = 8;                              // 没有数据时，假设建筑物有两层
    const mh = Math.min(h, +(p.min_height ?? p.render_min_height) || 0);

    out.push({
      geometry: g, polys, lng: sx, lat: sy, h, mh,
      d: Math.hypot((sx - centerlng) * kx, (sy - center.lat) * ky),   // 以中心点为基准的距离，单位：米
      rf: 0.5 * Math.sqrt(Math.max(4, (x1 - x0) * kx * (y1 - y0) * ky)), // 近似半径，单位：米
    });
  }
  out.sort((a, b) => a.d - b.d);
  return out.slice(0, 9000);
}

kx和ky分别表示该城市所在纬度处，经度和纬度每变化一度时高度的变化量，其数值分别为111320 * cos(lat)和110574。后续的所有计算都是以震中为基准，用米作为单位进行的。

关于这些数据有三点需要注意：首先，querySourceFeatures会为每个与其相交的图块返回一个建筑对象，因此需要对这些对象进行去重处理；使用四舍五入后的质心坐标再加上建筑的数量作为键值对，是一个不错的选择。

其次，在Mapbox Streets中，这些建筑对象没有固定的ID，这一点在步骤11中非常重要。

第三，你只能获取那些已经被加载到系统中的数据，这意味着在收集建筑信息之前，你需要先放大地图到显示建筑的位置，然后等待相应的图块数据被下载完成。

在HazardSim系统中，程序会分两步进行操作：首先会等待直到map.areTilesLoaded()的值为true，且建筑数量在连续检测中不再发生变化，之后才会执行collectBuildings函数。如果过早开始收集数据，你得到的可能只是位于第一个图块中的那十栋建筑。

步骤9：使拉伸效果看起来像建筑物

默认情况下，拉伸效果表现为一个单色的平面立方体形状，看起来就像一张简单的图表。其实有一个很简单的方法可以改善这个效果：fill-extrusion-pattern函数可以从样式文件中获取图像，并将其覆盖在建筑物的墙壁上；你还可以在程序运行时使用map.addImage函数为这些图像添加新的内容。

可以在画布上绘制一些立面纹理效果（例如，墙上排列成行的窗户，颜色稍作变化并加入一些随机噪声，这样就不会显得像是计算机生成的图像），将这些纹理注册到系统中，然后根据每栋建筑的高度为其选择合适的纹理：

function makeFacade(map, id, wall, glass) {
  const cv = document.createElement("canvas");
  cv.width = 128; cv.height = 64;
  const c = cv.getContext("2d");
  c.fillStyle = wall; c.fillRect(0, 0, 128, 64);

  // 绘制四层窗户的图案，每层有16个窗口
  for (let f = 0; f < 4; f++) {
    for (let x = 3; x + 6 <= 128; x += 12) {
      const lit = Math.random() < 0.07;
      c.fillStyle = lit ? "#c6b68c" : glass;
      c.fillRect(x, f * 16 + 3, 6, 9);
    }
  }
  // 添加颗粒质感
  const im = c.getImageData(0, 0, 128, 64), d = im.data;
  for (let i = 0; i < d.length; i += 4) {
    const n = (Math.random() - 0.5) * 12;
    d[i] += n; d[i + 1] += n; d[i + 2] += n;
  }
  map.addImage(id, im);
}

makeFacade(map, "fac-glass", "#566a78", "#3c4e5c");
makeFacade(map, "fac-light", "#cdc8be", "#364048");
makeFacade(map, "fac-brick", "#8a5a44", "#2e2a26");

接下来需要设置图层的属性。建筑的高度和基底宽度都是通过表达式计算得出的，这两个数值都会乘以一个名为k的字段状态值，其默认值为1。这个乘数实际上决定了最终的渲染效果，因此现在需要将其加入代码中：

const H  = ["coalesce", ["get", "height"], ["get", "render_height"], 8];
const MH = ["coalesce", ["get", "min_height"], ["get", "render_min_height"], 0];
const K  = ["coalesce", ["feature-state", "k"], 1];

function addBuildingLayers(map, source, sourceLayer, before) {
  const common = { type: "fill-extrusion", source, minzoom: 13 };
  if (sourceLayer) {
    common["source-layer"] = sourceLayer;
    common.filter = ["all", ["!=", ["get", "extrude"], "false"], ["!=", ["get", "underground"], "true"];
  }
  map.addLayer({
    ...common, id: "bld-walls",
    paint: {
      "fill-extrusion-pattern": ["case",
        [">=", H, 45], "fac-glass",
        [">>", H, 14], "fac-light",
        "fac-brick"],
      "fill-extrusion-height": ["*", H, K],
      "fill-extrusion-base":   ["*", MH, K],
      "fill-extrusion-opacity": 1,
    },
  }, before);
  // 在建筑物顶部添加一层较薄的拉伸效果，以形成屋顶的纹理
  map.addLayer({
    ...common, id: "bld-roof",
    paint: {
      "fill-extrusion-pattern": "fac-roof",
      "fill-extrusion-height": ["+", ["*", H, K], 0.05],
      "fill-extrusion-base":   ["max", ["*", MH, K], ["-", ["*", H, K], 0.45]],
      "fill-extrusion-opacity": 1,
    },
  }, before);
}

由于模式无法针对每个特征进行旋转或缩放，因此窗户的位置与楼层的对齐情况并不完全一致。但在当前设置的摄像机高度下，这种差异根本看不出来。

步骤10：确定哪些建筑物会倒塌

现在来看模拟部分，这部分内容相对简短。对于每一栋建筑物来说，都需要知道该位置处的震动强度，以及在这种震动强度下该类型建筑物发生倒塌的概率。

const DEPTH_KM = 10;

function mmiAt(M, rKm) {
  const R = Math.sqrt(rKm * rKm + DEPTH_KM * DEPTH_KM);
  const v = 1.5 * M - 3.0 * Math.log10(R) - 0.0018 * R + 1.6;
  return Math.max(0, Math.min(12, v));
}

const FRAG = [[6.5, 0], [7, 0.004], [8, 0.025], [9, 0.08], [10, 0.2], [11, 0.38], [12, 0.55]];
const QUALITY = { strict: 0.2, avg: 1, poor: 3 };

function collapseProbability(mmi) {
  if (mmi <= FRAG[0][0]) return 0;
  for (let i = 1; i < FRAG.length; i++) {
    if (mmi <= FRAG[i][0]) {
      const [m0, p0] = FRAG[i - 1], [m1, p1] = FRAG[i];
      return p0 + (p1 - p0) * (mmi - m0) / (m1 - m0);
    }
  }
  return FRAG[FRAG.length - 1][1];
}

function assignDamage(buildings, M, epicentre, quality) {
  const kx = 111320 * Math.cos(epicentre.lat * Math.PI / 180), ky = 110574;
  const q = QUALITY[quality] || 1;

  for (const b of buildings) {
    const rKm = Math.hypot((b.lng - epicentrelng) * kx, (b.lat - epicentre.lat) * ky) / 1000;
    const mmi = mmiAt(M, rKm);
    const hf = b.h < 5 ? 0.7 : b.h < 35 ? 1.15 : b.h < 70 ? 0.8 : 0.45;
    const pCollapse = Math.min(0.9, collapseProbability(mmi) * q * hf);
    const pDamage = mmi >= 6 ? Math.min(0.6, pCollapse * 2.5 + 0.004) : 0;

    const r = Math.random();
    b.cat = r < pCollapse ? 2 : r < pCollapse + pDamage ? 1 : 0;   // 2表示倒塌，1表示受损，0表示完好无损
    b.start = 0.8 + Math.random() * 2.8;                            // 地震发生后的秒数
  }
}

步骤11：实现层的折叠动画效果而不重新构建这些层

这就是步骤9中采用的机制所带来的实际效果。在每一帧中，对于每一个正在倒塌的建筑物，都会将其特征状态k设置为一个介于1与其最终残骸高度之间的数值。这样，建筑物的形状就会逐渐缩小，但不会重新构建任何层，也不会重新生成GeoJSON数据，更不会重新解析地图瓦片。

map.setFeatureState({ source: "bld", id: b.id }, { k: 0.62 });

不过这里有一个问题，也正是这个问题让我在这个项目中花费了最多的时间：setFeatureState方法需要一个特征ID，而Mapbox Streets提供的建筑物数据并没有这样的ID。如果直接使用这个方法来处理这些数据，那么根本不会产生任何效果。

解决办法是将你收集到的这些建筑物数据复制到自己创建的GeoJSON数据源中，并为这些建筑物指定一个由你自己控制的id，然后让系统将这个id设置为特征ID：

const feats = buildings.map((b, i) => {
  b.id = i;
  return { type: "Feature", properties: { i, height: b.h, min_height: b.mh }, geometry: b.geometry };
});

map.addSource("bld", {
  type: "geojson",
  data: { type: "FeatureCollection", features: feats },
  promoteId: "i",
  maxzoom: 16,
});
addBuildingLayers(map, "bld", null, firstSymbolLayer);

现在，setFeatureState({ source: "bld", id: 17 }, ...)这种方法就可以用来处理第17栋建筑物了。在GeoJSON数据源报告isSourceLoaded之前，请保留原有的地图瓦片层；之后再删除这些层，这样就不会有任何帧出现城市画面闪烁的情况了。

k这个数值随时间变化的变化方式，决定了建筑物的倒塌效果是真正的倒塌，还是只是看起来像建筑物滑入了地面。最终我设计了五种不同的变化模式，每种模式都是针对特定的建筑物随机选择的，并且其变化幅度会受到建筑物高度的限制：

// k0表示起始高度倍数（默认为1），kf表示最终残骸高度，u表示变化进度（0到1之间）
function collapseK(b, u) {
  const k0 = b.k0, sp = k0 - b.kf;
  switch (b.mode) {
    case 1: {   // 逐层倒塌：先快速下降一段时间，然后暂停
      const s = u * b.ns, i = Math.min(b.ns - 1, Math.floor(s));
      let f = s - i;
      f = f < 0.45 ? f * f / 0.45 * 1.1 : 0.55 + (f - 0.45) * 0.45 / 0.55;
      return k0 - sp * Math.min(1, (i + Math.min(1, f)) / b.ns);
    }
    case 2: {   // 缓慢倒塌：先从底层开始倒塌，下降一层后暂停，然后继续倒塌
      const sd = Math.min(sp * 0.5, 3.2 / b.h);
      if (u < 0.12) return k0 - sd * (u / 0.12) ** 2;
      if (u < 0.5) return k0 - sd;
      const v = (u - 0.5) / 0.5;
      return k0 - sd - (sp - sd) * v * v;
    }
    case 3: {   // 上部结构先倒塌
      const pt = sp * 0.45;
      if (u < 0.3) return k0 - pt * (u / 0.3) ** 2;
      if (u < 0.55) return k0 - pt;
      const v = (u - 0.55) / 0.45;
      return k0 - pt - (sp - pt) * v * v;
    }
    case 4:     // 横向移动：速度保持稳定，但逐渐加快
      return k0 - sp * Math.pow(u, 1.6);
    default:    // 像煎饼一样倒塌：速度不断加快
      return k0 - sp * u * u;
  }
}

对于那些虽然受损但仍然矗立在原处的建筑物（cat === 1），kf的值接近1，因此这些建筑物的倒塌过程会表现为短暂的下降：它们会明显下降半米然后停止移动。这个细节虽然很小，但却能让整条街道看起来都受到了影响，而不仅仅是一些建筑物发生了倒塌。

请将dt的值限制在40毫秒以内。如果动画播放过程中某个帧的渲染时间超过了25毫秒，那么正在下降中的建筑物就会突然停止移动，这种变化是肉眼可以立刻察觉到的。

步骤12：地图中的Three.js场景

对于远处的建筑物，可以使用缩放技术来处理它们的显示效果；而对于靠近相机的建筑物，则需要让观众能够看到它们结构分解的细节——比如塔楼倾斜后断裂成碎片、碎屑从建筑物底部飞溅出来，或者废墟散落在街道上。为了实现这些效果，就需要使用three.js这个库。

自定义图层会接收两个回调函数：onAdd在图层被添加到地图中时被调用，此时会传入地图对象及其GL渲染上下文；render则会在每一帧绘制时被调用，此时会传入相应的渲染上下文以及地图的投影矩阵。需要注意的是坐标系统的转换：我们应该使用以米为单位的坐标系，而不是墨卡托坐标系，因此需要构建一个转换矩阵，并在每一帧绘制时将其应用到地图的投影计算中：

function makeScene(map, origin, before) {
  const mc = mapboxgl.MercatorCoordinate.fromLngLat([originlng, origin.lat], 0);
  const s = mc.meterInMercatorCoordinateUnits();
  // 原点处向东/向北/向上的距离（以米为单位）→ 转换为墨卡托坐标系的距离。在墨卡托坐标系中，y轴的方向是向南的，因此需要乘以-s。
  const toMercator = new THREE.Matrix4()
    .makeTranslation(mc.x, mc.y, mc.z)
    .scale(new THREE.Vector3(s, -s, s));

  const scene = new THREE.Scene();
  const camera = new THREE.Camera();
  scene.add(new THREE.HemisphereLight(0xfff4e4, 0x3c342d, 0.55));
  scene.add(new THREE.AmbientLight(0xffffff, 0.1));
  const sun = new THREE.DirectionalLight(0xfff8ee, 0.8);
  sun.position.set(-0.45, -0.6, 1);
  scene.add(sun);

  let renderer = null;
  const layer = {
    id: "debris3d",
    type: "custom",
    renderingMode: "3d",
    onAdd(m, gl) {
      renderer = new THREE.WebGLRenderer({ canvas: m.getCanvas(), context: gl, antialias: true });
      renderer.autoClear = false;             // 地图已经绘制了这一帧，所以不需要重新清除屏幕
      renderer.localClippingEnabled = true;   // 这个设置用于让建筑物看起来像埋在地下一样
    },
    render(gl, matrix) {
      camera.projectionMatrix.fromArray(matrix).multiply(toMercator);
      renderer.resetState();                  // 由于地图的更新，GL渲染状态可能发生了变化
      renderer.render(scene, camera);
    },
  };
  map.addLayer(layer, before);
  return { scene, camera, layer };
}

在MapLibre中，render函数的第二个参数是一个对象，而矩阵数据存储在args.defaultProjectionData.mainMatrix中。其他方面都是一样的。

renderingMode: "3d"这个设置使得地图能够与你共享其深度缓冲区。例如，当一块混凝土结构位于某个凸出物之后时，这块混凝土就会被该凸出物遮挡而无法显示出来——这就是为什么需要将渲染结果直接放入地图的上下文中，而不是在另一个独立的画布上显示的原因。

当有成千上万的碎片需要被模拟时，就可以使用实例化技术。 HazardSim会维护一个InstancedMesh对象池，每种类型的碎片都会对应一个这样的对象（比如板块、墙块碎片、钢筋、岩石或屋顶结构），每个对象的容量是固定的。每当新的碎片被创建时，系统就会为这些对象分配相应的索引。

function makePool(scene, geometry, material, capacity) {
  const mesh = new THREE.InstancedMesh(geometry, material, capacity);
  // 提前为所有对象分配颜色缓冲区（具体实现可参考相关文档）
  mesh.instanceColor = new THREEInstancedBufferAttribute(new Float32Array(capacity * 3).fill(1), 3);
  mesh.count = 0;
  mesh.frustumCulled = false;                 // 地图使用的相机并非three.js默认的相机类型
  mesh.instanceMatrix.setUsage(THREE.DynamicDrawUsage);
  scene.addmesh);

  const dummy = new THREE.Object3D(), color = new THREE.Color();
  dummy.rotation.order = "ZXY";
  return {
    n: 0, capacity,
    add(p) {
      if (this.n >= capacity) return -1;
      const i = this.n++;
      this.set(i, p);
      color.setRGB(p.rgb[0], p.rgb[1], p.rgb[2]);
      mesh.setColorAt(i, color);
      mesh.instanceColor.needsUpdate = true;
      mesh.count = this.n;
      return i;
    },
    set(i, p) {
      dummy.position.set(p.x, p.y, p.z);
      dummy.rotation.set(p.tx, p.ty, p.yaw);
      dummy.scale.set(p.sx, p.sy, p.sz);
      dummy.updateMatrix();
      mesh.setMatrixAt(i, dummy.matrix);
      mesh.instanceMatrix.needsUpdate = true;
    },
  };
}

const concrete = new THREE.MeshStandardMaterial({ map: concreteTexture, roughness: 1 });
const pools = {
  slab: makePool(scene, new THREE.BoxGeometry(1, 1, 1), concrete, 30000),
  wall: makePool(scene, new THREE.BoxGeometry(1, 1, 1), facadeMaterial, 1500),
  bar: makePool(scene, new THREE.BoxGeometry(1, 1, 1), new THREE.MeshStandardMaterial({ color: 0x6e4c38, metalness: 0.25 }), 20000),
  rock: makePool(scene, rockGeometry(), concrete, 9000),
};

虽然有三万块板块需要被模拟，但实际上只需要进行一次渲染操作而已。

地图还需要你做一件事：只有当某些内容发生变化时，地图才会重新绘制。从地图的角度来看，你的场景中的元素发生变化并不等同于地图本身需要更新。因此，在每一帧中，当你移动了任何对象之后，都应该调用map.triggerRepaint()方法，这样才能确保你的网格模型能够在屏幕上正确显示出来。

步骤13：从建筑物轮廓到会倾倒的盒子

矢量贴图可以生成建筑物的轮廓多边形。但要让建筑物倾倒，就需要一个“盒子”模型——这个盒子需要有一个中心点、一个方向以及两个边长。这种有方向的边界盒可以用以下代码表示：

function orientedBox(ring, kx, ky, origin) {
  // 将经纬度坐标转换为局部米坐标
  const loc = ring.map(([lng, lat]) => [(lng - origin.lng) * kx, (lat - origin.lat) * ky]);

  // 最长的边决定了盒子的方向
  let ang = 0, best = 0;
  for (let i = 1; i < loc.length; i++) {
    const ex = loc[i][0] - loc[i - 1][0], ey = loc[i][1] - loc[i - 1][1];
    if (ex * ex + ey * ey > best) { best = ex * ex + ey * ey; ang = Math.atan2(ey, ex); }
  }

  // 将每个顶点投影到该轴及其法线上，从而得到盒子的边界范围
  const ca = Math.cos(ang), sa = Math.sin(angle);
  let u0 = Infinity, u1 = -Infinity, v0 = Infinity, v1 = -Infinity;
  for (const [x, y] of loc) {
    const u = x * ca + y * sa, v = -x * sa + y * ca;
    if (u < u0) u0 = u; if (u > u1) u1 = u;
    if (v < v0) v0 = v; if (v > v1) v1 = v;
  }
  const cu = (u0 + u1) / 2, cv = (v0 + v1) / 2;
  return {
    x: cu * ca - cv * sa, y: cu * sa + cv * ca,    // 中心点，以原点为基准的坐标
    ang,
    L1: Math.max(2, u1 - u0), L2: Math.max(2, v1 - v0),
  };
}