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

Вестник МГТУ, том 13, №4/2, 2010 г.

стр.1087-1091

УДК 519.65

Построение зональных карт уровенной поверхности океана по спутниковым данным на основе в-сплайн-интерполяции

А.Д. Гомонов

Лаборатория биоэкономического мониторинга и краткосрочного прогнозирования ПИНРО; Политехнический факультет МГТУ, кафедра высшей математики и программного обеспечения ЭВМ

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

Abstract. At present, fishing fleets put a considerable amount of efforts in exploring prospective fishing grounds. However, if any information on directions of currents in the given areas of the ocean is available, one could identify the areas perspective for fisheries with less expenditures. Regional maps of oceanic currents can be particularly drawn based on appropriate maps of sea-surface level. The paper reviews the method of construction of sea level maps, where satellite altimetry data are applied, as opposed to existing methods. This allows to perform necessary calculations using data for a shorter-term period and in a broader range of latitudes. Mathematically, the method is based on interpolation of data using the B-spline technique.

Ключевые слова: интерполирование, прикладная математика, сглаживающие сплайны, альтиметрия, базисные сплайны Keywords: approximation, interpolation, applied mathematics, smoothing splines, altimetry

1. Введение

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

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

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

Необходимо учитывать, что в альтиметрических данных присутствуют ошибки: в координатах измерений (чем ближе к полюсам, тем больше ошибка), неточность самих измерений, учитываемая специализированными поправками к альтиметрическому сигналу. Сигнал от альтиметра проходит через атмосферу (атмосферные поправки) и отражается от мгновенного уровня морской поверхности (поправки на состояние и уровень подстилающей поверхности). Самая большая поправка связана с приливами. Существует несколько моделей приливных поправок: CSR 3.0, FES 95.2, GOT 99. Эти модели убирают приливную составляющую из альтиметрических данных, однако они не всегда хорошо работают в прибрежных зонах и в динамически активных районах Мирового океана.

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

1087

Гомонов А.Д. Построение зональных карт океанских течений...

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

Следует отметить, что подобные модели уже реализованы с использованием различных математических алгоритмов (Романов, 2004). Однако все они работают в диапазоне широт ±72° и, кроме того, для их работы необходимо накапливать спутниковые данные, как минимум, в течение 10 суток. Эти "особенности" в значительной мере снижают достоверность восстановленных полей на текущий момент времени.

В статье предложен численный метод построения полей аномалий уровенной поверхности океана (с целью последующего расчета течений) на основе применения 5-сплайнов. Метод позволяет учитывать характер локальной изменчивости альтиметрии (используя значения производных в заданных точках), что позволяет уменьшить период накопления спутниковых данных и определять в режиме реального времени быстро изменяющиеся структуры в океане. Кроме того, появляется возможность построения полей альтиметрии в диапазоне широт ±81.5°.

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

Данные спутниковых альтиметров представляют собой набор измерений вдоль трека спутника с дискретностью примерно 5-7 км. Учитывая то, что данные поступают с различных спутников, каждый из которых имеет свою траекторию движения, имеется конечное число измерений в интересующей акватории, в пространственных точках, в общем случае не совпадающих с узлами регулярной прямоугольной сети. Поставим задачу: исходя из конечного числа точек измерений уровенной поверхности океана, расположенных в узлах нерегулярной сети, определить значение уровенной поверхности океана в узлах регулярной прямоугольной сети.

Математическую постановку задачи будем осуществлять при следующих предположениях:

1. Пусть в прямоугольной области Q = [a,b] х [c,d] = {(x,y)eR2 : a < x < b; c < y < d}, введена регулярная прямоугольная сеть A=AxxAy, где:

4c={x*/ a = xo <...<xN = b, xk = Xk-i+hx, k = 1,2,.. ,,N};

Ay={yk / c = yo <-..<уы = d, yk = yk-1+hy, k = 1,2,.. ,,M};

hx и hy- шаги сетки по x и по y соответственно. В результате область Q разбивается на прямоугольники Ар = {(xy)eR2 : xe [xx+1],ye [yy+1] i=1,2,...,N;p=1,2,...,M}.

2. Пусть в области Q представлено L точек Pk = (xk,yk), kel, в которых заданы числа fk -измеренные значения исследуемой функции f(x,y). Здесь I = {k / k=1,...,L} - конечное множество индексов. Точки Pk, kel расположены в области Q произвольным образом и в общем случае не совпадают с узлами регулярной прямоугольной сетки A=AxxAy. В некотором подмножестве точек Pk могут быть заданы значения fs, seI1cI, которые являются значениями производных исследуемой функции по направлению ls: f s= df9/s.

Необходимо определить приближенные значения функции f(x,y) в узлах регулярной сетки А.

Решение поставленной задачи осуществляется посредством аппроксимации неизвестной функции f(xy) функцией g(x,y), удовлетворяющей следующим условиям:

- в точках Pk=(xk,yk), kel значения g(x,y) должны минимально отличаться от fk - измеренных значений f(x,y) в указанных точках;

- в точках Ps, seI1 производные функции g(x,y) по направлению ls должны иметь минимальные отличия от f s - заданных значений производных по направлению функции f(xy) в указанных точках;

- функция g(x,y) должна обладать достаточными свойствами гладкости.

3. Метод восстановления полей аномалий уровенной поверхности океана в терминах В-сплайнов

Теория приближения функций начала свое развитие с разработки методов многочленной аппроксимации в трудах П.Л. Чебышева, К. Вейерштрасса, С.Н. Бернштейна, М.В. Келдыша и др. (Хеминг, 1968; Вершинин, 1988). Однако практические возможности применения классической многочленной интерполяции ограничены. Это связано, в частности, с тем, что, во-первых, с увеличением числа узлов интерполяции на регулярных сетках возникают вычислительные сложности при решении систем линейных уравнений для нахождения коэффициентов многочленов. Кроме того, при вычислениях возникает быстрое накопление погрешностей округлений (Вершинин, 1988). Во-вторых, интерполяционные многочлены высокой степени имеют осцилляционный характер вне узлов интерполяции и не отражают реально происходящие физические процессы. Поэтому методы многочленной интерполяции нецелесообразно использовать для моделирования уровенной поверхности Мирового океана и ее аномалий.

При построении искомой аппроксимации g(x,y) необходимо учитывать, как было отмечено ранее, что данные измерений имеют погрешности. Очевидно, что эти погрешности будут присутствовать и в смоделированной уровенной поверхности Мирового океана. Однако с точки зрения конечной цели -

1088

Вестник МГТУ, том 13, №4/2, 2010 г.

стр.1087-1091

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

В качестве метода построения аппроксимации g(x,y) будем использовать сплайновые методы сглаживания. Для решения задачи построения карт уровенной поверхности океана достаточно использовать сглаживающие сплайны третьей степени, так как при их применении обеспечивается непрерывность первой и второй производных моделирующей функции g(x,y), что вполне достаточно с практической точки зрения. Исходя из особенностей решаемой задачи, удобно использовать так называемые 5-сплайны. Доказано, что 5-сплайны определенной степени образуют базис в соответствующих функциональных пространствах (Завьялов и др., 1980), а значит, существует возможность представления других функций этих пространств через 5-сплайны. Следовательно, всякий кубический сплайн дефекта 1, заданный на сетке А, может быть представлен через кубические 5-сплайны:

5(x)

' 0, x< x,;

(x-x,)3/6h3, xe [ x,, x*i];

J (h3+3h2(x-x,+i) + 3h(x-x,+1)2 - 3(x-x,+1)3)/6h3, xe [ xm, x^];

(h3+3h2(x,+3-x) + 3h(x,+3-x)2 - 3(x,+3-x)3)/6h3, xe [ x,+2, xM]; (x,+4-x)3/6h3, xe [ x,+3, x,+4];

^ 0, x> x+4,

где x,+i—x, = h, V i, в следующем виде:

N+3 M+3

g(x,y)=E E ap5,(x)5j-(y). (1)

i=-3 J=-3

Здесь aj - числовые коэффициенты, а всевозможные произведения 5i(x)5J-(y) образуют базис в пространстве кубических сплайнов двух переменных дефекта 1 в прямоугольной области £2 .

Традиционно (Завьялов, 1980; Вершинин, 1988) задача сглаживания для одномерного случая определяется, как поиск функции g, минимизирующей функционал Ф(д)= YФ\(&) + Ф2(§) где:

b

Ф\(g) = l[g(x)]2dx, Ф2(g)='LOk[g(Xk) - ft]2.

a ke I

Ф1g) обеспечивает свойство максимальной гладкости (минимум потенциальной энергии гибкой линейки), а Ф2(g) характеризует степень приближения функции g к точкам измерений (ok - вес k-го измерения). Параметр сглаживания у> 0 является количественным показателем "компромисса" между точной интерполяцией через точки и сглаживанием.

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

N+3 M+3

Ф2^) = E Ok [fk - E Eaj5,(x)5j-(y)f. (2)

ke I , =-3 J=-3

Еще одно слагаемое, входящее в минимизирующий функционал, по аналогии с одномерным случаем, должно характеризовать гладкость аппроксимирующей функции. Для двумерной задачи (Смоляк, 1971; Harder, Desmarais, 1971), для этой цели было предложено использовать модель упругой пластины бесконечной протяженности, которая деформируется лишь изгибом, т.е. ее отклонения от исходного состояния задаются в конечном числе независимых точек, где приложены соответствующие точечные нагрузки (Ландау, Лифшиц, 1965).

Гибкая пластина является очевидным аналогом одномерного сплайна - гибкой линейки. Соответствующий функционал гладкости (определяемый при этом полной свободой энергии изогнутой тонкой упругой пластинки) имеет вид (Смоляк, 1971; Harder, 1971):

j [(d 2g/dx2)2+2(d 2g/dxdyf+(d 2g/dy2)2]dxdy ^ min. (3)

Тогда, используя (1) и (3), в нашем случае можно записать:

N+3 M+3 N+3 M+3

ФА^ШЕ E{aj(d25,(xk)/dx2)5j-(yk))}2+2E E{a1-(d51(xk)/dx)(d5j(yk)/dy))}2+

ke I ,=-3 J =-3 ,=-3 J=-3

N+3 M+3

+E E{aj 5,(xk)(d25j-(yk)/dy2)}2],

i=-3 J=-3

где yk - параметр сглаживания k-го измерения.

(4)

1089

Гомонов А.Д. Построение зональных карт океанских течений...

Будем искать кубическую сплайн-функцию вида (1), которая будет минимизировать функционал

Ф(я)= Ф\(я) + Фг(Я),

(5)

где Ф1(^) и ФгС?) определяются формулами (4) и (2) соответственно. Для минимизации такой функции можно применить метод сопряженных градиентов, который обладает рядом преимуществ перед другими методами (штрафов, покоординатного спуска и т.д.). Так, в отличие от указанных методов, метод сопряженных градиентов обеспечивает гарантированное нахождение точки минимума квадратичной функции за конечное число шагов, которое не превышает числа ее переменных. Другие методы сходятся лишь в пределе (Васильев, 2002).

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

z = kx exp(-x2 - y2). (6)

Функцию такого вида удобно взять, так как представляется, что она достаточно хорошо позволяет описывать некую модель подъема и опускания вод уровенной поверхности океана. Для Северо-Восточной части Атлантического океана коэффициент k примем равным 20, что соответствует наблюдаемому на практике уровню поднятия и опускания вод в этом регионе. Для моделирования ошибок измерений на функцию (6) наложим шумовую составляющую, формируемую случайным образом (рис. 1), при этом максимальный уровень "шума" установим на уровне 10 % от максимальных значений генерируемых синтетических данных.

При восстановлении сгенерированных по формуле (6) синтетических данных посредством минимизации функционала (5) для функции (1) было принято %=const=0,996 (параметр сглаживания k-го измерения в (4)) и ck=const=1 (ok - вес k-го измерения в (2)), для keI. В результате был получен график сплайн-поверхности (рис. 2).

Оценку качества восстановленных синтетических данных по формуле (6) с "шумовой" составляющей можно оценить по введенным нормам:

Rq = max | [(gfey)- z(x„y_,))/z(xI-,y_,)] | - равномерная взвешенная норма,

kx ky

Lq = [(kx+1)(ky+1)]-1(E 'L[(g(xi,y]) - z(xiy))/z(xiy)]2)1/2 - среднеквадратическая взвешенная норма,

i=1 j=1

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

Результаты расчетов этих норм для функции (6) с наложенной на нее "шумовой" составляющей для различных шагов сети по x и по y приведены в следующей таблице.

Шаг сети восстановленных данных по x и по y

Норма 0,5 0,25 0,1

Rq 3,51 2,21 1,93

Lq 0,91 0,76 0,54

Как видно из данных, приведенных в таблице, двумерные 5-сплайны сходятся к аппроксимируемым функциям даже с некоторой наложенной на нее "шумовой" составляющей на последовательности прямоугольных сеток A=AxxAy, о чем свидетельствует уменьшение норм Rq и Lq при возрастании количества узлов сети восстанавливаемых данных.

Рис. 1. График функции (6) с наложенной шумовой составляющей

Рис. 2. Результат применения сплайн-фильтрации к функции (6) с наложенной шумовой составляющей

1090

Вестник МГТУ, том 13, №4/2, 2010 г.

стр.1087-1091

Как было отмечено ранее, зачастую кроме данных fk о значениях функции f, известны и значения производных по направлению в некоторых точках области Д то есть, даны df/dls, seIb Другими словами, в некоторых точках заранее известны пространственные тенденции изменения аномалий уровенной поверхности океана. Для их учета необходимо классический функционал (5) дополнить третьим слагаемым Ф3(д), которое будет учитывать данные о производных по направлению.

В общем случае производную по направлению для функции g(xb x2,..., xm) от m переменных xb x2,...,xm можно записать как dg/dl=(dg/dx1)cos(a1)+(dg/dx2)cos(a2)+(dg/dxm)cos(am). В случае функции двух переменных производную по направлению можно записать в виде:

dg/dl = (dg/dx)cos(a)+(dg/dy)sin(a). (7)

С учетом (1) и (7) Ф3(д), характеризующая интегральную степень приближения производных g(x,y) по направлению ls к заданным значениям этих производных в точках Ps, seI1, может быть представлена в виде:

N+3 M+3

Фъ^)= Хв [df/d ls - X 'Zaij((dBi(xs)/dx)B/(ys)cos(as) + B(xs)(dBj (ys)/dy)sin(as))f, (8)

se I1 i=-3 j=-3

где в вес производной по направлению в точке Ps, se I1.

В результате ищем аппроксимирующую функцию g(x,y) вида (1), минимизирующую функционал Ф^)=Ф^)+Ф2^)+Ф3^), слагаемые которого определяются формулами (2), (4) и (8).

N+3 M+3 N+3 M+3

Wg) = X yk[Z Е{ар(d2 Bl(xk)/dx2)Bj(yk)}2 + 2Е Е {a1](dB(xk)/dx)(dBr(yk)/dy)}2 +

keI i=-3 j=-3 i=-3 j=-3

N+3 M+3 L N+3 M+3

+Е Е {ajBl(xk)(d2Br(yk)/dx2)}2] + Xokfk-X XajB,(xk)Bjr(yk)]2 +

i=-3 j=-3 k=1 i=-3 j=-3

N+3 M+3

+Xh [dfdls - X Xalj((dBl(xk)/dx)Bj-(yk)cos(ak)+ B(xk)(dBj (yk)/dy)sin(ak))]2.

seIi i=-3 j=-3

4. Выводы

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

Данный подход был реализован программно. Кроме полей альтиметрии, программный комплекс позволяет осуществлять восстановление полей температуры, содержания хлорофилла и солености, что значительно повысило надежность прогноза перспективных участков промысла. Апробация изложенной методики начата летом 2010 года на промысле в районе Северо-Западной Африки. Полученные предварительные результаты позволяют утверждать, что применение данной методики позволило повысить производительность работы рыболовных судов, участвующих в эксперименте, до 15-20 % по сравнению с другими рыболовными судами.

Литература

Harder R.L., Desmarais R.N. Interpolation using surface splines. J. of Aircraft, v.9, N 2, p.189-191, 1972. Васильев Ф.П. Методы оптимизации. М., Факториал Пресс, 824 с., 2002.

Вершинин В.В., Завьялов Ю.С., Павлов Н.П. Экстремальные свойства сплайнов и задача сглаживания.

Новосибирск, Наука, 101 с., 1988.

Завьялов Ю.С., Квасов Б.И., Мирошниченко В.Л. Методы сплайн-функций. М., Наука, 352 с., 1980. Ландау Л.Д., Лифшиц Е.М. Теория упругости. М., Наука, 204 с., 1965.

Романов А.А. Восстановление мезомасштабной изменчивости аномалий высоты поверхности океана по данным спутниковой альтиметрии. Дис. ... канд. физ.-мат. наук, М., 125 с., 2004.

Смоляк С.А. Сплайны и их применение. Экономика и математические методы, т.7, вып.3, с.419-431, 1971. Хеминг Р.В. Численные методы. М., Наука, 398 с., 1968.

Шатохин Б.М. Исследование экологической роли уровенной поверхности океана и бифуркационных механизмов формирования зон повышенной рыбопромысловой продуктивности. Тезисы докладов конф. "Математическое моделирование и информационные технологии в исследованиях

биоресурсов Мирового океана", Владивосток, c. 133-136, 2004.

1091