Владивосток, оползни и логистическая регрессия16.08.2021 03:01

27–28 августа 2019 года во Владивостоке и Приморском крае произошли массовые оползни. К счастью, обошлось без жертв. Однако, материальные потери оказались существенными: разбитые автомобили, перекрытые дороги, поврежденные здания и детские площадки. Оползни сошли в момент прохождения мощного циклона с обильными дождями. Мы робко предположили что «осадки виновны», распаковали методы классической статистики и приступили к исследованию.

TL; DR

Копии архивов исторических данных вместе с кодом статистического анализа на Python здесь. Почитать опубликованную статью в Pure and Applied Geophysics можно тут.

Последствия оползневой активности в г. Владивосток, 28–29 августа 2019 г. (a) ул. ТОбольская 11; (b) ул. Сафонова; © ул. Юмашева 20; (d) ул. Почтовая; (e, f) ул. Снеговая 64Последствия схода оползня в г. Владивосток, 28–29 августа 2019 г., на 5-й Террасной улице. Стрелки показывают направление потока.

Данные

Для того что бы исследовать зависимость оползней от выпавших осадков нам нужны две вещи: оползневые случаи и осадки с развитием во времени. С осадками всё относительно просто. Во Владивостоке действует одна из самых старинных станций (WMO Index 31960). Исторические метеорологические данные практически непрерывны с 01 января 1914 года. Архив наблюдений содержит ежедневные измерение температуры воздуха, количества осадков, температуры почвы на глубинах. Кстати, если вам нужны метеоданные для ваших исследований или проектов, то вот сайт с ними.

Набор данных (ката в новом редакторе нет, извините):

С оползневыми случаями всё немного сложнее. Готового набора нет. Нам пришлось очень муторно и долго изучать новостные сводки, отмечая даты с оползнями. Каждый найденный оползневой случай соответствует отдельному дню в ряду наблюдений как бинарная переменная. В итоге получилась вот такая табличка. Все остальные дни были отмечены как неоползневые. Ввиду отсутствия объективных данных об оползнях до 1946 года, ряд наблюдений осадков для последующего анализа был сформирован с 01 января 1946 года по 01 апреля 2020 года.

Хорошо! Но с осадками всё не так просто. Дело в том, что на почву влияет не только дождь, который идёт прямо сейчас (Daily Rainfall или DR). Общий «параметр» увлажнения грунта зависит и от предшествующих осадков. И здесь у нас два варианта для создания производных переменных. Первый — это кумулятивное количество осадков c начала года (Cumulative Precipitation или CP). Второй параметр — это предшествующие осадки в определённом временном окне (Antecedent Rainfall или AR). AR придумали наши коллеги из Южной Кореи. Интуитивно, осадки выпавшие два месяца назад должны меньше влиять на некое триггерное увлажнение грунта после которого возникает оползень. В оригинальной работе параметр называется Inter event time definition (IETD) и рассчитывается он каждый час, так как в Южной Корее метеорологические данные более детальные. Параметр IETD в настоящем исследовании может варьироваться от 24 до 46 часов. Это связано с тем, что исходные данные представляют количество осадков за одни сутки и точно определить параметр IETD в часах по имеющимся историческим метеоданным невозможно. В общем случае, безосадочный период фиксируется отсутствием осадков (Rainfall per day = 0) в предшествующие сутки. То есть, если дождь лил непрерывно 7 дней, то AR будет учитывать кумулятивное значение за этот период. Если был 1 «сухой» день до оползня, то AR=0.

Разведочный анализ (Exploratory data analysis, EDA)

Scientific Dataist’ы скажут нам, что данные нужно покрутить, повертеть. Ну что ж, приступим. Для начала взглянем на распределение осадков по месяцам.

Распределение осадков по месяцам с января 1946 по апрель 2020 для города Владивосток

Здесь мы видим наибольшее количество осадков с мая по октябрь, что несюрприз для муссонного климата Владивостока. Наибольшее количество дождей выпадает в августе. А та аномальная точка на 500+ миллиметров — это как раз наш август 2019-го.

Теперь посмотрим на данные сверху:

Оползневые случаи с января 1946 по апрель 2020 для города Владивосток

Здесь у нас всего 34 оползня и с явным отрывом лидирует июль и август. Корреляция с осадками на лицо. Стоит отметить, что минимальное количество осадков за сутки для оползневого случая равно 54 мм. Возможно, это значение является критическим для оползневых процессов в нашем городе.

Матрица ковариаций, рассчитанная для всех переменных в результирующем наборе данных, показывает высокую положительную корреляцию температурных переменных 0.45≤ρ≤0.99. С другой стороны, переменные, связанные с осадками, имеют слабую положительную корреляцию между собой 0.058≤ρ≤0.14. В данных наблюдается выраженная положительная корреляция количества осадков в сутки и оползневыми случаями ρ=0.43.

Ковариационная матрица по всему набору данных

Таким образом, количество осадков в сутки и производные кумулятивные значения являются хорошими кандидатами для построения регрессионной модели. Среди температурных показателей достаточно выбрать только одну переменную, например среднюю дневную температуру (в финальном варианте мы её исключили). Построим отдельную матрицу только для осадков:

Ковариоционная матрица для «осадочных» переменных

Ну и конечно, выпускники факультетов математической статистики скажут нам, что эта табличка нужна в разведочном анализе:

Сводная статистика по набору данных (кат тут не нужен, табличка ещё пригодится)

Логистическая регрессия

Итак, мы добрались до главного. Тут как то надо обосновывать почему именно логистическая регрессия. Целью исследования было построить некую модель для предсказаний оползневых случаев по метеоданным. Вместе с этим, мы хотели посмотреть как сильно влияют на оползни те или иные переменные. То есть нам важна значимость коэффициентов с оценкой (машинное обучение отменяется). И здесь, кроме logit есть ещё probit. Выбор между этими инструментами — это дело вкуса, если хотите. В геофизике как то принято использовать logit. Кроме этого, logit можно всегда интерпретировать как соотношение шансов, в общем, он теплее.

Запишем отношение шансов формально:

Логистическая регрессия для оползневых случаев. DR, AR, CP — наши «осадочные» переменные. Бетты — коэффициенты регрессии.

В этой формуле нет температурных осадков. Мы их убрали выполнив линейную регрессию по всем переменным из за низкой статистической значимости. То есть p-значение для двустороннего T-теста больше 10% для каждой температурной переменной. Построить регрессию по всем переменным и убедится в этом можно на основе этого скрипта. Кстати, мы использовали statsmodels. А любителям R рекомендую вот эту статью)

Теперь нам нужно обсудить очень важный момент. Если мы посмотрим на сводную статистику (табличка выше, блок вся дата), то увидим, что оползневых «дней» всего 33 из 27119 (строка Landslide). Кроме того, что данные жутко несбалансированы, это может привести к проблеме quasi-complete separation (варианты перевода приветствуются). Подробно эта проблема описана в работе статистов Albert & Anderson в 1984 году. Интуитивно, выполняя оптимизацию целевой функции, модель нам прочертит границу по 54 мм, которую мы и так видим на рисунке с оползневыми случаями. Но нам, конечно, надо точнее. Нужно знать, если за день выпало больше 54 мм, то какая вероятность того что оползни будут в городской черте. Для этого, мы просто отрежем данные ниже 54 мм (таблица выше, блок данные для регрессии).

Запишем этот ход конем формально:

Отрезали данные ниже порогового значения

где

Пороговое значение это 54 мм ± стандартное отклонение в суточных вариациях осадков

Теперь давайте посчитаем наши коэффициенты:

Здесь мы видим, что, во первых статистическая значимость теста максимального правдоподобия в норме, т.е. больше 95% для всех переменных. А во вторых интересна оценка «вклада» каждого осадочного параметра в прогнозируемое значение возникновения оползня. В случае линейной регрессии, мы можем просто посмотреть на сами параметры, какой больше тот и круче. В случае logit мы должны продифференцировать относительно каждого коэффициента и получить так называемые маржинальные значения (marginal effects). И на первом месте, по влиянию на вероятность оползня, у нас дневные осадки. Следом, на порядок ниже, идут предшествующие осадки. И самый слабый (но значимый!) — это накопленные осадки с начала года.

Круто! Один результат есть. Но нам ещё нужна предиктивная модель. Для этого нам нужно определить порог срабатывания нашей logit (P) модели. По умолчанию, для logit всё что больше 0.5 и есть кейс. Но мы можем немного поварьировать этот порог и посмотреть на значение наших метрик:

Сбалансированная точность, потому что всего 33 позитивных случая из 200

Наши друзья precision и recall. Достаточно распространённые метрики в классификации. Интуитивно, precision это сколько оползней мы вообще правильно обнаружили, а recall это сколько «пропусков» для положительных вариантов. TP — это правильно предсказанный случай, FP — неправильно, FN — неправильно не предсказанный.

F1 оценка всё учитывает. Будем ориентироваться на неё.

Ну и scikit learn нам в помощь:

Значение метрик оценки для классификации оползней при разных порогах (predicted threshold)

Методом перебора по всем 33 случаям мы определили оптимальный порог в 0.24 на наших данных. Здесь надо сказать, что мы не делили данные для раздельного обучения и валидации, так как их очень мало.

Матрица несоответствий при классификации оползней на обучающей выборке моделью Logit (P) (200 samples, 0.24 predicted threshold):

Производительность модели Logit (P) как классификатора оползневых случаев по обучающей выборке метеоданных (200 samples, 0.24 predicted threshold)

Ну и, наконец, давайте упростим нашу формулу, что бы вот подставил значения осадков и получил вероятность, потом посмотрел больше ли оно 0.24 и предсказал!

Финальная предиктивная модель оползней на основе данных об осадках

В заключение нашего статистического экскурса стоит сказать, что оползни зависят не только от осадков, но и от других физических параметров, которых у нас не было в распоряжении. Финальная модель, конечно, далека от совершенства. Стоит выполнять раздельную регрессию и валидацию и так далее. Но в качестве первого приближения вполне годится, так как ряд наблюдений это десятки лет.

© Habrahabr.ru