Диссертация (1145332), страница 13

Просмтор этого файла доступен только зарегистрированным пользователям. Но у нас супер быстрая регистрация: достаточно только электронной почты!

Текст из файла (страница 13)

(2.20)0(4)Так как функция elf (r, ω) ∈ Cb (G0 × Ω0 ), а функция τe(r, ω, t) абсолютнонепрерывна по переменной t, то к левой части равенства (2.20) применимаформула интегрирования по частям:85d(r,−ω)Zexp(−eτ (r, ω, t))elf (r − tω, ω)dt = f (r, ω)−0d(r,−ω)Zf (r − d(r, −ω)ω, ω) exp(−eτ (r, ω)) −exp(−eτ (r, ω, t))×0× −µe(r − tω)f (r − tω, ω) + µe(r − tω)f (r − tω, ω) dt.

(2.21)Учитывая граничное условие (2.2), из (2.21) получаемd(r,−ω)Zexp(−eτ (r, ω, t))F (r − tω, ω)dt =0f (r, ω) − h(r − d(r, −ω)ω, ω) exp(−eτ (r, ω)).e и Se непосредственно вытеОтсюда и из свойств интегральных операторов A(4)кает, что функция f (r, ω) из пространства Cb (G0 × Ω0 ) для любых (r, ω) ∈G0 × Ω0 удовлетворяет уравнению (2.18). Необходимость доказана.Достаточность.

Пусть функция f (r, ω) ∈ D(G0 × Ω0 ) удовлетворяетуравнению (2.18). Тогда для любых (ξ, ω) ∈ Γ− из уравнения (2.18), учитывая, что d(ξ + tω, −ω) = t, имеемf (ξ + tω, ω) = exp(−eτ (ξ, −ω, t))×Zth(ξ, ω) + exp(eτ (ξ, −ω, t0 ))F (ξ + t0 ω, ω)dt0  . (2.22)0Так как неопределенный интеграл от интегрируемой по t функции является абсолютно непрерывной функцией, а сумма и произведение абсолютно непрерывных функций есть также абсолютно непрерывная функция, тоиз (2.22) заключаем, что функция f (ξ + ωt, ω) абсолютно непрерывна на{ξ + tω : t ∈ [0, d(ξ, ω)]} для любых (ξ, ω) ∈ Γ− .

Непосредственно убежда-86емся, что при t = 0 функция f (ξ + ωt, ω) удовлетворяет краевому условию(2.2).Далее, подействовав оператором el на f (ξ + tω, ω) и воспользовавшись равенством (2.22), получим при почти всех t ∈ [0, d(ξ, ω)]elf (ξ + tω, ω) = −eµ(ξ + tω) exp(−eτ (ξ, −ω, t))×Zth(ξ, ω) + exp(eτ (ξ, −ω, t0 ))F (ξ + t0 ω, ω)dt0  +0µe(ξ + ωt) exp(−eτ (ξ, −ω, t))×Zth(ξ, ω) + exp(eτ (ξ, −ω, t0 ))F (ξ + t0 ω, ω)dt0  +0exp(−eτ (ξ, −ω, t)) exp(eτ (ξ, −ω, t))F (ξ + tω, ω) = F (ξ + tω, ω).Таким образом, функция f (r, ω) из класса D(G0 ×Ω0 ), являющаяся решениемуравнения (2.18), удовлетворяет всем условиям в определении краевой задачи(2.1), (2.2). Лемма доказана.Обозначимefe0 (r, ω) = h(r − d(r, −ω)ω, ω) exp(−eτ (r, ω)) + (AJ)(r,ω).Теперь докажем окончательное утверждение о корректности прямой задачи(2.1), (2.2).(4)Теорема 2.4.

Пусть eh(r, ω) ∈ K, J(r, ω) ∈ Cb (G0 × Ω) ∩ K и справедливыусловия (2.6), (2.7), тогда в конусе K решение задачи (2.1), (2.2) существует, единственно и представимо в виде ряда Нейманаf (r, ω) = fe0 (r, ω) +∞XeS)e n fe0 (r, ω),(A(2.23)n=1(4)сходящегося в норме пространства Cb (G0 × Ω0 ).Доказательство. По лемме 2.6. достаточно доказать разрешимость уравне-87ния (2.18). Заметим, что конус K является замкнутым подмножеством бана(4)хова пространства Cb (G0 × Ω0 ), следовательно K — полное пространство.eSe переводит K в K, функцияПоэтому достаточно доказать, что оператор AeSfe k < kf k для всех f ∈ K.fe0 (r, ω) ∈ K и kAe : K → K, то AeSe : K → K.Так как оператор Se : K → K и оператор A(4)Поскольку eh(r, ω) ∈ K, τe(r, ω) ∈ Cb (G0 × Ω) и J(r, ω) ∈ K ∩ C (G0 × Ω), тоbfe0 (r, ω) ∈ K.Используя неравенство (2.17), для всех f ∈ K имеемeSfe k4 = max k(AeSfe )i k = k(AeSfe )1 k ≤kA1≤i≤4λ 1 − e−µd kf1 k = λ 1 − e−µd kf k4 .И поскольку µ(r) ≥ µs (r), r ∈ G, то λ ≤ 1 иeSfe k4 ≤ λ 1 − e−µd kf k4 < kf k4 .kAСледовательно, утверждение теоремы вытекает из принципа сжимающих отображений.Отметим, что в случае µe(r) = µ(r) обоснование корректности прямой задачи совпадает с аналогичными рассуждениями, проведенными в §3.1.

Необходимость введения функции µe(r) вызвана потребностью обоснования вычислительного алгоритма решения прямой задачи, рассмотренного далее.Опишем алгоритм вычисления вектор-функции f , основанный на методеМонте-Карло. Пусть для любых r ∈ G выполняется неравенство µ(r) ≤ µ,где µ – некоторая константа. Тогда, полагая µe(r) = µ, по теореме 2.4. получаем представление решения в виде равномерно сходящегося ряда (2.23).Добавку (µ−µ(r))ϕ(r, ω) в выражении для интегрального оператора Se можнотрактовать как некоторое фиктивное рассеяние без изменения направленияраспространения фотона.

Данный метод носит название метода максимального сечения и позволяет более просто разыгрывать длину свободного пробега частицы даже в областях со сложной структурой. Подобный подход даетнеплохие результаты в случае, когда изменение коэффициента полного вза-88имодействия в среде невелико. Пусть m – учитываемое нами число членовряда Неймана, а n – число моделируемых траекторий, тогда приближенноезначение функции f (r, ω) можно вычислить по формуле:n1Xsi (r, ω),f (r, ω) ≈ f n (r, ω) =n i=1si (r, ω) = fe0 (r, ω) +jm YX1 − exp(−µd(ri,k−1 , −ω i,k−1 ))j=1 k=1µ(2.24)××(µ − µ(ri,k ) + µs (ri,k ))Q(ω i,k−1 , ω i,k )fe0 (ri,j , ω i,j ).Для генерации точек траекторий ri,k и направлений ω i,k , необходимых дляреализации алгоритма, используем начальные значения: ri,0 = r, ω i,0 = ω.На каждом последующем шаге вычисление точки ri,k осуществляется поформуле:ri,k = ri,k−1 − ω i,k−1 ti,k ,(2.25)где ti,k – независимая реализация случайной величины, распределенной на[0, d(ri,k−1 , −ω i,k−1 )] с плотностьюρ(t) =µ exp(−µt).1 − exp(µd(ri,k−1 , −ω i,k−1 ))(2.26)Далее разыгрываем реализацию αi,k равномерно распределенной на [0, 1] случайной величины.

Приαi,kµ − µ(ri,k )≥µ − µ(ri,k ) + µs (ri,k )(2.27)для определения величины ω i,k и матрицы Q(ω i,k−1 , ω i,k ) используем следующие расчетные формулы:qω1i,k=2 cos γ ,1 − νi,ki,k(2.28)ω2i,kq2 sin γ ,= 1 − νi,ki,k(2.29)89ω3i,k = νi,k ,(2.30)Q(ω i,k−1 , ω i,k ) = 4πP (ω i,k−1 , ω i,k ),(2.31)а приαi,kµ − µ(ri,k )<µ − µ(ri,k ) + µs (ri,k )следующие:ω i,k = ω i,k−1 ,Q(ω i,k−1 , ω i,k ) = E.(2.32)Здесь νi,k , γi,k есть независимые реализации равномерно распределенных насоответствующих промежутках [−1, 1] и [0, 2π) случайных величин, а E естьединичная матрица 4 × 4.Таким образом, алгоритм нахождения решения прямой задачи выглядитследующим образом.0. Пусть i = 1.1.

Берем i-ую траекторию и полагаем k = 1.2. Для k-ого звена в траектории разыгрываем длину свободного пробегаti,k с плотностью вероятности, задаваемой формулой (2.26).3. Вычисляем координаты новой точки столкновения ri,k по формуле (2.25)4. Разыгрываем реализацию αi,k равномерно распределенной на [0,1] случайной величины.5. Если выполняется условие (2.27), то новое направление ω i,k находитсяиз (2.28)-(2.30), а матрица Q(ω i,k−1 , ω i,k ) по формуле (2.31). Если жеусловие не выполняется, то согласно (2.32) направление не меняется, а в качестве матрицы Q(ω i,k−1 , ω i,k ) берется единичная.6.

Формируем слагаемые ряда Неймана и накапливаем сумму в (2.24).7. Если учтены не все члены ряда Неймана (при k < m), то переходим кпункту 2, увеличивая k на единицу.8. Если учтены не все траектории (при i < n), то переходим к пункту 1,увеличивая переменную i на единицу.Приведенный здесь алгоритм решения прямой задачи (2.1), (2.2) был программно реализован и далее использовался для нахождения выходящего излучения, необходимого для решения задачи компьютерной томографии.904.2.Численное решение задачи томографииРешение задачи томографии будет находиться с помощью алгоритма, основанного на теореме 2.2., согласно которой задача сводится к задаче обращения двумерного преобразования Радона от функции µd(r,ω)Ze(Rµ)(r,ω) :=µ(r + ωt)dt = gei (r, ω),(2.33)−d(r,−ω)гдеgei (r, ω) = ln[hi (r − d(r, −ω)ω, ω)][Hi (r + d(r, ω)ω, ω)](2.34)в любой горизонтальной плоскости {r = (r1 , r2 , r3 ) ∈ R3 : r3 = const}, имеющей общую точку с множеством G0 .

Эта задача имеет единственное решениев широком классе функций [88, 148].Из формул (2.33)-(2.34) непосредственно вытекает, что для нахожденияe(Rµ)(r,ω) можно использовать любые компоненты векторных функций h иH, имеющие ненулевые разрывы по угловой переменной. Это обстоятельствобудет использовано при проведении численных экспериментов, тестирующихалгоритм решения задачи томографии.Продемонстрируем работу алгоритма решения задачи томографии на следующем примере.6SПусть G0 =Gi , гдеi=1G1 = {r : (r1 − 0.3)2 + (r2 − 0.1)2 + r32 < 0.04},G2 = {r : (r1 + 0.3)2 + (r2 + 0.1)2 + r32 < 0.01},G3 = {r : (r1 − 0.3)2 + (r2 − 0.1)2 + r32 < 0.09}\G1 ,G4 = {r : r12 + r22 + r32 < 0.64}\{G1 ∪ G2 ∪ G3 },G5 = {r : 0.64 < r12 + r22 + r32 < 0.81},G6 = {r : 0.81 < r12 + r22 + r32 < 1},91и функция µ(r) является кусочно постоянной в области G, то есть µ(r) = µiпри r ∈ Gi .

В численном эксперименте были взяты следующие значениякоэффициентов: µ1 = µ4 = µ6 = 1, µ2 = µ3 = µ5 = 2, при 50% рассеяниив среде (µs (r) = 0.5µ(r)) и отсутствии внутренних источников излучения(J(r, ω) ≡ 0).Матрица P = P (ω, ω 0 ) определяется следующим образом [62]:P (ω, ω 0 ) = L(π − ψ)R(ω, ω 0 )L(−ψ 0 ),(2.35)где ψ – угол между плоскостью рассеяния и плоскостью, образованной векторами ω, e3 , а ψ 0 – угол между плоскостью рассеяния и плоскостью, образованной векторами ω 0 , e3 . Матрица поворота L(ψ) определяется следующимобразом:10000 cos 2ψ sin 2ψ 0L(ψ) = 0 − sin 2ψ cos 2ψ 0 .0001(2.36)Будем полагать, что рассеяние происходит по рэлеевскому закону, которыйописывается следующей угловой матрицей [62]:(ω · ω 0 )2 + 1 (ω · ω 0 )2 − 10(ω · ω 0 )2 − 1 (ω · ω 0 )2 + 1030R(ω, ω ) =16π 002ω · ω 000000.02ω · ω 0(2.37)Компоненты матрицы поворота определяются на основе следующих представлений:√√ν 0 1 − ν 2 − ν 1 − ν 02 cos(γ − γ 0 )pcos ψ =,1 − (ω · ω 0 )2psin ψ = 1 − cos2 ψ · sgn{sin(γ − γ 0 )}, (2.38)92√0√ν 1 − ν 02 − ν 1 − ν 2 cos(γ − γ 0 )p,cos ψ =1 − (ω · ω 0 )2psin ψ 0 = 1 − cos2 ψ 0 · sgn{sin(γ − γ 0 )}, (2.39)0где ν = cos θ, ν 0 = cos θ0 , а θ, γ и θ0 , γ 0 представляют собой сферические углыдля векторов ω и ω 0 :ω1 = cos γ sin θ, ω10 = cos γ 0 sin θ0 , ω2 = sin γ sin θ, ω20 = sin γ 0 sin θ0 ,ω3 = cos θ, ω30 = cos θ0 , γ, γ 0 ∈ [0, 2π),θ, θ0 ∈ [0, π].Были взяты следующие компоненты вектор-функции h(r, ω), соответствующей входящему излучению:(h1 =1,(ω3 ≥ 0,h2 =1.2, ω3 < 0,0.2, ω3 ≥ 0,1,ω3 < 0,(h3 =0.4, ω3 ≥ 0,0.1, ω3 < 0,h4 = 0.Восстановление функции µ(r) осуществлялось в плоскости r3 = 0.

Пересечение плоскости r3 = 0 и области G представляет собой круг радиуса 1.При восстановлении функции µ(r) использовалась параллельная схема сканирования [20]. Пустьr1 = ρ cos γ,r2 = ρ sin γ,ω = (− sin γ, cos γ, 0),ρ ∈ [−1, 1],γ ∈ [0, 2π),ω⊥ = (cos γ, sin γ, 0),ω · ω⊥ = 0,тогда равенство (2.33) можно переписать в виде√ 2Z1−ρµ(ρω⊥ + tω)dt = gi (ρ, γ),−√1−ρ2где gi (ρ, γ) = gei (r(ρ, γ), ω(γ)), r = ρω⊥ (γ), ω = (− sin γ, cos γ, 0). Таким образом, в сечении r3 = 0 мы получаем интегралы от следа функции µ по почти93всем прямым, проходящим через точки r = ρω⊥ (γ) в направлениях ω(γ),ρ ∈ [−1, 1], γ ∈ [0, 2π). То есть, задача определения функции µ(r) сводится кзадаче обращения преобразования Радона (Rµ)(ρ, γ) от функции µ(r).Для проведения вычислительных экспериментов было взято следующееразбиение множества [−1, 1] × [0, 2π), на котором определяется преобразование Радона (Rµ)(ρ, γ):ρl = −1 + l/60,Для нахожденияl = 0, 120,gi (ρl , γs )γs = sπ/90,s = 0, 179.необходимо вычислить величину скачка[H(r + d(r, ω)ω, ω)] в точках (rl,s , ω s ), где rl,s = ρl ω⊥ (γs ), ω s = ω(γs ).Вычисление функции H проводилось с помощью метода Монте-Карло.Вектор-функция аппроксимировалась суммой первых 21-го слагаемых рядаНеймана (2.23), который сходится как геометрическая прогрессия со знаменателем q = sup(µ−µ(r)+µs (r))/µ = 3/4.

Характеристики

Список файлов диссертации