Numerical solution of the inverse problem of identifying methane consumption rate in soil

Full Text

ВВЕДЕНИЕ

Основным (около 90 % от общего) стоком метана из атмосферы Земли является окисление гидроксил-радикалом в тропосфере. Вторым по важности (около 10 %) считается потребление метана микроорганизмами – метанотрофами, живущими в автоморфных (непереувлажнённых) почвах. Таким образом, изучение величин потребления (удельных потоков) СН₄, понимание процессов, обуславливающих его временную и пространственную динамику, а также моделирование потребления необходимы для построения обоснованных климатических прогнозов. Потребление метана в почве за счет окисления метанотрофными бактериями в автоморфных почвах – единственный известный биологический механизм стока для атмосферного метана (Prather, Holmes, 2017 [1]).

Вместе с эмпирическими методами для количественной оценки концентрации метана в почве применяется математическое моделирование. Существуют модели для оценки процессов поглощения метана почвой, такие как P96, R99 и C07. Модель P96 основывается на измененной версии закона Фика и оценивает диффузионный поток атмосферного CH₄ в почву. Модель R99 учитывает микробное окисление CH₄ в почве, устанавливая зависимость скорости окисления от параметров почвы, таких как температура, влажность и содержание азота. Модель C07 модифицирует этот подход, вводя аналитическое решение уравнения диффузии-реакции и скалярный модификатор для учета влияния влажности почвы и температуры на скорость окисления CH₄. Современная модель MeMo v1.0 (см. [2] и библиографию в этой работе), которую мы также используем ниже, является новой версией моделей по оценке глобального поглощения атмосферного метана почвой. В отличие от моделей R99 и C07, MeMo использует общее аналитическое решение одномерного уравнения диффузии-реакции. Данная модель позволяет учитывать поток метана из глубины почвы, который участвует в процессе метанотрофии в почве. Уравнение записывается в виде [2] (см. также стационарный случай в [1, 4, 3]):

$\begin{array}{c}{\partial }_{t}C={\partial }_z\left(a_z\left(z\right)C_z\right)-V\left(z\right)C=\mathrm{0,}C\left(0\right)=C_a,C_z{|}_{z=\beta }=0\left(z\in \left(\mathrm{0,}\beta \right)\right). \end{array}$

Здесь a_z – коэффициент диффузии CH₄ в порах почвы (м³ воздуха м^-1 почвы ч^-1), вычисляемый с помощью измеренных значений Φ (общая пористость) и Θ (объемная влажность, м³/м³), β – максимальное значение глубины, C – функция концентрации метана, C_a – измеренная концентрация CH₄ в атмосфере (мг /м³), V(z) – измеренная или найденная тем или иным способом скорость потребления метана при C=C_a (м^-3 воздуха м^-3 почвы ч^-1).

Таким образом, задача определения коэффициента V(z) – задача определения младшего коэффициента в эллиптическом или параболическом уравнении. Отметим, что задача определения младшего коэффициента в параболическом уравнении известна давно, и имеется большое количество работ, посвященных численному исследованию задачи. Основное внимание было уделено стационарному случаю, который наиболее часто использовался в приложениях. В стационарном случае уравнение и граничные условия на верхней z=0 и нижней (z=β) границах записываются таким образом:

$\begin{array}{c}{\partial }_z\left(D_z\left(z\right){\left(C_z\right)}_z\right)-V\left(z\right)C_z=\mathrm{0,}{C}_z\left(0\right)=C_a,{\left(C_z\right)}_z{\left.\text{}\right|}_{z=\beta }=0. \left(1\right)\end{array}$

Дополнительные условия имеют вид:

$\begin{array}{c}C\left(b_i\right)={\alpha }_i,i=\mathrm{1,2,}\dots ,r,b_i\in \left(\mathrm{0,}Z\right). \left(2\right)\end{array}$

Удельный поток (УП) метана Q (мг/м²) в час, прогнозируемый моделью, вычислялся по полученному решению:

$Q=D_z\left(\mathrm{0,005}\right)\frac{C_z\left(0.01\right)-C_z\left(0\right)}{\mathrm{0,01}}.$

Метан вносит значительный вклад в глобальное потепление, удерживая тепло в атмосфере и увеличивая температуру на поверхности планеты [5]. Количественные оценки выбросов метана в атмосферу в различных литературных источниках отличаются. Однако их можно классифицировать по двум основным категориям: антропогенные источники метана и природные источники метана. Антропогенные источники метана связаны непосредственно с деятельностью человека, к ним можно отнести: а) энергетическую и промышленную деятельность: добычу, транспортировку и переработку ископаемых; б) животноводство: ферментация пищи рогатым скотом способствует выделению метана; в) сельское хозяйство (рисовые плантации): бактериальное разложение органического материала на рисовых полях с выделением метана; г) полигоны твердых бытовых отходов: разложение органических отходов с выделением биогаза, состоящего преимущественно из метана. Суммарно антропогенные источники метана выделяют порядка 440,5 млн т в год [6, 7]. Существует множество практик, позволяющих уменьшить эмиссию метана, вызванную антропогенными источниками, некоторыми из них можно считать современные технологии утилизации газовых выбросов, понижение расхода топливного газа, биогазовые установки для переработки отходов животноводства и пищевой промышленности. Природные источники метана выделяют значительное количество метана в атмосферу. Основными природными источниками метана являются болотные системы, на них приходится до 30 % от общего объема выбросов метана в атмосферу [8], [9]. Сошлемся на работы, посвященные в том числе и использованию математических моделей в задачах описания эмиссии метана [10–17].

В работе [3] задача (1), (2) в стационарном случае решалась численно: функция концентрации метана определялась из уравнения (1) с помощью специальной сглаживающей функции, которая была получена при помощи регуляризации на основе интерполяции значений концентрации в известных точках. Мы рассмотрели тот же алгоритм и эту же задачу. Однако в отличие от работы [3] число замеров концентраций r у нас произвольно и может быть небольшим. В методе из работы [3] предполагалось, что концентрация метана замерена в достаточно большом количестве точек, что позволяет фактически точно с помощью интерполяции и с последующим сглаживанием с помощью аппроксимации получить функцию концентрации метана. Однако это не всегда возможно. Основные рассмотрения мы проводили в случае, когда замеры концентраций сделаны в достаточно небольшом количестве точек, и строили коэффициент поглощения в виде конечного отрезка ряда (в нашем случае – в виде разложения по базисным кусочно-линейным функциям метода конечных элементов). Была предложена другая, на наш взгляд, более совершенная схема решения задачи, которая сводилась к некоторому нелинейному уравнению и затем решалась методом наискорейшего спуска. Получены хорошие результаты, в частности, точности определения коэффициента поглощения.

Описание алгоритма

Рассмотрим задачу (1)–(2). Для численного решения используем метод конечных элементов. Пусть h=β/N – шаг по пространству и {φ_i } – набор базисных функций. Ищем приближенное решение уравнения (1) в виде C=u^N+c⁰

$u^N=\sum _{}_{k=1}^N{C}_k{\phi }_k\left(x\right),C_0=c^0,$

где C_i определяем из системы

(a(z) (u^N)_z ,φ_lz)+(V(z)u^N,φ_l )=-(c⁰V,φ_l ), l=1,2,…,N.

Здесь скобки обозначают скалярное произведение в L₂(G), т. е. (u,v)=∫₀^β u(z)v(z)dz. Базисные функции φ_i (i=1,2,…,N) определяются следующим образом:

${\phi }_i\left(z\right)=\left\{\begin{array}{c}\frac{\left(h\left(i+1\right)-z\right)}{h},x\in \left(hi,h\left(i+1\right)\right)\\ \frac{z-h\left(i-1\right))}{h},z\in \left(h\left(i-1\right),hi\right))\\ \mathrm{0,}z\notin \left(h\left(i-1\right),h\left(i+1\right)\right)\end{array}\right.,$

${\phi }_0\left(z\right)=\left\{\begin{array}{c}\frac{h-z}{h},z\in \left(\mathrm{0,}h\right))\\ \mathrm{0,}z\notin \left(\mathrm{0,}h\right)\end{array}\right.$

${\phi }_N\left(z\right)=\left\{\begin{array}{c}\frac{z-\left(N-1\right)h}{h},z\in \left(\left(N-1\right)h,Nh\right))\\ \mathrm{0,}z\notin \left(\left(N-1\right)h,Nh\right).\end{array}\right.$

Систему можно представить в матричном виде:

$\begin{array}{c}A\overrightarrow{C}=\overrightarrow{F}. \left(3\right)\end{array}$

Необходимо найти вектор C⃗, выразим данный вектор из системы (3):

$\begin{array}{c}\overrightarrow{C}=A^{-1}\overrightarrow{F},\end{array}$

где C⃗=(C₁,C₂,…,C_N ), а вектор F⃗ имеет следующий вид:

F⃗^T=(V(z)c⁰,φ₁ ),(V(z)c⁰,φ₂ ),…,(V(z)c⁰,φ_N-1 ),(V(z)c⁰,φ_N ).

Матрица A – матрица жесткости с элементами:

$\begin{array}{c}{\text{A}}_{\text{ij}}=\left(\text{a}\left(\text{z}\right){\text{φ}}_{\text{jz}},{\text{φ}}_{\text{iz}}\right)+\left(\text{V}\left(\text{z}\right){\text{φ}}_{\text{j}},{\text{φ}}_{\text{i}}\right),j=\mathrm{1,2,}\dots ,N,i=\mathrm{1,2,}\dots ,N.\end{array}$

Матрицу A можно представить в виде суммы двух матриц: A=A₀+A_V, где элементы этих матриц A₀ есть числа (a(z)φ_jz,φ_iz). Приближенно элементы матрицы A_0 можно записать в следующем виде:

$\begin{array}{c}\begin{array}{c}\begin{array}{c}{R}_{ij}=\left\{\begin{array}{c}\frac{1}{h^2}{\int }_{\left(i-1\right)h}^{\left(i+1\right)h}a\left(z\right)dz,при i=j\\ -\frac{1}{h^2}{\int }_{ih}^{\left(i+1\right)h}a\left(z\right)dz,при\left|i-j\right|=1\\ \mathrm{0,}при\left|i-j\right|\ge 2\end{array}\right.\\ i=\mathrm{1,2,}\dots ,N-\mathrm{1,}j=\mathrm{1,2,}\dots ,N,\end{array}\\ R_{NN}=\frac{1}{h^2}{\int }_{\left(N-1\right)h}^{Nh}a\left(z\right)dz.\end{array}\end{array}$

Элементы матрицы A_v вычисляются путем нахождения скалярного произведения:

$\begin{array}{c}\begin{array}{c}{A}_v=\left(V\left(z\right){\phi }_j,{\phi }_i\right)\\ i=\mathrm{1,2,}\dots ,N,j=\mathrm{1,2,}\dots ,N\end{array} \end{array}$

Обращая матрицу A, найдем вектор C⃗. Таким образом построим решение прямой задачи.

При решении обратной задачи ищем функцию V(z) в виде V(z)=∑^r_(j=1) B⃗_j ψ_j(z), где ψ_j – некоторый набор базисных функций, τ=β/r – шаг по пространству. С шагом τ построим кусочно-линейные функции ψ_i , используемые в методе конечных элементов. Обозначим узлы сетки через z_i=hi, zⁱ=τi. Пусть узлы z_ik k=1,2,…,r отвечают точкам замеров, т. е. z_ik≈z^k. Приближенно систему можно записать в виде:

-C⃗/c⁰=(A₀+A_v )^-1F₀ , F₀=(F₁,F₂,…,F_N )^T, F_i=(V(z),φ_i).

Обозначим правую часть через g(β⃗)=(A₀+A_v)^-1F₀. F⃗₀=((V(z),φ₁),…,(V(z),φ_N))^T. Если мы возьмем строчки этой системы с номерами i_k, k=1,2,…,r, то получим систему

$\begin{array}{c}\overrightarrow{\gamma }=P_{0}g\left(\overrightarrow{\beta }\right),{\gamma }_j=1-\frac{{\alpha }_j}{c_0}>\mathrm{0,}\left(4\right)\end{array}$

где P₀ g – вектор длины r с координатами g_i1,…,g_ir и индекс (i_j) – взятие координаты с индексом(i_j). Пусть $\overrightarrow{g}=\left(g_1,\dots ,g_N\right)$. Найдем производные ${\partial }_{{\beta }_l}\overrightarrow{g}\left(\overrightarrow{\beta }\right)$. Имеем $\left(A_0+A_V\right)\overrightarrow{g}=F_0$. Найдем частные производные по β_l

$\begin{array}{c}\frac{\partial }{\partial {\beta }_l}\left(A_0+A_v\right)\overrightarrow{g}\left(\overrightarrow{\beta }\right)=A_0{\overrightarrow{g}}_{{\beta }_l}+{\partial }_{{\beta }_l}{A}_v\overrightarrow{g}={\partial }_{{\beta }_l}{F}_0.\end{array}$

Элементы матрицы ∂_βl A_v и производные F_iβl имеют вид

${\partial }_{{\beta }_l}{a}_{ij}^v={\partial }_{{\beta }_l}\left(V{\phi }_j,{\phi }_i\right)=\left({\psi }_l{\phi }_j,{\phi }_i\right),иF_{i{\beta }_l}=\left({\psi }_l,{\phi }_i\right).$

Далее имеем

$\left(A_0+A_v\right){\overrightarrow{g}}_{{\beta }_l}=F_2={\partial }_{{\beta }_l}{F}_0-\left({\partial }_{{\beta }_l}{A}_v\right)\overrightarrow{g},$

или

g_kβl=((A₀+A_v )^-1F₂ )_ik .

Обозначим матрицу {g_kβl} через J(β⃗), а транспонированную матрицу через J^*. Ищем решение системы (4). Обозначим матрицу системы через A(β⃗), т. е. решаемая система имеет вид A(β⃗)=γ⃗. Итерационная процедура записывается в виде

β⃗_i+1=β⃗_i-τ_iJ'(β⃗_i), J (β⃗)=J^*(β⃗)(A(β⃗)-γ⃗),

где параметр τ_i определяется так:

τ_i = ∥ J' (β⃗_i) ∥²∥ JJ^*(A(β⃗_i )-γ⃗ ) ∥².

Здесь под нормой вектора понимаем квадратный корень из суммы квадратов координат. В качестве начального приближения стоит взять вектор β⃗₀ с положительными координатами, например, β⃗₀=(1,1,…,1).

РЕЗУЛЬТАТЫ И ОБСУЖДЕНИЕ

В этом разделе проанализируем результаты численных экспериментов, полученных в результате решения прямой и обратной задач в рамках экспериментов на ЭВМ. Характеристики испытуемой ЭВМ: процессор: Intel i5 9500; ОЗУ: 16,00 ГБ; тип системы: Windows 11. Численные эксперименты были проведены в программном комплексе MATLAB 2024a.

В результате расчётов было восстановлено исходное значение функции в точках потребления метана по глубине на основе распределения концентрации метана по глубине. В качестве входных данных были взяты значения функции в точках замеров. Значения одномерного массива метана были получены в ходе решения прямой задачи согласно системе (3), где входным вектором скорости потребления метана выступал некоторый экспериментальный одномерный массив данных [3], представленный на рисунке 1.

 

Рисунок 1. Форма тестового вектора скоростей для N = 100

 

В эксперименте использовались следующие глобальные параметры:

1. a(z)=1+z – функция молекулярной диффузии.
2. c₀=1 – значение концентрации метана на поверхности почвы.
3. N=100 – количество узлов сетки.

Для данного вектора скорости потребления метана были получены значения функции концентрации метана в точках, представленные на рисунке 2.

 

Рисунок 2. Значения функции концентрации метана C⃗(z)

 

Для нахождения точности полученных результатов были применены следующие оценки: MAE и MAPE.

1. $MAE=\frac{1}{N}\sum _{}_{i=1}^N|y_i-{\widehat{y}}_i|$ – средняя абсолютная ошибка: позволяет оценить качество полученного решения в абсолютных единицах, где y является эталонным решением, ŷ является полученным решением.
2. $MAPE=\frac{100}{N}\sum _{}_{i=1}^N\left|\frac{y_i-{\widehat{y}}_i}{y_i}\right|$ – средняя относительная ошибка: позволяет оценить качество полученного решения в процентном соотношении, где y является эталонным решением, ŷ является полученным решением.

В ходе решения обратной задачи был получен вектор скорости потребления метана, который в среднем имеет отклонение от исходного на 0,016768 (MAE). При этом полученный вектор скорости имеет ошибку 0,979798 %. Полученное решение обратной задачи представлено на рисунке 3.

 

Рисунок 3. Графики полученных решений в ходе обратной задачи

 

Покажем отношения восстановленного вектора скорости к исходному на рисунке 4.

 

Рисунок 4. Сопоставление скоростей потребления метана

 

Результаты восстановления значений вектора зависимости скорости потребления метана от глубины удовлетворительны. Ряд восстановленных значений примерно соответствует эталонным значениям, однако небольшое количество точек значимо отклоняется от эталонного значения, что может быть обусловлено используемыми базисными функциями. Для следующего эксперимента увеличим количество точек, для которых необходимо получить значения вектора концентрации C⃗. Количество замеренных точек r=100, количество точек, для которых необходимо получить значения концентрации, N = 200. Обратная задача решается для r точек. Прямая задача решается для N точек по полученному во время решения обратной задачи вектору B⃗. Результат решения представлен на рисунке 5.

 

Рисунок 5. Результаты решения для r = 100, N = 200

 

При интерполяции исходной функции с увеличением количества значений функции в два раза незначительно изменилась ошибка: MAE = 0,016834, MAPE = 1,299812 %, при этом значение функции зависимости концентрации метана от глубины в точках отклонилось от эталонного решения примерно на 2,280873 %.

Оценим сходимость алгоритма решения обратной задачи к глобальному минимуму, где для каждого эксперимента определим случайный одномерный массив B⃗, равномерно распределённый в диапазоне от 0 до 1 размерностью N=20, и установим минимальный шаг сходимости ϵ=10^-5. В ходе проведения серии экспериментов были получены оценки точности работы алгоритма решения обратной задачи: MAE для одномерного массива C⃗, MAE для вектора V⃗(z), MAPE для одномерного массива V⃗(z). Результаты численных экспериментов оценки сходимости алгоритма решения обратной задачи представлены в таблице 1.

 

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

MAE ($\overrightarrow{C}$)	MAE(V(z))	MAPE(V(z))
0,000602	0,263407	15,306
0,000756	0,252482	14,42567
0,000557	0,276899	16,81246
0,00078	0,283667	16,88795
0,00059	0,26605	15,61719
0,000628	0,288329	17,56511
0,000601	0,293131	18,00528
0,000525	0,260345	15,28794
0,000561	0,262915	15,3175
0,000536	0,281399	17,06343
0,000601	0,290857	18,03111
0,000584	0,295519	18,4257
0,000978	0,288084	17,23596
0,0007	0,306291	19,07516
0,000813	0,282567	16,76015
0,000676	0,269621	15,94521
0,000531	0,253104	14,26699
0,000479	0,278698	16,86203
0,00043	0,26466	15,54362
0,000501	0,260129	14,91595
Математическое ожидание
0,00062145	0,275908	16,46752
	Дисперсия
1,70475E-08	0,000233	1,866449

 

Приведенные в таблице 1 результаты показывают, что для случайно заданного стартового одномерного массива B⃗ ожидаемое решение будет иметь ошибку 16,46752 %, ожидаемое отклонение от данной ошибки не превысит 1,866449 %.

ЗАКЛЮЧЕНИЕ И ВЫВОДЫ

Численный алгоритм основан на методе наискорейшего градиентного спуска. Алгоритм решения обратной задачи стабильно сходится к точке глобального минимума с ожидаемым отклонением не более ≈ 1,9 % при минимальном размере шага алгоритма и малом количестве точек замера. В ходе ряда экспериментов была установлена сходимость алгоритма. Алгоритм можно считать устойчивым: при различных случайных значениях одномерного массива B⃗ алгоритм сходится к одному решению.

About the authors

Julia A. Tukmacheva

Yugra State University

Email: y_tukmacheva@ugrasu.ru

Postgraduate Student, Engineering School of Digital Technologies

Russian Federation, Khanty-Mansiysk

Sergey G. Pyatkov

Yugra State University

Author for correspondence.
Email: s_pyatkov@ugrasu.ru

Doctor of Physical and Mathematical Sciences, Professor Engineering School of Digital Technologies

Russian Federation, Khanty-Mansiysk

References

Supplementary files