Цифровая обработка многомерных сигналов the contourlet transform and its application to digital image restore

THE CONTOURLET TRANSFORM AND ITS APPLICATION TO DIGITAL IMAGE RESTORE

Sergeev E.

Yaroslavl State University

Currently, there are many different methods to restore noisy signals and images. Images are 2D-signals with the inherent geometric structure, which is the main feature of the visual information. Features are located along a smooth curve – contours, due to smooth surfaces of section of displayed objects. The use of separable or inseparable wavelet transform does not perform the processing of the images in the form of curved contours, boundaries, etc. In this paper we describe an algorithm for image reconstruction from noisy data based on contourlet transform. The new transform has the following properties [1-3]:

1) Multiresolution. The representation allow images to be successively approximated, from coarse to fine resolutions.

2) Localization. The basis elements in the representation localized in both the spatial and the frequency domains.

3) Critical sampling. The representation has form a basis, or a frame with small redundancy.

4) Directionality. The representation contains basis elements oriented at a variety of directions, much more than the few directions that are offered by separable wavelets.

5) Anisotropy. To capture smooth contours in images, the representation contains basis elements using a variety of elongated shapes with different aspect ratios.

Algorithm of image filtering is as follows:

1. 
   Computing of the contourlet transform images.
   
2. 
   Changing of the coefficients of contourlet transform for a particular rule. In this paper we used an algorithm based on a hard threshold processing transform coefficients.
   
3. 
   Computing of the inverse contourlet transform.
   

References

1. P.J. Burt and E. H. Adelson, “The Laplacian pyramid as a compact image code,” IEEE Trans. Comm., vol. COM-31, pp. 532-540, April 1983/

2. M.N. Do, M. Vetterli, “The contourlet transform: an efficient directional multiresolution image representation,” IEEE Trans. Image Proc., vol. 14, no. 12, December 2005.

3. R.H. Bamberger and M. J. T. Smith, “A filter bank for the directional decomposition of images: theory and design,” IEEE Trans. Signal Proc., vol. 40, no. 4, pp. 882-893, April 1992.



СИНТЕЗ БИОРТОГОНАЛЬНЫХ МНОГОМЕРНЫХ БАНКОВ ФИЛЬТРОВ С ПОМОЩЬЮ МЕТОДА ЛИФТИНГА

Чобану М.К.

Московский энергетический институт (технический университет)

Существуют различные методики синтеза многоскоростных систем для обработки многомерных (ММ) сигналов, включая применение преобразования МакКлеллана, структурные методы, полиномиальные методы. Методика лифтинга («подъема») является одним из полиномиальных подходов при синтезе таких систем. Она дает новый угол зрения для рассмотрения и изучения банков фильтров (БФ) и порождаемых ими вейвлетов и кратномасштабного анализа [1], [2]. Изначально причиной их создания было желание синтезировать меняющиеся во времени БФ с со свойством точного восстановления сигнала (ТВС), или так называемое второе поколение вейвлетов. Основная идея лифтинга заключается в простом соотношении между всеми БФ, имеющими один и тот же низкочастоный (или высокочастотный) фильтр [3]. Получаемые с помощью лифтинга реализации имеют лестничный вид.



Рис. 1.

В [1] предложен метод, основанный на лифтинге, для формирования БФ и вейвлетов для заданной решетки и заданного числа нулевых производных. Это построение включает в себя два шага подъема: предсказание и уточнение. Предсказывающий фильтр принадлежит к новому классу интерполяционных фильтров и синтезируется, используя специальный алгоритм для ММ полиномиальной интерполяции. Такая конструкция обладает многими преимуществами лифтинга, включая замещение при вычислении, возможность преобразования “целое число в целое число” (отсутствие ошибок квантования) и др.

Лифтинг-преобразование включает в себя три этапа: разделение (split), предсказание (predict) и обновление (update) (Рис. 1). Лифтинг-схему можно представить в виде классического 2-х канального банка фильтров со свойством точного восстановления сигнала.

По внешнему виду такое представление очень близко к полифазной реализации:



Полифазная матрица для банка синтеза равна присоединенной матрице по отношению к матрице (знак означает транспонирование, комплексное сопряжение и замену :



Исходя из этих матриц, можно найти традиционные НЧ и ВЧ фильтры анализа и синтеза , , , [4], [5]. Операторы P и U определяются так: пусть , тогда можно построить банк фильтров с первичными нулевыми моментами и N вторичными (дуальными) нулевыми моментами, если в качестве оператора предсказания выбрать интерполяционный фильтр порядка N, реализующий аппроксимацию в точке с координатой . Для шахматной матрицы децимации точка, в которой производится интерполяция (предсказание) отсчета, равна (см. Рис. 2).



Рис. 2. 2-D интерполяция различного порядка ([4])

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

Как показано в [5] должен быть многочленом с дробным сдвигом и иметь порядок (M – матрица децимации). В качестве такого многочлена можно выбрать сопряженный многочлен , деленный на 2. Тем самым основная задача синтеза банков фильтров с помощью метода лифтинга сводится к синтезу многочленов с заданными порядком и ММ дробным сдвигом. Если применить фильтр с дробным сдвигом к многочлену, дискретизированному на заданной решетке, получится тот же самый многочлен, но дискретизированный на исходной решетке со сдвигом . Фильтр является фильтром с ММ дробным сдвигом и имеет порядок , если (1)

Рассмотрим сдвиг многочлена , в -области. Тогда или

Следуя определению (1) результатом будет или

Тем самым должно выполняться равенство для всех . Верна следующая теорема (которая в [1] сформулирована только как предположение и доказана автором в [4]): фильтр является фильтром с дробным сдвигом и имеет порядок тогда и только тогда, когда его импульсная характеристика удовлетворяет условиям (1)

Данная теорема позволяет синтезировать ММ фильтры произвольной размерности и с произвольным ММ сдвигом и тем самым реализовать ММ многоскоростные системы с помощью метода лифтинга.

Пример. 2-D фильтры с переменным дробным сдвигом. Для синтеза многочлена с произвольным дробным сдвигом . применим данную теорему для порядка . Тогда импульсная характеристика будет равна



норма фильтров с дробным сдвигом

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

10 моментов

В этом случае ее норма равна



16 моментов

Если все 16 возможных моментов будут обеспечены, то и (2)

Оптимизированный случай (для 10 моментов)

Если использовать полученное ранее решение, можно попытаться улучшить его с точки зрения минимизации его нормы. Допустим, что заданы 10 моментов: (3)

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

После их подстановки новая норма будет равна

Все три зависимости норм показаны на Рис. 3. Можно заметить, что норма для случая, когда все 16 моментов заданы, для достаточно большой области значений больше, чем норма более простого и с меньшим числом моментов (10 вместо 16) решения, полученного с помощью ММ ряда Тейлора. Решение с минимизированной нормой дает лучший результат. Однако если сравнить c имульсной характеристикой оптимизированного решения (не приведена из-за ограниченности места), то можно заметить, что число вычислений (умножений, сложений) для решения, полученного с помощью ММ ряда Тейлора, значительно меньше числа вычислений для оптимизированного решения. При этом отличие их норм составляет не более 50%.



Рис. 3. l₂ норма для решения с 10 моментами, с 16 моментами и оптимизированное решение

(для 10 моментов)
Литература

1. 
   Kovaćevi J., Sweldens W. Wavelet families of increasing order in arbitrary dimensions // IEEE Trans. Image Proc.— 2000.— March. — Vol. 9, no. 3.— Pp. 480–496.
   
2. 
   Sweldens W. The lifting scheme: A construction of second generation wavelets // SIAM Journ. on Math. Anal. — 1998.— Mar.— Vol. 29, no. 2.— Pp. 511–546.
   
3. 
   Vetterli M., Herley C. Wavelets and filter banks: Theory and design // IEEE Trans. Acoust., Speech, and Signal Proc.— 1992.— Sep.— Vol. 40, no. 9.— Pp. 2207–223.
   
4. 
   Чобану М.К. Синтез основных элементов многомерных многоскоростных систем. Многомерные фильтры с дробным сдвигом // Сибирский журнал вычислительной математики. 2008. т.11, №2. с. 219-238.
   
5. 
   Tchobanou M. Design of basic elements of multidimensional multirate systems. Multidimensional filters with fractional shift// Numerical Analysis and Applications, 2008, v.1, 2, p.179-195.
   

SYNTHESIS OF BIORTHOGONAL MULTI-DIMENSIONAL FILTER BANKS BY LIFTING TECHNIQUE

Tchobanou M.

Moscow Power Engineering Institute (Technical University)

Lifting technique is a polynomial approach for the design of multi-dimensional multirate systems. A technique developed by the author is being applied in the paper for the synthesis of such multirate systems, based on the design of multi-dimensional filters with fractional MD delay. l₂ norm is given for designed filters, that is important for denoising applications.



МЕТОД СВЕРХРАЗРЕШЕНИЯ ИЗОБРАЖЕНИЙ, СЖАТЫХ С ПОТЕРЯМИ С ПОМОЩЬЮ ВЕЙВЛЕТ-ПРЕОБРАЗОВАНИЯ

Баранов М.В., Дука Д.С., Чобану М.К.

Московский энергетический институт (технический университет)

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

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

В качестве алгоритма сжатия в этой работе рассматривается SPIHT [2]. В нем используется очевидная взаимосвязь между высокочастотными полосами ДВП. Используя принцип прогрессивной побитовой передачи, он формирует выходной битовый поток, в котором сначала передаются старшие значащие биты низкочастотной полосы ДВП, затем высокочастотные области высоких уровней разложения и т.д. Т.о. основная погрешность работы алгоритма определяется отбрасыванием младших бит коэффициентов ДВП в высокочастотных областях. На величину ошибки влияет коэффициент сжатия, базис и количество уровней разложения ДВП. Число уровней разложения влияет на качество прямо пропорционально его величине. Оптимальный базис ДВП зависит от структуры изображения. Чаще всего применяется разделимый базис, в виду простоты его применения, однако неразделимый базис позволяет получить лучшие результаты. [4, 7] Неразделимые фильтры выделяют двумерные пространственные частоты, тогда как разделимые - горизонтальную и вертикальную составляющие. На практике это позволяет реализовывать более адаптированный к особенностям человеческого восприятия алгоритм сжатия.

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

Такое различие в алгоритмах сжатия, основанных на разделимом и неразделимом ДВП, приводит к различиям в формировании погрешности в изображении.


(b)

(а)

(с)

Рис. 1 (a) Верхние и нижние оценки коэффициентов в-т преобразования,(b) ИХ фильтра восстановления

(с) результат

Для простоты рассмотрим способ анализа ошибок вносимых сжатием разделимым ДВП. Пусть в матрице K содержатся количества непереданных бит каждого отсчета W – разделимого ДВП исходного изображения Ki. Также пусть ΔW обозначает матрицу погрешностей для W: (W - ΔW) ≤ W ≤ (W + ΔW).

ΔW(y,x)=(2^N-1) & (2^Ki(y,x)-1-1)|Ki – исходное изображение, N – длина разрядной сетки.

При обратном ДВП одного уровня разложения, в случае применения разделимых фильтров, пространственный отсчет Ki(x,y) получается в результате последовательной свертки коэффициентов ДВП, находящихся на строках и затем на столбцах, с ИХ фильтра синтеза.

Для простоты изложения рассмотрим одномерный случай, см. рис. (1).

Введем обозначения:

x –вектор исходных значений

h – импульсная характеристика фильтра восстановления

y = x*h вектор результата, получающийся в результате свертки

x^sup, x^inf – верхняя и нижняя оценки отсчетов вектора x с учетом погрешностей

y^sup,y^inf – искомые оценки результирующего вектора

h⁺ = 0.5 (h+|h|) – вектор с отсчетами вектора h, с нулями вместо отрицательных отсчетов

h ^– = 0.5 (h-|h|) – вектор с отсчетами вектора h, с нулями вместо положительных отсчетов

Тогда легко видеть, что

y^sup = x^sup * h⁺ - x^inf * h^–

y^inf = x^inf * h⁺- x^sup * h ^–

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

Далее, полученные оценки диапазонов [y^inf, y^sup] для различных кадров исходной сцены накладываются друг на друга на сетке координат высокого разрешения. Вследствие совпадения координат отсчетов разных кадров на сетке высокого разрешения, диапазоны этих точек пересекаются. Это в конечном счете приводит к уменьшению ошибки в этих точках. Далее для получения значений точек изображения высокого разрешения необходимо проводить интерполяцию на нерегулярной сетке координат [3]. Для этого использовался алгоритм основанный на ДВП – Multiresolutional Basis Fitting Reconstruction.

Для простоты рассмотрим его одномерный случай. Пусть имеется функция f(t), для которой необходимо получить M равномерно распределенных значений на диапазоне t = 0, 1…M-1. Дано P неравномерно распределенных значений этой функции в точках t = t₀, t₁…t_p-1. В обычной постановке проблемы P < M. Примем равномерно распределенные отсчеты с шагом 1 за пространство V₀ с ортогональным базисом шкалирующей функции φ(t). Далее применяя вейвлет разложение f(t) можно представить как сумму декомпозиций f_J(t) и g_J(t) принадлежащих пространствам V_J и W_J с базисными шкалирующими функциями φ_J,k (t) = φ(2 ^J t – k ) и вейвлет-функциями ψ_J,k (t) = ψ(2 ^J t – k ). Тогда разложения f_J и g_J по этим базисам являются вектора коэффициентов a_J,k и b_J,k.

f(t) = ∑a_J,k φ_J,k (t) + ∑∑b_J,k ψ_J,k (t)

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

f(t_i) = ∑a_J,k φ_J,k (t) + ∑∑b_J,k ψ_J,k (t), i = 0, … P-1

Таким образом, мы получили систему из P линейных уравнений, в векторной форме ее можно записать как f = G_J a_J + ∑ H_J b_J. Эта задача является плохо определенной, поэтому ее решение достигается только в оптимальном смысле. В [3] авторы сначала приближенно находят решение f ≈ G_J a_J, затем ищут f - G_J a_J ≈ H_J b_J.

Описанный выше алгоритм легко распространяется на двумерный случай. Также, он не привязан к разделимому ДВП и применяется для сверхразрешения изображений, сжатых алгоритмом сжатия на основе неразделимого ДВП.

Литература

1. 
   Xin L., “Improved Wavelet Decoding via Set Theoretic Estimation”, IEEE Trans. on Circuits and Systems for Video Technology, vol. 15, NO. 1, JANUARY 2005
   
2. 
   Said A., Pearlman W.A., “A New Fast and Efficient Image Codec Based on Set Partitioning in Hierarchical Trees”, IEEE Trans. on Circuits and Systems for Video Technology, vol. 6., p. 243-250, June 1996.
   
3. 
   Nguyen N., Milanfar P., “A wavelet-based interpolation-restoration method for superresolution (wavelet superresolution)” , Circuits systems signal process Vol. 19, No. 4, 2000, PP.321-338
   
4. 
   Чобану М.К. Иерархический алгоритм кодирования для неразделимых решеток и банков фильтров // Вычислительные технологии. 2007. №4. с. 106-119.
   
5. 
   Чобану М.К. Аналитический синтез многомерных многоскоростных систем // Успехи современной радиоэлектроники. 2007. №4. СС. 61-80.
   
6. 
   Чобану М.К., Батлук А.В. Исследование применения банков фильтров для сжатия изображений // Цифровая обработка сигналов. 2005. №4. сс.29-40.
   
7. 
   Чобану М.К., Черников А.В. Современный метод сжатия изображений на базе вейвлет-преобразования и иерархического алгоритма кодирования. // Цифровая обработка сигналов. 2005. №3. сс.40-59.
   

SUPER-RESOLUTION RECONSTRUCTION OF IMAGE SEQUENCES COMPRESSED BY DWT

Baranov M., Duka D., Tchobanou M.

Moscow Power Engineering Institute (Technical University)

This paper describes the methodology for super-resolution of image sequences compressed with both of separable and non-separable dwt-based techniques. Methodology of analysis of errors produced by such algorithms also was described.



МЕТОД УЛУЧШЕНИЯ ИЗОБРАЖЕНИЯ СЖАТОГО АЛГОРИТМОМ СЖАТИЯ С ПОТЕРЯМИ НА ОСНОВЕ ВЕЙВЛЕТ-ПРЕОБРАЗОВАНИЯ

Дука Д.С., Баранов М.В., Чобану М.К.

Московский энергетический институт (технический университет)

Последнее время практика применения вейвлетов для решения задач сжатия и обработки изображений и видео, являющихся нестационарными по своей природе, получает большое развитие. В этой области применение вейвлет-преобразования позволило достичь одновременного снижения сложности и повышения эффективности кодеков. В настоящее время разработаны международные стандарты по сжатию неподвижных изображений и видео, включая JPEG-2000, H.264/AVC, WMP10 и MPEG-4. Задача повышения эффективности вейвлетных кодеков становится все более актуальной. В отличие от алгоритмов сжатия, основанных на дискретном косинусном преобразовании, основной проблемой которых являются границы блоков, алгоритмы сжатия на ДВП размывают изображение и создают эффект ауры на резких границах на изображении.

В


Рис. 1: Зависимости узлов в пространственно-ориентированном дереве (SPIHT)
предлагаемом докладе рассматривается вопрос применения метода сверх-разрешения при сжатии изображений методом SPIHT [3; 4] (Set Partitioning In Hierarchical Trees, рис. 1). В этом методе используется очевидная взаимосвязь между высокочастотными полосами ДВП. Используя принцип прогрессивной побитовой передачи, он формирует выходной битовый поток, в котором сначала передаются старшие значащие биты низкочастотной полосы ДВП, затем высокочастотные области высоких уровней разложения и т.д. Т.о. основная погрешность работы алгоритма определяется отбрасыванием младших бит коэффициентов ДВП в высокочастотных областях. На величину ошибки влияет коэффициент сжатия, базис и количество уровней разложения ДВП. Число уровней разложения влияет на качество прямо пропорционально его величине. Оптимальный базис ДВП зависит от структуры изображения. Для определения диапазона допустимых значений точек разжатого изображения предлагается использовать следующий метод. Пусть в матрице K содержатся количества непереданых бит каждого отсчета W – ДВП исходного изображения Ki. Также пусть ΔW обозначает матрицу погрешностей для W:

(W - ΔW) ≤ W ≤ (W + ΔW) ΔW(y,x)=(2^N-1) & (2^Ki(y,x)-1-1)

Ki – исходное изображение, N – длина разрядной сетки.

При обратном ДВП одного уровня разложения, в случае применения разделимых фильтров, пространственный отсчет Ki(x,y) получается в результате последовательной свертки коэффициентов ДВП, находящихся на строках и затем на столбцах, с ИХ фильтра синтеза.

В статье [1] приводится способ пост-обработки разжатого изображения. Алгоритм, описанный в этой статье, относиться к классу POCS. Вводимые множества:

• 
  Множество изображений, вейвлет-преобразование которых лежит в пределах погрешностей, вызванных алгоритмом сжатия с потерями (). (1)
  
• 
  Множество кусочно-постоянных функций, оператор проекции на которое является решением дифференциального уравнения и приводится в [2]. Легко видеть, что кусочно-постоянные функции принадлежат пространству Бесова, т.к. в них допускается конечное число разрывов первого рода и сами по себе формируют выпуклое множество [1].
  
• 
  Множество сигналов ограниченных по амплитуде. (2)
  

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

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

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

Литература

[1] Li X. «Improved Wavelet Decoding via Set Theoretic Estimation», IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS FOR VIDEO TECHNOLOGY, VOL. 15, NO. 1, JANUARY 2005

[2] You Y., Xu W., Tannenbaum A., Kaveh M., «Behavioral analysis of anysotropic diffusion in image processing», IEEE Trans. Image Process. Wall. 5, No. 11, pp.1539-1553, Nov. 1996

[3] Чобану М.К., Черников А.В. «Современный метод сжатия изображений на базе вейвлет-преобразования и иерархического алгоритма кодирования». // Цифровая обработка сигналов. 2005. №3. cc. 40-59.

[4] Said A., Pearlman W.A., «A New Fast and Efficient Image Codec Based on Set Partitioning in Hierarchical Trees», IEEE Trans. on Circuits and Systems for Video Technology, vol. 6., p. 243-250, June 1996.
METHOD FOR IMPROVING LOSSY COMPRESSED IMAGE USING WAVELET-TRANSFORM

Duka D., Baranov M., Tchobanou M.

Moscow Power Engineering Institute (Technical University)

This paper describes technique for improving image quality on decoder side via a set theoretic approach at low bit rates. The decoding rule is project-onto-convex-sets, which exploits a priori information about the quantization and geometric constraints of the image source. Proposed approach utilizes the space domain samples error estimation.



Варианты реализации графического расчёта в алгоритме моделирования распространения нефтяного загрязнения для задач экологического мониторинга

Втюрин С.А.

Институт космических исследований РАН

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

Постановка общей задачи

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

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

Задача данного доклада

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

• 
  требуемым для их реализации упрощениям и допускам, исходной модели;
  
• 
  точности итогового прогноза;
  
• 
  сложности алгоритма и характера её зависимости от масштабов модели.
  

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

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

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

Принятая модель, учитывает механизмы вязко-гравитационного растекания и дрейфа (под воздействием течений и ветра) растекающегося объема вязкой несжимаемой жидкости – нефти и предполагает использование следующих допущений:

• 
  скорость вязкостно-гравитационного растекания равномерна, т.е. не изменяется в зависимости от уменьшения “толщины пленки”;
  
• 
  каждому объему нефти определенной марки соответствует свой предельный радиус или предельная площадь растекания[1];
  
• 
  Используемый для расчетов фрагмент цифровой карты с данными о загрязнении должен обеспечивать пространство для моделирования в интересующий период.
  
• 
  Исходя из характера задачи, учет любых возможных допусков и неточностей производится по схеме "на худший случай".
  
• 
  Границами области определения модели являются границы водоема, при этом в соответствии с предыдущим пунктом, не учитывается возможные потери вещества при соприкосновении с берегом или иными препятствиями.
  

Исходная математическая модель.

Математической основой расчета, использованной как основа данной модели, является формализация изменения радиуса вязкостно-гравитационного растекания нефтепродуктов по водной поверхности выполненная по формуле Букмастера[2]: (1), где - радиус растекания во времени, м; g = 9,81 м/с²; = 1-; - плотность нефтепродуктов, г/см³; - плотность воды, г/см³; V - объем вылива, м³; - кинематический коэффициент вязкости воды, , м²/с; a - температурный коэффициент; t - время в секундах (с).

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

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

Алгоритм прогнозного моделирования.

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

Рассмотрено три принципиально отличающихся подхода к реализации процедуры моделирования растекания и дрейфа:

• 
  Подробное физическое моделирование. По вышеупомянутым причинам данный вариант признан нецелесообразным и в дальнейшем не рассматривается.
  
• 
  Моделирование на основе растрового представления данных.
  
• 
  Моделирование на основе векторного расчёта перемещения границ нефтяного пятна.
  

Растровая модель

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

Построение изменений производится на основе изображения пятна, при помощи применения двумерного фильтра с маской, соответствующей размытию по Гауссу, дополненного пороговой функцией что, позволяет автоматически учесть сложные конфигурации нефтезагрязнений, применяя фактически стандартные алгоритмы фиксированной сложности O(n²), зависящей от величины растра. Учет препятствий и границ водоёма производится путём определения достижения условного цвета и соответствующей коррекцией при расчете времени шага.

Векторная модель

Векторный вариант расчёта эффективен для простых односвязных областей и предполагает автоматическое или автоматизированное оконтуривание области загрязнения и представления её контура в виде многоугольника с заданной точностью. На текущем шаге итерации каждая из вершин многоугольника может быть смещена на вектор, определяемый как сумма векторов растекания, направленного в данной точке по перпендикуляру к касательной, и компоненты дрейфа, вычисляемой по формуле (2) и направленной, в случае наличия подробных данных, согласно карте течений и направлению ветра. Такая модель позволяет точнее учитывать условия и особенности заданной местности, а также более точно рассчитывать изменение формы пятна. Учёт препятствий при этом производится аналогично растровому случаю. Таким образом сложность основной части алгоритма линейно зависит от количества точек на периметре O(n). Однако при моделировании загрязнений сложной формы возможно возникновение самопересечений кривой контура. Контроль возникновения, правильная интерпретация и исключение этих пересечений требует подключения вспомогательных алгоритмов, сложность которых достигает O(n²) от количества точек, задающих многоугольник, а реализация вычисления площади в случае областей сложной формы- O(n²) от величины области исходного растра.

Сравнение вариантов реализации

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

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

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

Результаты сравнения различных вариантов процедуры моделирования по различным параметрам представлены в таблице.

Параметр	Растровая модель	Векторная модель
Требуемые упрощения модели	Область загрязнения меняет форму только под действием границ и препятствий	Область загрязнения является односвязной
Точность предсказания	Более грубая оценка области, не полностью учитывает изменение формы области.	Обеспечивает более точный прогноз при доступности подробных исходных данных об условиях для данного района.
Вычислительная сложность алгоритма	Фиксированная O(n2), зависит от величины растра	Зависит от формы области, точности аппроксимации периметра, этапа итерации. От O(n) до O(n2)
Особенности	Масштабирование модели на большую область существенно увеличивает время работы или уменьшает точность.	Возможно взаимодействие с ГИС при помощи векторного представления. Легче поддаётся масштабированию.
Область применения	Более грубая оперативная оценка и прогноз области загрязнения, возможно в условиях недостатка исходных данных.	Более подходит для точного моделирования в ходе проектирования и исследования потенциальных источников загрязнения.