К теории определения поля перемещений на основе уравнения переноса в дискретном случае Текст научной статьи по специальности «Математика»

ВЕСТНИК САНКТ-ПЕТЕРБУРГСКОГО УНИВЕРСИТЕТА

Сер. 10. 2010. Вып. 3

УДК 517.9:519.6 Е.Д. Котина

К ТЕОРИИ ОПРЕДЕЛЕНИЯ ПОЛЯ ПЕРЕМЕЩЕНИЙ НА ОСНОВЕ УРАВНЕНИЯ ПЕРЕНОСА В ДИСКРЕТНОМ СЛУЧАЕ

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

До настоящего времени из всех математических моделей наибольшее применение в радионуклидной диагностике нашли линейные камерные модели, теория которых достаточно хорошо разработана. Они описывают кинетику индикатора системой линейных однородных дифференциальных уравнений первого порядка с постоянными коэффициентами [1]. Также используются циркуляционные модели на основе принципа Стюарта-Гамильтона и интегральных уравнений [2]. При изучении физиологических систем, в которых пространственные характеристики транспорта индикатора не менее важны, чем его временные параметры, эти модели не обеспечивают получение нужной диагностической информации, в таких случаях для построения моделей применяют дифференциальные уравнения в частных производных [2].

В данной статье предлагается новая математическая модель и разрабатывается математический аппарат для определения поля перемещений РФП на основе уравнения переноса для дискретного случая.

Учитывая дискретный характер измерений, представляется естественным использование дискретной системы для описания транспорта РФП.

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

ремещения РФП; к - номер кадра. Наша задача состоит в определении данного вектора. В этой связи рассмотрим распределение плотности РФП по кадрам и будем обозначать плотность (концентрацию) РФП р = р(к,х(к)).

Изменение плотности РФП вдоль траекторий системы (1) происходит в соответствии с уравнением переноса [3, 4]

Котина Елена Дмитриевна — кандидат физико-математических наук, доцент кафедры теории управления факультета прикладной математики-процессов управления Санкт-Петербургского государственного университета. Количество опубликованных работ: 56. Научные направления: математическое моделирование, численные методы, методы оптимизации. E-mail: [email protected].

© Е.Д. Котина, 2010

x(k +1) = x(k) + u(k, x(k)), k = 0, 1, 2, ...,

(і)

вектор пе-

p(k + 1,x(k + 1)) = A 1 (k,x(k))p(k,x(k)), k = 0,1, 2,... .

(2)

В (2) А(к,х(к)) - якобиан преобразования (1),

где

А(к, х(к))

ди(х(к)) ґ и1,х2

дх(к + 1)

дх(к)

Е +

ди(х(к))

дх(к)

и2,х1 и2,х2

дх(к)

, при использовании таких обозначений:

(3)

диі(к,х(к)) диі(к,х(к)) ди2(к,х(к)) ди2(к,х(к))

_ о 7ТТ — ^1,Ж2? о 7~Г\ — ^2,Х2 і о 7~Г\ — ^2,жі-

дх1(к)

дх^(к)

дх^(к)

дх1(к)

Для простоты дальнейших преобразований будем пока рассматривать только один шаг процесса (т. е. два соседних кадра):

х(1) = х(0) + и(0,х(0))

и обозначать

х(0) = х, и(0, х(0)) = и(х) = и, х(1) = х + и(х) = х + и, р(0, х(0)) = р(х), р(1, х(1)) = р(х + и(х)) = р(х + и).

Для определения вектора перемещения и = и(х) введем следующий функционал:

Фсіх,

(4)

Ш

где

Ф

Е +

Ф = Ф2 + а2 • П2, ди(х)

дх

р(х + и) — р(х),

^2 = (и1,х1)2 + (и1,х2)2 + (и2,х1)2 + (и2,х2)2,

а2 - весовой коэффициент.

Далее рассмотрим задачу минимизации функционала (4) по и. Он состоит из двух слагаемых, минимизация первого дает нам вектор перемещения, а второе является ре-гуляризирующим фактором, часто используемым при определении оптического потока [5, 6]. Перемещение и, доставляющее минимум функционалу (4), есть решение нашей задачи.

Заметим, что

1 п. ди

1 + шум +

Р + ди Е+ді

Запишем Ф в виде

Ф = 1+ р1(\\уи + в2

ди

дх

дх

р(х + и) — р(х),

где = р2 = 1. (Заметим, что если положить в = 0, то мы рассматриваем в (3) только линейную часть, если же в = в = 0, то получаем, что плотность вдоль траектории не меняется, т. е. имеем случай так называемого оптического потока [5, 6].)

Запишем уравнения Эйлера-Лагранжа, которым должна удовлетворять вектор-функция и, минимизирующая функционал (4):

дФ д дФ д дФ

О,

dui dxi dui,xi дх2 дп1}Х2

дФ д дФ д дФ

ди2 дхi ди2,х1 дх2 дщ,х2

После некоторых преобразований получаем

0.

т г, г, ди \ др(х + и) . „ . . . дФ

Ф ( 1 + /?idiv U (32 — )---------—-------------(/?1 + ^2^2,х2)р{х + U)~Q^T ~

/п п \тдр(х + и) дФ

— (/?1 + p2li2,x2)^-----д--------Ь p2U2,xlp{x + и)—---------h

дх1 дх2

+ Ф/32Ц2,Ж1 + ^ - a2div grad mi = 0, (5)

дх2

др(х + и) дФ др(х + и)

1 rJ2Ul 1 1 -

п ди \ др(х + и) \дФ

ф 1 + /?ldiv U + р2 ТГ~ ----------д------Ь p2Ui,x2Р(х + M) Tj-------Н ^р21Ч,х2

\ дх J ди2 дх1

I 2'-" 1,Х2Г\~ 1 У ^ 1 АГ^1,Х2 гл

ди2 дх1 дх1

........................дФ _ „ чдр( х + и) 2

- (/?i + fhuitXi)p{x + и)---------Ф(/?1 + p2Ui,xi)--------------a div grad и2 = 0. (6)

дх2 дх2

Далее для удобства записи введем следующие обозначения:

др(х + и)

дхi

dpxijx + и) dxj

= Pxi(х + и), i =1, 2, (7)

~ pxi,xj(х + u), i,j 1 2, (8)

если при дифференцировании не учитывается зависимость и от х.

Тогда с учетом (7) запишем

др(х + и)

= Рх1{Х + и){ 1 + и\х\} + рх 2(х + и)и2,х1, (9)

дх1

др(х + и) дх2

= Рх1 (х + и)и1,х2 + Рх2 (х + и)(1 + U2,x2), (10)

если при дифференцировании учитывается зависимость и от х.

Заметим, что в формулах (5), (6) и^ = рх1(х + и), г = 1,2, а _

ди* дх*

г = 1, 2, вычисляются по формулам (9, 10).

Обозначим

ди

7 = 7(и) = и + /?2 •

Тогда Ф = (1 + 7)р(х + и) — р(х), а уравнения (5), (6) принимают такой вид:

( дФ дФ \

Ф(1 + 7)/эж1(ж + "“)+(— ~о^г(@ 1 + /32^2,х2~) + р2и2,х1 ) р(х -\- и) —

- Ф(в + в2и2,х2)(Рх1(х + u)(l + U1,x1) + Рх2(х + u)u2,x1) +

+ Ф^2U2,x1 (Рх1 (х + u)u1,x2 + Рх2(х + u)(1 + U2,x2)) — a2div gradU1 = 0, (11)

дФ дФ

ф(1 + j)Px2(x + и) + ( -7^ihui,x2 - + P2Ul,xl Jp(x + u) +

+ Ф p2U1,x2(Px1 (х + u)(1 + U1,x1) + Px2 (х + u)u2,x1) —

— Ф(в1 + вU1 ,x1 )(px 1 (х + u)u1,x2 + Px2(х + u)(1+ U2,x2)) — a2div gradU2 = 0. (12)

Можно также записать уравнения (11), (12) в векторной форме:

((1 + y^rad^^ + u) + 0 grad р(х + и))Ф + р(х + u)0 grad Ф — а2Аи = 0,

где 0 = 0(и) = —&Д — в ( U2'x2 —U2’x1); Аи = (А^ = fdivgradU11-

V —U1,x2 U1,x1 J \Au2j ydiv grad U2 J

Отметим, что gradp^ + u) вычисляется следующим образом:

^ Px1 (х + U)(1 + U1,x1) + Px2 (х + U)u2,x1

grad р(х + u) = i , , , W1

' Рх1(х + u)U1,x2 + Px2 (х + u)(1+ U2,x2

в то время как

^ Рх 1 (х + и)

Решая систему (11), (12), получаем искомый вектор перемещения и.

Построение итерационной схемы. Для ее решения построим итерационную схе-

к 1,

му. Обозначим ик = (uk ,и%)т и запишем итерационное уравнение

((1 + yk)gradup(x + ик) + 0kgrad р(х + ик))Фк+1 +

+ р(х + ик)0кgrad Фк+1 — а2 Аик+1 = 0. (13)

Здесь и далее значок к означает, что вместо переменной и берем итерационную переменную ик, т. е. Y,t = y(u^, 0к = 0(ик), Фк = Ф(х,ик) и т. д. Заметим, что хотя здесь

мы написали р(х + ик), иногда, для краткости, будем также использовать обозначение рк = р(х + ик).

Определим ик+1 так:

к+1 к | j к

и = и + du ,

т. е. разделяем неизвестную итерационную переменную ик+1 на известную ик и неизвестное приращение dик.

Рассмотрим, как будет вычисляться Фк+1:

Фк+1 = (1 + y к+1)р(х + ик+1) — р(х). (14)

Для этого рассмотрим сначала y^1:

д(ик + dU)

Y 1 = в^^(ик + du) + в2 = e1div ик + ^div^u^ + в2

дх

дик д^ик)

+

дх дх

(15)

Вычислив определитель (15) записать в виде

ди2 d(duk)

Yk+1 = ^idiv uk + в2

dx dx

duk

dx

и сделав некоторые преобразования, можно d(duk)

+ eidiv(du ) + в2

dx

+ uk (uk ,duk),

где

u (u ,du )= в((u 2,x2, -u2,xi)grad(du1) + (uljxi, -ul x2)grad(du2)).

Нетрудно увидеть, что получилось следующее:

7 к+1(ик+1) = 7к(ик)+ 7к (¿ик)+ик (ик,Сик). Возвращаемся к формуле (14), переписав ее таким образом:

Фк+1 = (р(х + ик+1) — р(х))(1 + 7 к+1(ик+1)) + р(х)7к+1 (ик+1).

(16)

Напомним, что

р(х + ик+1) = р(х + ик + Змк).

Далее с учетом вышеизложенного и используя разложение р(х + ик + ¿ик) по ¿ик, перепишем (16):

Фк+1 = (р(х + ик + ¿ик) — р(х))(1 + 7к+1(ик+1)) + р(х)7 к+1(ик+1) =

= (p(x + uk) - p(x) +

dp(x + uk) k

dx

duk + o(||duk||))(1 + y2(uk) + Yk(duk) + uk(uk, duk)) +

+ p(x)(Yk(uk) + Yk(duk) + uk(uk, duk)) = (p(x + uk) - p(x))(1 + 7k(uk)) + p(x)7k(uk) + dp(x+v2) k

+

dx

(1 + Yk(uk))duk + p(x+uk)(7k(duk) + uk(uk, duk)) + o(||d

^ др(х + ик) к к к

Отметим, что здесь -------—------= (рЖ1(ж + и ), рж2(ж + и )) = (gradир{х + и )) .

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

фк+i = фк + ,1 + 7^dP(x + uk)duk + / + uk\tfjldiv(duk\ + ш*(иъ duk\\

дх

заметим, что здесь и далее yк = Y(ик).

Теперь вычислим grad Фк+1:

grad Ф2+1 = grad Ф2 +

grad

(1+ Y2)

др(х + ик) дх

du2 + p(x + uk)(e1 div(du2) + u2(u2, du2))

(17)

Выполнив преобразования в формуле (17), будем вычислять градиент по схеме:

grad Фк+1 = grad Ф2 + |(grad + ^ + ^Як|ймк+

+ (1 + Y2 )(pxi(x + uk )grad(duk) + px2(x + uk )grad(du2)) + (^idiv(duk) +

+ u2(u2, du2))grad p(x + u2) + p(x + uk)grad(eidiv(duk) + u2(u2, du2)),

здесь Hк вания.

pkxixj = Рxi,xj(х + ик), i,j = 1, 2, * - знак транспониро-

Подставляя полученные выражения для Фк+1, grad Фк+1 в итерационное уравнение (13), получаем

[(1 + Yk)((1 + Yk)graduрк + 0кgrad рк)(graduРk)* +

+ pk0k(grad Yk(gradupk)* + PkHk)\duk + рк0к((1 + Yk)(pХ1g™d(duk) +

+ pХ2grad(duk)) + (e1div(duk) + wk)grad рк + pkgrad(в1div(duk) + wk)) —

— о?Aduk + Фк[(1 + Yk)gradupk + 0к grad рк] + рк0кgrad Фк — а2Аик = 0, (18)

Таким образом, находя duk, мы определяем uk+1 и повторяем снова итерационный процесс. В итоге получили схему для построения итерационного процесса с целью нахождения поля вектора перемещения.

Заметим, что для определения duk , т. е. решения линейной системы (18), можно также использовать итерационный процесс. Следует отметить, что для решения этой задачи могут быть применены различные итерационные схемы, в частности метод релаксации, который показал себя достаточно хорошо при решении задачи определения оптического потока [7-9]. Наша задача сводится к нахождению оптического потока при @1 = р2 = 0.

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

Литература

1. Беллман Р. Математические методы в медицине/ пер. с англ. А. Л. Асаченкова, Н. А. Шально-вой; под ред. Л. Н. Белых. М.: Мир, 1987. 200 с.

2. Радионуклидная диагностика / под ред. Ф. Н. Лясса. М.: Медицина, 1983. 304 с.

3. Овсянников Д. А. Моделирование и оптимизация динамики пучков заряженных частиц. Л.: Изд-во Ленингр. ун-та, 1990. 312 с.

4. Kotina E.D. Discrete optimization problem in beam dynamics // Nuclear Instruments and Methods in Physics Research. 2006. Vol. A 558. P. 292-294.

5. Horn B. K. P., Schaunck B. G. Determing of Optical Flow // Artificial Intellegence. 1981. Vol. 17. P. 185-203.

6. Papenberg N., Bruhn A., Brox T. et al. Highly Accurate Optic Flow Computation with Theoretically Justified Warping // Intern. J. of Computer Vision. 2006. Vol. 67, N 2. P. 141-158.

7. Young D. M. Iterative Methods for Solving Partial Difference Equations of Elliptic Type: Ph. D. Thesis. Cambridge, Mass.: Harvard University, 1950. 74 p.

8. Young D. M. Iterative Solution of Large Linear Systems. New York: Academic Press, 1971. 563 p.

9. Фаддеев Д. К., Фаддеева В. Н. Вычислительные методы линейной алгебры. М.; Л.: Физматгиз, 1963. 734 с.

Статья рекомендована к печати проф. Д. А. Овсянниковым.

Статья принята к печати 1 апреля 2010 г.