Слушая землетрясения: трассировка лучей сквозь планету

S-волна, дошедшая до границы ядро–мантия, закончилась. Внешнее ядро жидкое, у жидкости нет сопротивления сдвигу, а поперечной волне нужно что-то сдвигать — так что на глубине 2891,5 км фаза заканчивается. По ту сторону крупного землетрясения это оставляет сейсмическую теневую зону размером с континент, и это одно из наблюдений, которые рассказали нам, что середина планеты расплавлена, — за десятилетия до того, как её сумели смоделировать.

Эта тень появляется на экране Tremor всякий раз, когда вы выбираете крупное событие, и то, как она туда попадает, — весь замысел в миниатюре. Никакой код тень не рисует. Никакого порога по угловому расстоянию в рендерере нет нигде. Мы трассируем лучи через скоростную модель, спрашиваем каждый, где он находится в текущий момент воспроизведения, и записываем ответы в заранее выделенный буфер.

Лучи, умершие на ядре, пишут NaN. Линейный рендерер обрывает на NaN серию и начинает новую. Пропуски и есть теневые зоны.

Модель времён пробега ak135 в 22 узлах

Всё, что дальше по потоку, стоит на одной таблице. ak135 — стандартная одномерная референсная модель Земли: для каждой глубины продольная и поперечная скорости, усреднённые по всей планете; достаточно приблизительная, чтобы описать её одним абзацем, и достаточно хорошая, чтобы сейсмологи до сих пор определяли по ней координаты настоящих землетрясений. Наша версия — 22 строки [depthKm, vp, vs], интерполируемые кусочно-линейно по радиусу, а не по глубине, и отходящие от опубликованной модели примерно на 0,1 км/с. Исходник называет это визуальным классом и не лукавит: это инструмент рисования, а не таблица времён пробега ak135, по которой вы стали бы выверять момент вступления фазы.

Хитрость таблицы в том, как в ней закодированы разрывы: глубина, встречающаяся дважды, и есть разрыв. Построитель проходит по списку узлов, режет его на непрерывные сегменты везде, где глубина повторяется, и записывает скачок. Больше ничто в движке не знает границы по именам: Мохо, 410, 660, граница ядро–мантия и граница внутреннего ядра — это одна структура данных, и добавить ещё одну значит добавить две строки.

Жидкое внешнее ядро — не особый случай в движке. Это ноль в таблице.

Этот ноль делает колоссальную работу. Когда S-луч доходит до границы ядро–мантия, шаг преломления смотрит поперечную скорость по ту сторону, не находит, во что преломляться, и заканчивает луч — трассировщик так и не узнаёт, что упёрся в ядро. Ниже границы внутреннего ядра поперечная скорость возвращается, потому что внутреннее ядро твёрдое; трассировщик с удовольствием повёл бы луч и там, если бы поперечная волна вообще могла туда прийти. Ни одна не приходит, и ассерт проверяет, что ни одна и не пришла.

Как изогнуть 360 лучей сквозь планету

Сферически симметричная Земля дарит вам целое измерение: луч, вышедший из очага, никогда не покидает плоскость, в которой лежат очаг и центр планеты. Поэтому трассировщик двумерный — он ведёт каждый луч по этой большекруговой плоскости и хранит его как радиус, угловое расстояние и накопленное время: три float на сохранённую выборку, в плоском Float32Array, по одной выборке на каждые шесть шагов.

Интегрирование — классическая лучевая теория: при показателе преломления n = 1/v уравнение луча d/ds(n·d̂) = ∇n говорит, что луч загибается в сторону большего показателя, то есть к более медленной среде. В радиально симметричной среде градиент направлен вдоль радиуса, поэтому весь изгиб сворачивается в дюжину строк скалярной арифметики, выполняемых раз в восемь километров шага.

if (r > 0.5) {
  const dndr = -velocityGradientAt(r, phase) / (v * v);
  const k = (dndr * DS_KM) / r;
  const n = 1 / v;
  const wx = n * dx + k * x;
  const wz = n * dz + k * z;
  const wl = Math.hypot(wx, wz);
  if (wl > 0) {
    dx = wx / wl;
    dz = wz / wl;
  }
}

Прямые шаги хороши ровно до встречи с разрывом, где приближение сказывается немедленно: пересеките границу ядро–мантия на пару километров позже — и угол преломления неверен, а луч уходит туда, где ему делать нечего. Поэтому пересечения решаются точно. stepToRadius() находит корень квадратного уравнения, ставящего луч ровно на границу, шаг обрезается там же, и закон Снеллиуса применяется в локальном базисе «касательная—нормаль»: sin θ₂ = sin θ₁ · v₂/v₁. Если результат больше единицы, прошедшего луча нет и направление зеркалится — полное внутреннее отражение, и именно так появляются отражённые от ядра фазы, о которых никто не просил. Иначе преломлённое направление пересобирается и сдвигается на 0,02 км в сторону, чтобы следующий шаг не обнаружил то же пересечение заново.

Триста шестьдесят лучей веером расходятся из гипоцентра с углами выхода от 0,25° до 179,75°, глубина очага зажата в 1–700 км. Полный веер стоит 10–20 мс на фазу — в главном потоке, без Web Worker, и оправдать это можно лишь тем, что он считается раз на выбор события, а не раз на кадр.

Фазы, которых мы не писали

У сейсмологии богатый алфавит фаз — P, S, PKP, PKIKP, PcP, ScS — и соблазнительный способ его нарисовать это прописать каждую фазу отдельным случаем: вот эту дугу для ядерной фазы, прямую P за таким-то углом спрятать. В Tremor нет ничего подобного. Трассировщик знает скорость, градиенты, границы и закон Снеллиуса; каждая фаза на картинке следует из этих четырёх вещей.

Веер S останавливается у ядра, потому что поперечная скорость там нулевая. Веер P делает кое-что поинтереснее: на входе во внешнее ядро продольная скорость падает с 13,66 до 8,0 км/с, поэтому лучи резко преломляются вниз, прямые приходы на поверхность заканчиваются после 104°, а прошедшая через ядро энергия выходит наружу гораздо дальше по кругу как PKP. Обе кромки теневой зоны P получаются из преломления, а не из if — и ровно поэтому их надо проверять.

Поэтому движок несёт исполняемую самопроверку, работающую в разработке, в том же духе, что и ассерт про МКС в Orrery: физическому движку, который никто не оценит на глаз, нужен тест, который может.

• Время пробега вдоль каждого луча растёт монотонно. (Луч, идущий назад во времени, означает сломанный интегратор.)
• Прямая P на 60° приходит около 600 с, с допуском ±60 с; на 90° — около 780 с с допуском ±80 с.
• Прямая S на 60° приходит около 1090 с, с допуском ±110 с.
• Последний прямой приход P на поверхность попадает в диапазон 90°–118° — достаточно широкий, чтобы стерпеть модель визуального класса, и достаточно узкий, чтобы поймать сломанное преломление.
• Приходы PKP существуют дальше 130°.
• Ни один S-луч никогда не обнаруживается глубже чем на 5 км ниже границы ядро–мантия.

Одна изохрона на кадр, ноль аллокаций

Прослеженный луч — это путь в пространстве и во времени; фронт волны — срез через него. sampleIsochron() двоичным поиском находит в выборках каждого луча интервал, содержащий текущее время воспроизведения, и интерполирует из него одну позицию — по паре «радиус — расстояние» на луч, в заранее выделенный массив. Лучи, не дожившие до этого времени или умершие на ядре, пишут (NaN, NaN).

Фронты рисуются как LineSegments, а не как полилиния, и этот выбор несущий. Полилиния соединила бы последнюю живую выборку перед разрывом с первой после него, приварив хорду прямо через тень. Сегменты так не делают: NaN заканчивает серию, следующая корректная пара начинает новую, и тень вырезается из самой геометрии. Светящаяся лента собрана из независимых квадов по той же причине — общая полоса приварила бы разрыв, который линии только что оставили.

Ничто здесь не пересобирается на каждом кадре: 720 слотов выборок, 2880 вершин линий и 1440 квадов ленты выделяются один раз, а работу делает setDrawRange. Каждый фронт рисуется дважды, зеркально на ±дельту, потому что фронт волны в сферически симметричной Земле — это поверхность вращения вокруг оси очага, а разрез открывает две её полуплоскости. Фронты держатся 120 с модельного времени и гаснут с возрастом — P от почти белого через угольно-красный, S от бледно-мятного через бирюзу, — так что насыщенное воспроизведение читается как последовательность, а не как смазанное пятно.

Разрезанная Земля из геометрии, а не из отсечения

Чтобы смотреть, как фронт идёт сквозь недра, планету надо вскрыть, и очевидный способ — плоскости отсечения — для оболочек не годится: фрагментная работа на каждом пикселе каждого слоя, а в итоге ещё и полые грани. Вместо этого оболочки строятся неполными. SphereGeometry с phiStart = π и phiLength = 1.5π — это оболочка в 270° с недостающим квадрантом, повёрнутая к камере по умолчанию. Дыра сидит в самом меше, поэтому не стоит ничего.

Один слой этим трюком воспользоваться не может. Береговые линии и градусная сетка обязаны вращаться — в этом их смысл, — поэтому они единственное в сцене, что честно отсекается на GPU, двумя плоскостями с включённым clipIntersection. Данные береговых линий декодируются из модуля, упакованного нами для Orrery, и заново разбиваются слерпом по большим кругам примерно каждые 2°, чтобы длинные хорды не провисали под поверхность и не пропадали внутри сферы.

Дальше сцена живёт в двух системах отсчёта, и это та часть, которой мы тихо гордимся больше всего. Система разреза не двигается никогда: клин, слои, волосяные линии границ и полярная ось статичны, поэтому мировые плоскости отсечения — константы, а графика граней разреза запекается в один индексированный меш и один вызов отрисовки. Географическая система — береговые линии, сетка, весь каталог — доворачивается слерпом так, чтобы выбранный гипоцентр ехал по полярной оси +Y, ребру, общему для обеих открытых граней. Любая плоскость через ось очага — допустимое большекруговое сечение сферически симметричной Земли, а здесь две из них по совпадению ещё и стены, которые видно. Разрез честен, какое бы событие вы ни выбрали, — и мгновенен, вовсе без анимации, для тех, кто попросил уменьшить движение.

Четыре тысячи землетрясений в одном меше

Каталог приходит из двух сводных фидов USGS — прошедшие сутки любой магнитуды плюс M4.5 и выше за прошедший месяц — через прокси, который ужимает каждое событие до позиционного кортежа и сам владеет своим кешированием. Ревалидация на уровне маршрута закешировала бы 503 с тем же удовольствием, что и данные, и на Orrery это однажды заморозило офлайн-состояние на шесть часов. Здесь TTL в 90 секунд привязан к более старой из двух отметок времени, так что упавший фид продолжают перезапрашивать, пока живой продолжает отдавать. Устаревшее лучше пустого.

На экране весь каталог — до 4000 событий — это один InstancedMesh из октаэдров. Каждый инстанс сидит на радиусе (6371 − depth) / 1000 единиц сцены: на своей настоящей гипоцентральной глубине, а не приколот к поверхности, поэтому зона субдукции читается как плита, наклонно уходящая в мантию. Масштаб логарифмичен по магнитуде, 0.012 + 0.02 · 1.9^(mag − 4), зажат с обоих концов: сверху — чтобы M8 осталась бусиной, а не луной, снизу — чтобы микроземлетрясения всё-таки существовали. Цвет — это возраст: от горячего оранжевого к угольному за первые 24 часа, дальше от угольного к сланцевому на тридцатидневном хвосте. История прошедшей недели видна ещё до того, как вы что-нибудь нажмёте.

Озвучка землетрясений из самописного декодера miniSEED

Сейсмограммы открыты. Сырое движение грунта может забрать кто угодно из сервисов FDSN dataselect, а станции IU сети Global Seismographic Network пишут непрерывно десятилетиями. Поэтому последнее, что Tremor делает с событием, — предлагает вам его проиграть: настоящая озвучка землетрясения, а не синтезатор, который её изображает.

Сам браузер этого не может по двум независимым причинам. Данные приходят как сжатый по Стейму miniSEED, разностно кодированный двоичный формат, о котором ни один браузер не слышал, а сервис не отдаёт CORS-заголовков, так что страница не смогла бы даже получить эти байты, умей она их читать. Обе половины живут на нашем сервере: прокси забирает, самописный декодер распаковывает, и до клиента доезжает нормализованный буфер Int16.

Выбор станции — маленькая задача ранжирования с физическим ограничением внутри. Tremor знает 31 широкополосную станцию IU и хочет сорок минут записи вертикальной компоненты, но ближайшая не значит лучшая: ближе 2° прибор зашкаливает на крупном событии, дальше 90° ядро начинает затенять P. Поэтому в ранжировании ведёт телесейсмический диапазон, расстояние решает ничьи, три кандидата пробуются по очереди, а запись моложе 45 минут отвергается сразу — архивам нужно время, и клиент предупреждает этот отказ, вместо того чтобы дать вам его обнаружить.

Декодер miniSEED — примерно настолько низкий уровень, насколько веб-код вообще опускается: DataView, чтение в big-endian и заголовок blockette-1000, сообщающий, с какой кодировкой мы имеем дело. Стейм упаковывает отсчёты в 64-байтовые кадры из шестнадцати 32-битных слов, первое из которых — карта полубайтов: шестнадцать двухбитных кодов, описывающих упаковку остальных пятнадцати. Steim-1 предлагает одну, две или четыре разности на слово; Steim-2 добавляет подкоды, укладывающие в одно слово пять 6-битных, шесть 5-битных или семь 4-битных разностей, плюс 30-битный слот для выбросов. Первое и второе слова нулевого кадра вовсе не разности, а константы интегрирования X0 и Xn, что делает поток самопроверяющимся. Реализованы обе кодировки; там, где тайминг записей расходится больше чем на 1,5 отсчёта, мы считаем пропуск и всё равно склеиваем, потому что для озвучки непрерывность важнее точности.

И вот мы дошли до части, в которой нет кода. Сейсмическая энергия живёт около 0,05–10 Гц; слух начинается примерно с 20 Гц. Учебниковое решение — питч-шифтер, и мы не написали его ни одной строчки, потому что Web Audio API позволяет объявить частоту дискретизации буфера. Сорок отсчётов в секунду, объявленные как 40 × 150 = 6000 Гц, играют в 150 раз быстрее и в 150 раз выше — внутри допустимого окна частот буфера от 3000 до 768 000 Гц — и сажают сейсмический диапазон на 7,5–1500 Гц. Сорок минут превращаются в шестнадцать секунд. Это не преобразование данных; это сами данные, прочитанные быстро.

Динамический диапазон — задача потруднее. Телесейсмическое вступление P может быть в сотню раз больше идущей за ним коды, и при линейном проигрывании посетитель слышит один щелчок и тишину. Поэтому клиент компандирует с показателем около 0,6 — в браузере, а не на сервере, потому что дело сервера честные отсчёты, а формовка должна жить там, где посетитель слышит, как её применяют. Перед передачей у записи снимается среднее (та же SNZO стоит примерно на 1500 отсчётов от нуля) и делается нормировка по пику. Спарклайн рисует пик на корзину, а не RMS, по той же причине: приход — это всплеск, а усреднение ровно та операция, которая его спрячет.

Сам плеер намеренно маленький: один GainNode, рампы по 30 мс, запланированные на шаг вперёд, аудиоконтекст, создаваемый прямо в обработчике клика, чтобы политике автовоспроизведения не к чему было придраться, и остановка проигрывания, когда вкладка уходит. Мьют у Tremor не свой — он читает общесайтовый мьют синтезатора, разделяемый со всеми играми Лаборатории. Заглушите игру, и планета тоже замолчит.

Чип в углу экрана говорит времена пробега визуального класса, по той же причине, по которой в Orrery он говорит распространение визуального класса. Одномерная модель без станционных поправок и без трёхмерной структуры ставит фронт туда, где его ждёт глаз, а не туда, где его записал бы прибор. Честно для глаза, не для прибора; никогда не используйте это ни для чего, кроме изумления. Сказать это прямо не стоит ничего, а всей остальной картинке покупает доверие.

Честные пробелы: трассировка лучей идёт в главном потоке, а у страницы пока нет автоматических сквозных тестов — только ассерты времени разработки, которые проверяют физику, а не кнопки. И то, и другое в списке.