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

Full Text

Введение

В последние годы растущий интерес к проблемам усвоения данных наблюдений в математических моделях геофизической гидродинамики обусловлен значительным прогрессом в создании все более мощных вычислительных комплексов, развитием новых спутниковых и контактных измерительных технологий, а также новых методов и численных алгоритмов. Необходимость синтеза данных измерений и данных моделирования породила сложнейшие научные и технологические проблемы, связанные с усовершенствованием вычислительных методов моделирования и анализа природных явлений. Методы усвоения данных позволяют связать модельные расчеты с реальными данными для построения или уточнения неизвестных входных данных или параметров и повышения точности прогноза (см. [1]–[9]).

При рассмотрении задач вариационного усвоения данных наблюдений важную роль играет анализ чувствительности оптимального решения и его функционалов по отношению к входным данным (см. [10]–[17]). В [18] проведено исследование чувствительности функционалов от оптимального решения задачи вариационного усвоения данных с целью восстановления потоков тепла на поверхности моря при использовании ковариационных матриц ошибок данных наблюдений и ошибок начального приближения (бэкграунда). Настоящая работа обобщает результаты работы [18] на случай задачи вариационного усвоения данных при одновременном восстановлении начального состояния и потока тепла для модели термодинамики моря. Исследована чувствительность функционалов от оптимального решения по отношению к входным данным о потоке тепла, доказана теорема о представлении градиента функционала через гессиан функции стоимости и приведены результаты численных экспериментов для модели динамики Черного моря, разработанной в ИВМ РАН.

1. Задача вариационного усвоения данных для модели термодинамики моря

Рассмотрим задачу термодинамики моря в виде (см. [19], [20])

$T_t+(\overline{U},Grad)T-Div({\stackrel{\wedge }{a}}_T\cdot GradT)=f_T$ в $D\times \left(\mathrm{0,}\overline{t}\right),$

$T=T_0$ при $t=0$ в $D$

$-{\nu }_T\frac{\partial T}{\partial z}=Q$ на ${\Gamma }_S\times \left(\mathrm{0,}\overline{t}\right),$ $\frac{\partial T}{\partial N_T}=0$на ${\Gamma }_{w,c}\times \left(\mathrm{0,}\overline{t}\right),$

${\overline{U}}_n^{\left(-\right)}T+\frac{\partial T}{\partial N_T}={\overline{U}}_n^{\left(-\right)}{d}_T+Q_T$ на ${\Gamma }_{w,op}\times \left(\mathrm{0,}\overline{t}\right),$

$\frac{\partial T}{\partial N_T}=0$ на ${\Gamma }_H\times \left(\mathrm{0,}\overline{t}\right),$ (1.1)

где $T=T(x,y,z,t)$ — неизвестная функция температуры, $t\in \left(\mathrm{0,}\overline{t}\right)$, $(x,y,z)\in D=\Omega \times \left(\mathrm{0,}H\right)$, $\Omega \subset R^2$, $H=H(x,y)$ — функция рельефа дна, $Q=Q(x,y,t)$ — суммарный приток тепла, $\overline{U}=(u,v,w)$, ${\stackrel{\wedge }{a}}_T=$ diag $\left({\left(a_T\right)}_{ii}\right)$, ${\left(a_T\right)}_{11}={\left(a_T\right)}_{22}={\mu }_T$, ${\left(a_T\right)}_{33}={\nu }_T$, $f_T=f_T(x,y,z,t)$ — заданные функции. Скорости $u,v,w$ зависят в общем случае от пространства и времени и считаются заданными в данной постановке, а коэффициенты ${\mu }_T,{\nu }_T$ предполагаются зависящими только от пространственных переменных на рассматриваемом интервале по времени. Граница области $\Gamma \equiv \partial D$ представляется как объединение четырех непересекающихся частей ${\Gamma }_S$, ${\Gamma }_{w,op}$, ${\Gamma }_{w,c}$, ${\Gamma }_H$, где ${\Gamma }_S=\Omega$ (невозмущенная поверхность моря), ${\Gamma }_{w,op}$ — жидкая (открытая) часть вертикальной боковой границы, ${\Gamma }_{w,c}$ — твердая часть вертикальной боковой границы, ${\Gamma }_H$ — дно моря. Другие обозначения и детальное описание постановки задачи можно найти в работах [21], [22], [23].

Запишем задачу (1.1) в форме операторного уравнения в ${\left(W_2^1\right(D\left)\right)}^{*}$:

$T_t+LT=F+BQ$ для п.в. $t\in \left(\mathrm{0,}\overline{t}\right),$

$T=T_0$ при $t=\mathrm{0,}$(1.2)

где равенство понимается в обобщенном смысле, а именно,

$(T_t,\widehat{T})+(LT,\widehat{T})=F\left(\widehat{T}\right)+(BQ,\widehat{T})\forall \widehat{T}\in W_2^1\left(D\right),$ (1.3)

а операторы L, $F,B$ определяются интегральными соотношениями

$(LT,\widehat{T})\equiv \underset{D}{\int }(-TDiv(\overline{U}\widehat{T}\left)\right)dD+\underset{{\Gamma }_{w,op}}{\int }{\overline{U}}_n^{(+)}T\widehat{T}d\Gamma +\underset{D}{\int }{\widehat{a}}_{T}Grad\left(T\right)\times Grad\left(\widehat{T}\right)dD,$

$F\left(\widehat{T}\right)=\underset{{\Gamma }_{w,op}}{\int }(Q_T+{\overline{U}}_n^{(-)}{d}_T)\widehat{T}d\Gamma +\underset{D}{\int }{f}_T\widehat{T}dD,$

$(T_t,\widehat{T})=\underset{D}{\int }{T}_t\widehat{T}dD,(BQ,\widehat{T})=\underset{\Omega }{\int }Q\widehat{T}{|}_{z=0}d\Omega ,$

при этом функции ${\widehat{a}}_T$, $Q_T$, $f_T$, $Q$ таковы, что равенство (1.3) имеет смысл.

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

Рассмотрим задачу об усвоении данных о температуре поверхности моря, следуя [22], [24], [18]. Предположим, что в задаче (1.1) функции $T_0\in L_2\left(D\right)$ и $Q\in L_2(\Omega \times (\mathrm{0,}\overline{t}\left)\right)$ не известны. Пусть задана функция данных наблюдений $T_{obs}(x,y,t)$ на $\overline{\Omega }\equiv \Omega \cup \partial \Omega$ при $t\in \left(\mathrm{0,}\overline{t}\right)$, которая по своему физическому смыслу есть приближение к функции поверхностной температуры на W, т.е. к $T{|}_{z=0}$. Считаем, что $T_{obs}\in L_2(\Omega \times (\mathrm{0,}\overline{t}\left)\right)$. Допускается случай, когда $T_{obs}$ имеется лишь на некотором подмножестве из $\Omega \times \left(\mathrm{0,}\overline{t}\right)$, характеристическую функцию которого обозначим через $m_0$.

Будем предполагать, что данные наблюдений $T_{obs}$ заданы с ошибками, а именно,

$T_{obs}=m_0{T}^t{|}_{z=0}+{\xi }_{obs},$

где $T^t$ — точное решение задачи (1.1) при некотором $Q=Q^t$, а ${\xi }_{obs}\in Y_{obs}=L_2(\Omega \times (\mathrm{0,}\overline{t}\left)\right)$ рассматривается как ошибка наблюдений в пространстве наблюдений $Y_{obs}$. Предполагается, что ошибки ${\xi }_{obs}$ случайные и они распределены по нормальному закону (гауссовские) с нулевым математическим ожиданием и ковариационным оператором $R\cdot =E\left[\right(\cdot ,{\xi }_{obs}\left){\xi }_{obs}\right]$, $R:Y_{obs}\to Y_{obs}$, где E — математическое ожидание. В дальнейшем мы будем предполагать, что $R$ положительно определен и, значит, обратим.

Рассмотрим следующую задачу вариационного усвоения данных: найти T, $T_0$ и $Q$, такие что

$T_t+LT=F+BQ,$

$T=T_0$ при $t=\mathrm{0,}$ (1.4)

$J\left(T_0,Q\right)=\underset{T_0,Q}{\inf }J\left(T_{\mathrm{0,}}Q\right),$

где

$\begin{array}{c}J(T_0,Q)=\frac{\alpha }{2}\underset{D}{\int }|T_0-T^{\left(0\right)}{|}^{2}dD+\frac{1}{2}\underset{0}{\overset{\overline{t}}{\int }}\underset{\Omega }{\int }(Q-Q^{\left(0\right)}\left){\mathcal{B}}^{-1}\right(Q-Q^{\left(0\right)})d\Omega dt+\\ +\frac{1}{2}\underset{0}{\overset{\overline{t}}{\int }}\underset{\Omega }{\int }(m_{0}T{|}_{z=0}-T_{obs})\times R^{-1}(m_{0}T{|}_{z=0}-T_{obs})d\Omega dt,\end{array}$

$T^{\left(0\right)}=T^{\left(0\right)}(x,y,z)$, $Q^{\left(0\right)}=Q^{\left(0\right)}(x,y,t)$ — заданные функции, $\alpha =$const$>0$, $\mathcal{B}:Y_{obs}\to Y_{obs}$ — ковариационный оператор ошибок бэкграунда. Функция $Q^{\left(0\right)}$ обычно выбирается в качестве начального приближения для неизвестного потока $Q$ (так называемый бэкграунд или фоновый поток). Цель вариационного усвоения данных — используя $T^{\left(0\right)},Q^{\left(0\right)}$, найти лучшие оценки для $T_0,Q$, согласованные с решением модели и наблюдениями для дальнейшего моделирования и прогноза.

Слагаемые с параметром a и с весовым оператором ${\mathcal{B}}^{-1}$ играют роль регуляризации по Тихонову (см. [25]), они считаются заданными при рассмотрении задачи. Если оператор ${\mathcal{B}}^{-1}$ положительно определен, то при $\alpha >0$ поставленная задача вариационного усвоения данных имеет единственное решение. Существование оптимального решения следует из классических результатов теории экстремальных задач (см. [2]), так как можно показать, что решение задачи (1.2) непрерывно зависит от начального состояния $T_0$ и потока Q (имеют место априорные оценки в соответствующих функциональных пространствах).

Необходимое условие оптимальности $gradJ=0$, которое определяет решение сформулированной задачи вариационного усвоения данных, приводит к системе оптимальности (см. [2]):

$T_t+LT=F+BQ,t\in \left(\mathrm{0,}\overline{t}\right),$

$T=T_0$ при $t=\mathrm{0,}$ (1.5)

$-{\left(T^{*}\right)}_t+L^{*}{T}^{*}=B{R}^{-1}{m}_0(B^{*}T-T_{obs}),t\in \left(\mathrm{0,}\overline{t}\right),$ 

$T^{\ast }=0$ при $t=\overline{t}$ (1.6)

$\alpha (T_0-T^{\left(0\right)})+T^{*}{|}_{t=0}=0$ в $D,$ (1.7)

${\mathcal{B}}^{-1}(Q-Q^{\left(0\right)})+B^{*}{T}^{*}=0$ на $\Omega \times \left(\mathrm{0,}\overline{t}\right),$ (1.8)

где $L^{*},B^{*}$ — операторы, сопряженные к $L,B$ соответственно. Последние два уравнения в этой системе означают, что градиент функции стоимости $J(T_0,Q)$ по отношению к начальному состоянию $T_0$ и потоку тепла Q обращается в нуль. Неизвестными в (1.5)–(1.8) являются функции $T,T^{*},T_0$ и Q.

Отметим, что если $T_{obs}$ — случайная гауссовская величина, то и оптимальное решение будет случайной величиной, распределенной по Гауссу, в силу того что оптимальное решение линейно зависит от $T_{obs}$, согласно системе оптимальности (1.5)–(1.8).

2. Чувствительность к входным данным о потоке тепла

Рассмотрим функцию $G(T,T_0,Q)$, зависящую от $T,T_0,Q$, которая предполагается вещественнозначной и может рассматриваться как функционал на $X=L_2(D\times (\mathrm{0,}\overline{t}\left)\right)\times L_2\left(D\right)\times L_2(\Omega \times (\mathrm{0,}\overline{t}\left)\right)$. Нас интересует чувствительность функционала $G(T,T_0,Q)$ к входным данным о потоке тепла $Q^{\left(0\right)}$ при условии, что $T,T_0,Q$ получены после вариационного усвоения из системы оптимальности (1.5)–(1.8). Как известно из [1], [26], чувствительность функционала определяется градиентом по $Q^{\left(0\right)}$, который является производной Гато:

$\frac{dG}{d{Q}^{\left(0\right)}}=\frac{\partial G}{\partial T}\frac{\partial T}{\partial Q^{\left(0\right)}}+\frac{\partial G}{\partial T_0}\frac{\partial T_0}{\partial Q^{\left(0\right)}}+\frac{\partial G}{\partial Q}\frac{\partial Q}{\partial Q^{\left(0\right)}}.$ (2.1)

Обозначим через $\delta Q^{\left(0\right)}$ вариацию функции $Q^{\left(0\right)}$. Из (1.5)–(1.8) выводим систему оптимальности для вариаций $\delta T,\delta T^{*},\delta T_0,\delta Q$:

$\delta T_t+L\delta T=B\delta Q,t\in \left(\mathrm{0,}\overline{t}\right),$

$\delta T=\delta T_0$ при $t=\mathrm{0,}$ (2.2)

$\begin{array}{c}-{\left(\delta T^{*}\right)}_t+L^{*}\delta T^{*}=B{R}^{-1}{m}_0{B}^{*}\delta T,\\ t\in \left(\mathrm{0,}\overline{t}\right),\end{array}$

$\delta T^{\ast }=0$ при $t=\overline{t},$ (2.3)

$\alpha \delta T_0+\delta T^{\ast }{|}_{t=0}=0$ в $D,$ (2.4)

${\mathcal{B}}^{-1}(\delta Q-\delta Q^{\left(0\right)})+B^{*}\delta T^{*}=0$ на $\Omega \times \left(\mathrm{0,}\overline{t}\right).$ (2.5)

Отметим, что данные наблюдений $T_{obs}$ и функция $T^{\left(0\right)}$ уже не входят в систему (2.2)–(2.5), в отличие от (1.5)–(1.8). Нетрудно видеть, что система (2.2)–(2.5) эквивалентна следующей задаче оптимального управления для определения $\delta T,\delta T_0,\delta Q$ таких, что

$\delta T_t=B\delta Q,t\in \left(\mathrm{0,}\overline{t}\right),$

$\delta T=\delta T_0$ при $t=\mathrm{0,}$

$S\left(\delta T_0,\delta Q\right)=\underset{T_0,Q}{\inf }S\left(T_0,Q\right),$ (2.6)

где

$S(\delta T_0,\delta Q)=\frac{\alpha }{2}\underset{D}{\int }|\delta T_0{|}^{2}dD+\frac{1}{2}\underset{0}{\overset{\overline{t}}{\int }}\underset{\Omega }{\int }(\delta Q-\delta Q^{\left(0\right)}\left){\mathcal{B}}^{-1}\right(\delta Q-\delta Q^{\left(0\right)})d\Omega dt+\frac{1}{2}\underset{0}{\overset{\overline{t}}{\int }}\underset{\Omega }{\int }{m}_0\delta T{|}_{z=0}{R}^{-1}\delta T{|}_{z=0}d\Omega dt.$ (2.7)

Справедлива следующая

Лемма 1. Гессиан $\mathcal{H}$ функционала (2.7) определяется на $\Theta ={(\xi ,\eta )}^T,\xi \in L_2\left(D\right),\eta \in L_2(\Omega \times (\mathrm{0,}\overline{t}\left)\right)$ последовательным решением задач

${\psi }_t+L\psi =B\eta ,t\in \left(\mathrm{0,}\overline{t}\right),$

$\psi =\xi$ при $t=\mathrm{0,}$ (2.8)

$-{\left({\psi }^{\ast }\right)}_t+L^{\ast }{\psi }^{\ast }=B{R}^{-1}{m}_0{B}^{\ast }\psi ,t\in \left(\mathrm{0,}\overline{t}\right),$

${\psi }^{\ast }=0$ при $t=\overline{t},$ (2.9)

^{$\mathcal{H}\Theta ={(\alpha \xi +{\psi }^{*}{|}_{t=0},{\mathcal{B}}^{-1}\eta +B^{*}{\psi }^{*})}^T.$} (2.10)

Доказательство. Согласно системе оптимальности (2.2)–(2.5) градиент функционала (2.7) определяется по формуле

^{$gradS={(\alpha \delta T_0+\delta T^{*}{|}_{t=0},{\mathcal{B}}^{-1}(\delta Q-\delta Q^{\left(0\right)})+B^{*}\delta T^{*})}^T,$} 

где $\delta T^{*}$ — решение сопряженной задачи (2.3). Продифференцируем последнюю формулу еще раз по $\delta T_0,\delta Q$, чтобы получить правило действия гессиана

$\mathcal{H}\Theta ={(\alpha \xi +{\psi }^{*}{|}_{t=0},{\mathcal{B}}^{-1}\eta +B^{*}{\psi }^{*})}^T,$

где ^{$\Theta ={(\xi ,\eta )}^T$}, — вариация ^{${(\delta T_0,\delta Q)}^T$}, а ${\psi }^{*}$ — решение сопряженной задачи (2.9), которая есть не что иное, как продифференцированная задача (2.3). При этом $\psi$ — решение задачи (2.8), которая получена из (2.2) дифференцированием по $\delta T_0,\delta Q$. Лемма доказана.

Введем оператор $\mathcal{P}:L_2(\Omega \times (\mathrm{0,}\overline{t}\left)\right)\to L_2\left(D\right)\times L_2(\Omega \times (\mathrm{0,}\overline{t}\left)\right)$, действующий на функции $g\in L_2(\Omega \times (\mathrm{0,}\overline{t}\left)\right)$ по формуле

$\mathcal{P}g={\left(\mathrm{0,}{\mathcal{B}}^{-1}g\right)}^T.$ (2.11)

Используя (2.8)–(2.10), нетрудно видеть (см. [27]), что система (2.2)–(2.5) эквивалентна уравнению для вариации оптимального решения $\delta \Theta ={(\delta T_0,\delta Q)}^T:$

$\mathcal{H}\delta \Theta =\mathcal{P}\delta Q^{\left(0\right)}.$ (2.12)

Гессиан $\mathcal{H}$ действует в $L_2\left(D\right)\times L_2(\Omega \times (\mathrm{0,}\overline{t}\left)\right)$ с областью определения $D\left(\mathcal{H}\right)=L_2\left(D\right)\times L_2(\Omega \times (\mathrm{0,}\overline{t}\left)\right)$, он ограничен, самосопряжен и неотрицательно определен. Если ${\mathcal{B}}^{-1}$ положительно определен, то $\mathcal{H}$ положительно определен при $\alpha >0$, поскольку $(\mathcal{H}\Theta ,\Theta )\ge \alpha (\xi ,\xi )+({\mathcal{B}}^{-1}\eta ,\eta )$. В последнем случае уравнение (2.12) имеет единственное решение

$\delta \Theta ={\mathcal{H}}^{-1}\mathcal{P}\delta Q^{\left(0\right)}.$ (2.13)

Формула (2.13) дает в явном виде выражение для вариаций оптимального решения $\delta T_0,\delta Q$ через вариацию функции начального приближения (бэкграунда) $\delta Q^{\left(0\right)}$. Следуя [18], уравнение вида (2.13) может быть положено в основу исследования чувствительности оптимального решения и его функционалов к ошибкам входных данных о потоке тепла.

Справедлива следующая

Теорема 1. Градиент функционала $G(T,T_0,Q)$ по $Q^{\left(0\right)}$ имеет вид

$\frac{dG}{d{Q}^{\left(0\right)}}={\mathcal{P}}^{*}{\mathcal{H}}^{-1}\mathcal{F},$ (2.14)

где

^{$\mathcal{F}={\left({\varphi }^{\ast }{|}_{t=0}+\frac{\partial G}{\partial T_0},B^{\ast }{\phi }^{\ast }+\frac{\partial G}{\partial Q}\right)}^T,$},(2.15)

${\mathcal{P}}^{*}$ — оператор, сопряженный к $\mathcal{P}$, $\mathcal{H}$ — гессиан, определенный формулами (2.8)–(2.10), а ${\varphi }^{*}$ — решение сопряженной задачи

$-{\left({\varphi }^{*}\right)}_t+L^{*}{\varphi }^{*}=\frac{\partial G}{\partial T},t\in \left(\mathrm{0,}\overline{t}\right),$

${\varphi }^{\ast }=0$ при $t=\overline{t},$ (2.16)

Доказательство. Рассмотрим значение градиента (1) на вариации $\delta Q^{\left(0\right)}$:

${\left(\frac{dG}{d{Q}^{\left(0\right)}},\delta Q^{\left(0\right)}\right)}_{Y_{obs}}={\left(\frac{\partial G}{\partial T},\delta T\right)}_Y+{\left(\frac{\partial G}{\partial T_0},\delta T_0\right)}_{L_2\left(D\right)}+{\left(\frac{\partial G}{\partial Q},\delta Q\right)}_{Y_{obs}},$ (2.17)

где $\delta Q^{\left(0\right)}$ — вариация функции $Q^{\left(0\right)}$, $\delta T=\frac{\partial T}{\partial Q^{\left(0\right)}}\delta Q^{\left(0\right)}$, $\delta T_0=\frac{\partial T_0}{\partial Q^{\left(0\right)}}\delta Q^{\left(0\right)}$, $\delta Q=\frac{\partial Q}{\partial Q^{\left(0\right)}}\delta Q^{\left(0\right)}$ — решения системы (2.2)–(2.5), $Y=L_2(D\times (\mathrm{0,}\overline{t}\left)\right)$.

Задача (2.16) является сопряженной по отношению к (2.2), поэтому в силу соотношения сопряженности

${\left(\frac{\partial G}{\partial T},\delta T\right)}_Y+{({\varphi }^{*}{|}_{t=0},\delta T_0)}_{L_2\left(D\right)}={({\varphi }^{*},B\delta Q)}_Y={(B^{*}{\varphi }^{*},\delta Q)}_{Y_{obs}}+{({\varphi }^{*}{|}_{t=0},\delta T_0)}_{L_2\left(D\right)}.$ (2.18)

Из (2.17), (2.18) получаем

${\left(\frac{dG}{d{Q}^{\left(0\right)}},\delta Q^{\left(0\right)}\right)}_{Y_{obs}}={({\varphi }^{*}{|}_{t=0}+\frac{\partial G}{\partial T_0},\delta T_0)}_{L_2\left(D\right)}++{\left(B^{*}{\varphi }^{*}+\frac{\partial G}{\partial Q},\delta Q\right)}_{Y_{obs}}=(\mathcal{F},\delta \Theta ),$ (2.19)

где ^{$\delta \Theta ={(\delta T_0,\delta Q)}^T$}, $F$ определяется по формуле (2.15), а $(\cdot ,\cdot )$ обозначает скалярное произведение в $L_2\left(D\right)\times L_2(\Omega \times (\mathrm{0,}\overline{t}\left)\right)$.

Уравнение для $\delta \Theta$ определяется формулой (2.13), отсюда

$(\mathcal{F},\delta \Theta )=(\mathcal{F},{\mathcal{H}}^{-1}\mathcal{P}\delta Q^{\left(0\right)})={({\mathcal{P}}^{*}{\mathcal{H}}^{-1}\mathcal{F},\delta Q^{\left(0\right)})}_{Y_{obs}}.$ (2.20)

Таким образом, из (2.17)–(2.20) заключаем, что

${\left(\frac{dG}{d{Q}^{\left(0\right)}},\delta Q^{\left(0\right)}\right)}_{Y_{obs}}={({\mathcal{P}}^{*}{\mathcal{H}}^{-1}\mathcal{F},\delta Q^{\left(0\right)})}_{Y_{obs}},$ (2.21)

откуда следует утверждение теоремы.

Пусть $\Theta ={\mathcal{H}}^{-1}\mathcal{F}$, тогда из (2.14) имеем

$\frac{dG}{d{Q}^{\left(0\right)}}={\mathcal{P}}^{*}\Theta .$(2.22)

Нетрудно убедиться в том, что сопряженный оператор ${\mathcal{P}}^{*}:L_2\left(D\right)\times L_2(\Omega \times (\mathrm{0,}\overline{t}\left)\right)\to L_2(\Omega \times (\mathrm{0,}\overline{t}\left)\right)$ определяется на ^{$\Theta ={(\xi ,\eta )}^T$}, $\xi \in L_2\left(D\right),\eta \in L_2(\Omega \times (\mathrm{0,}\overline{t}\left)\right)$ по формуле

${\mathcal{P}}^{*}\Theta ={\mathcal{B}}^{-1}\eta .$ (2.23)

Тогда

$\frac{dG}{d{Q}^{\left(0\right)}}={\mathcal{P}}^{*}\Theta ={\mathcal{B}}^{-1}\eta .$ (2.24)

Из (2.22)–(2.24) и (2.14)–(2.16) заключаем, что тем самым доказана следующая

Теорема 2. Градиент $dG/d{Q}^{\left(0\right)}$ определяется последовательным выполнением следующих шагов:

1) решить сопряженную задачу

$-{\left({\varphi }^{*}\right)}_t+L^{*}{\varphi }^{*}=\frac{\partial G}{\partial T},t\in \left(\mathrm{0,}\overline{t}\right),$

${\varphi }^{\ast }=0$ при $t=\overline{t},$ (2.25)

полагая

^{$\mathcal{F}={\left({\varphi }^{*}{|}_{t=0}+\frac{\partial G}{\partial T_0},B^{*}{\varphi }^{*}+\frac{\partial G}{\partial Q}\right)}^T;$}

2) найти ^{$\Theta ={(\xi ,\eta )}^T$} как решение уравнения с гессианом

$\mathcal{H}\Theta =\mathcal{F};$

3) вычислить градиент функционала по формуле

$\frac{dG}{d{Q}^{\left(0\right)}}={\mathcal{B}}^{-1}\eta .$ (2.26)

В численных примерах, следуя [18], рассматривались функционалы вида

$G(T,T_0,Q)=\underset{0}{\overset{\overline{t}}{\int }}dt\underset{\Omega }{\int }{F}^{*}(x,y,t)T(x,y\mathrm{,0,}t)d\Omega ,$ (2.27)

где $F^{*}(x,y,t)$ — весовая функция, связанная с полем температуры на поверхности $z=0$. Для определения средней температуры в избранной акватории океана $\omega$ при $z=0$ в интервале $t_1-\tau \le t\le t_1$ в качестве $F^{*}$ выбирается функция

$F^{*}(x,y,t)=\left\{\begin{array}{cc}1/\left(\tau mes\omega \right),& если(x,y)\in \omega ,t_1-\tau \le t\le t_1,\\ 0& в противном случае,\end{array}\right.$ (2.28)

где mes w означает площадь района w. В этом случае функционал (2.27) представляется в виде

$G(T,T_0,Q)=\frac{1}{\tau }\underset{t_1-\tau }{\overset{t_1}{\int }}dt\times \left(\frac{1}{mes\omega }\underset{\omega }{\int }T(x,y\mathrm{,0,}t)d\Omega \right).$ (2.29)

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

$G(T,T_0,Q)=\underset{0}{\overset{\overline{t}}{\int }}(B{F}^{*},T)dt={\left(B{F}^{*},T\right)}_Y,Y=\left(L_{2}D\times \left(\mathrm{0,}\overline{t}\right)\right).$

В силу

${(\frac{\partial G}{\partial T},\delta T)}_Y={(B{F}^{*},\delta T)}_Y$

производные от $G$ по $T,T_0,Q$, входящие в (2.25), определяются по формулам

$\frac{\partial G}{\partial T}=B{F}^{*},\frac{\partial G}{\partial T_0}=\frac{\partial G}{\partial Q}=0.$ (2.30)

Алгоритм (2.25), (2.26) с учетом конкретного вида производных (2.30) использовался при численных расчетах для оценки чувствительности функционалов, связанных с температурой после вариационного усвоения, по отношению к изменениям функции входных данных о потоке тепла $Q^{\left(0\right)}$.

3. Результаты численных экспериментов

Для численных экспериментов использовалась трехмерная численная модель гидротермодинамики Черного моря, разработанная в ИВМ РАН на основе метода расщепления (см. [23]) и дополненная процедурой усвоения температуры поверхности моря (ТПМ) для восстановления начального состояния $T_0$ и тепловых потоков Q с учетом ковариационных матриц ошибок наблюдений и ошибок бэкграунда.

Параметры рассматриваемой области Черного моря и ее географические координаты можно описать следующим образом: s-сетка $306\times 200\times 27$ (широта, долгота и глубина соответственно). Первая точка “сетки C” (см. [28]) имеет координаты ${26.65}^{\circ }$ E и ${40.15}^{\circ }$ N. Шаги сетки по x и y постоянны и равны 0.05 и 0.036 градуса соответственно. Шаг по времени равен Dt = 2.5 мин.

В качестве данных наблюдений $T_{obs}$ в данном эксперименте были выбраны данные ТПМ спектрометра MODIS на спутнике Aqua, предоставленные спутниковой службой ‘See the Sea’, входящей в состав ЦКП “ИКИ Мониторинг” (см. [29]). Данные ТПМ были пересчитаны на сетку численной модели (см. [30]).

Значения среднего климатического теплового потока $Q^{\left(0\right)}$ использовались в процедуре усвоения данных в качестве начального приближения (бэкграунда). Поток тепла на поверхности моря рассчитан по данным реанализа Era 5 за период с 1979 по 2021 г. (см. [31]). По данным за 1979–2021 гг. рассчитаны средние значения и дисперсии теплового потока по дневным и ночным данным для каждого дня года. Полученные дисперсии представляют собой диагональные элементы ковариационной матрицы ошибок бэкграунда B. Аналогичным образом (см. [32]) на основе данных ТПМ сервиса Copernicus (data.marine.copernicus.eu) за период с 1982 по 2019 г. рассчитывались элементы ковариационной матрицы ошибок данных наблюдений R.

Используя модель гидротермодинамики Черного моря, дополненную процедурой усвоения температуры поверхности $T_{obs}$, были проведены численные расчеты на исследуемой акватории. Алгоритм усвоения работал лишь в некоторые моменты времени $t_k$, при этом $t_{k+1}=\overline{t}=t_k+\Delta t$. При реализации процедуры усвоения на одном шаге по времени $(t_k,t_{k+1})$ рассматривалась система вида (1.5)–(1.8) с диагональными ковариационными матрицами $\mathcal{B}$ и $R$. На каждом из таких шагов был проведен расчет чувствительности рассматриваемого функционала (2.27) к входным данным о потоке тепла согласно алгоритму (2.25), (2.26).

На фиг. 1 представлены диагональные элементы матриц B (фиг. 1а) и R (фиг. 1б) на 27 февраля 2021 г. 10 ч 35 мин. На фиг. 2a приведены данные наблюдений, полученные со спутника, на момент расчета чувствительности функционала от решения задачи. Пример рассчитанного потока тепла на 27 февраля 2021 г. приведен на фиг. 2б.

Фиг. 1. Значения диагональных элементов матриц B (а) и R (б) в численном эксперименте.

Фиг. 2. Данные, используемые при расчетах: (а) — температура поверхности моря, °C; (б) — поток тепла на поверхности, Вт/м².

Разработанный алгоритм (2.25), (2.26) позволяет определять области наибольшего градиента функции отклика как средней по интервалу температуры поверхности океана для избранного района, полученной в результате процедуры усвоения. Результаты расчета градиента функционала от решения задачи по отношению к входным данным о потоке тепла согласно (2.25), (2.26) для 27 февраля 2021 г. показаны на фиг. 3а. Из фиг. 3а следует, что наибольшая чувствительность функционала к ошибкам входным данных о потоке тепла с максимальными значениями градиента порядка $1.0\times {10}^{-6}$ прослеживается у побережья Черного моря, особенно в южной части моря. Наименьшая чувствительность к ошибкам входным данных наблюдается в центральных районах Черного моря в области с наибольшей глубиной, и здесь градиент составляет $0.5\times {10}^{-9}$, что на 3 порядка меньше максимального значения.

Фиг. 3. Градиент функционала G(T,T₀,Q): (а) — dG/dQ⁽⁰⁾; (б) — dG/dT_obs.

Результаты расчета градиента функционала от решения задачи по отношению к ошибкам данным наблюдений согласно алгоритму из [24] для 27 февраля 2021 г. приведены на фиг. 3б. Из фиг. 3б следует, что наиболее мелкая западная часть Черного моря наиболее чувствительна к этим ошибкам c максимальным значением градиента $5.0\times {10}^{-3}$. Отметим, что южная и восточная части побережья Черного моря также чувствительны к ошибкам данных наблюдений. Центральная часть Черного моря, характеризующаяся наибольшей глубиной, согласно проведенным расчетам, является, наоборот, наименее чувствительной областью. Здесь значения градиента функционала менее $5.0\times {10}^{-5}$, т.е. на 2 порядка меньше, чем в областях с большей чувтствительностью.

Из фиг. 3а и 3б следует, что в целом чувствительность рассматриваемого функционала к данным наблюдений ТПМ на 4 порядка больше, чем чувствительность к входным данным о потоке тепла.

Эти результаты подтверждаются прямыми вычислениями функционала $G(T,T_0,Q)$ в соответствии с (29), полученного после вариационного усвоения, путем введения возмущений в данные наблюдений $T_{obs}$ и в данные о потоке тепла $Q^{\left(0\right)}$, следуя работе [24].

На фиг. 4 представлены результаты применения процедуры усвоения данных для модели динамики Черного моря на 27 февраля 2021 г. 10 ч 35 мин. На фиг. 4а приведена разность между ТПМ, полученной при расчете по модели без использования усвоения, и данными наблюдений, а на фиг. 4б с использованием процедуры усвоения данных наблюдений. Согласно этим результатам можно сделать вывод о том, что использование усвоения данных позволяет получать более близкие к данным наблюдений поля ТПМ, что несомненно улучшает как реанализ расчетных полей, так и прогностические свойства модели.

Фиг. 4. Отклонение рассчитанных значений ТПМ от данных наблюдений, °C: (а) — для модели без усвоения и данных наблюдений; (б) — для модели с усвоением и данными наблюдений.

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

Заключение

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

Авторы благодарны рецензенту за полезные замечания, которые позволили улучшить представление результатов и изложение материала в статье.

About the authors

В. П. Шутяев

Author for correspondence.
Email: victor.shutyaev@mail.ru
Russian Federation, 119333 Москва, ул. Губкина, 8

Е. И. Пармузин

Email: e.parmuzin@inm.ras.ru
Russian Federation, 119333 Москва, ул. Губкина, 8

References

Supplementary files

Supplementary Files

Action

1. JATS XML

Download

2. Fig. 1. Values of diagonal elements of matrices B (a) and R (b) in a numerical experiment.

Download (27KB)

Indexing metadata

3. Fig. 2. Data used in calculations: (a) — sea surface temperature, °C; (b) — heat flux on the surface, W/m2.

Download (26KB)

Indexing metadata

4. Fig. 3. Gradient of the functional G(T,T0,Q): (a) — dG/dQ(0); (b) — dG/dTobs.

Download (31KB)

Indexing metadata

5. Fig. 4. Deviation of calculated SST values from observational data, °C: (a) - for the model without assimilation and observational data; (b) - for a model with assimilation and observational data.

Download (30KB)

Indexing metadata