Method of processing data from gamma-spectrometric field measurements in geophysics using the concept of stationary intervals

Sergey Kataev; Sophia Kataeva; Andrey Miller; Evgeniy Ermolin

doi:doi:

Home / Journals / Vestnik AGGE / Issue 4 / Method of processing data from gamma-spectrometric field measurements in geophysics using the concept of stationary intervals

Method of processing data from gamma-spectrometric field measurements in geophysics using the concept of stationary intervals

Submit manuscript Download PDF
Text

To cite

Citations:

Method of processing data from gamma-spectrometric field measurements in geophysics using the concept of stationary intervals

Journal: VESTNIK AGGE № 4 , 2025

Rubrics: ENGINEERING AND MINING GEOPHYSICS

VAK Russia 1.2.2

Sergey Kataev ¹

Sophia Kataeva ²

Andrey Miller ³

Evgeniy Ermolin ⁴

Author and publication information

Authors:

1. TSPU (kafedra fiziki i metodiki obucheniya fizike, professor)
from 01.01.1995 to 01.01.2025
Tomsk, Russian Federation

2. National Research Tomsk State University

Tomsk, Tomsk, Russian Federation

3. VNIIGIS

Saint-Petersburg, St. Petersburg, Russian Federation

4. GM-Service

Saint-Petersburg, St. Petersburg, Russian Federation

Type:

Article

Pages:

from 93 to 104

Status:

Published

Received:

01.12.2025

Accepted:

29.12.2025

Published:

29.12.2025

Subject area:

VAK Russia 1.2.2

Language:

Russian

Keywords:

gamma spectrometer, gamma mapping, graphs, random events flow, method for detecting structures

Abstract and keywords

Abstract:
The features of processing field gamma-spectroscopic data by method of revealing structures using the notion of stationarity interval in various geophysical studies are considered.

Keywords:
gamma spectrometer, gamma mapping, graphs, random events flow, method for detecting structures

Text

Text (PDF): Read Download

Введение

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

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

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

Подход к обработке гамма-спектроскопических данных

Можно попытаться улучшить результаты обработки данных измерений, если предположить, что существуют интервалы (пространственные или временные), в пределах которых параметры гамма-поля не меняются. Геологически это будет соответствовать пластам определенных пород. Назовем такие интервалы интервалами стационарности, то есть интервалами, в пределах которых случайный процесс поступления гамма-квантов в детектор можно считать стационарным. Но статистическая ошибка измерений при таком подходе должна быть меньше, чем разница между отдельными пластами, границы которых мы пытаемся зафиксировать. Введение понятия интервала стационарности (ИС) можно считать оправданным, учитывая, что в каждой точке измерения поток радиоактивных частиц можно считать пуассоновским, исходя из условий его формирования. Такой подход был нами применен при анализе скважинных измерений методом гамма-каротажа [1], а также при наземных исследованиях [2] с использованием одного из алгоритмов метода выявления структур (МВС) [3]. А именно алгоритм APFLOW [4], который ранее нами применялся для определения параметров случайного потока событий, являющегося стандартным объектом в теории массового обслуживания. Данный алгоритм позволяет выявлять структуру потока, т. е. представлять исходный массив данных в виде последовательности интервалов, каждый из которых характеризуется своим набором параметров. Измеренные количества гамма-квантов в секунду при погружении датчика в скважину либо при перемещении его по поверхности представляют собой ряд случайных чисел $n_1,n_2,…,n_m.$ В каротажных исследованиях в каждой точке датчик принимает поток радиоактивных частиц, идущий из горизонтального слоя, и этот поток, исходя из физики процесса, можно считать пуассоновским. Т. е. если проводить измерения в одной точке определенное (не малое) время, то разности времен срабатывания датчика будут подчиняться пуассоновскому распределению. Это утверждение будет справедливым для любой точки, находящейся на произвольной глубине в скважине. Причем число срабатываний датчика за единицу времени в точках, относящихся к одному и тому же слою, будет примерно одинаковым, так как, проводя измерения в разных точках одного слоя, мы фиксируем интенсивности $n_1,n_2,…,n_m$ фактически одного и того же пуассоновского потока, полученные для разных его участков. Учитывая свойства пуассоновского потока, а именно стационарность и однородность, можно утверждать, что эти интенсивности должны быть приблизительно одинаковыми. А это означает, что плотность распределения значений $n_1^j,n_2^j,…,n_{m_j}^j$ , где $j$ — номер слоя, не будет являться пуассоновским распределением. И ясно, что чем больше время, за которое определяются интенсивности $n_1^j,n_2^j,…,n_{m_j}^j$ , тем более близкими друг к другу должны быть эти значения и тем меньше должна быть дисперсия этого ряда. Эти рассуждения позволяют сформулировать критерий определения границ слоев, основанный на идее минимума дисперсии для $n_1^j,n_2^j,…,n_{m_j}^j$ , относящихся к одному слою.

Исходя из вышесказанного, задача определения границ слоев при каротажных измерениях может быть сформулирована следующим образом: найти номера отсчетов $m_{j}, j = 1, . . ., k$ , где $k$ — количество слоев, разделяющих различные «интервалы стационарности». Мы сохраняем это принятое в теории массового обслуживания понятие, понимая под ИС последовательность идущих подряд номеров ряда, интенсивности которых примерно одинаковы, а суммарная дисперсия выделенных таким образом интервалов — минимальна. Приведенные рассуждения полностью применимы и в случае перемещения датчика по поверхности.

Метод нахождения интервалов стационарности

Для нахождения интервалов стационарности мы используем МВС [3]. Для этого необходимо и достаточно построить матрицу отношений между объектами. В качестве объектов в задаче выступают номера отсчетов $i=1,...,m,$ а в качестве признаков — накопленные значения интенсивности $t_i, i=1,...,m.$ Для того чтобы поток был упорядоченным, значения строились рекуррентно, т. е.

$t_i=t_{i-1}+n_i.$ (1)

При этом, как видно из (1), разность между двумя соседними значениями равна просто $n_i$ . В качестве характеристики отношения между объектами $i$ и $j$ удобно выбрать «расстояние» между ними, т. е. величину $d_{ij},$ , определяемую по (2):

(2)

Определенные таким образом элементы образуют верхнетреугольную матрицу $D=[d_{ij} ]_{m×m}$ . Тогда формализация МВС для данной задачи может быть изложена в следующей редакции. Сопоставим каждому i-му событию потока, $i=1,m$ , вершину с тем же номером i некоторого графа G(V,Е), где V — множество вершин графа мощности N, а Е — множество пар вершин, т. е. множество ребер графа. Стоимостью каждого ребра графа, связывающего i-ю и j-ю (i < j) вершины, будем считать величину $d_{ij}$ , т. е. соответствующий элемент матрицы D. Мощность множества ребер графа составит $m(m–1)/2.$

Далее, поскольку МВС отыскивает на графе заранее заданные структуры, надо определить, что в данной задаче является структурой. Определение структуры можно сформулировать исходя из понятия интервала стационарности. А именно, если группа событий с номерами $i_1^p,i_2^p,...,i_k^p$ составляет участок стационарности потока с интенсивностью $λ=λ_p$ , то это означает, что:

а) среднее значение стоимостей всех ребер, связывающих между собой вершины с номерами $i_1^p,i_2^p,...,i_k^p$ , примерно равно интенсивности $λ_p$ :

$\frac{2}{k(k-1)} ∑_{(i,j)∈I}d_{ij} ≅λ_p$ ,

где I — множество всевозможных упорядоченных по номерам пар из $\{i_1^p,i_2^p,...,i_k^p\}$ ;

б) номера $i_j^p$ , j=1,k идут подряд, т. е. $i_j^p=i_{j-1}^p+1$

в) все $d_{ij}, i,j∈I$ принадлежат в основном некоторому интервалу $(λ_p–\Delta λ, λ_p+\Delta λ)$ .

Необходимо отметить, что оценка интенсивности, полученная по всей группе событий, т. е. $(i_1^p,i_k^p)$ -я компонента матрицы D, должна давать наиболее точную оценку $λ_p$ . Изложенные соображения определяют идею алгоритма выделения участков стационарности потока. В самом деле, при выборе удачного интервала $(λ_p–\Delta λ, λ_p+\Delta λ)$ и удалении из графа G всех ребер, стоимость которых не попадает в указанный интервал, в графе выделятся компоненты связности G_k=(V_k,E_k), соответствующие участкам стационарности, с примерно одинаковой интенсивностью (G_k является подграфом графа G). Разумно ожидать, что участку стационарности в наибольшей степени отвечает такой подграф, в котором, во-первых, максимальное число ребер, а во-вторых, каждая вершина должна быть соединена с наибольшим числом вершин, т. е. степени вершин должны быть примерно одинаковыми. После выделения нескольких ИС матрица D очищается от «шума» путем удаления тех элементов, номера которых относятся к разным ИС.

На рисунке 1 схематично показана структура матрицы D для данных, распределенных по трем ИС. Элементы матрицы, помеченные красным, относятся к одному из трех ИС, а желтым — образуют «шум».

Рисунок 1. Схематическое изображение структуры матрицы оценок D

Для понимания того, как выбирать интервал сканирования $(λ_p -∆λ, λ_p +∆λ)$ и для оценки качества найденных ИС полезно использовать огибающую гистограммы функции плотности распределения элементов матрицы D.

В качестве иллюстрации применения МВС (алгоритма APFLOW) рассмотрим имитирующие гамма-каротаж данные модельного эксперимента [1], представляющие собой регистрируемое датчиком число импульсов в секунду (интенсивность). На рисунке 2 помещены исходные данные, а на рисунке 3 показаны гистограммы исходной матрицы D_in, построенной для этих данных по (2), и результирующей матрицы D_fin после проведения процедуры очистки от шума.

Рисунок 2. Исходные данные модельного эксперимента

Рисунок 3. Гистограммы исходной матрицы D_in (пунктирная линия) и полностью очищенной матрицы D_fin (сплошная линия)

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

Можно полагать, что ИС, имеющие приблизительно одинаковые средние значение интенсивности, относятся к одному и тому же состоянию потока. Поэтому удобно конкретный поток данных характеризовать набором, состоящим из найденного числа ИС, распределенных по M состояниям. Так, поток данных приведенного примера содержит пять ИС, распределенных по четырем состояниям.

Таким образом, алгоритм APFLOW дает возможность, выявляя ИС, получить структуру потока данных — упорядоченный набор номеров членов ряда $m_1 , m_2 , …$ , которые определяют границы ИС.

Результаты работы предлагаемого подхода демонстрируется следующим примером.

Пример

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

На территории Сибири была произведена пешеходная съемка с использованием гамма-спектрометра GM-GAMMA-1000 с объемом кристалла около 1000 см³ и приемником GNSS. Всего было пройдено около 60 трасс с общим количеством точек порядка 150 000. В идеале в каждой пространственной точке результатом измерения гамма-спектрометром является энергетический спектр, который показывает, сколько гамма-квантов было зафиксировано в определенном интервале энергии за единицу времени. Основные источники радиации — это изотопы атомов K, U и Th. В реальности при больших объемах проводимых исследований получаемый спектр в отдельной пространственной точке не дает возможности определять аномальные содержания указанных радиоактивных элементов в силу небольшого числа поступивших импульсов в единицу времени. По этой причине, как и в проделанных нами ранее исследованиях [2], для последующего анализа исходные данные группировались в пять интегральных каналов: k₁, k₂, k₃, k₄, k₅. Количество импульсов в каждом интегральном канале получалось суммированием по идущим подряд аппаратурным каналам (от 1 до 4084). Причем границы каналов k₂, k₃, k₄, k₅ подбирались таким образом, чтобы выделенные области энергии включали в себя основные пики, соответствующие гамма-излучению атомов Cs, K, U и Th соответственно.

Рисунок 4. Спектр модельных источников гамма-излучения, нормированный на время измерения, в логарифмическом масштабе по вертикальной оси

Таким образом, каждой пространственной точке (x, y) исследуемых профилей присваивался упорядоченный набор из пяти чисел. Учитывая, что более объективным показателем распределения числа импульсов по интегральным каналам k₁ — k₅ является не абсолютное значение — количество импульсов в данном интегральном канале, а его относительная величина — доля этих импульсов относительно полного числа импульсов, то дальнейшее исследование проводилось именно для долей каналов калия (k₃), урана (k₄) и тория (k₅) относительно полного числа импульсов (формула 3).

$dol_i (x,y)=\frac{Imp_{k_i} (x,y)}{∑_{j=1}^5Imp_{k_j} (x,y)}, i=3,4,5,$

где $Imp_{k_i} (x,y)$ — суммарное число импульсов интегрального канала $k_i$ .

Результаты

Было проведено разбиение на ИС каждой трассы на указанной территории для данных по каналу K₄ (уран) c целью получения представления о распределении на исследуемой территории этого изотопа и нахождения возможных аномальных областей. Было выделено порядка 1000 ИС. В качестве примера на рисунке 5 приводятся исходная и очищенная гистограммы для одной из трасс (tr20). Видно, что алгоритм APFLOW выделил на этом профиле пять ИС с интенсивностями от 1,6 до 3,9 (в соответствии с максимумами на красной кривой).

Рисунок 5. Гистограмма исходной (пунктирная линия) и полностью очищенной (сплошная линия) матрицы D для трассы tr20

Исходя из анализа распределения интенсивностей, полученных ИС и целей исследования, было выделено три обобщенных класса со значениями относительной интенсивности 0–2,2; 2,2–3,15; 3,15–5,0. Таким образом, в каждом классе будет содержаться явно выраженный пик относительной интенсивности излучения в окне урана по результатам работы алгоритма APFLOW.

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

Из рисунка 6 видно, что выделяются две области, относящиеся к 3-у классу, с повышенным содержанием доли урана по спектру. Суммарный спектр точек, располагающихся в этих областях, также показывают наличие урановой аномалии. На рисунке 7 показаны сглаженные спектры, построенные по 3000 точек, входящих в каждый из трех выбранных классов. Доля урана в области E = 1760 кэВ для точек класса 3 практически в два раза превышает значения для класса 2 и в три раза фоновые значения (класс 1).

Таким образом, можно утверждать, что аномалии не носят аппаратурный характер и связаны именно с изменением спектра регистрируемого излучения. То, что аномалии оказались приурочены к определенным профилям, можно объяснить тем, что первоначально направление профилей было выбрано вдоль геологических границ.

Также аномалии не связаны с рельефом или с общим радиоактивным фоном. Квадрат коэффициента корреляции между номером класса и высотой по GPS (данные приемника в спектрометре) составляет 0,0192, между номером класса и МЭД — 0,0199.

Точки, относящиеся к третьему классу, лежат на линии тектонического контакта (красная линия на геологической основе) и линии несогласного залегания пород (черная линия). Точки первого и второго класса можно отнести к рядовым (фоновым) значениям по содержанию урана.

Рисунок 6. Распределение точек исследуемой территории по классам

Рисунок 7. Спектры для трех классов по излучению в окне урана: зеленый — первый класс, красный — второй класс, черный — третий класс. Стрелка указывает на положение уранового пика при E = 1,76 МэВ

Выводы

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

Благодарности

Авторы благодарят сотрудников фирмы «ООО «ДЖИ М Сервис» за предоставленные геофизические и геологические данные.

References

1. Kataev S. On a New Method of Detecting Layer Boundaries at the Borehole Measurements by Gamma Logging / Kataev S., Kataeva S., Miller A. // Near Surface Geoscience Conference & Exhibition 2019, September 2019, The Hague, Netherlands. — P. 8–12. — https://doi.org/10.3997/2214-4609.202152102.

2. Kataev S. G. Primenenie statisticheskikh metodov obrabotki spektra estestvennogo gamma-izlucheniya na primere s`emok na poligone tverdykh bytovykh otkhodov v Smolenske / Kataev S. G., Kataeva S. S., Miller A. A., Ermolin E. Yu. // Engineering Geophysics 2024: proceedings. — Kazan, 2024 (in Russ.).

3. Kataev S. G. Algorithmic and software of analytical and noninvasive control of parameters of the environment based on a method of detection of structures. The thesis. Doctor of technical science 05.11.13. National research Tomsk Polytechnic University. — Tomsk, 2012 (in Russ.).

4. Bekkerman E. N., Kataev S. G., Kataeva S. S. Heuristic approximation method for a random flow of events by an mc-flow with arbitrary number of states / Bekkerman E. N., Kataev S. G., Kataeva S. S. // Automation and Remote Control. — 2013. — Vol. 74. — № 9. — P. 1449–1459.

Submit manuscript Download PDF
Text JATS XML