CC BY
42
\ УПРАВЛЕНИЕ В МЕДИЦИНЕ И БИОЛОГИИ
УДК 519.6; 621.391
АЛГОРИТМЫ ОБРАБОТКИ СИГНАЛОВ В МЕДИЦИНСКОЙ ДИАГНОСТИКЕ С ИСПОЛЬЗОВАНИЕМ ОПОРНОГО СЛУЧАЙНОГО ПРОЦЕССА
Э. П. Тихонов,
канд. техн. наук, доцент Санкт-Петербургский государственный электротехнический университет (ГЭТУ)
В настоящее время для решения сложных задач медицинской диагностики, в том числе редукционной проблемы Рэлея, например в цифровой рентгеноскопии и томографии, широко используются математические методы решения обратных задач, реализуемые программно на современных компьютерах. В качестве альтернативы в работе предлагается метод измерения с опорным случайным сигналом, который позволяет решать аналогичные задачи на аппаратно-программном уровне с достижением, по мнению автора,лучших метрологических характеристик.
Mathematical methods of the decision reverse problems are applied to the decision of problems of medical diagnostics. Alternative technologies in work are offered the method of measurement with bearing casual signal. The method examined in work allows to solve similar problems at a achievement the best metrological characteristics
Исходные сведения и постановка задачи
Основу современных средств медицинской диагностики составляют бесконтактные методы измерения, использующие эффекты взаимодействия зондирующих физических полей (сигналов) с исследуемым объектом. К зондирующим физическим полям относятся, например, ультразвуковые, световые, электромагнитные, рентгеновские, радиоизо-топные поля [1]. В зависимости от вида зондирующего сигнала (ЗС) изменяются закономерности его взаимодействия с исследуемым объектом и соответственно методы, алгоритмы и устройства диагностики, но сохраняется общее положительное свойство —получение информации бесконтактным способом. Некоторые трудности возникают при выявлении степени изменения ЗС в зависимости от патологического состояния диагностируемых объектов, в том числе внутренних органов тела человека. Указанные трудности вызваны следующими причинами:
— нелинейным характером взаимодействия ЗС с исследуемым объектом;
— низкой чувствительностью ЗС к изменению свойств и состояния диагностируемого объекта;
— особенностью преобразования в датчиках энергии ЗС в электрическую энергию;
— низким уровнем отношения мощности сигнала к мощности помехи.
Особую роль в медицине играют бесконтактные методы диагностики на основе использования в качестве ЗС рентгеновского излучения [1,2], причем прогрессирующее положение занимают цифровая рентгенография и томография [2, 3].
В качестве параметра исследуемого локализованного в пространстве и во времени объекта или его части в цифровой рентгенографии и томографии рассматривается изменение плотности в дискретных подобластях Цу объекта, на которые условно разбивается тем или иным способом наблюдаемая область или весь объект V исследования; при этом выполняется требование
л т
^ = Е 2Х-
(=1 7=1
Таким образом, как бы осуществляется объемная дискретизация исследуемого объекта и тем самым дискретизация его физического параметра-плотности, причем информация, которую несет рентгеновский ЗС об измеряемой плотности, проходя через подобласть Ц,, усредняется по объему с торцами, ограниченными площадями Бу.. и высотой (толщиной) djj. Центры площадей Бу имеют
координаты х, и уу-. В результате взаимодействия рентгеновского излучения с биологическим веществом в указанных объемах происходит ослабление исходной интенсивности ЗС в соответствии с законом Вера.
В современных цифровых рентгеновских установках, в отличие от томографических, зондирующий луч формируется специальным коллиматором и попадает на линейку, состоящую из конечного числа детекторов (рис. 1), которые перемещаются одновременно с рентгеновской трубкой с постоянной скоростью только в одной плоскости в направлении, например, оси у [2, 3]. В этом случае направление рентгеновского луча совпадает с осью г и, согласно упомянутому закону, наблюдаемая интенсивность выходного сигнала от каждой подобласти Цу изменяется в соответствии с законом [1, 2]
- |р{г)с1г
1(х,у) = 10(х,у)е
где /0(х, у), /(х, у) — функции, описывающие интенсивности рентгеновского луча, проходящего через центры Бу.. на входе и выходе исследуемой подобласти; р(г) — усредненная по площади Эу искомая функция плотности.
В томографах помимо поступательного перемещения с постоянной скоростью источники и приемники (детекторы) рентгеновского ЗС поворачиваются дискретно на угол ср. Причем сканирование подобным образом выполняется полностью для всего исследуемого объекта.
Функция плотности р(г) зависит от исследуемого объекта и содержит информацию о свойствах объекта, выраженных в равенстве
р(г) = 14г)рП,
где |!(7) — массовый коэффициент ослабления рентгеновского излучения; р(г) — плотность вещества, подвергающегося рентгеновскому облучению.
В медицинской диагностике наибольшую плотность имеют костные ткани, массовый коэффици-
Мас II
■ Рис. 1. Схема цифровой рентгеновской установки
ент ослабления которых в 68 раз превышает аналогичный параметр мягких тканей, поэтому в костных тканях рентгеновский луч имеет наибольшее поглощение [1]. Как будет показано ниже, этот факт имеет решающее значение для предлагаемого и исследуемого в работе алгоритма измерения плотности.
В цифровых рентгеновских установках и томографах происходит накопление информации за время экспозиции Т об измеряемой плотности по интенсивности /(х, у), которая зависит также от четвертой координаты — времени, т. е.
г
Я(х, у, Т) - |/(х, у, 0с/£,
о
где Я(х, у, Т) — результат накопления информации на выходе исследуемой подобласти за время Г.
В реальных цифровых рентгеновских установках накопление информации осуществляется в темпе с преобразованием энергии рентгеновских лучей в электрические сигналы. В зависимости от типа детектора электрический сигнал принимает тот или иной вид. Например, при использовании в качестве детекторов пропорциональных многопроволочных камер рентгеновские кванты преобразуются за счет ударной ионизации в электрические импульсы, которые после нормализации поступают в каждом канале на вход накапливающих двоичных счетчиков. На счетчиках осуществляется преобразование за время экспозиции Т средней частоты следования случайного потока импульсов, пропорционального интенсивности рентгеновского излучения и имеющего пуассоновское распределение вероятностей, в цифровой двоичный код. Далее этот код поступает на ЭВМ, где осуществляются дальнейшая обработка полученной информации и отображение результатов обработки на дисплее или других носителях информации. При этом интеграл для определения Я(х, у, Т) преобразуется в сумму, а результаты измерения дискретизируются (квантуются) не только в пространстве, а и во времени с интервалом дискретизации, равным Т.
Если же в качестве датчиков используются фотодиоды со слоем люминофора для преобразования рентгеновского излучения в электрический сигнал, то возможны и другие формы представления электрического сигнала, например, в виде реализации непрерывного случайного процесса. Количественные условия формирования реализации случайного процесса из импульсного исходного сигнала описаны в работе [4], а информация об интенсивности рентгеновского излучения содержится в его математическом ожидании. В этом случае целесообразно исходный случайный процесс представить в виде модели случайного квазистацио-нарного процесса, для которого принимается гипотеза о постоянстве математического ожидания и других его вероятностных характеристик в течение времени экспозиции.
Если же за время экспозиции принятая гипотеза не выполняется, то следствием этого является возникновение дополнительной, т. е. динамической
составляющей погрешности измерения интенсивности рентгеновского излучения. Этот недостаток присущ всем подобным «накопительным» статистическим алгоритмам, усредняющим исходную информацию за фиксированное время Т. Таким образом, фиксация времени Т, как правило, приводит к возникновению динамической составляющей погрешности при изменении плотности в течение времени экспозиции Т.
Для рентгеновской томографии удобно пользоваться компактной формулой для непрерывного представления функции плотности в виде [1, 2]
- |р(*. у)с&'
/(/, ф) = /0(/, Ф) е 1-(,*ф) , (1)
где / — координата детектора; ср — угол поворота источника и детектора; /(/, ф), /0(/, ф) — интенсивности рентгеновский ЗС на входе и выходе исследуемого объекта; р(х, у) — плотность исследуемой подобласти объекта, находящейся в поле зрения луча ЗС, с координатами х и у; при этом луч ЗС распространяется по прямой линии вида
хсоб(ф) + у э1п(ф) = 1.
При записи уравнения (1) в форме
0(/, Ф)= |р(х,у)с/з (2)
Щ. ф)
получают так называемое преобразование, или уравнение Радона, где
е(/, ф)= — !п[/(/, ф)//0(/, ф)].
Решение данного уравнения позволяет найти искомую функцию поглощения, непосредственно связанную простой функциональной зависимостью с плотностью исследуемого объекта. Учитывая связь функции поглощения с функцией плотности, далее по тексту функция 0(/, ф) будет называться либо функцией поглощения, либо функцией плотности.
Трудности и методы решения уравнения Радона хорошо известны [2,3], так как оно относится к классу так называемых некорректно поставленных задач, при которых погрешности измерения оказывают существенное влияние на результаты решения. Несмотря на успешное применение известных математических методов решения некорректных задач, вопросы, связанные с повышением достоверности получения исходной информации, входящей в правую часть уравнения, по результатам измерения остаются актуальными. Это связано с тем, что качество решения данного уравнения Радона в значительной степени зависит от погрешности представления левой части уравнения (2), которую можно получить только в результате измерения. Попытка непосредственного уменьшения погрешности измерения путем программного решения обратной задачи уже применительно к результатам измерения, связанным математически с измеряемым параметром посредством уравнения Фредгольма I рода, приводит к определенным
положительным результатам [1]. Однако представляется, что возможно эффективное решение проблемы Рэлея, но уже с привлечением аппаратно-программных средств.
В данной статье на основе метода измерения со случайным опорным сигналом, предложенного и разработанного в работах [5,6]с использованием подхода, рассмотренного в работах [7, 8] для бесконтактного измерения паросодержания, предлагается итерационный алгоритм измерения функции поглощения 0(/, ф) и исследуются его свойства. Как показала практика, применение алгоритма с опорным сигналом уменьшает не только инструментальные погрешности, но и динамическую составляющую погрешности измерения, и непосредственно в процессе измерения линеаризирует левую часть уравнения (2).
Итерационный алгоритм с опорным сигналом
Суть предлагаемого алгоритма состоит в следующем. Используя информацию о ЗС на выходе исследуемого объекта при максимальном (или минимальном) поглощении рентгеновского излучения, в измерительный алгоритм дополнительно вводится случайный опорный сигнал, представляющий собой пуассоновский процесс [7, 8]. Благодаря введению такого случайного опорного процесса удается связать интенсивность рентгеновского излучения на выходе исследуемого объекта с соответствующим цифровым эквивалентом на выходе измерительного устройства посредствам линейной градуировочной характеристики. Преимущества подобного метода в отличие от чисто математического способа линеаризации заключается в следующем. Опорный сигнал может формироваться в дополнительном канале, у которого максимальная или минимальная плотность нормированы путем соответствующего поглощения рентгеновского (или другого, например, ультразвукового) излучения введенными в измерительный канал фантомами. Причем параметры фантома могут подстраиваться по параметрам исследуемого объекта в процессе калибровки. Поскольку данный опорный канал непосредственно входит в измерительную линейку, то он позволяет исключить или, по крайней мере, значительно ослабить целый ряд влияющих на результат измерения отрицательных факторов и осуществить линеаризацию градуировочной характеристики, связывающей результат измерения с измеряемым параметром — искомой плотностью [7].
В математической форме данный алгоритм для каждого канала при фиксированном векторе 2 = = (/, ф) имеет вид
/с/[(/7 + 1)Д*] = к,{пМ)~
- у,- (Л/(пМ) ехр [-0/ (г) - а к, (лДГ)] --Л/о(лЛОехр[-р0/тах +ак1{пМ)]}, (3)
где к,(пМ) и ki[(n+^)At] — цифровой эквивалент искомой функции плотности на п-м и п + 1-м так-
тах итерации для /'-го канала (результат измерения искомой плотности); Af—временной интервал, через который опрашивается выход /'-го датчика в/'-м канале; г), — случайная пуассоновская величина, соответствующая появлению электрического импульса в момент опроса /'-го канала; г\ю — случайная пуассоновская величина, соответствующая появлению электрического импульса в опорном канале, совпадающем с моментом опроса /'-го канала; у/ — заданная постоянная величина или постоянная уравновешивания для /'-го канала; 0,(z) — искомая функция плотности; 0/тах —максимальное значение плотности, которое может принимать функция плотности при данном рентгеновском обследовании для /'-го канала, сформированное на фантоме; а и (3 — заданные постоянные.
Для рассматриваемого алгоритма дополнительно выполняются следующие условия и ограничения:
• /7=1, 2..А/, гдеN — число итераций;/'= 1, 2,
т, гдет — число каналов;
• 0 < 0/(z) < 0/ max, 0 < к/(nAt) < 0/ max, Y/ = 0/ max/^0> где N0 — заданное положительное число, соответствующее числу квантования максимального значения измеряемой плотности;
• начальное значение/c,(nAf) равно нулю при числе итераций п, равном нулю, а для каждой новой позиции линейки датчиков сохраняется предыдущее значение к/{ ).
Отметим, что постоянную у, целесообразно выбирать у = const независимо от /'-го канала из условия
max 0/max _ i_________
Условия функционирования рассматриваемого алгоритма можно изменить, если связать момент изменения искомого параметра к{пМ) с моментом появления импульса на выходе детектора /'-го канала с одновременной проверкой наличия импульса в опорном канале.
При исследовании подобных итерационных алгоритмов необходимо проверить следующие три условия, определяющие возможность использования данных алгоритмов для измерения искомых параметров. Эти условия определяются выполнением следующих требований:
• сходимостью в установленном смысле измеряемого параметра к{пМ) к его установившемуся значению;
• числом тактов, определяющих время сходимости или переходный процесс от начального к установившемуся (измеренному) значению параметра;
• конечностью финальной дисперсии или дисперсией, характеризующей флуктуацию измеряемого параметра относительно установившегося значения и тем самым характеристику случайной составляющей погрешности измерения при числе тактов п ->
Если за время наблюдения выполняются условия стационарности случайного потока, то можно установить [5, 6], что искомый параметр в алгорит-
ме (2) при п > N (где N некоторое положительное число) и надлежащем выборе постоянной у стремится к такому значению кю, при котором выполняется в среднем равенство
{п/ (nAt) exp [- 0,- (z) - а /с,о (лДГ)]} =
= {n/O(nAf)exp[-(30max + aki0(nAt)]}, где мл {п/ {nAt) exp [- 0/ (z) - а /С/о {nAt)]}=
= Ij exp[— 0/(z) - akjo(nAt)] и мч {п/о (nAt) exP [- P emax + a ki0(лАГ)]} =
= // о exp[-(30/max + aki0(nAt)],
т. e. выполняется сходимость по распределению [9].
Если обеспечить идентичность интенсивностей /,• и /ю, то при п > N выполняется в среднем условие
«/oW=Pe,'ma*~e,'(Z>. (4)
2а
которое при соответствующем выборе постоянных а и (3 определяет искомую градуировочную характеристику измерительного устройства, функционирующего в соответствии с исследуемым алгоритмом. Например, для а = 0,5 и (3 = 1 получаем в среднем равенство
^io(z) = 0/max _0/(^)i
где z = nAt.
Алгоритм (3) можно представить в виде следующей модификации:
к, [(п +1) At]= к, {nAt) -
-y/{ri/(nAf)exp[- 0,-(z) + akj{nAt)]~
_ Л/о (nAt) exp [-a kj (лА£)]}. (5)
Для данного алгоритма при тех же условиях и а = 0,5 выполняется в среднем равенство
ki0(nAt) = Qi0(z). (6)
Попутно заметим, что, выбирая соответствующим образом постоянную а, можно достичь эффекта усиления результатов измерения и тем самым повысить чувствительность измерительного устройства, так как
— = (2 а) ~1,
de
и для а < 0,5 следует 8 > 1.
Анализ погрешности
Особенность алгоритмов (3) и (5) состоит в том, что по заданному значению у можно подобрать такое значение е(у) > 0, при котором для всех п > N выполняются условия:
к j (nAt)-
Р 0/max 0/(лДҐ)
2а
>є(у) =0;
Р{|/С;(ЛА0- 0/(пЛГ)| > е(у)}= О,
где Р{...} — вероятность выполнения соответствующих неравенств.
Это значит, что при постоянном конечном значении параметра у можно гарантировать сходимость по вероятности /<7(лД?) только в некоторую е(у)-ок-рестность искомого значения/с0/(лА?). Следовательно, е(у)-окрестность определяет абсолютное значение методической погрешности измерения искомого параметра посредством рассматриваемых алгоритмов, которая фактически эквивалентна погрешности оценки математического ожидания, возникающей из-за конечности объема выборки.
Величину е(у) можно оценить как для асимптотического случая, т. е. когда л стремится к°°, так и при конечном числе тактов л, если воспользоваться, например, неравенством Чебышева
к і (nAt)
2a
>E(y) <
D\к і (nAt) -
2a
є2(У)
(7)
Для того чтобы установить меру отклонения цифрового эквивалента или результата измерения /с,(лД£) при конечном л от искомого значения /со/(лД0, необходимо найти дисперсию, которая характеризует методическую случайную составляющую погрешности измерения искомого параметра при л < А/0.
Переходный процесс, характеризующий сходимость в среднем /с,(лд?) к /с0/(лДГ), и его дисперсия для алгоритма (5) в первом приближении определяются по следующим формулам;
кі [лД£] = (1 - vA,,-)n [ki(0)-k0i] + k0j]
D{kj (nAt)} ~ D-, [1-(1-vX/)2n],
(8)
(9)
где X; - Yi If, kj(nAt) — значение измеряемого параметра на n-м такте итерации; k0j— искомое установившееся значение измеряемого параметра; v = At exp (-6(z)); D- — финальная дисперсия или ее предельное значение,
D; =
2aexp(-6/(z)) - ya (1 + 2 I, At)
(Ю)
Подставляя значение у и пренебрегая вторым слагаемым в знаменателе, получаем следующую оценку финальной дисперсии:
(11;
max 0/ exp max 0,- (z)
D* < —--------—----------
2N§a
Пусть абсолютное значение погрешности е(у) составляет долю v < 1 от max 0/, тогда формула (11) приобретает вид 1
D, <
exp I max 0/ (z)
2/S/oav
В соответствии с неравенством (7) оценка вероятности непревышения абсолютной погрешностью Д(у) измерения искомой плотности заданной величины е(у) равна
Р{|Д(у)|< £(у)} =•
exp max 0/ (z)
2Л/оау
Полученное значение вероятности характеризует фактически доверительный интервал искомой методической погрешности, надежность которого определяется значением указанной вероятности [9]. Эта вероятность обычно задается априори, и по ее значению и по известному максимуму измеряемого параметра тахв/ из выражения (12) можно найти параметр Л/0.
Формулы (8) — (11) показывают, что при надлежащем выборе параметра Л/0, т. е. таком его значении, при котором выполняется условиеX < 1 для всех /', условие сходимости текущего значения измеряемого параметра /с(лД?) к его установившемуся значению выполняется. При этом сходимость алгоритма обеспечивается при конечной финальной дисперсии О,-, значение которой оценивается в соответствии с формулами (10) и (11) и полностью контролируется параметром/\/0. Заметим, что полученные формулы для оценки параметра А/0 дают завышенный результат из-за использования при выводе этих формул неравенства Чебышева.
Анализ динамического режима измерения
Преимуществом данного итерационного алгоритма является также тот факт, что его работоспособность не только сохраняется при изменении за время экспозиции искомой функции плотности, но и позволяет повысить точность измерения по сравнению с известным алгоритмом обычного усреднения. Этот факт является важным хотя бы потому, что диагностируется живой объект, полную неподвижность которого за время экспозиции обеспечить невозможно.
Действительно, пусть искомая функция плотности за наблюдаемый отрезок времени изменяется в первом приближении линейно
0(лд 0 « 0(О)+0'(О)лД?,
тогда для алгоритма (1) при (3 = 1 для искомого значения параметра получим
к[(п + 1)Дґ]= (1-vy/
МО)-
-0(0)
2а
. в max - 8(0) - вЩШ - Г3 [1 - (1 - уу)"]}
2а
Из последнего выражения следует, что при достаточно большом числе тактов л искомый параметр для линейно изменяющейся функции плотности запаздывает относительно истинного значения на время Т3 - Af/v у, которое не влияет на результат решения уравнения Радона.
Положительная особенность предлагаемого алгоритма проявляется прежде всего в высокой помехоустойчивости относительно аппаратурных искажений, обусловленных рассеянием рентгеновского луча, взаимным влиянием детекторов и нестабильностью коэффициента регистрации е излучения детектором. Погрешность, возникающая из-за аппаратурных искажений, можно представить в виде дополнительного к основному излучению для данного канала фонового излучения /ф(/, ф). Это фоновое излучение при обычном накоплении (интегрировании) приводит к возникновению дополнительной составляющей погрешности, так как алгоритм (1) преобразуется к виду
к, [(п +1) At] = к, {nAt) — у {в-| [r\j {nAt) exp [- 0/ (z)] + +г|ф {nAt)] exp [-akj {nAt)] --e2[ti/o (пАЦ + цф{пА^] exp [-(30/max + akj (nAt)]},
где e-, и e2 — коэффициенты регистрации детекторами в /-м и опорном каналах; г)ф(nAt) — фоновое излучение.
Нетрудно показать [8], что максимальная относительная дополнительная погрешность, возникающая из-за неидентичности каналов, определяется по формуле
= б.
1-U
+ ^Ф Л/о
где 5е — относительная вариация коэффициента регистрации, 8е = (е2 - е-,) /е-,.
Так как 1 /0, тах« 1 и Лф« Л/о. то 8тах = 5е, и при хорошей технологии выполнения детекторов эта погрешность может достигать очень малых значений. Заметим, что данный способ исключения дополнительных аппаратурных погрешностей при томографии более эффективен, чем рассматриваемый, например, в работе [2] численный метод решения проблемы Рэлея, которая предшествует основному решению уравнения Фредгольма I рода.
Результаты моделирования
Анализируемый подход был использован в реальных условиях для бесконтактного измерения ра-диоизотопным методом плотности двухфазной среды [7, 8], который в целом хорошо себя зарекомендовал. Техническое решение радиоизотопного измерения с использованием опорного сигнала было предложено, например, в работах [10, 11] и ре-
ализовано в измерительно-вычислительных комплексах, которые успешно эксплуатировались в Сибирском энергетическом институте СО АН (г. Иркутск) и в филиале ВНИИАЭС — ЭНИС (г. Электро-горск).
Однако для наглядного представления процесса измерения применительно к рассматриваемому случаю, с учетом того, что теоретические результаты по исследованию переходных процессов получены при определенных допущениях и ограничениях, целесообразно провести более углубленное исследование на основе имитационного моделирования.
В соответствии с программой, разработанной в математическом пакете «Статистика» (см. Приложение), на рис. 2, 3 приведены результаты моделирования алгоритма для измерения значений функции плотности И. Значения Я соответственно равны 0,05; 0,10; 0,15; 0,20; 0,25 и за время измерения не меняются. В программе установлены следующие значения параметров: у = 0,001; объем выборки или число итераций N - 8000; интенсивность / = 0,01; искомое значение функции плотности изменяется от 0,05 до 0,25 с шагом дискретизации 0,05. На рис. 2 приведены финальные значения указанных на рис. 1 параметров, т. е. без учета данных о переходном процессе, результатов измерения.
Графики, приведенные на рис. 2, имеют следующие значения:
Значения Я 0,05 0,1 0,15 0,2 0,25
MEAN 0,05422 0,10109 0,15079 0,19882 0,23521
(среднее) MEDIAN 0,05400 0,10200 0,15100 0,20000 0,23200
(медиана) SD (с. к. о.) 0,00665 0,00560 0,00589 0,00951 0,00745
Для того чтобы оценить полученные теоретически формулы для оценки свойств рассматриваемых алгоритмов, на рис. 4-8 приведены результаты сравнения аналитических расчетов и имитационного моделирования алгоритма (3) в соответствии с программой,приведенной в Приложении.
Как следует из рис. 2-8, полученные в результате моделирования значения функции плотности практически совпадают с их априорно заданными
0,30 0,25 0,20 0,15 0,10 0,05 0,00
0 800 1600 2400 3200 4000 4800 5600 6400 7200 8000
■ Рис. 2. Результаты моделирования алгоритма для измерения значений функции плотности
а)
0,26
0,24
0,22
0,20
0,18
0,16
0,14
0,12
0,10
0,08
0,06
0,04
0,02
4
Кг
^0 = 0,25
К
е = 0,2
1 1
^\0 = 0,15
1 V
'0 = 0,1
i ►_ _____ (
^0 = 0,05
Среднее значение
Медиана
б) С. к. о.
0,0155
0,0145
0,0135
0,0125
0,0115
0,0105
0,0095
0,0085
0,0075
0,0065
0,0055
0,0045
0,0035
0,0025
0 = 0,25
0 = 0,2
= 0,15 У
0,1
0 = 0,05
С. к. о.
Рис. 3. Средние значения и медианы измеряемых параметров, полученные по результатам моделирования без учета переходных процессов (а) и финальные среднеквадратические отклонения (с. к. о.),характеризующие случайную составляющую погрешности измерения искомой функции плотности для этих же измеряемых параметров и при тех же условиях (б)
■ Рис. 4. Гоафики, иллюстрирующие результаты расчета (1) и машинного эксперимента (2) по определению с. к. о. случайной составляющей погрешности измерения функции плотности (искомая плотность равна 0,05)
О 1000 2000 3000 4000 5000 6000
■ Рис. 5. Расчетные (1), (2) и экспериментальные (3) и (4) графики для различных интенсивностей (для 1, 2, 3, 4 они соответственно равны 0,01, 0,02, 0,02, 0,03) случайных пуассоновских импульсов, иллюстрирующие с учетом переходного процесса изменение с. к. о. измеряемой плотности, равной 0,05.
Ш Рис. 6. Графики, иллюстрирующие результаты расчета среднего значения (1) измеряемого параметра О, равного 0,05, и машинного эксперимента (2) с учетом переходного процесса, а также значения 25 % (3) и 75% (4) квантилей, определяющие случайный разброс измеряемого параметра, возникающего из-за случайного характера зондирующего рентгеновского сигнала
■ Рис. 7. Графики, иллюстрирующие с учетом переходного процесса усредненные данные машинного эксперимента для одного и того же значения измеряемого параметра, равного
0,05 при различных интенсивностях I случайных пуассоновских импульсов ЗС:
------ / = 0,02;---- / = 0,03;. - / = 0,04
■ Рис. 8. Графики, полученные при различных постоянных усреднениях (у) с учетом переходного процесса, результатов измерения изменяющегося по синусоидальному закону искомого параметра (а — шаг равен 0,03; б, в — шаг равен 0,04; на рис. в значение измеряемого параметра увеличено на 0,6 по сравнению с его значениями, указанными на рис. а и б):
------входной синусоидальный сигнал плотности;
----— измеренный синусоидальный сигнал плотности
значениями, что подтверждает работоспособность предложенных алгоритмов для измерения функции плотности с использованием ЗС в целях медицинской диагностики.
Заключение
В отличие от известных математических методов решения обратных задач, например, для повышения точности результатов измерения, в предлагаемом подходе фактически решается обратная
задача на аппаратно-программном уровне. При этом решение задачи осуществляется путем аппаратного выполнения функции сравнения результатов преобразования как образцовой, так и измеряемой физической величины идентичными операторами.
Действительно, в момент уравновешивания выполняется равенство
А{-х - к) = А*{к-хтах),
или
А[-х + к) = А* (-к),
где А и А — идентичные операторы преобразования в измерительном и опорном каналах соответственно; х — значение измеряемого параметра из интервала (0,хтах); к — уравновешиваемая физическая величина; хтах — максимальное значение образцовой меры, идентичность которой максимальному значению измеряемой величины устанавливается априори.
Из исходного уравнения сразу же следует решение
-х-к = А~1 ■А*{к-хтах).
Если операторы А~А*,то Л-1 - А ~ 1 и х= (хтах --к).
Имея программный доступ к регулировке операторов А и А*, достаточно добиться их идентификации в процессе проведения калибровки непосредственно перед проведением измерений в одной точке измеряемой зависимости (в точке максимума или минимума). Тем самым достигается практически полная идентификация операторов А и А* за счет устранения составляющих инструментальной погрешности, связанных с дрейфом и старением элементов блоков и узлов, в том числе источника ЗС, образующих измерительную цепь.
Для уменьшения случайной составляющей погрешности, возникающей из-за ведения дополнительного канала со случайным опорным сигналом, можно осуществить раздельно операцию накопления информации по опорному каналу в процессе калибровки непосредственно перед проведением операции измерения и реализовать, например, алгоритм (5) в виде
к/ [(п +1) ДГ ] = к, (лД£) -
- У/ {п/ (пМ) ехр [- 0/ {г) + а к, (лДГ) ] -
-Л/о ехр [-а/с,- (лД£)] },
где Лю — результат усреднения информации, поступающей по опорному каналу после калибровки.
Таким образом, в рассматриваемой статье предлагается аппаратно-программный метод решения проблемы Рэлея, возникающей при использовании различных зондирующих сигналов в бесконтактных измерениях, в том числе в цифровой рентгенографии и томографии.
Приложение
Основная программа для моделирования процесса измерения искомой плотности в соответствии с алгоритмом (2) имеет вид:
^oг ]:=1 Ю М с!о Data(1,j):=0;
Тог \:=4 ^ М do 1Ъг к=1 1о N do
Data(\+^,^)■=Data{\,])+D*({0.5*{^+S\gn{Rnd(L)-
-01*ехр(-р1+0.5*0а1а(и)))))-(0.5*(1+31дп(рпс1(1_1)-
-□1*ехр(-0.5*0а1а(у))))));
Программа для моделирования изменяющейся за время экспозиции функции плотности:
0:=0.001; 1М:=7000; М:=10; 1_:=10; 01:=9.9; В1 :=0.1; 1_1 :=10;
]:=1 Ш М с!о Оа1а(1 ,]):=0;
]:=2 1о М с1о йэг \■=^ Ш N с1о
Оа1аО+1,П:=Оа1а(и)+ОТ((0.5*(1+8|дп(Рпс1(1)--О1*ехр(-Р1*(0.5+8т(0.2*4*Р1*!/1\1))+0.5*Оа1а(и)))))--(0.5*(1+81дп(Рпс)(1_1 )-01 *ехр(-0.5*0а1а(и))))));
к=1 1о N с!о Оа1а(1+1,1 ):=Р1 *(0.5+зт(0.2*4*Р1*1/[\1));
В программах введены следующие параметры и постоянные: у = О = 0,001; число итераций N = 6000 или (в зависимости от рис. 8000); объем выборки М = 200 или (в зависимости от рис. 8000); L = 10 — параметр, характеризующий равномерное распределение вероятностей, используемое для имитации пуассоновского потока квантов рентгеновского излучения; D^ = 9,8 (или другие указанные на рисунках значения) — параметр, определяющий интенсивность потока квантов рентгеновского излучения; Ж =0,05 (или другие значения, указанные на рисунках) — значение измеряемой функции плотности.
Кроме того, помимо данной программы, использовались встроенные в математический пакет «СТАТИСТИКА» операторы и программы для вычисления соответствующих вероятностных характеристик результатов эксперимента.
Литература
1. Ремизов А. Н. Медицинская и биологическая физика: Учебник для вузов. — М.: Высш. шк., 1999. — 616 с.
2. Сизиков В. С. Математические методы обработки результатов измерений: Учебник для вузов. — СПб: Политехника, 2001. — 240 с.
3. Троицкий И. Н. Статистическая теория томографии. — М.: Радио и связь, 1989. — 240 с.
4. Рытов С. М. Введение в статистическую радиофизику. Ч. 1. Случайные процессы. — М.: Наука, 1976.— 494 с.
5. Тихонов Э. П. Измерения с опорным случайным процессом // Метрология. — 1985. — №10. — С. 20-29.
6. Тихонов Э. П., Селиванова М. П. Основные принципы, области применения и перспективы развития измерений с опорным случайным процессом // Измерения, контроль, автоматизация. — 1990. — № 3. — С. 3-6.
7. Селиванова М. П., Тихонов Э. П. Сравнительный анализ и выбор оптимального алгоритма радиоизо-топного измерения плотности // Автометрия. — 1990. — №3. — С. 69-72.
8. Селиванова М. П., Тихонов Э. П. Метод измерения плотности двухфазных сред в динамике // Автоматика. — 1989. — № 5. — С. 79-83.
9. Ширяев А. Н. Вероятность. — М.: Наука, 1980. — 570 с.
Ю.А.с. ви 1384001 А1, СССР, МКИ й01 N 23/02 /
Радиоизотопное устройство измерения плотности в динамике/ Э. П. Т и х о н о в, А. Н. Г у с е в.
11.А.С. Эи 1517534 А1, СССР, МКИ 001 N 23/02 /
Радиоизотопное устройство измерения плотности в динамике / Э. П. Т и х о н о в, Г. Г. Ж и в и л о в.
