Доклады Российской академии наук. Математика, информатика, процессы управления, 2020, T. 494, № 1, стр. 17-20

ДВУХЭТАПНЫЙ МЕТОД РЕШЕНИЯ СИСТЕМ НЕЛИНЕЙНЫХ УРАВНЕНИЙ И ЕГО ПРИЛОЖЕНИЕ К ОБРАТНОЙ ЗАДАЧЕ ЗОНДИРОВАНИЯ АТМОСФЕРЫ

Член-корреспондент РАН В. В. Васин 12*, Г. Г. Скорик 12**

1 Институт математики и механики им. Н.Н. Красовского
Екатеринбург, Россия

2 Уральский федеральный университет им. Б.Н. Ельцина
Екатеринбург, Россия

* E-mail: vasin@imm.uran.ru
** E-mail: skorik@imm.uran.ru

Поступила в редакцию 04.06.2020
После доработки 04.06.2020
Принята к публикации 16.07.2020

Полный текст (PDF)

Аннотация

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

Ключевые слова: нелинейная система, регуляризация, итерационный процесс, инфракрасный спектр, зондирование атмосферы, восстановление HDO

ВВЕДЕНИЕ

Обратные задачи теплового зондирования атмосферы, связанные с определением концентрации изотопомеров парниковых газов и водяного пара по спектрам инфракрасного излучения, измеренным сенсорами спутникового или наземного базирования, как правило, сводятся к решению системы нелинейных уравнений

(1)
$P(z) = f,\quad z \in {{\mathbb{R}}^{n}},\quad f \in {{\mathbb{R}}^{m}}.$

При этом размерность m существенно больше, чем n, что приводит к переопределенности системы (1). В работе [1] был предложен метод решения системы нелинейных уравнений (1), в котором можно выделить следующие этапы построения приближенного решения:

тихоновская регуляризация системы (1)

(2)
$inf\left\{ {{{{\left\| {P(z) - {{f}_{\delta }}} \right\|}}^{2}} + \alpha \left\langle {B(z - \bar {z}),\;z - \bar {z}} \right\rangle {\text{:}}\;z \in {{\mathbb{R}}^{n}}} \right\},$
где B – симметричная положительно определенная матрица, $\bar {z}$ – некоторый фиксированный вектор пространства ${{\mathbb{R}}^{n}}$, $\left\| {f - {{f}_{\delta }}} \right\| \leqslant \delta $;

в предположении, что оператор P дважды непрерывно дифференцируем, квадратичная аппроксимация целевой функции $\Psi (z)$ в задаче (2) и привлечение итерационного prox-метода для решения задачи минимизации с полученной квадратичной функцией, т.е.

(3)
$\begin{gathered} {{z}^{{k + 1}}} = arg\,min\left\{ {\Psi ({{z}^{k}}) + \Psi {\text{'}}({{z}^{k}})(z - {{z}^{k}}){{ + }_{{_{{_{{_{{_{{_{{_{{}}}}}}}}}}}}}}}} \right. \\ + \;\frac{1}{2}\langle \Psi ''({{z}^{k}})(z - {{z}^{k}}),z - {{z}^{k}}\rangle + \left. {\frac{1}{2}\beta {\text{||}}z - {{z}^{k}}{\text{|}}{{{\text{|}}}^{2}}{\text{:}}\;z \in {{\mathbb{R}}^{n}}} \right\}. \\ \end{gathered} $

Для численного нахождения очередного приближения в (3) привлекался метод наискорейшего спуска.

Как показано в [2], при условии разрешимости задачи в смысле наименьших квадратов

(4)
$inf\{ {{\left\| {P(z) - f} \right\|}^{2}}{\text{:}}\;z \in {{\mathbb{R}}^{n}}\} $
регуляризованная задача (2) имеет решение zα и при $\alpha (\delta ) \to 0$, ${{\delta }^{2}}{\text{/}}\alpha (\delta ) \to 0$, δ → 0 все предельные точки последовательности ${\text{\{ }}{{z}_{{\alpha (\delta )}}}{\text{\} }}$ являются решениями задачи (4), ближайшими к вектору $\bar {z}$, т.е. $\bar {z}$-нормальными квазирешениями. Кроме того, была установлена сходимость процесса (3) к некоторому регуляризованному решению uα и его {zα}-фейеровость.

Предположим, что множества решений зада-чи (2) и уравнения

(5)
$\Phi (z) + \alpha B(z - \bar {z}) = 0,$
где $\Phi (z) = P\,{\text{'}}{{(z)}^{{\text{т}}}}(P(z) - {{f}_{\delta }})$, совпадают при всех данных fδ, $\left\| {f - {{f}_{\delta }}} \right\| \leqslant \delta $, δ, α.

Запишем процесс (3) в эквивалентном виде

(6)
$\begin{gathered} {{z}^{{k + 1}}} = {{z}^{k}} - {{[\Phi {\text{'}}({{z}^{k}}) + \alpha B + \beta I]}^{{ - 1}}} \times \\ \times \;[\Phi ({{z}^{k}}) + \alpha B(z - \bar {z})]. \\ \end{gathered} $

При построении итераций согласно (6) наиболее трудоемкой операцией является вычисление и обращение симметричной матрицы

$\begin{gathered} A({{z}^{k}}) = \Phi {\text{'}}({{z}^{k}}) + \alpha B + \beta I = \frac{1}{2}({\text{||}}P({{z}^{k}}) - {{f}_{\delta }}{\text{||}}_{{{{\mathbb{R}}^{n}}}}^{2} + \\ + \;\alpha \langle B({{z}^{k}} - \bar {z}),{{z}^{k}} - \bar {z}\rangle + \beta {\text{||}}{{z}^{k}} - \bar {z}{\text{|}}{{{\text{|}}}^{2}})'' \\ \end{gathered} $
в каждой итерационной точке ${{z}^{k}}$. Поэтому естественно рассмотреть модифицированный вариант метода (6)
(7)
$\begin{gathered} {{z}^{{k + 1}}} = {{z}^{k}} - \gamma {{[\Phi {\text{'}}({{z}^{0}}) + \alpha B + \beta I]}^{{ - 1}}} \times \\ \times \;[\Phi ({{z}^{k}}) + \alpha B(z - \bar {z})], \\ \end{gathered} $
в котором матрица вторых производных A(z) вычисляется один раз в точке z0 в течение всего процесса итераций, а также дополнительно введен параметр шага γ. Считаем, что точка z0, в которой вычисляется матрица A(z), совпадает со стартовой точкой процесса (7).

Заметим, что при формировании как основного процесса (3), так и его модифицированного варианта (7) наряду с квадратичной аппроксимацией (см. например [3, гл. 5 §10]) применяется prox-метод, что позволяет не только получить сходимость итераций, но и оценить скорость сходимости.

1. ИССЛЕДОВАНИЕ СХОДИМОСТИ МОДИФИЦИРОВАННОГО МЕТОДА

Предположим, что при заданном уровне погрешности δ правой части уравнения (1) выбран параметр $\alpha (\delta )$, при котором ${{z}_{{\alpha (\delta )}}}$ аппроксимирует с приемлемой точностью некоторое $\bar {z}$-нормальное квазирешение $\hat {z}$ задачи (4). Введем обозначения:

$\begin{gathered} \mathcal{D} = (\Phi {\text{'}}({{z}^{0}}) + \alpha B + \beta I), \\ F(z) = {{\mathcal{D}}^{{ - 1}}}\left[ {\Phi (z) + \alpha B(z - \bar {z})} \right], \\ {{\mathcal{D}}_{{\text{1}}}} = (\Phi {\text{'}}({{z}^{0}}) + \alpha B),\quad m = \mathop {inf}\limits_{\left\| z \right\| = 1} \left\langle {{{\mathcal{D}}_{1}}z,z} \right\rangle , \\ M = \mathop {sup}\limits_{\left\| z \right\| = 1} \left\langle {{{\mathcal{D}}_{1}}z,z} \right\rangle . \\ \end{gathered} $

Установим сильную фейеровость оператора шага T процесса (7) и погрешность аппроксимации регуляризованного решения ${{z}_{{\alpha (\delta )}}}$ итерационным процессом (7).

Теорема 1. Пусть выполнены условия

$\begin{gathered} \left\| {\Phi ({{z}_{1}}) - \Phi ({{z}_{2}})} \right\| \leqslant {{N}_{1}}\left\| {{{z}_{1}} - {{z}_{2}}} \right\|, \\ \left\| {\Phi {\text{'}}({{z}_{1}}) - \Phi {\text{'}}({{z}_{2}})} \right\| \leqslant {{N}_{2}}\left\| {{{z}_{1}} - {{z}_{2}}} \right\| \\ \end{gathered} $
в шаре $S({{z}^{0}},d)$, содержащем ${{z}_{{\alpha (\delta )}}}$, $\hat {z}$, а также

${\text{||}}{{z}_{{\alpha (\delta )}}} - {{z}^{0}}{\text{||}} \leqslant {{r}^{0}},\quad {{r}^{0}} = m{\text{/}}5{{N}_{2}},\quad m > 0.$

Тогда при $\gamma < m(m + \beta ){\text{/}}{{({{N}_{1}} + \alpha \left\| B \right\|)}^{2}}$ оператор T есть сильно ${\text{\{ }}{{z}_{{\alpha (\delta )}}}{\text{\} }}$-фейеровский в шаре $S({{z}^{0}},{{r}^{0}})$ и zk сходится к ${{z}_{{\alpha (\delta )}}}$, а при $\gamma = m(m + \beta ){\text{/}}(2{{({{N}_{1}} + \alpha \left\| B \right\|)}^{2}})$ справедлива оценка

(8)
$\begin{gathered} {\text{||}}z_{{\alpha (\delta )}}^{k} - {{z}_{{\alpha (\delta )}}}{\text{||}} \leqslant {{q}^{k}}{\text{||}}{{z}^{0}} - {{z}_{{\alpha (\delta )}}}{\text{||}}, \\ q = \sqrt {1 - {{m}^{2}}{\text{/}}{{{({{N}_{1}} + \alpha \left\| B \right\|)}}^{2}}} , \\ \end{gathered} $
где z0– начальная точка в методе (7).

Доказательство. Используя введенные обозначения, условия теоремы и самосопряженность оператора ${{\mathcal{D}}_{1}}$, имеем для $z \in S({{z}_{{\alpha (\delta )}}},{{r}_{0}})$ соотношения

$\begin{gathered} \left\langle {F(z),z - {{z}_{\alpha }}} \right\rangle = \left\langle {F(z),z - {{z}_{\alpha }}} \right\rangle - \left\langle {F({{z}_{\alpha }}),z - {{z}_{\alpha }}} \right\rangle = \\ = \;\left\langle {{{\mathcal{D}}^{{ - 1}}}\left[ {\Phi (z) - \Phi ({{z}_{\alpha }}) + \alpha B(z - {{z}_{\alpha }})} \right],z - {{z}_{\alpha }}} \right\rangle = \\ = \;\left\langle {{{\mathcal{D}}^{{ - 1}}}} \right.\left[ {\Phi (z)} \right. - (\Phi (z) + \Phi {\text{'}}(z)({{z}_{\alpha }} - z) + \xi ) + \\ + \;\left. {\left. {\alpha B(z - {{z}_{\alpha }})} \right],z - {{z}_{\alpha }}} \right\rangle = \\ = \left\langle {{{\mathcal{D}}^{{ - 1}}}} \right.\left[ {\Phi {\text{'}}(z)} \right.(z - {{z}_{\alpha }}) + \alpha B(z - {{z}_{\alpha }}) + \\ \end{gathered} $
$\begin{gathered} + \;(\Phi {\text{'}}({{z}^{0}})(z - {{z}_{\alpha }}) - \left. {\left. {\Phi {\text{'}}(z)(z - {{z}_{\alpha }})) - \xi } \right],z - {{z}_{\alpha }}} \right\rangle \geqslant \\ \geqslant \;\frac{m}{{m + \beta }}{{\left\| {z - {{z}_{\alpha }}} \right\|}^{2}} - \frac{1}{{m + \beta }} \times \\ \times \;\left[ {{{N}_{2}}{{{\left\| {z - {{z}_{\alpha }}} \right\|}}^{2}}{\text{||}}z - {{z}^{0}}{\text{||}} + \frac{1}{2}{{N}_{2}}{{{\left\| {z - {{z}_{\alpha }}} \right\|}}^{3}}} \right] \geqslant \\ \geqslant \;\frac{m}{{m + \beta }}{{\left\| {z - {{z}_{\alpha }}} \right\|}^{2}} - \frac{1}{{m + \beta }}{{N}_{2}} \times \\ \times \;\left( {\left\| {z - {{z}_{\alpha }}} \right\| + \,{\text{||}}{{z}_{\alpha }} - {{z}^{0}}{\text{||}} + \frac{1}{2}\left\| {z - {{z}_{\alpha }}} \right\|} \right){{\left\| {z - {{z}_{\alpha }}} \right\|}^{2}}. \\ \end{gathered} $

Подставляя в соотношения ${{z}_{\alpha }} = {{z}_{{\alpha (\delta )}}}$ и учитывая в правой части последнего неравенства $z\, \in \,S({{z}_{{\alpha (\delta )}}}$, r0), ${\text{||}}z - {{z}^{0}}{\text{||}} \leqslant {\text{||}}z - {{z}_{{\alpha (\delta )}}}{\text{||}} + \,{\text{||}}{{z}_{{\alpha (\delta )}}} - {{z}^{0}}{\text{||}} \leqslant 2r$, $r = m{\text{/}}5{{N}_{2}}$, окончательно получаем оценку снизу:

(9)
$\left\langle {F(z),z - {{z}_{{\alpha (\delta )}}}} \right\rangle \geqslant \frac{m}{{2(m + \beta )}}{{\left\| {z - {{z}_{{\alpha (\delta )}}}} \right\|}^{2}}.$

В силу условия Липшица для $\Phi (z)$, справедлива оценка сверху для ${{\left\| {F(z)} \right\|}^{2}}$:

$\begin{gathered} {{\left\| {F(z)} \right\|}^{2}} = {{\left\| {F(z) - F({{z}_{{\alpha (\delta )}}})} \right\|}^{2}} = \\ = \;{{\left\| {{{\mathcal{D}}^{{ - 1}}}\left[ {\Phi (z) - \Phi ({{z}_{{\alpha (\delta )}}}) + \alpha B(z - {{z}_{\alpha }})} \right]} \right\|}^{2}} \leqslant \\ \end{gathered} $
(10)
$ \leqslant \frac{1}{{{{{(m + \beta )}}^{2}}}}{{({{N}_{1}} + \alpha {\text{||}}B{\text{||}})}^{2}}{\text{||}}z - {{z}_{{\alpha (\delta )}}}{\text{|}}{{{\text{|}}}^{2}}.$

Объединяя (9) и (10), приходим к оценке

(11)
${{\left\| {F(z)} \right\|}^{2}} \leqslant \frac{{2{{{({{N}_{1}} + \alpha \left\| B \right\|)}}^{2}}}}{{m(m + \beta )}}\left\langle {F(z),z - {{z}_{\alpha }}(\delta )} \right\rangle .$

Сильная {uα}-фейеровость оператора шага T в процессе (7) означает выполнение неравенства (см. определение 1.5 в [4])

(12)
${{\left\| {T(z) - {{z}_{{\alpha (\delta )}}}} \right\|}^{2}} - {{\left\| {z - {{z}_{{\alpha (\delta )}}}} \right\|}^{2}} + \nu {{\left\| {T(z) - z} \right\|}^{2}} \leqslant 0$
для некоторого ν > 0, что эквивалентно неравенству

(13)
${{\left\| {F(z)} \right\|}^{2}} \leqslant \frac{2}{{\gamma (1 + \nu )}}\left\langle {F(z),z - {{z}_{{\alpha (\delta )}}}} \right\rangle .$

Сравнивая неравенства (11) и (13), убеждаемся, что при $\gamma < m(m + \beta ){\text{/}}{{({{N}_{1}} + \alpha \left\| B \right\|)}^{2}}$ оператор T является сильно ${\text{\{ }}{{z}_{{\alpha (\delta )}}}{\text{\} }}$-фейеровским, что влечет сходимость zk к ${{z}_{{\alpha (\delta )}}}$ [4].

Если теперь подставить оценки (8) и (9) в соотношение

(14)
$\begin{gathered} {{\left\| {T(z) - {{z}_{\alpha }}} \right\|}^{2}} = \\ = {{\left\| {z - {{z}_{\alpha }}} \right\|}^{2}} - 2\gamma \left\langle {F(z),z - {{z}_{\alpha }}} \right\rangle + {{\gamma }^{2}}{{\left\| {F(z)} \right\|}^{2}} \\ \end{gathered} $
и проминимизировать правую часть полученного неравенства по γ, то получим значение γ = m(m + + β)/ $2{{({{N}_{1}} + \alpha \left\| B \right\|)}^{2}}$, при котором при z = zk (14) перейдет в оценку (8).

Замечание 1. В работе [2] для основного процесса (6) доказана сходимость итераций к регуляризованному решению и свойство фейеровости оператора шага. Переход к модифицированному варианту (7), введение дополнительного параметра γ и использование иной схемы доказательства позволили установить сильную фейеровость оператора T, что влечет не только сходимость, но и оценку скорости сходимости итераций.

2. ОБРАТНАЯ ЗАДАЧА ЗОНДИРОВАНИЯ АТМОСФЕРЫ

В качестве примера реализации представленного метода выбрана задача восстановления высотных профилей $\delta {{D}_{{{\text{HDO}}}}}(h)$ относительной концентрации молекул HDO в водяном паре по инфракрасным спектрам пропускания солнечного света. Количественные данные об относительной концентрации тяжелой воды являются важным маркером скрытых трендов в гидрологических процессах, происходящих в атмосфере.

Относительная концентрация $\delta {{D}_{{{\text{HDO}}}}}$ определяется формулой

$\delta {{D}_{{{\text{HDO}}}}} = \left( {\frac{R}{{{{R}_{0}}}} - 1} \right) \times 1000\% ,\quad R = \frac{{{{N}_{{{\text{HDO}}}}}}}{{{{N}_{{{{{\text{H}}}_{{\text{2}}}}{\text{O}}}}}}},$
где ${{N}_{{{\text{HDO}}}}}$ и ${{N}_{{{{{\text{H}}}_{{\text{2}}}}{\text{O}}}}}$ – соответственно число молекул HDO и H2O в единице объема; ${{R}_{0}} = 3.1069 \times {{10}^{{ - 4}}}$ – среднее значение R в океанической воде. Величина $\delta {{D}_{{{\text{HDO}}}}}$ в атмосфере лежит в пределах – 1000–0‰.

Исследуемая задача сводится к решению переопределенной системы нелинейных уравнений

(15)
$P(z) = {{\bar {P}}_{\delta }},$
где ${{\bar {P}}_{\delta }}$ – измеренный с ошибкой $\delta $ вектор значений ИК-спектра солнечного света, на сетке частот из микроокон, содержащих выбранные спектральные линии HDO/H2O. Спектр измерялся фурье-интерферометром FTIR (Fourier Transform InfraRed “Bruker” IFS 125M), установленным на Уральской атмосферной станции, расположенной на территории Коуровской астрономической обсерватории УрФУ, находящейся в окрестности Екатеринбурга. Вектор $z$ составлен из значений $ln{{N}_{{{{{\text{H}}}_{{\text{2}}}}{\text{O}}}}}({{h}_{i}})$ и $\delta {{D}_{{{\text{HDO}}}}}({{h}_{i}})$ на заданной сетке высот {hi}, P – нелинейный оператор, вычисляющий модельный спектр по заданному z. Оператор P строится на основе данных, полученных в результате работы программы FIRE-ARMS, описание которой содержится в работе [5], решающей прямую задачу вычисления модельного спектра по известным профилям концентрации газов, давления и температуры, полученных из реанализа метеоданных на момент измерения спектра. Построение нелинейного оператора $P(z)$ подробно описано в [2].

Решение задачи (15) сильно чувствительно к ошибкам измерения спектра, поэтому вначале применяется тихоновская регуляризация со стабилизатором, полученным на основе метода максимального правдоподобия,

(16)
$min\left\{ {{{{\left\| {P(z) - {{{\bar {P}}}_{\delta }}} \right\|}}^{2}} + \alpha \langle {{C}^{{ - 1}}}(z - \bar {z}),z - \bar {z}\rangle {\text{:}}\;z \in {{R}^{n}}} \right\},$
где $C = (K + \varepsilon I)$, $\bar {z}$ и K – среднее значение и ковариационная матрица априорного распределения вектора z, $\varepsilon $ – малый положительный параметр, выбираемый из условия положительной определенности матрицы $C$.

Для аппроксимации решения zα задачи (16) использовался итерационный процесс (7) с γ = 1 и несколько видоизмененный, что позволило исключить обращение матрицы C,

(17)
$\begin{gathered} {{z}^{{k + 1}}} = {{z}^{k}} - {{[C\Phi {\text{'}}({{z}^{0}}) + \alpha I + \beta C]}^{{ - 1}}} \times \\ \times \;(C\Phi ({{z}^{k}}) + \alpha ({{z}^{k}} - \bar {z})), \\ \end{gathered} $
где $\Phi (z) = P\,{\text{'}}{{(z)}^{{\text{т}}}}(P(z) - {{\bar {P}}_{\delta }})$, $\beta $ – параметр prox-метода. 

Для сравнения, задача минимизации (16) решалась также с помощью немодифицированного итерационного процесса с вычислением матрицы вторых производных на каждом шаге итерации

(18)
$\begin{gathered} {{z}^{{k + 1}}} = {{z}^{k}} - {{[C\Phi {\text{'}}({{z}^{k}}) + \alpha I + \beta C]}^{{ - 1}}} \times \\ \times \;(C\Phi ({{z}^{k}}) + \alpha ({{z}^{k}} - \bar {z})). \\ \end{gathered} $

Метод (17) был протестирован на множестве наборов (около 260) измеренных спектральных данных, описание которых содержится в [6]. Также было произведено сравнение времени обработки спектров методами (17) и (18). Для данной задачи новый метод дает прирост производительности на несколько процентов.

Замечание 2. Незначительная разница по затратам машинного времени при использовании модифицированного (17) и немодифицированного (18) процессов в описанном эксперименте обусловлена невысокой размерностью задачи (dim(z) = = 100). Однако в задачах высокой размерности преимущество модифицированных методов по времени является существенным (см. [7]).

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

  1. Skorik G.G. // Evrasian Journal Mathematical and Computer Applications. 2018. V. 6. Iss. 1. P. 56–64.

  2. Skorik G.G., Vasin V.V. // Evrasian J. Math. and Comp. Appl. 2019. V. 7. Iss. 2. P. 79–88.

  3. Васильев Ф.П. Методы оптимизации. Кн. 1 М.: МЦНМО, 2011. 619 с.

  4. Vasin V.V., Eremin I.I. Operators and Iterative Processes of Fejér Type. Theory and Application. B., N.Y.: Walter de Gruyter, 2009. 155 p.

  5. Gribanov K.G., Zakharov V.I., Tashkun S.A. et al. // Quantitative Spectroscopy and Radiative Transfer. 2001. V. 68. Iss. 4. P. 435–451.

  6. Skorik G.G., Vasin V.V., Gribanov K.G. et al. // Doklady Earth Sciences. 2014. V. 454. Iss. 2. P. 208–212.

  7. Васин В.В., Скурыдина А.Ф. // Тр. Ин-та математики и механики УрО РАН. 2019. Т. 23. № 1. С. 57–74.

Дополнительные материалы отсутствуют.

Инструменты

Доклады Российской академии наук. Математика, информатика, процессы управления