Обработка данных гамма-каротажа

«Обработка» — это очень громко сказано. Я напишу как я просто «поигрался» с данными, которые снимал в течение 30 минут с Модуля Гамма-Каротажа (МГК). Мы немного посчитаем и построим графики.

«Свежеиспечённый» Модуль Гамма-Каротажа (МГК) перед вводом в эксплуатацию требует проверки. При испытании мне показалось, что после включения прибора его характеристики «плывут» какое-то время. В народе этот непродолжительный период называется обычно «прогревом». Что, в общем-то, близко к истине, ведь рассеиваемая модулем мощность составляет примерно 2 Вт — это вполне ощутимое количество, если учесть, что эта мощность рассеивается на плате размером 112 х 27 мм².

Я собрал установку, подключил Регистратор Каротажных Сигналов (РКС) к компу. Можно было, конечно, направить поток данных сразу в файл, но в пакете, приходящем от модуля МГК много ненужных мне данных, поэтому я написал предельно простую прогу на Python.

Данные от РКС поступают в комп через USB-порт в виде ASCII-строк. Один пакет — одна строка. Внутри строки данные разделены пробелами. Всё очень просто — прога должна принять строку, разбить ее на слова, и вывести в стандартный поток только данные по гамма-каротажу и данные по напряжению питания. Ну, я еще в этот вывод добавил дату-время. Дело в том, что я запустил прогон на выходные, а остановлю его в Пятницу. Установка будет «молотить» чуть более двух суток. Поэтому, если вдруг что-то пойдет не так, то хотелось бы знать когда это «не так» случилось.

Вот эта страшная супер-пупер прога:

#!/usr/bin/env python
#coding:utf-8


import serial
import time

for line in serial.Serial('/dev/ttyUSB0', 115200, timeout = 10):
  try:
    #print line
    l = line.split()
    V = 30.0 / 700 * int(l[2])
    tm = time.localtime(time.time())
    print "%02d:%02d:%02d:%02d - %s %5.1f %s" % (tm.tm_mday, tm.tm_hour, tm.tm_min, tm.tm_sec, l[3], V, l[1])
  except IndexError:
    pass

К стати, прога дурная — через 12 часов работы выкинуна исключение и вывалилась. Подробности см. в самом конце.

К сожалению, я не успел подготовить комп для удаленного доступа. Дело было в Пятницу. Вечером. Я был уставшим. Помниться, работая в Спектроне, я мог мониторить свои установки в любое время суток и из любой точки — хоть из дома, хоть из Канады — Интернет позволяет сократить расстояния, а Линукс обеспечивает безопасность и легкость управления. Ну, не важно! В общем, на данный момент единственное решение было писать данные в файл и не выпендриваться с удаленным доступом.

Так вот, 30-минутное начало этого файла я скопировал себе на флешку, и сейчас мы с ним поиграемся.

Во-первых, полученные данные в файле выглядят так:

20:17:47:07 - 29  33.1
20:17:47:08 - 29  31.6
20:17:47:09 - 28  30.8
20:17:47:10 - 28  30.6
20:17:47:11 - 30  30.6
20:17:47:12 - 25  30.7
20:17:47:13 - 30  30.6
20:17:47:14 - 35  30.6
20:17:47:15 - 21  30.6
20:17:47:16 - 44  30.6
...

и далее чуть более 1800 таких же строк.

В первой колонке указаны дата-время в формате:

число_месяца:час:минута:секунда

Во второй колонке — находится дефис ‘-‘. Только не спрашивайте, зачем я его туда втолкнул! — Не знаю! Вечер. Пятница. Усталость…

В третьей колонке находится собственно сами значения от гамма-модуля.

Затем идут два пробела. Это еще одна моя плюха, немного погодя поймём почему это есть плюха.

Далее идет значение напряжения (в Вольтах) на входе модуля. Оно нам тоже понадобится для построения графика.

Чтобы было удобно работать, я решил разбить этот файл на два независимых файла — в одном будут содержатся данные по гамма-каротажу, в другом — данные по напряжению питания. Для этого я выполнил две консольные команды:

$ cut -d' ' -f3 gamma0.txt > gamma
$ cut -d' ' -f5 gamma0.txt > uline

, здесь gamma0.txt, gamma и uline — это имя входного файла, и имена файлов с гамма- и элктрическими данными.

Консольная команда cut «вырезает» из входного файла заданный столбец и выводит его в стандартный вывод. Для гамма-данных это третий столбец (-f3), а для данных напряжения питания — это пятый столбец (-f5).

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

Опция -d’ ‘ — задает разделитель столбцов. В нашем случае это есть пробел.

В результате работы этих двух команд в директории появляются два файла с таким содержимым:

файл gamma:

29
29
28
28
30
25
30
35
21
44
...

и файл uline:

33.1
31.6
30.8
30.6
30.6
30.7
30.6
30.6
30.6
30.6
...

Теперь займемся построением графиков. Запустим gnuplot и в нем выполним команду для построения графика — сначала гамма:

gnuplot> plot 'gamma' u 1 w l

Но, вот, не задача — по оси абсцисс график начинает не от нуля. Не проблема! Зададим диапазон в ручную:

gnuplot> set yrange[0:50]

и перерисуем заново:

gnuplot> replot

Отлично! Но хотелось бы еще видеть линии сетки. Сделаем еще пару движений руками:

gnuplot> set grid
gnuplot> replot

Ну, вот, кажись получился шедевр. Теперь можно скриншотить и публиковать:

gamma

Однако, на графике присутствует много шумов, из-за них сложно понять — есть ли дрейф или нет. Вроде как нет. Но полной уверенности тоже нет. Давайте сгладим график, точнее «сгладим» данные.

Для этого я написал еще одну простецкую прогу:

#!/usr/bin/env python
#coding: utf-8

import sys

def smooth(X, n) :
  Y = []
  i = 0
  while i <= (len(X) - n) :
    dX = X[i:(i + n)]
    Y.append(float(sum(dX)) / n)
    i += 1
  return Y

#############################################
X = []
for x in open(sys.argv[1]).readlines() :
  X.append(int(x))

#print len(X), min(X), max(X)

Y = smooth(X, int(sys.argv[2]))

for y in Y :
  print y

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

Далее для диагностики я добавил строку для распечатки длины списка, а так же минимального и максимального значений в списке. Сейчас эта строка закомментирована, так уже не нужна.

Затем вызывается функция smooth, которая возвращает список сглаженных значений. И наконец вызывается цикл для распечатки значений списка.

Функция smooth принимает два параметра — список исходных значений и длину, по которой этот список усредняется.

Я назвал прогу незатейливо — smooth.py. Прога вызывается так:

$ ./smooth.py gamma 100

, здесь gamma — имя входного файла, 100 — это длина усреднения.

Поскольку, прога вывод осуществляет в стандартный поток вывода, а нам желательно строить график на основе данных их файла, то вывод проги следует перенаправить в файл, например, так:

$ ./smooth.py gamma 100 > gamma-100

, здесь gamma-100 — имя выходного файла

В общем, для тех, кто давно в Линуксе, нет ничего заумного и не понятного.

Я произвел усреднение по 100 и по 500 значений. И вот что получил:

gamma-100

gamma-500

Ну, что можно сказать? — Да, дрейфа по гамма-измерениям нет. Посмотрим, что теперь нам скажет измерение напряжения питания.

Снова вызываем gnuplot и выполняем следующие команды:

gnuplot> set yrange[0:50]
gnuplot> set grid
gnuplot> plot 'uline' u 1 w l

В результате имеем следующий график:

uline

Хорошо видно, что в начале, примерно на протяжение 3-4 минут напряжение «шумит» и немного падает. После чего уже стоит мертво на протяжении всего оставшегося времени.

Вообще, все рисунки кликабельны, можете посмотреть их «поближе».

Я решил более детально рассмотреть начальную область графика. Для этого можно было в gnuplot выполнить пару команд по масштабированию по X и по Y, но я поступил проще. На графике можно придавить правую кнопочку мышки и обвести интересующую вас область. График автоматически заново перерисуется:

uline_starting

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

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

Что бы построить этот график придется написать еще одну прогу, которая должна посчитать статистику на основе входных данных от модуля МГК. Прога тоже не особо заумная:

#!/usr/bin/env python
#coding: utf-8

import sys

def density(X) :
  D = {}
  for x in X :
    if D.has_key(x) :
      D[x] += 1
    else :
      D[x] = 1
  return D

def correct0(D) :
  maxkey = max(D.keys())
  for k in range(1, maxkey) :
    if not D.has_key(k) :
      D[k] = 0
  return D

X = []
for x in open(sys.argv[1]).readlines() :
  X.append(int(x))
  
#print len(X), min(X), max(X)

D = density(X)
D = correct0(D)

#print len(D)

for k in D.keys() :
  print k, D[k]

Запуск проги простой — ей нужно указать имя входного файла и перенаправит ее вывод в файл:

$ ./density.py gamma > gd

Теперь снова вызываем gnuplot и выполняем в нем три команды:

gnuplot> set xrange[0:60]
gnuplot> set grid
gnuplot> plot 'gd' u 1:2 w boxes

В результате получаем вот такой нафиг-график:

gamma-density
Ну, что? Ждем понедельника. В понедельник получим результаты за двое суток наблюдений. А на сегодня — всё! Приятных вам выходных!

UPDATE 23.06.2014-13:27

Питоновская супер-пупер прога, отработав 12 часов, решила устроить саботаж. А может решила, что работать в Субботу в 6 утра — это ниже ее достоинства. В общем, выкинула исключение

Traceback (most recent call last):
  File "./show_gamma.py", line 8, in <module>
    for line in serial.Serial('/dev/ttyUSB0', 115200, timeout = 10):
  File "/usr/lib/python2.7/dist-packages/serial/serialposix.py", line 456, in read
    raise SerialException('device reports readiness to read but returned no data (device disconnected?)')
serial.serialutil.SerialException: device reports readiness to read but returned no data (device disconnected?)

и вывалилась в консоль.

Не знаю, что ей не понравилось! Вроде бы и timeout был с хорошим-таким запасом — 10 секунд, и аппаратура (модуль МГК и регистратор) работают безупречно. Не знаю! В общем, вместо ожидаемых 60 часов записи, у меня получилось всего 12 часов. Ну, плюс еще небольшая получасовая запись сегодня утром.

Графики для этой 12-часовой записи:

Необработанный график количества гамма-импульсов по времени:
gamma-12h
График гамма-импульсов сглаженный окном в 100 точек:

gamma-12h-100

График гамма-импульсов сглаженный окном в 500 точек:
gamma-12h-500

График плотности вероятности:
gamma-12h-density

Хочу обратить внимание на плавность графика. Количество данных (43 тысячи замеров против 1800 в предыдущем графике) — несомненно дает о себе знать.

График изменения напряжения питания я решил не делать. Визуально во всех измерениях было 30.0 В. Чтобы убедиться, что это так, в интерактивном режиме Питона я ввел следующую команду:

>>> for u in open("uline").readlines() :
>>>   if not u == '30.0\n' :
>>>     print 'Not 30.0 V'

При выполнении Питон ничего не вывел, следовательно, во всех записях присутствует напряжение «30.0».

UPDATE 23.06.2014-19:08

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

Среднее значение, оно же — средне-арифметическое, оно же — математическое ожидание, оно же — «средняя температура по больнице», вычисляется совсем просто. Нужно просуммировать все элементы последовательности и полученное значение разделить на количество элементов.

Допустим, X — это последовательность. На языке Си — это массив значений. На языке Python — это список. Нам нужно найти сумму его членов. В Питоне это делается крайне просто:

M = sum(X) / len(X)

Однако, проблема в том, что в нашем случае X — это список целых значений. Его сумма — целое значение. И длина списка — тоже целое. При делении целого на целое получим целое, а это не совсем то, что на бы хотелось. Поэтому сумму преобразуем в число с плавающей точкой:

M = float(sum(X)) / len(X)

Для статистики ни доверительный интервал (разность между наибольшим и наименьшим значениями), ни среднее отклонение — неудобны. Вместо них используется понятие «дисперсия» и «стандартное отклонение». Стандартное отклонение численно равно квадратному корню из дисперсии. Дисперсия же определяется как сумма квадратов всех отклонений от среднего, поделенная на длину последовательности. Звучит страшно, но в Питоне выражается одной строкой:

V = sum([(x - M) ** 2 for x in X]) / len(X)

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

То же самое со статистическими расчетами. Если У вас есть среднее значение, то это полная фигня! Недаром появилась поговорка: «средняя температура по палате». А уж как правительство любит разводить свой народ на тему «средняя зарплата» по стране — это просто классика!

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

Я не буду приводит всю прогу целиком, напишу только одну функцию, которая на основе входных данных (списка) вычисляет среднее и дисперсию.

def get_MV(X) :
  M = float(sum(X)) / len(X)
  V = sum([(x - M) ** 2 for x in X]) / len(X)
  return (M, V)

При вызове функции в качестве аргумента указываем список значений. Функция возвращает кортеж из двух чисел.

Так вот, получасовые замеры гамма-фона дают следующие значения:

длина последовательности = 1851
min значение = 11
max значение= 48
среднее = 27.285251,
стандартное отклонение = 5.405611

12-часовые измерения:

длина последовательности = 43062
min значение = 8
max значение= 56
среднее = 27.982676
стандартное отклонение = 5.503473

Полуторачасовые измерения утром в Понедельник:

длина последовательности = 4897
min значение = 12
max значение= 48
среднее = 27.400245
стандартное отклонение = 5.390860

Как видно, повторяемость результатов очень хорошая.

Advertisements

2 responses to “Обработка данных гамма-каротажа

  1. Спасибо за очередную статью!

    Почему бы сразу не формировать строку в читаемом виде сразу в микроконтроллере и не отправлять её по RS232?

    Стоит добавить, что gnuplot умеет сам читать числа из файла, если они разделены пробелами. Например, если хотим построить второй столбец из файла и столбцы разделены пробелами, то пользуемся командой:

    plot «datafile.dat» using 2 with lines

    а для шестого столбца:
    plot «datafile.dat» using 6 with lines

    И вручную файл командой cut можно не делить.

    Для построения гистораммы и нахождения минимумов и максимов и СКО и т.п. обработки данных лучше использовать систему Octave. У неё Matlab-совместимый синтаксис. Это open-source клон Matlab.

    Она может строить гистограмму в одну команду:
    hist(data,Nintervals)
    Для СКО, среднего и сглаживания также есть готовые функции.

    Python также по-моему может всё это делать в одну функцию через SciPy, NumPy.

    • А Вам спасибо за Ваш комментарий!

      >> Почему бы сразу не формировать строку в читаемом виде сразу в микроконтроллере и не отправлять её по RS232?
      А я, собственно, так и делаю. Только тут есть одно обстоятельство.

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

      Во вторых, МЭК формирует большие пакеты, по 800 с лишним байт. Внутри пакета находится так же и информация от МГК, которая имеет объем всего два байта. Пакеты бинарные, не текстовые!

      Регистратор каротажных сигналов (РКС) принимает эти пакеты, проверяет их на валидность, выполняет еще кое-какие функции и — ВНИМАНИЕ! — в текстовом виде отправляет информацию в комп по USB. То есть из двоичного формата пакеты преобразуются в текстовый. При этом объем в байтах у пакетов возрастает примерно до 1200-1500 байт.

      Но мне для тестирования МГК вся эта информация не нужна. Поэтому я ее на приеме сразу откидываю. Беру только данные по гамма-каротажу, номер пакета и напряжение (см. текст проги).

      Так вот, проблема заключается в том, что перепрограммировать РКС ради тестирования МГК — это как-то не правильно. Правильнее написать прогу на стороне компа, которая принимает пакеты полностью, но на выход пишет только нужные данные, остальное (99%) просто откидывает.

      Надеюсь, вроде понятно объяснил.

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

      Тем не менее, делать, как Вы предлагаете, тоже можно. Главное — получить результат. И каждый сделает так, как ему удобнее и как он знает.

      И еще. Мне придется признаться — я не знаю ни Octave, ни MatLab, не умею ими пользоваться.

      Про SciPy и NumPy — я знаю, но, к сожалению, не имею практики использования.

      Времени просто катастрофически не хватает. А тут еще события на Украине…

Добавить комментарий

Заполните поля или щелкните по значку, чтобы оставить свой комментарий:

Логотип WordPress.com

Для комментария используется ваша учётная запись WordPress.com. Выход / Изменить )

Фотография Twitter

Для комментария используется ваша учётная запись Twitter. Выход / Изменить )

Фотография Facebook

Для комментария используется ваша учётная запись Facebook. Выход / Изменить )

Google+ photo

Для комментария используется ваша учётная запись Google+. Выход / Изменить )

Connecting to %s