О РЕШЕНИИ ЗАДАЧИ СИНТЕЗА МАГНИТНОГО ПОЛЯ
В МР-ТОМОГРАФЕ МЕТОДОМ РЕГУЛЯРИЗАЦИИ С ОГРАНИЧЕНИЯМИ И МЕТОДОМ АППРОКСИМАЦИИ
5-ФУНКЦИЕЙ С.В. Иванов, Н.Г. Рущенко, В.С. Сизиков, Д.Ю. Соколов, Е.В. Хомутникова
Среди многочисленных задач магнитно-резонансной томографии (МР-томографии) [1-7], [8, с. 33-62] рассмотрим задачу синтеза магнитного поля в МР-томографе, уже рассмотренную в ряде работ (1, 2, 4, 7-12 и др.).
Формулировка задачи
В МР-томографии исключительно важной является задача создания высокооднородного магнитного поля внутри МР-томографа. Относительная неоднородность напряженности Н основного (статического) поля в пределах рабочего
объема томографа должна быть АН / Н ~10-5 -10-6 [4-7]. Это позволит с высоким разрешением решать основную задачу МР-томографии - задачу реконструкции МР-изображений [5], [6, р. 270-274], [8, с. 45-51].
Как правило, задача создания высокооднородного магнитного поля решается путем введения соленоидальных компенсирующих (корректирующих) катушек различного порядка [1, с. 274], [7]. При этом параметры катушек определяются из условия компенсации членов ряда Тейлора, в который разложено магнитное поле. Однако это довольно сложная техническая (а также математическая) задача.
Рассмотрим другой подход к решению задачи создания высокооднородного поля в МР-томографе. Этот подход основан на расчете такого закона распределения тока вдоль обмотки катушки томографа, что он обеспечит высокую однородность магнитного поля внутри катушки. Та и другая задача может ставиться как прямая, или задача анализа (расчет магнитного поля по распределению тока), и обратная, или задача синтеза (расчет распределения тока и конфигурации катушек по заданному магнитному полю) [10, с. 17]. Обратная задача является значительно более сложной, чем прямая. При этом рассматривается локальный метод синтеза (к нему относится метод компенсирующих катушек [1,7]) и интегральный метод синтеза на оси катушки [2, 9, 11, 12] и во всем объеме [7].
Задача синтеза
Рассмотрим задачу синтеза магнитного поля внутри катушки МР-томографа [4], [8, с. 55-61], [9], [10, с. 17], [11, 12], причем сначала рассмотрим ее простейший вариант - синтез на оси цилиндрической катушки с бесконечно тонкой обмоткой.
Одним из первых задачу интегрального синтеза магнитного поля сформулировал с учетом ее некорректности и решил с применением метода регуляризации Тихонова К. Адамиак [9]. Эту задачу рассмотрел также Л. Луганский [2] и др. [11,12]. В данной работе дается дальнейшее развитие этой задачи.
Итак, рассмотрим следующую обратную задачу МР-томографии: определить распределение плотности тока J(а) вдоль бесконечно тонкой обмотки цилиндрической катушки по заданному магнитному полю (напряженности) Н(г) на ее оси. Это - задача интегрального синтеза магнитного поля на оси катушки МР-томографа [4]. Здесь а -
координата вдоль обмотки катушки, z - координата вдоль ее оси. Рассмотрим случай H(z ) = const.
На рис. 1 представлена цилиндрическая катушка с бесконечно тонкой обмоткой, где a е [-/, l] - расстояние вдоль обмотки катушки, z е [-/, l] - расстояние от центра катушки вдоль ее оси, l - полудлина катушки, H(z) = H = const, z е [-/, l] - заданная напряженность магнитного поля на оси катушки, J(a), a е [-l, l] - искомое распределение тока вдоль обмотки катушки, R - радиус катушки.
J (a)
/_ a
н ( z) • / / ^R
-l 0 l z
Рис. 1. Цилиндрическая катушка с бесконечно тонкой обмоткой
Согласно [1, с. 37], обмотка в случае J (a) ^ const называется соленоидом с переменной плотностью тока или с неоднородным распределением тока. Реализовать неоднородное распределение тока можно, например, путем разбиения соленоида на несколько отдельных секций [1, с. 37]. Другие способы (создание изолированных витков с индивидуальной подводкой тока через единый источник или создание соленоида с единым, намотанным на цилиндр проводом, но имеющим на каждом витке свое сопротивление) предложены в работах [11, 12].
Вывод интегрального уравнения задачи синтеза
Как известно [1, с. 23], [7, 9], [10, с. 20], [13, с. 330-333], [14, с. 327], напряженность магнитного поля на оси тонкого витка (кругового тока) равна (см. рис. 2)
3 К2
Н(г) = к ,_ _=■, -х< z <х,
V
R2 + ( z - a )2
где к - коэффициент пропорциональности, 3 - ток в витке, К - радиус витка, z -координата точки, в которой вычисляется Н, а - z-координата центра витка (ось z направлена перпендикулярно плоскости витка).
R^ н (z)
0 а 1 z z
Рис. 2. Напряженность на оси тонкого витка
Для упрощения дальнейших записей будем полагать к = 1. Тогда
H ( z ) =
JR2
[2 + (z - a )2
— x < z < X.
(1)
Далее полагаем, что имеется бесконечное множество витков, намотанных на цилиндр радиуса Я и полудлины I и имеющих плотность тока J = J(а). В этом случае, интегрируя (1) по виткам, получим следующее выражение для суммарной напряженности магнитного поля в точке с координатой г на оси цилиндра с бесконечно тонкой обмоткой (см. рис. 1).
H (z)=/
J (a) R da
V[R2 + (z — a)2
— X < z < X .
(2)
Вычисление поля H (z) по заданному току J (a), согласно (2), есть прямая задача. Сделаем важное замечание: формула (2) справедлива как для точек внутри цилиндра ( | z | < l), так и вне его ( | z | > l). Выполним краткий анализ прямой задачи. Если J(a) = J = const , то интегрирование (2) дает:
H ( z ) =
z+l
z — l
■Jr2 + (z +1 )2 д/R2 + (z — l )2
J , — x < z < x .
(3)
Частные случаи (3) - поле в центре катушки h (0) = 2 J + (R /1 )2
, поле на краю
катушки H(l) = J+ (R/2l)2 . Если R = l, то H(0) = 2 JV2, H(l) = 2 JV5, т.е.
напряженность в центре катушки в д/5/2 «1.58 раз больше, чем на ее краях. Если
2 3
z —^ x, то [13, с. 333] H(z) — 2 JR l/z , т.е. вне катушки поле H(z) падает как
— 3
~ z . Такая асимптотика справедлива и для тонкого витка (см. (1)), а также будет иметь место и для цилиндрической катушки при J(a) ^ const .
Если l >> R (бесконечно длинная катушка), то H(0) = 2 J, H(l) = J .
Если, например, l = 10R, то H(0) = 1.9901 J, H(R) = 1.9898 J, т.е.
| H (R ) — H (0)|/H (0) = 1.5-10 —4. Видим, что в случае длинной катушки в принципе
можно получить неоднородность поля порядка 10—4 в пределах | z | < R, однако такую катушку технически трудно изготовить.
Вышеприведенный краткий анализ прямой задачи показывает, что при J(a) = const поле H(z) падает от центра катушки к ее краям, а вне катушки
_ 3
H(z ) ~ | z | при | z | —X.
Перейдем к рассмотрению обратной задачи. Запишем (2) в виде:
R
-l^2 + ( z — a )2
-J ( a ) da = H ( z ),
— l < z < l.
(4)
Соотношение (4) есть интегральное уравнение Фредгольма I рода, где H(z)-задаваемая правая часть (напряженность магнитного поля на оси катушки), например, H(z) = H = const, а J(a) - искомая функция (распределение тока вдоль бесконечно тонкой обмотки цилиндрической катушки МР-томографа). Задача решения уравнения (4) является некорректной [8, с. 178].
Используем безразмерные переменные: s = ^R, х = z/R, ^ = l/R . Тогда (4) можно записать в виде
l
s0
J K ( x, s ) J ( s ) ds = H ( x ),
- s0 ^ x ^ s0,
s0
где
K (x, s) =
VR+(x - s)2 ]3
(5)
(6)
- ядро интегрального уравнения.
Решение уравнения (5) и техническая реализация решения позволят в принципе создать на оси катушки МР-томографа поляризующее поле с заданным законом изменения его напряженности, в частности, H(x) = H = const .
Решение уравнения методом регуляризации с ограничениями
Рассмотрим вопрос о решении уравнения (5).
Вначале выполним качественный анализ. Заметим, что запись (5)-(6) справедлива не только для x е [-s0,s0], но и для x £ [-s0,s0], т.е. для x е (-да, да).
Пусть H(x) = const при x е [-s0, s0]. Тогда при x £ [-s0, s0] поле H(x) ^ const и
будет монотонно падать с ростом | x |, переходя при | x да в асимптотику
H(x) ~| x |-3 (см. рис. 3).
—о 0 ¿0 х
Рис. 3. Поведение напряженности поля Н (х) внутри и вне катушки
Видим, что функция Н(х) имеет разрывные производные Н'(х), Н"(х), ... при х = ± . Отсюда следует, что искомая функция 3(¿) уравнения (5) должна быть сингулярной при ^ = -¿0 + 0 и ^ = ¿0 -0, а при | 51>¿0 ток 3(5) = 0, как это следует из физической постановки задачи. Кроме того, 3(-¿) = 3(¿) , если Н(-х) = Н(х) .
Таким образом, функция 3(¿) должна быть сингулярной и симметричной, т.е. иметь качественный вид, изображенный на рис. 4.
Теперь перейдем к количественному алгоритму решения.
Уравнение (5) в принципе может быть решено аналитически. Однако аналитического решения получить не удается. Поэтому рассмотрим вопрос о численном решении уравнения (5). Учитывая некорректность уравнения (5), воспользуемся методом регуляризации Тихонова [8, 15].
В работе [9] также использовался метод регуляризации Тихонова. При этом, как показало численное решение, при малых значениях параметра регуляризации а регуляризованное решение 3а (¿) имеет большие знакопеременные флуктуации. В то же время параметр регуляризации а должен быть действительно очень мал (порядка
1
10 5 -10 10). Это обусловлено тем, что погрешность задания Н(х) отсутствует (есть только погрешность численного алгоритма). Кроме того, регуляризованное решение Аа (з) должно иметь очень большой диапазон значений (из-за сингулярности
точного решения). При этом решение Аа (з) должно быть монотонно возрастающим от центра катушки (з = 0) к ее краям (з = ± з0 ).
А(з)
—0 0 з0 з
Рис. 4. Принципиальный вид функции А (з)
Чтобы регуляризованное решение Аа (з) было таковым, а именно,
неотрицательным и монотонно возрастающим от центра катушки к ее краям, нужно воспользоваться методом регуляризации Тихонова с ограничениями на решение [15, с. 118].
Ограничения в данном случае состоят в том, что решение Аа (з) ищется на множестве неотрицательных монотонно невозрастающих функций [15, с. 118]. Однако такому условию удовлетворяет лишь левая половина решения Аа(з), т.е. при з е (-з0,0]. Поэтому модифицируем уравнение (5), учитывая симметрию решения: А(-з) = А(з) . Для этого запишем (5) в виде:
Н (х)=Й
0
= 1
А (з) йз
(х - з)2]; А (з) йз
з0 I
А (з) йз
Л
1 + (х - з)"
и
I
1 + (х - з) • А (-з) йз
1 + (х + з) '
Поскольку А ( з ) = А (-з ), т.е. распределение тока симметрично, то Н ( х ) = Н ( - х ), т. е. магнитное поле Н(х) также симметрично. Справедливо и обратное: если Н ( х ) = Н ( - х ), то А ( з ) = А (-з ).
В результате уравнение (5) может быть записано в виде
| Я( х, з ) А ( з ) йз = Н ( х ),
_з0
где новое ядро равно Я( х, з) =
1
- з0 < х < 0,
1
Я
(7)
(8)
[1 + (х + з)2 ]3 д/[1 + (х - з)2 ]3 Итак, А (з) будем искать в левом полупространстве путем решения уравнения (7), т.е. при з е (-з0,0], после чего в правом полупространстве А(з) = А(-з), з е [0,з0). При
этом дополнительно будем учитывать, что функция 3(5), 5 е [—,0], является неотрицательной монотонно невозрастающей (типа изображенной на рис. 5).
J(s)
-0
0 s
Рис. 5. Ветвь функции J(s)
Для получения численного решения Ja (s) методом регуляризации Тихонова на
множестве неотрицательных монотонно невозрастающих функций можно воспользоваться программой PTIPR на Фортране [15, с. 118, 174]. В этой программе выполняется минимизация сглаживающего функционала методом проекции сопряженных градиентов на множество неотрицательных монотонно невозрастающих функций при каждом значении параметра регуляризации а > 0 .
В [11, 12] приведены результаты численного решения уравнения (7) при H (х) = const методом регуляризации Тихонова с ограничениями на решение.
Регуляризованное решение Ja (s) при s ^ -s0 + 0 имеет крутое нарастание (тем
больше, чем меньше а ), а в остальной области практически Ja (s) = const. Технически
это означает, что катушка должна состоять из соленоида с однородным распределением тока, за исключением нескольких крайних витков, в которых ток идет по нарастающей к краю катушки.
Однако, как видно из работ [11, 12], напряженность магнитного поля
не получается высокооднородной ни при каких значениях а.
К более однородному полю Н(х), как следует из [11, 12], приводит решение задачи синтеза путем приближения 3а (5) 5 -функцией. Рассмотрим этот вопрос подробнее.
Положим, что искомое распределение тока 3(5) представляется суммой двух 5 -функций, точнее, полу-5 -функций:
0
Решение путем приближения к 5 -функции
J(s) = 51 /2 (s + s0) + 51/2 (s - s0) (см. рис. 6). Тогда в соответствии с (5) и (6),
(9)
2LV[1 + (х + s0)2]3 V[1 + (х-s0)2]3 J. (10)
1
+
1
А А(з)
-1-
-з0 0 з0 з
Рис. 6. Представление функции ^в) в виде суммы полу- 5 -функций
На рис. 7. представлена Н(х) для ряда значений з0, вычисленная согласно (10). Видим, что при з0 = 0.5 поле Н(х) имеет для такой простой катушки степень однородности | АН | / Н ~ 10-4 для | х |< 0.4з0.
-1.1
1.1
Рис. 7. Напряженность Н(х) вдоль оси Данный результат можно улучшить, если представить распределение тока А (з) в
виде
А(з) = А0 • [51/2 (з + з0) + 51/2 (з - з0)] + А
(11)
Подбирая з0 и отношение А / А0 можно с помощью (11) получить однородность
| АН |/ Н ~10-4 для | х |< 0.4з0.
В работах [11, 12] А(з) представлен суммой "усеченных" 5 -функций:
|А0, з = ±з0,
А (з) =
А, | з |< з0
(12)
где А0 и А - конечные константы, причем А0 >> А. В этом случае, подбирая А0 / А, также можно добиться | АН |/Н ~ 10 4 -И0 5 для | х |< 0.4з0.
X
Можно сделать следующий вывод. Если распределение тока J(s) вдоль обмотки катушки МР-томографа сделать однородным, за исключением крайних витков (т. е. типа (9), (11), или (13)), то это позволит повысить однородность магнитного поля H (х) на оси катушки. Данный способ (приближение к 5 -функции), хотя и является более простым, чем метод регуляризации Тихонова с ограничениями на решение, но приводит к более точному результату (к более однородному полю), т. е. дает более адекватное математическое описание задачи и, как следствие, более точное решение.
Литература
1. Монтгомери Д.Б. Получение сильных магнитных полей с помощью соленоидов. М.: Мир, 1971.
2. Луганский Л.Б. Расчет соленоида с заданной формой магнитного поля // Журн. техн. физики. 1985. Т. 55. Вып. 7. С. 1263-1271.
3. Тихонов А.Н., Арсенин В.Я., Тимонов А.А. Математические задачи компьютерной томографии. М.: Наука, 1987.
4. Лухвич А.А., Чурило В.Р. Источники поляризующего магнитного поля и его градиентов для я.м.р.-томографии (обзор) // ПТЭ. 1987. № 5. С. 172-173.
5. Лич М. Получение ЯМР-изображений с пространственной локализацией // Физика визуализации изображений в медицине. / Под ред. С. Уэбба. М.: Мир, 1991. Т. 2. С. 105-231.
6. Cho Z.H., Jones J.P., Singh M. Foundations of medical imaging. NY: Wiley, 1993.
7. Галайдин П.А., Замятин А.И., Иванов В.А. Расчет и проектирование электромагнитных систем магниторезонансных томографов. / Уч. пособие. СПб: ИТМО, 1998.
8. Сизиков В.С. Математические методы обработки результатов измерений. СПб: Политехника, 2001 (Электронный вариант: Сизиков В.С. Устойчивые методы обработки результатов измерений. Эл. учебник. http://de.ifmo.ru/library/sizikov.pdf) .
9. Adamiak K. Method of the magnetic field synthesis on the axis of cylinder solenoid // Appl. Phys. 1978. Vol. 16. P. 417-423.
10. Афанасьев Ю.В., Студенцов Н.В., Хорев В.Н. и др. Средства измерений параметров магнитного поля. Л.: Энергия, 1979.
11. Сизиков В.С., Марусина М.Я., Иванов С.В., Колобухова Т.Б., Николаев Д.Б., Соколов Д.Ю., Хомутникова Е.В. Прямая и обратная задачи синтеза магнитного поля в ЯМР-томографе. // Научно-техн. вестник СПбГИТМО(ТУ). 2001. Вып. 3. С. 209-214.
12. Сизиков В.С., Ахмадулин Р.И., Николаев Д.Б. Синтез магнитного поля вдоль оси катушки ЯМР-томографа. // Изв. вузов. Приборостроение. 2002. Т. 45, № 1. С. 5257.
13. Фриш С.Э., Тиморева А.В. Курс общей физики. Т. 2. МЛ: ГИТТЛ, 1952.
14. Дружкин Л.А. Задачи теории поля. М.: МИРГЭ, 1964.
15. 15.Тихонов А.Н., Гончарский А.В., Степанов В.В., Ягола А.Г. Численные методы решения некорректных задач. М.: Наука, 1990.