Теоретические основы вычислительной диагностики

14. Представление одномерно неоднородной среды с непрерывным изменением параметров, как предела слоистой среды

На очень высоких частотах все среды (диэлектрические и проводящие) имеют диэлектрическую и магнитную проницаемости, близкие к постоянным вакуума ₀, ₀ [6]. В соотношении для коэффициента прохождения слоистой структуры (11.4) частичные коэффициенты прохождения с ростом частоты стремятся к 1. Основной вклад в общий коэффициент прохождения определяется экспоненциальным множителем . Для оценки влияния волновые числа _j запишем в следующем виде:

. (14.1)

Из этого соотношения видно, что второе слагаемое дает конечный вклад, независящий от частоты.

Кроме того предел отношения синуса угла падения к синусу угла преломления приближается к единице:

, (14.2)

и угол преломления приближается к углу падения _j=₀. То есть излучение распространяется в плоскослоистой среде прямолинейно без преломления. Рассмотренные предельные соотношения справедливы для коротковолнового ультрафиолетового и рентгеновского излучений.

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

, (14.3)

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

- единичный вектор направления распространения луча.

В соответствии с соотношениями для составляющих поля плоской волны (5.9) (5.10) получаем для высоких частот:

, (14.4)

где - мнимая часть волнового числа среды,

- координата, отсчитываемая от некоторой начальной точки на луче (=0),

I₀=SE₀²/W₀ - интенсивность в начальной точке луча.

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

, (14.5)

где - скорость света в вакууме (310⁸ м/с).

В высокочастотном приближении, таким образом, интенсивность излучения не зависит от частоты поля, а затухание определяется проводимостью на отдельных участках траектории. Формула (14.5) естественно обобщается на общий случай неоднородных сред с произвольной пространственной зависимостью проводимости от координат:

, (14.6)

где f(z) = (z)/c, L - прямая распространения луча, по которой производится интегрирование.

15. Основное уравнение компьютерной томографии

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

Пусть задана функция двух координат . Интеграл этой функции по произвольной прямой определяется положением этой прямой. Опустим перпендикуляр на прямую на начала координат:

, ( 15.1)

где - единичный вектор, направленный вдоль перпендикуляра к прямой, а s - одномерная координата, определяющая положение прямой относительно начала координат. Значение интеграла по прямой функции от функции является функцией аргумента :

, (15.2)

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

. (15.3)

Решение этой задачи было получено Иоганом Радоном в 1917 году и само преобразование , называется преобразованием Радона. Зафиксируем положение вектора и будем рассматривать функцию скалярного аргумента s

. (15.4)

В задачах компьютерной томографии изучаемый объект обычно находится в ограниченной области и переменная s меняется в ограниченных пределах s_min < s < s_max, сама функция по физическому смыслу также ограничена. Поэтому существует интеграл:

, (15.5)

который является преобразованием Фурье от функции .

Сама искомая функция также ограничена и ее аргументы меняются в ограниченных пределах и поэтому существует и ее преобразование Фурье по переменным x и y:

, (15.6)

где ds=d_xd_y - элемент площади плоскости xoy,

- вектор пространственных частот.

Между спектральными функциями и существует однозначная связь. Для ее установления на плоскости пространственных частот _x, _y рассмотрим значения на прямой :

. (15.7)

Для интегрирования по плоскости S введем переменные и , входящие в соотношение (15.7), так что :

.

С учетом (15.2) и (15.5) получаем:

. (15.8)

Полученное соотношение является математической формулировкой теоремы о центральном сечении.

16. Решения уравнения компьютерной томографии

Мы теперь подготовлены к тому, чтобы найти формулу обращения. Функция может быть найдена с помощью обратного преобразования Фурье по известной спектральной функции :

, (16.1)

где ds=d_xd_y - элемент площади плоскости _x, _y.

Для того, чтобы воспользоваться соотношением (15.8) интегрирование в (16.1) будем проводить в полярных координатах:

,

так что _x=cos, _y=sin, ds=dd.

Полярные координаты на плоскости пространственных переменных введем следующим образом:

.

Тогда (16.1) приобретет вид:

. (16.2)

Полученная формула и является искомой формулой обращения.

Однако полученная форма решения неудовлетворительна по причине явной избыточности требуемой информации. Функция Радона симметрична:

(16.3)

и для 2 (-S) мы просто повторяем измерения. Преобразование Фурье также обладает свойством симметрии

R_-(-)= R(), (16.4)

в чем нетрудно убедиться непосредственно с помощью (15.6). Поэтому в (16.2) можно интегрирование по проводить в интервале (0, ):

. (16.5)

Формула (16.5) лежит в основе так называемого Фурье-алгоритма восстановления. Используя теорию обобщенных функций можно перейти от тройного интегрального преобразования к двойному. Для этого рассмотрим преобразование Гильберта:

, (16.6)

где несобственный интеграл понимается в смысле главного значения Коши. Преобразование Фурье от Hh(s) по переменной s дается формулой:

, (16.7)

где h() - преобразование Фурье от h(s).

Обращаясь к (16.5) усматриваем, что

, (16.8)

что означает преобразование Фурье от преобразования Гильберта первой производной функции Радона по переменной s. С учетом этого формулу (16.5) можно представить в виде:

. (16.9)

Во внутреннем интеграле сделаем замену переменной t=s+q, после чего его главное значение Коши может быть вычислено с помощью приема, сводящего эту операцию к обычному интегралу:

.

Используя эту формулу и меняя порядок интегрирования, получим:

.

Интеграл от второго слагаемого с учетом соотношения симметрии функции Радона (16.3) преобразуется следующим образом:

.

Таким образом

, (16.10)

где - усредненная функция Радона со смещением.

Вновь используя свойство симметрии (16.3) можно доказать четность по аргументу q:

=

и, соответственно, нечетность первой производной . После чего получаем окончательную формулу в том виде, в котором она впервые была подучена самим Радоном:

. (16.11)

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

, (16.12)

ставящее в соответствие данной точке интеграл от функции f вдоль прямой, перпендикулярной и проходящей через эту точку. Это преобразование называется обратной проекцией f в точку вдоль прямой ортогональной вектору , или двойственным оператором. Интеграл от (16.12) по всем возможным направлениям, усредненный по углу называется суммарной обратной проекции f в точке . За счет свойства четности функции интегрирование можно осуществлять в интервале .

, (16.13)

где та же функция, что и в (16.11).

Аналогично определим оператор, называемый потенциалом Рисса I, действующий на функцию по переменной s. Он задается своим преобразованием Фурье по закону

, (16.14)

где - преобразование Фурье от . Нетрудно видеть, что

. (16.15)

Потенциал Рисса может быть введен и для исходной функции f:

.

Интегрирование в полярных координатах дает:

.

По теореме о центральном сечении (проекционной теореме) преобразование Фурье искомой функции можно заменить на :

.

Пользуясь свойством четности переходим к интегрированию в интервале (-):

.

Внутренний интеграл есть потенциал Рисса от в степени -1:

,

а последняя формула в свою очередь определяет усредненную обратную проекцию функции :

.

Применяя, наконец, к обеим частям оператор I^- получаем окончательную общую формулу обращения преобразования Радона для двумерной задачи:

. (16.16)

Рассмотрим частные случаи полученной формулы. Для =0 получим:

, (16.17)

что возвращает нас с учетом определений потенциала Рисса и двойственного оператора к формуле (16.5). Другой вид приобретает формула обращения при =1:

. (16.18)

Здесь потенциал Рисса в отличие от (16.17) действует на функцию двух переменных. Поэтому для преобразования Фурье от имеем:

, (16.19)

где - преобразование Фурье от :

.

В развернутом виде формула обращения (16.17) приобретает следующий вид:

. (16.20)

17. Численная реализация алгоритмов восстановления и понятие о некорректных задачах

Существование различных формул обращения в какой-то мере объясняет и наличие множества численных алгоритмов восстановления, которые при наличии точных значений на непрерывном множестве точек , s (0, S_min s S_msx) и возможности точного вычисления интегралов (например, аналитически) позволяют решить задачу. Однако, в реальности, когда известны приближенные значения на дискретном множестве точек , получение достоверного решения наталкивается на принципиальные трудности, которые можно разделить на проблему дискретизации и разрешения и проблему некорректности задачи восстановления.

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

. (17.1)

Ширина спектра этой функции составляет . Для более общего вида функции ширина спектра составляет b=a. В качестве меры ширины спектра будем для удобства использовать сторону квадрата a (ширина спектра по каждой переменной):

где через обозначен квадрат на плоскости S.

Возможность восстановления функции по ее значениям на дискретном множестве точек регламентируется теоремой отсчетов:

Теорема. Пусть f - функция с шириной спектра a, h - шаг дискретизации, выбранный из условия ha. Тогда f однозначно определяется величинами , где - двухмерный вектор с целочисленными проекциями и допускает представление

(17.2)

сходящееся в среднем на S. При этом преобразование Фурье:

(17.3)

сходится в среднем в .

Условие ha называется условием Найквиста. Ряд (17.1) представляет сумму однотипных функций со сдвигом аргумента на целочисленное число шагов вдоль координат x и y. Это означает, что разложение (17.1) не может представлять функцию с размерами деталей вдоль каждой координаты меньше 2h. Если условие Найсквита не выполняется, выбору значений называют недостаточной. При этом могут порождаться как ложные детали размером не более удвоенного шага выборки, так и глобальные искажения изображения.

Восстановление по ее выборке - это последний этап всей процедуры восстановления. До этого мы должны надежно восстановить функцию Радона по ее значениям на дискретном множестве узлов и параметров s. Условия восстановления функции Радона удобно сформировать для , определенной на шаре единичного радиуса (чего всегда можно добиться выбором масштаба пространственных переменных). Если имеет ширину спектра а, то такую же ширину имеет спектр и по переменной s. Выборка будет достаточной, если выполняется условие Найквиста:

s_i = i/q (-q i q) q a/.

Необходимое количество направлений р по угловой координате _i (0 i p-1) также связано с шириной спектра. Функцию с шириной спектра а<р можно надежно восстановить по различным направлениям из интервала [0,]. Можно таким образом утверждать, что условия дискретизации и разрешения исследованы достаточно для практического восстановления алгоритмов восстановления. Основной трудностью при этом является то, что функции в компьютерной томографии обычно отличны от нуля в области конечных размеров (финитны), а, следовательно, обладают неограниченным спектром. Поэтому необходимо предварительно апроксимировать также функции функциями с ограниченным спектром. Обычно это удается сделать с удовлетворительной для практики точностью.

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

.

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

(х, у) = 0, тогда и только тогда, когда х = у;

(х, у) = (у, х);

(х, у) + (у, z) (х, z).

Конкретный вид метрики в каждом из рассматриваемых множеств F и R может быть свой, поэтому целесообразно метрики обозначить соответственно _f и _R. Множества функций F и R с введенными на них метриками называют метрическими пространствами. Задача нахождения решения f F по “исходным данным” R R называется корректно поставленной на паре метрических пространств (F, R), если выполнены условия:

Для каждой функции R R существует решение f F;

Решение определяется однозначно;

Задача устойчива на пространствах (F, R), т.е. для любого >0, можно указать такое ()>0, что из неравенства _R(R₁, R₂)() следует _f(f₁, f₂), где , .

Задачи, не удовлетворяющие хотя бы одному из перечисленных условий, называются некорректно поставленными. Важно отметить, что определение некорректно поставленной задачи относится к конкретному выбору метрик в пространствах F и R. При выборе других метрик та же задача может оказаться корректно поставленной.

Практически некорректность задачи компьютерной томографии будет проявляться в том, что если пользоваться точными формулами обращения, тот даже при небольших искажениях в исходных данных решение может существенно отличаться от искомой функции f(). Для получения решения, близкого к искомой функции, необходимо использовать регуляризованные формулы обращения. Рассмотрим, например, формулу (9.20), восстанавливающую f() по F₀() - преобразованию Фурье суммарной обратной проекции. Если исходные данные искажены некоторой погрешностью n(, s)

(17.4)

Обратную проекцию от добавки n(, s) обозначим F_оn(), а ее преобразование Фурье соответственно F_оn(). После подстановки в формулу (9.20) получим:

. (17.5)

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

Регуляризованная формула обращения имеет вид:

. (17.6)

Пользуясь (17.6) мы ограничиваем так же и спектр искомой функции, но сознательно идем на этот для получения устойчивого решения. Мы рассмотрели один из методов регуляризации. Более полное представление о методах регуляризации можно получить из работ [9] [7].

Диагностика в средне- и длинноволновой области

В случае, когда размеры и характерные детали изучаемого объекта становятся соизмеримыми с длиной волны, модель компьютерной томографии может в лучшем случае использоваться как грубое приближение к точной модели, поскольку все большую роль начинают играть отклонения от оптических законов дифракционные явления. В пределе, когда длина волны оказывается либо больше размеров объекта структура полей становится статической. В рассматриваемой области длин волн получить в явном виде формулы обращения не удается и в большинстве случаев приходится довольствоваться модельно прямой задачи. Для линейных сред и монохроматического зондирующего поля эта задача имеет следующую постановку. Пусть исследуемый объект занимает область пространства V_i с границей s, на которой определена внешняя нормаль . Электрофизические параметры объекта описываются магнитной проницаемостью _i и комплексной диэлектрической проницаемость _i как функциями пространственных координат. В частности, V_i может быть представлена в виде совокупности кусочно-однородных областей. Зондирующее поле внешних источников вызывает появление рассеянного объектом поля во внешней области и внутреннего поля .Поле внешних источников считается известным, а поля , подлежат определению. Во внешней области V_е, дополняющей V_i до полного пространства, поле удовлетворяет однородным уравнениям Максвелла:

(17.7)

в отличие от поля , удовлетворяющего соответственно неоднородными уравнениями

(17.8)

где - известная плотность зондирующих электрических токов. (В некоторых задачах во второе уравнение удобно вводить плотность магнитных токов ). Поле считаем невозмущенным, каким бывает поле в отсутствие объекта V_i. Т.е. уравнения (11.2) формально выполняются во всем пространстве V = V_e V_i и имеют явное решение:

(17.9)

где - векторный потенциал токов . Векторный потенциал вычисляется по известной формуле:

, (17.10)

где - функция Грина однородного пространства ( - расстояние от точки интегрирования до точки наблюдения ), область размещения источников зондирующего поля (W V_e).

Внутреннее поле удовлетворяет однородным уравнениям Максвелла в области V_i:

(17.11)

Предельные значения суммарного поля внешнего пространства и внутреннего поля соответственно извне и изнутри должны быть связаны граничными условиями сопряжения:

(17.12)

Задача нахождения известных полей , по заданному внешнему полю и известным V_i, , , _е, _е называется прямой задачей и сводится к решению уравнений (17.7) (17.11) при граничных условиях (17.12) и условиях излучения, налагаемым на внешнее поле [10].

, (17.13)

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

Задача вычислительной диагностики относится к обратным задачам. Известным является рассеянное поле (полное поле ) на некотором множестве точек на поверхности S_e. По его значениям необходимо восстановить параметры _i, _i объекта. Видно, что искомые функции _i, _i сложным образом входят в систему дифференциальных уравнений поля, решаемых также при достаточно непростых условиях (17.12) (17.13). Недавние исследования этой задачи показали, что её решение неоднозначно. Так что к трудностям ее решения, связанным с неустойчивостью, добавляется проблема неоднозначности. Становится очевидной принципиальная необходимость в априорной информации об электрофизических свойствах и геометрических характеристиках исследуемого объекта, обеспечивающей единственность решения задачи.

Список литературы

Даутов О.Ш. Теоретические основы вычислительной диагностики . Учебное пособие. (Электронная версия). Казань, КГТУ им. А.Н.Туполева, 2008 г., - 97 с.

Календер В Компьютерная томография: основы, техника, качество изображений и области клинического использования. М., Техносфера, 2006 г., - 343 с.

Рис У.Г. Основы дистанционного зондирования. М., Техносфера, 2006 г. - 335 с.

Гонсалес Р., Вудс Р. Цифровая обработка изображений. - М., Техносфнера. -2005 г.- 312 с..

5. Компьютеры и суперкомпьютеры в биологии. Под ред. В.Д. Лахно и М.Н.Устинина. - Москва - Ижевск: Институт компьютерных исследований, 2002 г., - 528 с.

6. Тихонов А.Н., Арсенин В.П., Тимонов А.А. Математические задачи компьютерной томографии. М.: Наука. Гл. ред. физ.-мат. лит., 1987. - Проблемы науки и технического прогресса. - 160 с., ил.

7. Дж. А. Стрэттон. Теория электромагнетизма. Пер. с англ. М.С. Рабиновича, В.М. Харитонова. Под.ред. С.М. Рытова. ОГИЗ, ГИТ.-ТЛ. М.,-Л.- 1948. -539 с.

8. А.В. Ржанов, К.К. Свиташев, А.И. Семененко, Л.В. Семененко, В.К. Соколов. Основы эллипсометрии. Новосибирск, “Наука”, 1978. - 424 с.

9. К.К. Свиташев, А.И. Семененко, Л.В. Семененко. Оптический поляризационный метод исследования поверхности проводников // Некоторые проблемы физики и химии поверхности полупроводников. Сб. ст. под ред. А.В. Ржанова. “Наука”, Сиб. отделение, 1972. -с.114-180.

10. Г.Т. Марков,А.Ф. Чаплин. Возбуждение электромагнитных волн. М.-Л.: Энергия, 1967. - 376 с.

11. Л.Д. Ландау, Е.М. Лившиц. Электродинамика сплошных сред. ГИТ-ТЛ, М., 1957. -532 с., ил.

12. Наттерер Ф. Математические аспекты компьютерной томографии. Пер. с англ. - М.: Мир, 1990. -288 с.: ил. ISBN 5-03-

13. Троицкий И.Н. Статическая теория тотмограии. -М.: Радио и связь, 1989. -240 с.: ил. -ISBN 5-256-00182-5.

14. А.Н. Тихонов, В.Я. Арсенин. Методы решения некорректных задач. Главная ред. физ.-мат. лит. из-ва “Наука”, -М., 1974. -224с., ил.

15. Хенл Х., Мауэ А., Вестпфаль К. Теория дифракции. Из-во “Мир”, 1964.

Размещено на Allbest.ru