Материал: Симонов Томографические измерителные информационные системы 2011

−∞ −∞

∞

где Zα (r ) = ( 1 )3 ∫ ν 2 J0 ( ν2 (r )d ν 4 ) , r = (х2 + у2 )12 , 2π 1+ α(δ) ν 1+ ν

0

ν= (ν12 + ν22 )12 , J0 (u) – функция Бесселя нулевого порядка. Далее необходимо выбрать параметр регуляризации α(δ) так,

чтобы μα(δ) → μ0 при δ → 0 .

Рассмотрим некоторые моменты численной реализации регулирующего алгоритма, например, с выбором α по невязке. Заметим, что выражение (2.130) представляет собой свертку функции Zα (r )

с функцией Sδ (x, y) . Поэтому при проведении вычислений можно воспользоваться теоремой о свертке, т. е. вместо (2.130) использо-

образ регуляризированного решения μα (х, у) .

При численной реализации (2.131) основными операциями являются дискретное прямое и обратное преобразования Фурье.

Остановимся на вопросе вычисления функции Sδ (x, y) , имеющей вид

Поскольку измерения функции Pδ (l, θ) проводятся для конечного числа углов θ [0, π] и в узлах lK = −R0 + (K −1) l ( l = 2R0 (N −1), K =1,2,..., N ), то R0 – радиус исследуемого объек-

та для параллельной схемы сканирования. При замене (2.132) какойлибо квадратурной формулой надо осуществлять интерполяцию

функции Pδ (l, θ) . Обычно используют метод интерполяции по бли-

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

От этих недостатков свободны интерполяционные сплайны n-го порядка, т. е. функции θδ(n) (x, y; θ) , имеющие непрерывные про-

изводные по х, у до n-го порядка и совпадающие с Pδ (lK , θ) для каждого фиксированного θ [0, π] в узлах lK .

Аппроксимируя интеграл (2.132) какой-либо квадратурной формулой и заменяя функцию Pδ (l, θ) интерполяционным сплайном

Qδ(n) (x, y; θ) n-го порядка, получим формулу для приближенного нахождения функции Sδ (x, y) :

i=1

где a j – коэффициент квадратурной формулы.

Можно ли построить регуляризирующие алгоритмы для уравнения (2.129), в которых не используются итерационные методы? На этот вопрос имеется положительный ответ [11]. Такие алгоритмы можно построить, основываясь на дополнительной информации об искомом решении (2.129).

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

значения μ0 (x, y) имеет вид

172

где Ci – круг радиуса Ri , расположенный внутри кругов C1 C0 с радиусами R1 и R0 соответственно. В примере были взяты следующие параметры модели: R0 = 120 мм, R1 = 110 мм, C0 = 1000 Hv (Hv – величины Хаунсфилда), C1 = 0 Hv и x0 = x1 = 0 , y0 = y1 = 0 ,

173

На рис. 2.24 показана модель исследуемого объекта и график функции μ0 (0, у) .

Рис. 2.24. Модель исследуемого объекта (а) и график функции μ0 (0, у) (б)

Под «тонкой структурой» функции μ0 (х, у) будем понимать функцию μ0 (х, у) с областью определения, включающей круги Ci

при i = 2, 3, 6, 7, 8, 9, 10, 11.

Для указанной модели μ0 (х, у) исследуемого объекта вычисля-

лись интегралы Р(l, θ) = ∫ μ0 (x, y)dl , где L(l, θ) – прямая ли-

L(l, θ)

ния с нормальными координатами (l, θ) из семейства прямых, инвариантного относительно вращения вокруг центра круга C0 , в ко-

торые вводилась случайная ошибка по формуле

Р(l, θ) = P(l, θ) + δ X P(l, θ)C , θ (0,π),

где X – равномерное распределение чисел в диапазоне [0, 1]. Уровень погрешности δ брался равным 0,005, что соответствует

средним значениям компьютерного томографа РКТ-01.

174

Для построения приближенного решения μ(х, у) основного уравнения КТ

использовался регуляризирующий оператор Rgδ вида (2.124′). При

использовалась дополнительной информации об искомом решении

μ0 (х, у) .

Из табл. 2.1 видно, что

μ0 (х, у) = μ(х, у) + μ1 (x, y) ,

где

Вследствие непрерывности оператора T в (2.134) имеем

P(l, θ) = P (l, θ) + P1 (l, θ) θ [0, π],

функция P(l, θ) известна, функция P1 (l, θ) мала по сравнению с функцией P (l, θ) .

Это типично, например, для томографии головного мозга человека, когда функция μ(x, y) описывает распределение линейного коэффициента ослабления рентгеновского излучения в костных тканях человека, функция μ1 (x, y) – распределение этого коэффициента в веществе мозга. Так как костные ткани головы имеют постоянный коэффициент ослабления, то функция μ(x, y) всегда известна и много больше функции μ1 (x, y) .

Это позволяет приближенное решение искать в виде