Пикселная технология дискретизации акватории Мирового океана Текст научной статьи по специальности «Математика»

Вычислительные технологии

Том 14, № 5, 2009

Пикселная технология дискретизации акватории Мирового океана*

И. В. КИРЕЕВ, С.Ф. ПяТАЕВ Учреждение Российской академии наук Институт вычислительного моделирования СО РАН, Красноярск, Россия, e-mail: kiv@icm.krasn.ru, psf@icm.krasn.ru

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

Ключевые слова: триангуляция, сетка, компьютерная графика, батиметрия.

Введение

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

Методам построения сеток посвящено большое количество работ, обзор которых можно найти в [1-5], тем не менее дискретизация некоторых двумерных многосвязных областей до сих пор вызывает определенные трудности. Это относится и к задаче сеточного разбиения акватории Мирового океана П, которая до сих пор далека от своего полного решения. Достаточно сказать, что этой проблеме посвящены полностью два номера тома 58 журнала Ocean Dynamics за декабрь 2008 года (Springer Berlin / Heidelberg).

Сложность решаемой проблемы в первую очередь связана с задачей описания области, которая в общем случае рассматривалась еще Декартом и известна в литературе [6] как обратная задача аналитической геометрии. Очевидно, что описать область П, т.е. задать рельеф дна и контуры границ непосредственно по навигационным картам, сложно из-за обширности П. Для этой цели естественно использовать батиметрические базы данных серии ETOPO, которые можно свободно получить на интернет-сайте Национального центра геофизических данных США (http://www.ngdc.noaa.gov/mgg/global/ global.html).

* Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант № 08-01-00621-а "Вычислительные технологии для расчета течений несжимаемой жидкости").

© ИВТ СО РАН, 2009.

1. Постановка задачи

Цель данной работы — создание алгоритма дискретизации максимально возможной связной части П акватории П по заданной минимальной глубине к > 0 и заданным наименьшему и наибольшему шагам сетки со сгущением сетки около береговой линии. Элементами дискретизации должны быть квадраты или равнобедренные прямоугольные треугольники. Необходимая батиметрическая информация должна быть получена только из анализа данных базы ЕТОРО.

2. Дискретизация акватории Мирового океана

ETOPO представляют собой целочисленные базы данных превышений земной поверхности относительно уровня моря с определенным, в зависимости от версии, шагом по долготе Л £ [0°, 360°) и широте ^ £ [0°, 180°]; значение Л = 180° соответствует Гринвичскому меридиану, ^ = 0° — Северному, а ^ = 180° — Южному полюсам Земли. Если (Л, не являются координатами точки базы ETOPO, то соответствующее значение превышения можно аппроксимировать, например, с помощью линейной лагранжевой интерполяции на наименьшем квадрате, который содержит рассматриваемую точку и вершины которого являются элементами базы. Таким образом, в каждой точке M поверхности Земли определена величина d(M) как функция Л и задающая превышение (в метрах), отрицательное для точки, расположенной ниже уровня моря, и неотрицательное в противном случае.

Неравенство d(M) < — h задает несвязную область содержащую не только акваторию П (П0 = П), но и водные бассейны, с ней не связанные, такие как оз. Байкал, Каспийское море и т. п. Кроме того, в могут входить отдельные участки континентов, расположенные ниже уровня моря на h метров, например, часть территории Голландии. В итоге задача описания акватории П сводится к выделению П из Для ее решения воспользуемся методами вычислительной геометрии [7].

Представим всю поверхность Земли, развернутую в прямоугольник

P = [0°, 360°] х [0°, 180°],

в виде объединения 2N0 х N0 квадратов Q°0j0, i0 = 0,... , 2N0 — 1, j0 = 0,... , N0 — 1, со стороной Д0 = 180°/N0. При дискретизации эти квадраты будут соответствовать максимальным по размеру элементам сетки.

Введем понятие уровня вложенности квадратов: уровень вложенности максимальных по размеру квадратов Qj полагаем равным 0. Далее каждый из них разбиваем на четыре квадрата Qj1j1 со стороной Д1 = 2-1Д0 и приписываем им уровень вложенности 1; здесь I1 = (i0,i1), J1 = (jo, ji), ii, ji = 0,1. Задавая некоторое натуральное число L, называемое далее максимальным уровнем вложенности сетки, можно рекурсивно построить последовательность покрытий прямоугольника P квадратами Ql j , l = 1,... , L, со стороной Дг = 2-1Д0. Уровень вложенности этих квадратов полагаем равным l, и если Qlj г С QI-11Jг 1, то мультииндексы I, Ji определяются через "предыдущие" 1г-1, J^-1 следующим образом: I = (1г_ь i), J^ = J—, j), i, j = 0,1; I0 = (i0), J0 = (jo). Области значений мультииндексов 1г, J^ есть соответственно Dl = {0,1,... , 2N0 — 1} х {0,1}г, DJ = {0,1,... , N0 — 1} х {0,1}г, а целочисленные координаты ¿г, ji задают смещение левой нижней вершины квадрата QIj относительно аналогичной вершины квадрата

о1-1

QI-1J г-1.

В качестве примера на рис. 1, а изображен квадрат , в котором квадраты с номерами 1-4 имеют следующие уровни вложенности и мультииндексы:

I = 1, Ii = (i0,1), Ji = (j0,0);

l = 2, I2 = (io, 1, 0), J2 = (j0,1,1);

I = 3, I3 = (io, 1, 0,0), J3 = (jo, 1,0, 0);

I = 4, I4 = (io, 1, 0,0,1), J4 = (jo, 1, 0,1,1).

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

В конечном итоге возникает максимальное (по количеству квадратов) покрытие прямоугольника P квадратами QjpLjL, которые будем интерпретировать ниже как пикселы. Напомним [8], что пиксел (pixel) — наименьший неделимый элемент растрового изображения, характеризуемый прямоугольной формой с размерами, определяющими разрешение графического объекта. При дискретизации часть этих пиксел-квадратов будут соответствовать минимальным элементам сетки, расположенным вблизи береговой линии.

Точку M £ P назовем "мокрой", если глубина в ней не меньше h метров, т. е. выполняется неравенство d(M) < — h, и назовем "сухой" в противном случае. Пиксел назовем "мокрым", если все его вершины — "мокрые", и назовем "сухим", если все его вершины — "сухие". Все оставшиеся пикселы, через которые проходит граница акватории, будем называть "прибрежными". На рис. 2, а приведен схематичный пример акватории,

2

4

3

1

Рис. 1. Построение вложенных квадратов в квадрате до с дроблением в 5 раз (Ь = 5)

Рис. 2. Прямоугольник Pa до заливки (а), после 1-заливки (б), после 2-заливки (в).

а

а

в

расположенной в некотором прямоугольнике. Здесь □ — "сухой" пиксел, 0 — "мокрый", а пикселы иного вида — "прибрежные".

Для построения дискретного аналога акватории Мирового океана необходимо выделить из множества всех пикселов прямоугольника P связное объединение "мокрых" пиксел из nh. Аналогичные задачи существуют в компьютерной графике и связаны с закрашиванием или заполнением областей экрана монитора [8, 9]. Решение этих задач осуществляется так называемыми алгоритмами заливки.

Из всего многообразия алгоритмов заполнения двумерных областей остановимся на итеративной реализации простого алгоритма заливки гранично-определенной 4-связной области [9], в котором доступ к соседним пикселам осуществляется по четырем направлениям — горизонтально влево и вправо и вертикально вверх и вниз. Иллюстрация работы этого алгоритма применительно к акватории, изображенной на рис. 2, а, приведена на рис. 2, б, где Ш - "залитый мокрый" пиксел. Этот алгоритм заливки может порождать связные цепочки шириной в 1 пиксел, как на рис. 2, б. Поскольку линейный размер пиксела совпадает с минимальным шагом строящейся сетки, все пикселы таких цепочек являются граничными. Поэтому имеет смысл несколько изменить упомянутый выше алгоритм заливки.

Зададимся некоторым натуральным числом K. Квадрат, состоящий из K2 пикселов, будем называть K-курсором. Выделяя любой курсор из расчетной области, рассмотрим новый, передвигая выделенный курсор на 1 пиксел вдоль одной из координатных осей через границу ранее залитой части двумерной области аналогично тому, как происходит обычная заливка 4-связной области. Если новый курсор состоит только из "мокрых" пикселов, то пополняем залитую область новыми пикселами из этого K-курсора. Будем назвать эту процедуру K-заливкой. Результат работы 1-заливки области на рис. 2, а совпадает с обычной заливкой и изображен на рис. 2, б, а результат 2-заливки той же области изображен на рис. 2, в, где 0 — "2-залитый мокрый" пиксел.

После того, как проведена K-заливка акватории Мирового океана и тем самым выделено связное множество "мокрых" пикселов П^, С П, перейдем к построению дискретного аналога Dh акватории В простейшем варианте в качестве элементов сетки можно рассматривать пиксел-квадраты, содержащиеся в П^. Однако это приводит к существенному увеличению числа узлов сетки, поскольку в подобластях, которые достаточно удалены от границ расчетной области и в которых, как правило, градиенты искомых функций малы, целесообразнее использовать более крупные конечные элементы. Таким образом, необходимо построить покрытие П набором квадратов так, чтобы линейные размеры двух соседних квадратов были бы равны между собой или различались в 2 раза. Это условие позволяет сгущать в нужных подобластях сетку и гарантирует плавность перехода от разреженных участков к участкам сгущения.

Определим множества П^, l = 0,... , L — 1, как объединение квадратов уровня l, целиком содержащихся в П£:

П ={ U QJ : QJ С П£ VI £DI, VJi £Dj| , l = 0,1,... , L — 1. U j J

Здесь условие QIj г С П^ означает, что все содержащиеся в QIjг пиксел-квадраты содержатся в П£. Очевидно, что имеют место включения

П С П С ... С С П С П+1 С ... С П£ С ПЛ.

Тогда дискретный Г? аналог акватории П^, описывается набором мультииндексов ^ = 1£ю где Г? = {(1о, 1о) : ^оЛ С ; Я? = {(I, : С П^П^}, а I =

Построенное таким образом разбиение области может содержать соседние квадраты, уровни вложенности которых различаются больше чем на 1 (на рис. 1, а, например, квадраты 1 и 3, 2 и 4). Во избежание такого резкого перехода проводится итерационное сглаживание сетки с контролем максимальности уровней вложенности граничных квадратов. Схематично результат работы процедуры сглаживания сетки на рис. 1, а представлен на рис. 1, б. Обычно для выполнения этих требований достаточно до 10 итераций сглаживания. После завершения процесса сглаживания узлы сетки нумеруются. Одновременно устанавливается принадлежность узлов и ребер границе и выделяются номера узлов, лежащих на серединах сторон квадратов. Возможно также нарушение условия периодичности по А, как, например, показано на рис. 3, а, где квадрат I, имеющий соседями квадраты II и III, не имеет узла на правом ребре А = 360°. В этом случае на середине правой стороны квадрата I (рис. 3, б) наносится узел 2 и он отождествляется с узлом 1, имеющим координату А = 0°.

Географические координаты (А, узлов сетки удобно вычислять через координаты (А(1ьлг), <^(1г,лг)) левой нижней вершины квадрата уровня /, связанные с его цело-

численными координатами I = (¿0, ¿1,... ¿г), ^ = (^0, ... ^г) соотношениями

г г

А(1г,лг) = До X] ¿к/2к, = До£ ^/2к. к=0 к=0

При этом географические координаты правой верхней вершины квадрата ^^ будут равны (А(1ьлг) + Дг, <^(1г,лг) + Дг), где Дг = Д0/2г. Вычисление по найденным координатам (А, глубины в каждом узле завершает построение сетки с квадратными элементами.

360° 0°

-О

-О

360° 0°

О

С>

о-

II

III

о

I

-о

а

о-

о-

II

III

а б

Рис. 3. Введение дополнительных узлов на Гринвичском меридиане при нарушении условия периодичности

1

2

1

I

а б в г

Рис. 4. Шаблоны разбиения квадрата на треугольники: а — узлами являются только вершины квадрата; б — середина одной стороны квадрата является узлом сетки; в — середины двух сторон квадрата являются узлами сетки; г — середины трех сторон квадрата являются узлами сетки

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

Для построения такой триангуляции в квадратах более низкого уровня вложенности, соседствующих с квадратами более высокого уровня, дополнительно вводим узлы в центрах квадратов и триангулируем их так, что узлы на серединах сторон становятся вершинами прямых углов треугольников. Для новых точек вычисляются географические координаты и глубины. Эта триангуляция представляет собой сетку, большая часть элементов которой — равные равнобедренные прямоугольные треугольники. Вблизи границы дQh расчетной области происходит сгущение сетки подобными треугольными элементами с коэффициентами подобия (\/2)-fc, k = 1,..., 2L, по отношению к треугольнику с катетом Ao.

Программный комплекс, реализующий описанный алгоритм, выполнен на языке программирования Фортран-95 без использования операторов, специфичных для каких-либо операционных систем или сред разработки, и может применяться как на ПК с операционной системой DOS или MS Windows-95 и выше, так и на вычислительных комплексах под управлением ЮНИКС-системы. Комплекс применялся для дискретизации акватории Мирового океана П с использованием батиметрической базы ETOPO2 версии 2001 года (http://www.ngdc.noaa.gov/mgg/global/relief/ETOPO2/ETOPO2-2001/), процедуры 2-заливки при максимальном уровне вложенности L = 3, основном шаге сетки Ao = 4° и при нескольких значениях минимальной глубины h.

На рис. 5 изображена дискретизация акватории Мирового океана квадратными конечными элементами при h = 0 м. На рис. 6 и 7 показаны триангуляции области при минимальных глубинах h = 0 и 30 м соответственно.

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

Рис. 5. Дискретизация акватории Мирового океана квадратными конечными элементами: Н = 0 м, До = 4°, Ь = 3; число элементов 26063, число узлов 14716

Рис. 6. Триангуляция области 0^: Н = 0 м, Д0 = 4°, Ь = 3; число элементов 65777, число узлов 36834

Рис. 7. Триангуляция области 0^: Н = 30 м, Д0 = 4°, Ь = 3; число элементов 70365, число

узлов 39046

этих систем может содержать в расчетной области особые точки, обусловленные спецификой географической системы координат. Для исключения сингулярности в таких системах достаточно перейти в новую сферическую систему координат, полюсные особенности которой располагаются на суше. Например, Северный полюс N (рис. 8) можно перенести на территорию Гренландии так, что Южный полюс по-прежнему будет расположен на территории Антарктиды. Для этого достаточно поместить "новый Северный полюс" в точку N' со старыми географическими координатами Ар = 44°51'30'' з.д., рр = 78°8'30'' с.ш. Если обозначить через А' € [0°, 360°] и р' € [-90°, 90°] координаты точки в новой сферической системе координат, старые сферические координаты которой равны А и р, то в качестве формул перехода можно принять соотношения, приведенные на рис. 8.

вш р' = сов(Л + Лр)сов р сов рр + вш р вш рр;

сов р' = ^Д — вш2 р' ;

вш Л' = вш(Л + Лр)сов р/ сов р';

сов Л' = [сов(Л + Лр) сов р вш рр — сов рр вш р]/ сов р'; Л' (—Лр, рр) = —Лр; р'(—Лр, рр) = 90°

Рис. 8. Перемещение полюсной особенности на территорию Гренландии

Рис. 9. Триангуляция акватории Мирового океана: Н = 0 м, До = 4°, Ь = 3; число элементов 66714, число узлов 37284

Рис. 10. Триангуляция акватории Мирового океана: Н = 200 м, До = 4°, Ь = 3; число элементов 69480, число узлов 38650

На рис. 9 и 10 представлены триангуляции акватории Мирового океана в новой системе координат Л', р' с глубинами от 0 и 200 м соответственно.

Детализация сетки акватории Мирового океана определяется минимальным шагом Al и зависит от конкретной задачи. В случае больших значений Al часть островов и заливов, размеры которых соизмеримы с Al, на сетке могут не отразиться. Так, на приведенных рисунках острова Курильской гряды не проявились (за исключением двух) из-за достаточно большого Al = 0.5°, обусловленного визуализацией сетки. При необходимости степень детализации можно повысить, уменьшая Д0 и(или) увеличивая значения параметра уровня вложенности L.

Заключение

В задачах численного моделирования динамики Мирового океана границу dтруд-но установить точно из-за суточных и сезонных колебаний уровня Мирового океана и шага базы ETOPO. Вместе с тем, как показано в работе [12], аппроксимация границы береговой линии по возможности без резких угловых изломов, важная при моделировании прибрежных процессов, незначительно влияет при расчетах крупномасштабных процессов в Мировом океане. Поэтому предложенная дискретизация области оказывается удовлетворительной для большого класса задач, а применение описанной выше пикселной технологии значительно ускоряет процесс формирования сетки. Так, с приведенными в работе значениями параметров Д0, L и h при использовании базы ETOPO2 время дискретизации акватории Мирового океана на ПК класса Pentium IV с тактовой частотой примерно 1.5 GHz составляло не более 10 с, включая считывание батиметрических данных и запись вычисленных параметров сетки. В случае необходимости более точного описания границы расчетной области можно провести граничную коррекцию [3] построенной сетки. В заключение отметим, что триангуляция акватории Мирового океана по предложенному алгоритму была использована в работах [13-15] для моделирования гравитационных и поверхностных волн в Мировом океане.

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

[1] Frey P.J., Geotge P.L. Mesh Generation Application to Finite Elements. HERMES Science, 2000. 817 p.

[2] Глллнин М.П., Щеглов И.А. Разработка и реализация алгоритмов трехмерной триангуляции сложных пространственных областей: итерационные методы. М., 2006. (Препр. ИПМ им. М.В. Келдыша. РАН. № 9.)

[3] Галанин М.П., Щеглов И.А. Разработка и реализация алгоритмов трехмерной триангуляции сложных пространственных областей: прямые методы. М., 2006. (Препр. ИПМ им. М.В. Келдыша РАН. № 10.)

[4] НЕмировский Ю.В., Пятаев С.Ф. Автоматизированная триангуляция многосвязных областей со сгущением и разрежением узлов // Вычисл. технологии. 2000. Т. 5, № 2. С. 82-91.

[5] Методы римановой геометрии в задачах построения разностных сеток / Ю.И. Шокин, В.Д. ЛисЕйкин, А.С. Лебедев, Н.Т. Данаев, И.А. Китаева. Новосибирск: Наука, 2005. 255 с.

[6] Рвачев В.Л. Теория R-функций и некоторые ее приложения. Киев: Наукова думка, 1982. 552 с.

[7] Препарата Ф., Шеймос М. Вычислительная геометрия: Пер. с англ. М.: Мир, 1989. 478 с.

[8] Роджерс Д. Алгоритмические основы машинной графики: Пер. с англ. М.: Мир, 1989.

[9] Вельтмандер П.В. Машинная графика: Учеб. пособие в 3 кн. Кн. 2. Основные алгоритмы. Новосибирск: Изд-во НГУ, 1997.

[10] Owen S.J. A Survey of Unstructured Mesh Generation Technology // Proc. of 7th International Meshing Roundtable. Dearborn. Michigan. 1998. Р. 239-269.

[11] Гилл А. Динамика атмосферы и океана. Пер. с англ. М.: Мир, 1986. Т. 2.

[12] Марчук Г.И., Залесный В.Б., Кузин В.И. О методах конечных разностей и конечных элементов в задаче глобальной ветровой циркуляции океана // Изв. АН СССР. Физика атмосферы и океана. 1975. Т. 11, № 12. С. 1294-1300.

[13] Kamenshchikoy L.P., Karepoya E.D., Shaiduroy V.V. Simulation of surface waves in basins by the finite element method // Russian J. Numer. Anal. Math. Modelling. 2006. Vol. 21, № 4. P. 305-320.

[14] Kamenshchikoy L.P., Karepoya E.D., Shaiduroy V.V. Numerical solution of the boundary problem for shallow water equations for modelling surface waves in World ocean by finite elements methods // Proc. of Fourth International Conf. Finite Difference Methods: Theory and Applications (FDM:T&A'06). Bulgaria: Rousse. 2007. P. 227-233.

[15] Agoshkoy V.I., Kamenshchikoy L.P., Karepoya E.D., Shaiduroy V.V. Numerical solution of some direct and inverse mathematical problems for tidal flows // Comp. Science and High Perf. Computing III, E. Krause et al. (Eds). NNFM 101. Berlin Heidelberg: SpringerVerlag, 2008. P. 31-43.

Поступила в редакцию 19 марта 2009 г., в переработанном виде —13 мая 2009 г.