Кристаллография, 2019, T. 64, № 2, стр. 173-183

ВВЕДЕНИЕ

Как известно [1–6], дифракционная рентгеновская топография (ДРТ) является высокочувствительным неразрушающим инструментом диагностики различных дефектов кристаллической решетки, таких как полосы роста, границы зерен, дефекты упаковки, отдельные включения и дислокации различного типа. Все эти дефекты приводят к изменению положений отдельных атомов кристаллической структуры относительно их правильного положения. При этом в схемах лауэвской и/или брэгговской дифракции реальная структура кристалла изучается по его 2D-проекциям. Как правило, для интерпретации и анализа экспериментальных ДРТ-изображений дефектов на 2D-проекциях они сопоставляются с 2D-проекциями, рассчитанными на основе компьютерного моделирования изображений дефектов с использованием численных методов решения динамических уравнений Такаги–Топена [1–3, 6–8].

В последние 20 лет в структурном анализе реальных кристаллов широко применяется метод дифракционной рентгеновской топо-томографии (ДРТТ). В методе ДРТТ кристалл-образец поворачивается вокруг оси, перпендикулярной семейству отражающих плоскостей кристалла (ось вращения выбирается вдоль вектора дифракции $h$). Поэтому при различных значениях соответствующего угла поворота получается набор дифракционных 2D-проекций, каждая из которых отвечает разным ориентациям плоскости дифракционного рассеяния рентгеновского излучения по отношению к собственной системе координат кристалла-образца (рис. 1).

Первые ДРТТ-эксперименты были успешно выполнены на синхротронных источниках рентгеновского излучения ESRF [9] и SPring-8 [10], а затем и в России на лабораторных установках с использованием характеристического рентгеновского излучения [11, 12]. При этом для компьютерной 3D-реконструкции дефекта в кристалле по соответствующему набору 2D-проекций ДРТ использовались различные модификации алгоритмов поглощающей томографии [13].

Особый интерес, напрямую связанный с применением метода ДРТТ, представляет идея компьютерной реконструкции пространственного положения дефектов в кристалле, а также, что особенно важно, локальных полей смещений (напряжений) вокруг дефектов в кристаллических веществах.

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

В этой связи представляются важными нахождение и использование аппроксимационных аналитических решений уравнений Такаги–Топена, которые позволили бы с достаточной точностью описать особенности контраста дефектов на ДРТ-топограммах, связанные с формированием контраста в той или иной области кристалла вблизи дефекта, что стало бы ключевым моментом для решения обратной задачи ДРТТ, в частности компьютерной 3D-реконструкции функции поля смещений вокруг дефектов.

Настоящая работа посвящена построению аппроксимационного аналитического решения уравнений Такаги–Топена, на базе которого строится решение обратной задачи ДРТТ.

Затем с использованием полукинематического приближения для амплитуды дифрагированной σ-поляризованной волны E_h(r) построена теория 3D-реконструкции функции поля упругих статических смещений u(r – r₀) точечного дефекта кулоновского типа (r₀ – радиус-вектор положения точечного дефекта в кристалле).

На примере одной 2D-проекции с вектором дифракции h = [$\bar {2}20$] для кристалла Si(111) с применением итерационных алгоритмов моделирования отжига [14] и квазиньютоновского спуска [15, 16] найдено решение обратной задачи ДРТТ, а именно, проведена компьютерная реконструкция трехмерной функции поля смещений точечного дефекта, $f({\mathbf{r}} - {{{\mathbf{r}}}_{0}}) = {\mathbf{h}} \cdot {\mathbf{u}}({\mathbf{r}} - {{{\mathbf{r}}}_{0}})$.

Отметим, что первая попытка трехмерной реконструкции статического поля точечного дефекта в кристалле была предпринята в [17], где использовался алгоритм модифицированного алгебраического метода (алгоритм SART [11]), широко применяемого в поглощающей рентгеновской томографии. К сожалению, оказалось, что трехмерные решения для функции поля смещений точечного дефекта, полученные с использованием алгоритма SART, сходятся только при ограниченном количестве разбиений области кристалла вокруг точечного дефекта (не более чем 5 × 5 × 5 = 125 вокселей).

2D-ИЗОБРАЖЕНИЕ ТОЧЕЧНОГО ДЕФЕКТА В КРИСТАЛЛЕ В УСЛОВИЯХ НАКЛОННОЙ СИММЕТРИЧНОЙ ЛАУЭ-ДИФРАКЦИИ. ПОЛУКИНЕМАТИЧЕСКОЕ ПРИБЛИЖЕНИЕ

Целью данного раздела является построение решения динамических уравнений Такаги–Топена в полукинематическом приближении, которое будет использовано в качестве базового для трехмерной реконструкции поля статических смещений по наклонным 2D-проекциям ДРТ.

Как известно [1–3], и это подтверждается неоднократными численными расчетами [16], прямое изображение дефекта в кристалле, если оно есть, обусловлено межветвевым рассеянием блоховских рентгеновских волн в сильно искаженной области кристалла непосредственно вблизи дефекта, что можно интерпретировать как кинематическое рассеяние блоховских волн, распространяющихся в совершенном кристалле вдали от центральной области дефекта.

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

В рассматриваемом случае распространение рентгеновского волнового поля внутри искаженного кристалла в условиях симметричной двух-волновой лауэ-дифракции описывается уравнениями Такаги–Топена [7, 8]:

вместе с граничными условиями на входной поверхности кристалла z = 0 в виде

Граничные условия (2) можно непосредственно учесть, модифицируя уравнения (1), в частности добавляя в правую часть первого уравнения дополнительный член $ - \frac{{2i}}{k}{\delta }({{s}_{0}} + {{s}_{h}})$.

Ипользуя известные подстановки ${{E}_{0}}\, \to \,{{E}_{0}}{{{\text{e}}}^{{ik{{{\chi }}_{0}}\frac{{{{s}_{0}} + {{s}_{h}}}}{2}}}}$, ${{E}_{h}} \to {{E}_{h}}{{{\text{e}}}^{{ik{{{\chi }}_{0}}\frac{{{{s}_{0}} + {{s}_{h}}}}{2}}}}$ и сохраняя для амплитуд проходящей и дифрагированной волн $\quad{{E}_{0}}({{s}_{0}},{{s}_{h}})$, ${{E}_{h}}({{s}_{0}},{{s}_{h}})$ прежние обозначения, легко показать, что амплитуда дифрагированной волны ${{E}_{h}}({{s}_{0}},{{s}_{h}})$ удовлетворяет неоднородному дифференциальному уравнению гиперболического типа в частных производных второго порядка

Здесь ${{E}_{0}}({{s}_{0}},{{s}_{h}})$ и $\quad{{E}_{h}}({{s}_{0}},{{s}_{h}})$ – амплитуды проходящей и дифрагированной волн в среде, зависящие от координат ${{s}_{0}},\;{{s}_{{h\quad}}}$ в плоскости рассеяния; $\frac{\partial }{{\partial {{s}_{0}}}}$ и $\frac{\partial }{{\partial {{s}_{h}}}}$ – производные вдоль направлений проходящей и дифрагированной волн, удовлетворяющих точному условию Брэгга; θ_B – угол Брэгга; ${{{\chi }}_{0}}$, ${{{\chi }}_{h}}$, ${{{\chi }}_{{\bar {h}}}}$ – нулевая, h и $\bar {h}$ компоненты Фурье диэлектрической проницаемости совершенного кристалла соответственно; k = 2π/λ, λ – длина волны падающего излучения; ${\mathbf{u}}({\mathbf{r}})$ – вектор поля смещений, описывающий искажение идеальной кристаллической структуры; ${{{\mathbf{k}}}_{0}}$ – волновой вектор плоской волны, падающей на кристалл-образец в точном брэгговском положении; ${{{\mathbf{s}}}_{0}} = \frac{{{{{\mathbf{k}}}_{0}}}}{{{{k}_{0}}}}$, ${{{\mathbf{s}}}_{h}} = \frac{{{{{\mathbf{k}}}_{h}}}}{{{{k}_{h}}}}$; $h$ – вектор дифракции.

Интересует поле упругих статических смещений точечного дефекта кулоновского типа u(r – r₀), которое имеет следующий вид:

где r₀ – радиус-вектор положения точечного дефекта.

В качестве модельного кристалла выбрана плоскопараллельная пластина Si(111), вектор дифракции h = [$\bar {2}20$], длина волны падающего рентгеновского характеристического излучения λ = = 0.071 нм. Использовалась лабораторная рентгеновская трубка с молибденовым анодом, энергия излучения 17.48 КэВ; Λ – длина экстинкции, равная 36.287 ммк, угол Брэгга θ_B = 10.65°.

Отметим, что второй член в правой части первого из двух уравнений Такаги–Топена отвечает плоской волне ${{E}_{0}}({\mathbf{r}}) = {{{\text{e}}}^{{i{{{\mathbf{k}}}_{0}}{\mathbf{r}}}}}$, падающей на входную поверхность кристалла z = 0.

На рис. 2 показана дискретная треугольная сетка, построенная на векторах ${{{\mathbf{s}}}_{0}}$, ${{{\mathbf{s}}}_{h}}$ с одинаковым шагом $p$, которая может быть использована для компьютерного моделирования и анализа наклонных 2D-проекций ДРТ на основе численного решения уравнений Такаги–Топена (1) с граничными условиями (2) в переменных S, X [11].

Отметим, что S, X – координаты косоугольной системы координат SXY (ось Y перпендикулярна плоскости (SX)), связаны с координатами s₀, s_h линейными соотношениями

Уравнение (3) можно записать в интегральной форме, вводя в рассмотрение функцию Грина, описывающую распространение дифрагированного излучения в идеальном кристалле от точечного источника в треугольной области (треугольник Бормана) с вершиной в точке характеристиками вдоль направлений ${{{\mathbf{s}}}_{0}} = \frac{{{{{\mathbf{k}}}_{0}}}}{k}$, ${{{\mathbf{s}}}_{h}} = \frac{{{{{\mathbf{k}}}_{h}}}}{k}$ ${\text{(}}k = {{k}_{0}} = {\text{ }}{{k}_{h}})$, а именно:

где J₀(u) – функция Бесселя нулевого порядка действительного аргумента u.

Дальнейшее рассмотрение компьютерной реконструкции 3D-поля статических смещений точечного дефекта (4) основывается на утверждении, что в области сильных искажений вблизи дефекта рассеяние носит кинематический характер, что обусловлено межветвевым рассеянием блоховских волн вблизи дефекта [1–6]. При этом волновое поле, формирующее так называемое “прямое” изображение дефекта на 2D-проекциях ДРТ, распространяется вдоль направления поля дифрагированной волны ${\text{ }}{{{\mathbf{s}}}_{h}} = \frac{{{{{\mathbf{k}}}_{h}}}}{k}$.

Математически это означает, что для описания прямого ДРТ-изображения дефекта соответствующее выражение для амплитуды дифрагированной волны может быть найдено в результате интегрирования по частям второго члена в квадратных скобках в правой части интегрального уравнения (5).

Кроме того, ограничиваясь только первым членом двойного интеграла в правой части (5), находим

Далее, следуя предположению, что прямое изображение дефекта на 2D-проекциях ДРТ отвечает кинематическому рассеянию проходящей волны в области сильных искажений вблизи дефекта, в правой части интегрального уравнения (6) амплитуду проходящей волны в первом приближении можно заменить на соответствующее выражение для идеального кристалла, т.е. ${{E}_{0}}({{s}_{0}},{\text{ }}s_{0}^{'}) \to E_{0}^{{(id)}}({{s}_{0}},{\text{ }}s_{h}^{'}),$ где $E_{0}^{{(id)}}({{s}_{0}},s_{h}^{'})$ есть не что иное, как решение для амплитуды проходящей волны в идеальном кристалле соответственно: Таким образом, для компьютерной 3D-реконструкции поля статических искажений вблизи точечного дефекта (4) для амплитуды дифрагированной волны будем использовать выражение

Формула (7) для амплитуды $\quad{{E}_{h}}({{s}_{0}},{{s}_{h}})$ описывает кинематическое рассеяние проходящей волны $E_{0}^{{\left( {i{\text{d}}} \right)}}({{s}_{0}},s_{h}^{'})$ сильно искаженной области вблизи дефекта, которая представляет собой фазовый объект. При этом волна ${{E}_{h}}({{s}_{0}},{{s}_{h}})$ распространяется вдоль направления дифрагированной волны ${{{\mathbf{s}}}_{h}} = \frac{{{{{\mathbf{k}}}_{h}}}}{k}.$ Таким образом, формула (7), полученная в полукинематическом приближении, будет использоваться для описания 2D-проекций ДРТ в той части, которая формируется за счет рассеяния рентгеновской волны в сильно искаженной области кристалла вблизи дефекта.

Для проведения численных расчетов выберем начало системы координат ${{s}_{0}}$, ${{s}_{h}}$ в левой верхней вершине трапеции (рис. 2) и введем косоугольную аффинную систему координат XYS (ср. с рис. 1), в которой ось Х направлена вдоль вектора дифракции h, а ось S вдоль вектора ${{s}_{h}}$, ось Y перпендикулярна плоскости XS.

При этом связь между координатами ${{s}_{0}}$, ${{s}_{h}}$ и X, S задается соотношениями

Кроме соотношений (8) потребуются соотношения между координатами ${{s}_{0}}$, ${{s}_{h}}$ и x, Z, а именно:

Имея в виду компьютерную 3D-реконструкцию функции упругих статических смещений точечного дефекта $f({\mathbf{r}} - {{{\mathbf{r}}}_{0}})$ по данным наблюдаемых 2D-проекций, ${{I}_{h}}({{s}_{0}},{{s}_{h}}) = {{\left| {{{E}_{h}}({{s}_{0}},{{s}_{h}})} \right|}^{2}},{\;}$ запишем формулу (7) в переменных X, y, z (ось вращения проходим через точку r₀ = (x₀, y₀, z₀)) где Ф – угол вращения образца вокруг вектора дифракции h, функция $f({\mathbf{r}} - {{{\mathbf{r}}}_{0}})$ определяется следующим выражением (ср. (4)):

Отметим, что координаты x, y, z “жестко” связаны с декартовыми координатами в кристалле.

В формуле (11) в явной форме учтена зависимость функции поля смещений точечного дефекта (4) в переменных X, y, z от угловых параметров θ_B, Φ.

Теоретическая формула (10) для амплитуды дифрагированной волны ${{E}_{h}}(X,y,z;{{{\theta }}_{{\text{B}}}},{\Phi })$ наряду с формулой (11) для функции поля смещений точечного дефекта $f({\mathbf{r}} - {{{\mathbf{r}}}_{0}})$ в области сильных искажений является основным выражением для решения обратной задачи ДРТТ.

Имея в виду компьютерную 3D-реконструкцию функции $f({\mathbf{r}} - {{{\mathbf{r}}}_{0}})$ на основе минимизации целевой χ²-функции вида

будем использовать итерационные алгоритмы моделирования отжига [14] и квазиньютоновского спуска [15, 16], адаптированные применительно к решению обратной задачи ДРТТ.

В (12) $\quad{{I}_{{h,obs}}}\{ X,y,z;{{{\theta }}_{{\text{B}}}},{\Phi }\} $ есть модельная (наблюдаемая) ДРТ-проекция, отвечающая значению угла наклона Ф, а ${{I}_{{h,calc}}}\{ X,y,z;{{{\theta }}_{{\text{B}}}},{\Phi }\} $ рассчитывается по формуле (10) в соответствии с выбранным итерационным алгоритмом, начиная с некоторой стартовой функции поля смещений ${{f}_{{in}}}({\mathbf{r}} - {{{\mathbf{r}}}_{0}})$.

КОМПЬЮТЕРНАЯ 3D-РЕКОНСТРУКЦИЯ ФУНКЦИИ СМЕЩЕНИЙ ТОЧЕЧНОГО ДЕФЕКТА НА ОСНОВЕ ПРИМЕНЕНИЯ МЕТОДОВ СИМУЛИРОВАННОГО ОТЖИГА И КВАЗИНЬЮТОНОВСКОГО СПУСКА

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

Стартуя с заданной в качестве начального приближения модели и варьируя ее параметры псевдослучайным образом, программа SA работает до достижения наилучшей сходимости расчетной модели с наблюдаемыми данными ДРТ-проекции дефекта (минимум целевой χ²-функции (12)).

Будем различать два состояния системы – функции поля смещений f(r – r₀): текущее и пробное значения функции $f(r\quad--\quad{{r}_{0}})$. В итерационном процессе текущее состояние, становясь пробным, служит в дальнейшем в качестве нового текущего состояния системы.

Псевдослучайное изменение системы применяется к текущему состоянию, которое затем становится пробной функцией $f(r - {{r}_{0}})$. Вероятность принятия пробной функции в качестве новой текущей определяется значением больцмановской функции exp(–Δχ²/T), где Δχ² – изменение целевой функции, T – параметр, называемый температурой отжига. Выбор нового текущего состояния системы существенно зависит от значения T. Легко видеть, что более высокая температура означает более высокую вероятность exp(–Δχ²/T)), это позволяет принять пробную модель в качестве новой текущей даже в случае худшего пробного приближения в сравнении с текущим, когда Δχ² > 0. Но если Δχ² < 0, то пробная модель всегда принимается в качестве новой текущей. В начале итерационной минимизации целевой χ²-функции температура T выбирается достаточно высокой, так что частота принятия худшего состояния системы в качестве текущего в 10–100 раз превышает частоту обновления лучших состояний системы (другими словами, изменения системы носят почти случайный характер). Это заставляет программу SA “блуждать” по пространству поиска минимальных значений функции системы f. В принципе, программа SA может преодолевать локальные минимумы целевой χ²-функции, которые являются одним из главных препятствий для многих других нелинейных методов минимизации целевой χ²-функции, например, для градиентных методов спуска [15, 16].

В итерационном процессе минимизации целевой χ²-функции температура T_n + 1 на следующем после предыдущего n-шага (внешний цикл по n) снижается монотонно как T_n + 1 = FT_n, где коэффициент отжига F равен 0.9–0.95, а значение целевой χ²-функции уменьшается аналогично тому, как снижается внутренняя энергия системы в процессе понижения температуры (так называемая машина Больцмана). По этой причине программа-алгоритм SA получила название “симулированного отжига” (алгоритм Метрополиса [18]). В табл. 1 представлены детали протокола итерационной программы SA согласно [14, 18].

Для простоты все последующие результаты расчетов приводятся в безразмерных координатах. Положение точечного дефекта задается радиус-вектором r₀, r₀ = nt/2, n – внутренняя нормаль к входной поверхности кристалла z = 0. Толщина t модельного кристалла Si(111) выбрана такой, что поглощением рентгеновского излучения в образце можно пренебречь.

На первом этапе решения обратной задачи ДРТТ программа SA была протестирована на примере минимизации целевой χ²-функции в случае одной 2D-проекции, когда в формуле (12) остается одно слагаемое с Φ = 0, а именно:

где ${{I}_{{h,{\;}obs}}}\{ X,y,z;{{{\theta }}_{{\text{B}}}},0\} $ наблюдаемая (модельная) 2D-проекция, рассчитанная на основе формулы (10) с учетом теоретической функции смещений (11) при значении масштабирующего коэффициента G = 1, а интенсивность рассчитывается по формуле (10) с пробными функциями смещений в процессе итерационной минимизации целевой χ²-функции (13).

Отметим, что масштабирующий коэффициент G в формуле (11) можно выбрать равным единице за счет соответствующего выбора шага дискретной пространственной сетки координат в кристалле, на которой задается функция поля смещений $f(r - {{r}_{0}})$.

Кроме того, преследуя цель оценить степень сходимости итерационного процесса минимизации целевой χ²-функции к правильной (теоретической) функции f_obs(r – r₀) (см. (11)), будем использовать контрольный параметр (CP), который определяется как

и представляет собой среднестатистическую оценку относительного отклонения текущего решения от истинного решения.

Для случая наблюдаемой 2D-проекции ДРТ (рис. 3) результаты компьютерной 3D-реконструкции функции смещений $f({\mathbf{r}} - {{{\mathbf{r}}}_{0}})$ для различных пространственных сеток, на которых она определена, а также ее различных линейных комбинаций в качестве стартовых функций с индексами убывания {p_i} в интервале значений i = 1–4 представлены в табл. 2.

Диапазон поиска каждого индекса убывания {p_i}, i = 1–4, в части действия алгоритма SA (внутренний цикл, табл. 1) задавался с равной вероятностью в интервале чисел 0.0–3.0.

Как видно из табл. 2, в случае одного индекса убывания {p_i}, i = 1, в качестве стартового значения теоретическое решение для функции смещений $f({\mathbf{r}} - {{{\mathbf{r}}}_{0}})$ с индексом убывания p = 1.5 достигается с точностью, по крайней мере CP = 10^–6, в то время как для большего, чем единица, числа индексов убывания {p_i}, i = 1–4, теоретическое решение достигается с точностью порядка CP = 10^–2.

На практике представляет интерес 3D-реконструкция функции поля смещений $f({\mathbf{r}} - {{{\mathbf{r}}}_{0}})$ деффекта в кристалле без необходимости его описания в аналитической форме. В отличие от случая, описанного выше, будем рассматривать функцию $f({\mathbf{r}} - {{{\mathbf{r}}}_{0}})$ в каждом узле дискретной пространственной сетки как искомый параметр. При этом будут использованы только свойства симметрии функции $f({\mathbf{r}} - {{{\mathbf{r}}}_{0}})$ по координатам x – x₀, y – y₀, z – z₀, а также “включено” требование монотонного убывания с увеличением расстояний |y – y₀| и/или |z – z₀|.

Наблюдаемая 2D-проекция ДРТ от тонкого (непоглощающего) кристалла Si(111) (рис. 3) рассчитана на основе формулы (10) с учетом (11) на пространственной сетке, состоящей из 15 × 15 × × 15 узлов. Видно, что проекция является симметричной относительно координаты y и в направлении распространения дифрагированной волны сдвигается от центра вдоль координаты x как целое на величину t/2 × tg θ_B (в данном случае этот сдвиг равен 4/3 в безразмерных координатах x и z).

Входные данные и результаты компьютерной 3D-реконструкции функции $f({\mathbf{r}} - {{{\mathbf{r}}}_{0}})$ по данным наблюдаемой 2D-проекции ДРТ представлены в табл. 3. В первом столбце указаны пространственные сетки, на которых проводятся расчеты, жирным шрифтом (здесь и далее) выделены номера плоскостей вдоль оси z, в которых значения функции $f({\mathbf{r}} - {{{\mathbf{r}}}_{0}})$ варьируются, в то время как в остальных плоскостях они фиксированы в соответствии с теоретическим индексом убывания p = = 1.5.

Из табл. 3 видно, что в случае пространственной сетки {15, 15, 15} удается снизить значение целевой χ²-функции более чем на 6 порядков, в то время как контрольный параметр СР растет до 0.118 против стартового значения, равного 0.0658.

Вероятно, это может происходить по двум причинам, первая из которых – неоднозначность решения обратной задачи ДРТТ. Такой вывод следует, в частности, из того факта, что в случае перехода к пространственным сеткам с меньшим числом узлов по толщине кристалла контрольный параметр СР уменьшается.

Вторая причина заключается в эффективности работы самого алгоритма SA для минимизации целевой χ²-функции. Расчеты показывают, что, если вместо минимизации целевой χ²-функции (13) применяется алгоритм SA для минимизации CP-функции, согласно (14), значение параметра CP уменьшается только примерно в 2 раза, например, для пространственной сетки {15, 15, 15} параметр CP изменяется от 0.0658 до 0.0371, в то время как его значение должно было бы снизиться до нуля в предположении, что алгоритм SA работает эффективно.

Моделирование на основе алгоритма квазиньютоновского спуска. Для сравнения с результатами, полученными с использованием алгоритма SA, был применен для минимизации целевой χ²-функции (13) алгоритм квазиньютоновского спуска по схеме Левенберга–Маркварда. Был использован программный код NL2SNO, имеющийся в открытом доступе (подробнее в [15, 16, 19], вариант NL2SOL с расчетом градиентов по методу конечных разностей).

Результаты минимизации целевой χ²-функции (13) c использованием алгоритма квазиньютоновского спуска для различных пространственных сеток и различных стартовых значений индекса убывания р_ini функции $f({\mathbf{r}} - {{{\mathbf{r}}}_{0}})$ представлены в табл. 4. Обозначения для идентификации пространственных сеток приняты аналогично тому, как это сделано в табл. 3. Видно, что для двух пространственных сеток {15, 15, 1–6 | 7–9 | 10–15} и {15, 15, 1–5 | 6–10 | 11–15} удается достичь значений контрольного параметра CP ∼ 5 × 10^–5 при конечных значениях целевой χ²-функции порядка 10^–22.

Отметим, что в случае пространственной сетки {15, 15, 15} и различных стартовых значениях индекса убывания р_ini контрольный параметр CP практически не уменьшается по сравнению со стартовым CP, а в некоторых случаях даже растет.

Как показывают расчеты, последнее обстоятельство можно преодолеть, если последовательно использовать итерационную схему, корректируя при каждой итерации значения функции $f({\mathbf{r}} - {{{\mathbf{r}}}_{0}})$ в периферийной области от центра точечного дефекта. Детали расчетов с использованием корректировки значений функции $f({\mathbf{r}} - {{{\mathbf{r}}}_{0}})$ периферийных областях пространственной сетки выходят за рамки данной работы и являются отдельной темой будущего исследования.

На рис. 4 приведены сечения трехмерной функции $f({\mathbf{r}} - {{{\mathbf{r}}}_{0}})$ в плоскостях z = const для значений z = 6, 8, 10 (z = 8 – центральное сечение).

Верхние рисунки 4а, 4б, 4в – сечения теоретической 3D-функции $f({\mathbf{r}} - {{{\mathbf{r}}}_{0}})$, индекс убывания р = 1.5, нижние – сечения 3D-функции $f({\mathbf{r}} - {{{\mathbf{r}}}_{0}})$ как результат компьютерной 3D-реконструкции с использованием метода квазиньютоновского спуска. Наилучшее согласие соответствующих изображений имеет место для центрального сечения z = 8. В принципе, это может означать, что наилучшим образом удается восстановить 3D-функцию $f({\mathbf{r}} - {{{\mathbf{r}}}_{0}})$ в области непосредственно вблизи центра дефекта.

ЗАКЛЮЧЕНИЕ

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

Для решения обратной задачи ДРТТ были использованы итерационные алгоритмы моделирования отжига и квазиньютоновского спуска, приведены результаты компьютерной 3D-реконструкции функции поля смещений точечного дефекта $f({\mathbf{r}} - {{{\mathbf{r}}}_{0}})$. Показано, что использованные в работе алгоритмы при определенных ограничениях на класс функций, на котором ищется функция $f({\mathbf{r}} - {{{\mathbf{r}}}_{0}})$, работают в случае одной 2D-проекции ДРТ, что позволяет восстановить теоретическую 3D-функцию поля смещений $f({\mathbf{r}} - {{{\mathbf{r}}}_{0}})$, заданную как в аналитической форме, так и в численном виде в узлах пространственной сетки кристаллической пластины.

Проведенные предварительные расчеты показали, что алгоритм квазиньютоновского спуска допускает определенную возможность улучшить работу итерационной схемы, корректируя при каждой итерации значения 3D-функции $f({\mathbf{r}} - {{{\mathbf{r}}}_{0}})$ в узлах пространственной сетки в периферийной области на некотором расстоянии от центра точечного дефекта.

Представляется весьма вероятным, что сходимость процесса минимизации целевой χ²-функции с использованием алгоритма квазиньютоновского спуска также должна улучшиться в случае использования набора наблюдаемых наклонных 2D-проекций (общая формула (12)), что является немаловажным обстоятельством для фильтрации шумов и практической обработки 2D-проекций в условиях ДРТТ.

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

Работа выполнена при поддержке Министерства науки и высшего образования в рамках выполнения работ по Государственному заданию ФНИЦ “Кристаллография и фотоника” РАН в области развития методов исследования структуры с помощью рентгеновского и синхротронного излучений.