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

В связи со сложностью задачи, нам потребуются серьезные статистические методы, такие, как полиспектральный анализ. Что интересно, такой анализатор у нас уже есть… в голове. Это легко подтвердить тем, что мы способны различать так называемый «малиновый звон» колоколов — этот эффект не проявляется на спектре, зато отлично виден на биспектре. Большинство людей отличает колокола с малиновым звоном, для этого даже не требуется наличие музыкального слуха. Опытный геолог, занимающийся визуальной дешифровкой космоснимков, способен вручную выделить на них элементы, сопутствующие различным погребенным геологическим структурам. Мы же, как обычно, воспользуемся вычислительными методами и построим 3D геологические модели для автоматизированного анализа.

Введение

Благодаря доступности платформы Google Earth Engine мы получили возможность статистического анализа огромного массива спутниковых данных и это изменило все. В самом деле, анализируя один оптический или радарный снимок, где интересующая нас геологическая картина замаскирована помехами как самого снимка, прошедшего многостадийную обработку, так и толщей случайных пород, мы мало что можем сделать для улучшения результатов — спектральный синтез гравики высокого разрешения и последующее решение обратной задачи только добавляют помехи. Зато в случае с ансамблем в сотню таких радарных снимков мы можем эффективно удалить помехи (в силу центральной предельной теоремы ожидаем их гауссовово распределение) и изучать ранее скрытое от нас геологическое строение. Не будем здесь углубляться в фундаментальную радиофизику, отмечу лишь, что спектры высших порядков от гауссова сигнала равны нулю — значит, используя биспектральный анализ, мы можем исключить как помехи спутниковой съемки, так и влияние толщи случайно распределенных осадочных пород. Попросту говоря, это эффективный метод выделения не-гауссовых составляющих сложных процессов и сигналов. В комментариях к предыдущей статье меня спрашивали про применимость (сканированных) карт для подобных задач — очевидно, что они полностью бесполезны для статистического анализа. Для специалистов добавлю, что триспектры (и спектры высших порядков) тоже подходят, но и ошибка вычисления моментов высших порядков растет, так что даже с биспектром это представляет непростую техническую задачу. Кстати, фрактальные геологические объекты такие, как рудные жилы и зоны, выделяются в биспектре геологической модели — поскольку фрактальность подразумевает связь масштабов, это означает связь пространственных частот и автоматически выделяется на полиспектрах. Также для специалистов замечу, что для меня очень полезной оказалась идея вычисления фрактальности через производные дробных порядков.

Мы уже использовали амплитудный биспектральный анализ для Индонезии — модели вертикальных аномалий получены именно этим методом. Но для Индонезии это не так важно, поскольку там для поиска достаточно модели плотности, полученной методом инверсии детального гравитационного поля, синтезированного на основе глобальной гравитационной модели низкого разрешения, рельефа ALOS 30 м и радарных снимков Sentinel-1 10 м. Как показано в предыдущей статье, построенная модель машинного обучения обеспечивает высокую точность, а дополнительный анализ позволяет нам непосредственно увидеть геологическое строение и понять логику работы нашей модели машинного обучения. С другой стороны, детальность биспектральной модели выше за счет снижения уровня шумов (это и визуально заметно на представленных картинках), а это важно всегда. По биспектру плотности можно восстановить исходную модель плотности с меньшим уровнем шумов, но на практике я этим не пользуюсь (в том числе, потому, что у нас нет способа оценить точность такой модели).

Выбираем территорию

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

Для целей статьи нам подойдет относительно небольшой участок (чтобы удобно было рассматривать модель), где априори известно наличие золота и других минералов. Почему я говорю про наличие других минералов? Для первичного анализа нам необходимо точно знать, где золото есть, и где его нет, а для этого как нельзя лучше подойдет месторождение любых других ископаемых — очевидно, что добывая, скажем, бокситы (алюминий), нельзя было бы не заметить наличие там же золота. Для качественной классификации оптимальный вариант подготовки данных заключается в наложении карты известных месторождений на выходы руды по нашей модели и выбора для классификации их пересечения, но для предварительного анализа мы ограничимся контурами известных месторождений, полагая все воксели модели внутри контуров золотоносных участков продуктивными и все воксели внутри не золотоносных участков не продуктивными. Если у нас слишком малое количество известных вокселей в одном из классов (продуктивном или не продуктивном), то можно несколько уменьшить размер продуктивных контуров или увеличить размер не продуктивных (ни в коем случае не наоборот, по понятным причинам).

Итак, возьмем на геологической карте СССР известное золоторудное месторождение в районе Егорьевского рудного узла на реке Суенга, рядом с которым находится месторождение молибдена. Смотрите региональные карты на сайте ВСЕГЕИ Материалы электронного издания ГГК-200/2, Материалы по листу N-45-XIII При наличии геологической модели региона таких карт вполне достаточно для наших целей. Для выбранного месторождения на нашей модели мы легко находим соответствующую рудоконтролирующую структуру подковообразной формы:

Здесь участок внутренней границы структуры определяет залегание золота и соответствующий ему внешний — молибдена, цвета контуров золоторудного и молибденового участков (показаны кольцами на поверхности) совпадают с соответствующими рудными структурами на биспектральной модели (вертикальных аномалий плотности). Как видим, контуры с карты полезных ископаемых как-то странно сдвинуты от структуры на нашей модели. Эта загадка легко разрешается, если добавить к модели еще и рельеф:

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

Структурный анализ

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

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

Из открытых данных известно, что рудные горизонты на этой территории залегают на глубинах 20-60 и 80 метров, и как раз для этих интервалов глубин мы нашли ожидаемые нарушения геологических структур.

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

Таким образом, полученные графики подтверждают все наши структурно-геологические предпосылки и свидетельствуют о верном выборе продуктивного и непродуктивного участков.

Построим соответствующие графики и для биспектра нашей модели плотности (модели вертикальных аномалий):

Как видим, у нас есть все инструменты для ручного анализа моделей и получения прогноза для золоторудных участков. Но мы не будем останавливаться на достигнутом и автоматизируем такой анализ с помощью методов машинного обучения.

Анализ с помощью машинного обучения (классификатора)

Для площадного анализа воспользуемся классификатором с гауссовым ядром, поскольку рудные области обладают типичным геологическим строением в плоскости поверхности. Дело в том, что мы ищем не отдельные (фрактальные) рудные жилы, а состоящие из групп таких жил геологические объекты. Как мы уже говорили, геологи не умеют работать непосредственно с фрактальными объектами, да и прекрасно без этого обходятся, поскольку в масштабе геологического исследования работают со скоплениями фрактальных структур, подчиняющихся гауссовой статистике. Итак, несмотря на фрактальное внутреннее строение рудонесущих структур, в плоскости выхода группы рудных жил наблюдается максимальная концентрация минерала в центре этой группы, а по радиальной координате плотность выходов жил и концентрации спадают приблизительно по гауссиане. При этом, крупные жильные выходы окружены более мелкими, образуя фрактальный рисунок гауссиан на поверхности. И на бумажных картах можно увидеть такую картину вокруг крупных месторождений, а при машинном анализе мы выделяем намного больше разномасштабных рудных зон. Для визуализации полученного прогноза мы исключаем мелкие объекты, при этом сохраняя их для последующей обработки такой, как вычисление индекса фрактальности золотоносных участков. Группируя найденные прогнозные участки по этому индексу, и сравнивая их с известными месторождениями, можно также сделать прогноз запасов руды и не только.

Таким образом, для нашей задачи отлично работает классификатор SVC с ядром RBF. Хотя мы делаем классификацию поверхности территории, но при этом анализируем трехмерную модель плотности, где сами выходы руды могут находиться как на поверхности, так и под ней. Ничто не мешает построенный классификатор скорректировать по вертикальной координате, получив 3D карту, но на практике этого не требуется, поскольку и геологи и заказчики предпочитают работать с обычными картами на поверхности и отметками глубин на них.

Заключение

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

И несколько слов о том, о чем мы не поговорили. Модели плотности всем хороши, кроме того, что они ничего нам не говорят о характеристиках вещества. В то же время, если совместить модели по данным гравиметрии и магнитометрии, то разница между ними расскажет многое о химических процессах и составе пород. Можно ли извлечь соответствующую информацию со спутниковых снимков? Да, с помощью фильтра Габора. Дело в том, что, как минимум в региональном масштабе, магнитные аномалии связаны с нарушениями строения и линеаментами (выделенными направлениями простирания геологических структур и нарушений). Соответственно, ориентированные системы нарушений и аномалии отлично выделяются фильтром Габора. Итого, глубинные магнитные аномалии приводят к эффектам на поверхности, выделение которых асимметричным ориентированным фильтром и последующий анализ позволяют косвенно выявить исходные магнитные аномалии и сделать соответствующие выводы. Как раз в Сибири наблюдаются сложные системы разломов, которые плохо поддаются линеаментному анализу. На рассматриваемой нами территории Егорьевского рудного узла можно выделить три системы разнонаправленных разломов, присутствующих на всей территории.