Recovery of Hydrochemical Parameters of the Streeter – Phelps – Shishkin Model Using the Nelder – Mead Method

全文:

ВВЕДЕНИЕ

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

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

В качестве математической модели для описания распространения загрязнения используется модель Стритера – Фелпса, модифицированная в работах [1, 2]. В [1] её предложено именовать моделью Стритера – Фелпса – Шишкина.

Эта модель позволяет описывать распространение загрязнений в малых реках – реках, где длина исследуемого участка значительно превышает размеры поперечного сечения, что приводит к предположению о равномерной концентрации загрязнения по всему сечению.

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

В [3] данная модель была протестирована in vivo в Или-Балхашском бассейне, что подтвердило её применимость для моделирования реальных условий.

Для упрощения обозначений далее будем считать, что площадь сечения реки равна 1, скорость течения также равна 1, благодаря чему время x численно равно расстоянию вдоль реки. Кроме того, загрязняющее вещество будет измеряться в кислородном эквиваленте, т. е. одна условная единица вещества требует для разложения одну условную единицу кислорода.

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

Модель описывается системой дифференциальных уравнений

$$\begin{gathered} \frac{d{y}_1}{dt}=-k{y}_1{y}_2+q\delta \left(x-x_0\right), \\ \frac{d{y}_2}{dt}=-k{y}_1{y}_2+a\left(Y-y_2\right), \end{gathered}$$ (1)

где х – время и одновременно расстояние, функции y₁ (x) и y₂ (x) – уровень загрязнения и кислорода соответственно. Параметр k представляет собой коэффициент скорости окисления, который зависит как от свойств загрязняющего вещества, так и от погодных условий. В работе [4] описаны экспериментальные подходы к его определению. Параметр Y обозначает максимальную концентрацию кислорода в чистой воде, зависящую от температуры воды. Параметр a – коэффициент аэрации, отражает скорость поступления кислорода в воду из атмосферы и также определяется погодными условиями. Слагаемое q∙δ(x–x₀), где δ – дельта-функция, определяет точечный сброс загрязняющего вещества в точке x₀ с интенсивностью q.

Получить аналитическое решение этих уравнений сложно, однако возможно использование численных методов в MATLAB. Ниже приведен код, задающий систему уравнений:

%% Диф. уравнение типа Стритера–Фелпса

function dy = sf(y,k,a, Y)

dy = zeros(2,1); % пока пустая матрица с производными

dy(1) = - k*y(1)*y(2); % производная загрязнения

dy(2) = - k*y(1)*y(2) + a * (Y - y(2); % производная кислорода

end

Для моделирования дельта-функции удобно разделить исследуемый участок реки, идущий от начальной точки (x=0) до конечной (x=l), на два отрезка: до x₀ и после x₀. Численное решение системы дифференциальных уравнений выполняется методом Рунге – Кутты (в MATLAB – функция ode45), после чего результаты для двух участков объединяются.

%% Решает уравнение sf со сбросом в x0 интенсивности q

function [x,y] = SolveSF(x0,q,k,a ,Y,l)

f = @(x,y) sf(y,k,a, Y) ;% Сокращенное имя для удобства вызова

[x,y] = ode45(f , [x0 , l] , [ q , Y] ); % Решение на участке после x0

x= [0 ; x0-0.1 ; x] ;% Добавление чистой воды до точки сброса

y= [0 , Y ; 0 , Y ; y] ;% Сборка значений концентраций

end

На рисунке 1 представлен результат моделирования процесса разложения загрязняющего вещества после точечного сброса в точке x₀ = 100 с концентрацией q = 0,5 и гидрохимическими параметрами a = 0,01, Y = 1, k = 0,01. Параметры реки и загрязнения выбраны именно такими для обеспечения наглядности графиков.

 

Рисунок 1. Распределение уровня кислорода и загрязнения при x₀ = 100, q = 0,5, a = 0,01, Y = 1, k = 0,01

 

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

Нижний график показывает уровень загрязнения. До точки х₀ = 100 он равен нулю, в точке х₀ происходит резкое увеличение концентрации до q = 0,5, после чего загрязнение постепенно снижается.

 

Рисунок 2. Распределение уровня кислорода при q = 0,5, q = 0,6 и q = 0,7

 

Изменяя параметры загрязнения и коэффициент аэрации, можно заметить, что распределение кислорода зависит от интенсивности сброса, скорости реакции и коэффициента аэрации следующим образом: увеличение концентрации загрязнения вполне предсказуемо приводит к падению уровня кислорода.

Увеличение коэффициента окисления закономерно усиливает крутизну графика на начальном этапе, а также за счёт более быстрого разложения загрязняющего вещества увеличивает крутизну на конечном участке.

Рост коэффициента аэрации почти не оказывает влияния на начальный участок графика, но ускоряет восстановление концентрации кислорода до нормального уровня на конечном этапе.

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

Особый интерес представляет возможность выявления места сброса x₀, что даст шанс быстро обнаружить и ликвидировать источник загрязнения. Параметр k, отражающий скорость реакции, вероятно, может помочь в определении типа загрязняющего вещества.

 

Рисунок 3. Распределение уровня кислорода при k = 0,01, k = 0,02 и k = 0,03

 

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

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

 

Рисунок 4. Распределение уровня кислорода при a = 0,01, a = 0,02 и a = 0,03

 

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

Математическая постановка сводится к решению обратной задачи для системы дифференциальных уравнений (1).

Итак, предположим, что известны данные замера уровня кислорода в нескольких точках вдоль реки (см. рисунок 5).

 

Рисунок 5. Данные замеров уровня кислорода

 

Задача состоит в том, чтобы определить параметры x₀, q, k, a, Y, обеспечивающие наилучшее совпадение смоделированного распределения кислорода с экспериментальными данными.

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

Определяется функция (далее – errfun) с таким алгоритмом:

1. На вход этой функции поступают параметры x₀, q, k, a, Y и «данные измерения» уровня кислорода V в точках, определенных в массиве P.
2. Численно решается система уравнений (1) с начальными условиями x₀, q, как было описано ранее. По найденному решению определяется уровень кислорода в заданных точках (P), соответствующий этим параметрам.
3. Функция возвращает сумму квадратов разностей между измеренными значениями кислорода (V) и вычисленными значениями (W) из предыдущего шага.

Пример реализации этой функции приведён ниже:

%% Вычисление квадрата отклонения

function e = errfun(P, V, x0, q, k, a, Y, l)

[x, y] = SolveSF(x0,q,k,a,Y,l);% Решение дифференциального уравнения

W = interp1(x, y(:,2), P);% Интерполяция значений кислорода в точках P

e = sum((V - W).^2);% Сумма квадратов отклонений

end

С помощью функции fminsearch, реализующей метод Нелдера – Мида, подбираются параметры x₀, q, k, a, Y, минимизирующие значение errfun, что позволяет приблизительно определить характеристики загрязнения.

Пример кода для решения обратной задачи:

%% Решение обратной задачи

er = @(z) errfun(P, V, z(1), z(2), z(3), z(4), z(5), l);

% Упрощение вызова errfun

rez = fminsearch(@(z) er(z), [x1, q1, k1, a1, Y1]);

% Запуск метода Н-М

% от начальных значений

% x1 , q1 , k1, a1 , Y1

x2 = rez(1); q2 = rez(2); k2 = rez(3);

a2 = rez(4); Y2 = rez(5);

% Извлечение результатов

В описанном ниже числовом эксперименте значения уровня кислорода V в точках P, использованные как «результаты измерений», соответствуют параметрам x₀ = 100, q = 0,5, k = 0,01, a = 0,01, Y = 1.

 

Рисунок 6. Распределение уровня кислорода при правильном решении обратной задачи

 

На рисунке 7 представлено решение обратной задачи для модельного примера. Для метода Нелдера – Мида в качестве стартовых значений параметров были установлены: x₁ = 50, q₁ = 0,7, k₁ = 0,02, Y₁ = 0,9, a₁ = 0,02.

 

Рисунок 7. Модельный пример и результат поиска после 145 итераций

 

Алгоритм завершил работу, когда значение целевой функции (errfun) стало меньше 10^-6, для этого потребовалось 145 итераций. При этом были получены следующие значения искомых параметров: x₀ = 100,1446, q = 0,5440, k = 0,0093, a = 0,0108, Y = 1,0002, что вполне приемлемо для практического применения.

При дальнейшем увеличении числа итераций до 328 были достигнуты более точные значения: x₀ = 100,000038936103, q = 0,50000970878046, k = 0,00999982269511033, a = 0,0100001902307945, Y = 0,999999960029293.

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

Модифицированная модель Стритера – Фелпса – Шишкина, описывающая распространение органических загрязнений в малых реках, позволяет численно моделировать распределение уровня кислорода. Анализ влияния параметров загрязнения (интенсивности сброса загрязняющего вещества, коэффициентов окисления и аэрации) показал их определяющую роль в формировании распределения уровня кислорода вдоль реки, что открывает возможности для оценки характеристик загрязнения по форме графиков уровня кислорода.

Численное решение обратной задачи методом Нелдера – Мида на основе данных замера уровня кислорода в достаточно большом количестве мест вдоль реки позволяет определять ключевые параметры загрязнения с приемлемой точностью.

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

Практическая значимость работы заключается в возможности создания системы быстрого выявления и ликвидации источников загрязнения малых рек, что способствует повышению экологической безопасности водных объектов.

Код MATLAB, использованный в вычислениях, можно найти по ссылке: https://disk.yandex.ru/d/C9S_9TtEIo8Bfw.

作者简介

Sergey Semenov

Yugra State University

编辑信件的主要联系方式.
Email: ssp@ugrasu.ru

Candidate of Physics and Mathematics, Associate Professor at the School of Digital Technologies Engineering

俄罗斯联邦, Khanty-Mansiysk

Maria Kurkina

Khanty-Mansiysk State Medical Academy

Email: mavi@inbox.ru

Candidate of Physics and Mathematics, Associate Professor at the Department of Physiology and Sports Medicine

俄罗斯联邦, Khanty-Mansiysk

Anton Finogenov

Yugra State University

Email: a_finogenov@ugrasu.ru

Candidate of Physics and Mathematics, Associate Professor at the School of Digital Technologies Engineering

俄罗斯联邦, Khanty-Mansiysk

参考

补充文件