Методика восстановления карт аномалий морской поверхности методом Бш-сплайнов.
Романов А. А. (romulas@ncmc.ru)
Всероссийский научно-исследовательский и проектно-конструкторский институт экономики, информации и автоматизированных систем управления рыбного хозяйства
(ФГУП "ВНИЭРХ")
Введение
В настоящий момент одной из основных задач спутниковой альтиметрии становится восстановление по данным альтиметрических спутников мезомасштабной изменчивости высоты морской поверхности для определенных регионов океана. Действительно, с появлением современных систем (Jason-1 и Envisat), прообразом которых стал спутник Topex/Poseidon подобная задача, наконец, может быть решена. Однако на пути ее решения возникает целый ряд проблем методического характера.
Пространственный и временной масштаб мезомасштабной изменчивости находятся в пределах 50-100 км и 10-100 дней на средних широтах [1 - 3]. Поскольку, при разработке альтиметрических миссий необходимо достигать компромисса между распределением подспутниковых треков на поверхности и длинной изомаршрутного цикла, к сожалению, на данный момент, невозможно восстанавливать и высоко, и низкоэнергетичные мезомасштабные структуры, независимо от используемого метода, используя данные только одной альтиметрической миссии. Достаточно трудно достичь необходимого пространственного разрешения, в случае проведения пространственно-временной интерполяции исходных данных [4, 5]. Например, для того чтобы разрешить структуру на поверхности океана размером 50 км на средних широтах, нужно использовать неприемлемо длинный изомаршутный цикл альтиметрической миссии длительностью около 2-х месяцев, что, несомненно, сильно усложнит восстановление мезомасштабной изменчивости структур на поверхности океана по времени. Для шаблона распределения наземных треков спутника Topex/Poseidon характер невязки восстановления является сильно неоднородным, из-за большого расстояния между треками по долготе, но из-за короткого 10-дневного изомаршрутного цикла спутника Topex/Poseidon, невязки восстановления карт с использованием его данных относительно однородны по времени. Характер распределения ошибок для шаблонов данных спутников ERS и Geosat более сложный: меньше в некоторые моменты времени и для некоторых пространственных точек благодаря более близко расположенным трекам спутников, однако неоднородны в пространстве из-за достаточно длинных изомаршрутных циклов (35 дней ERS, 17 дней Geosat) [6].
Возможным подходом для решения указанных проблем становится комплексирование данных от различных альтиметрических миссий. При этом возникает очередная сложность. Дело в том, что данные, полученные от различных альтиметров, изначально не являются пригодными для совместного использования. Причин несовместимости несколько: использование различных поправок к альтиметрическому сигналу при расчете аномалий, различные инструментальные погрешности альтиметров, а главное использование разных орбитальных моделей и, следовательно, разных ошибок определения орбиты спутников. На практике, подобные погрешности приводят к чрезмерно большим разницам восстановленных аномалий высот поверхности океана в так называемых кроссоверах, точках пересечения различных подспутниковых треков как одной, так и разных систем. Однако, с появлением спутника Topex/Poseidon с его небольшой орбитальной ошибкой, по сравнению с существовавшими на тот момент системами, появилась возможность «откалибровать» данные других спутниковых систем,
а именно Geosat и ERS-2 и уменьшить погрешность, вносимую неточной орбитой, в расчет аномалии морской поверхности в несколько раз. Методика подобных преобразований подробно описана в работах [7 - 10].
Следующая сложность - проблема пространственно временной интерполяции альтиметрических данных. Необходимо создать метод, который бы позволил провести подобную операцию с хорошей точностью и минимальными затратами по времени, поскольку использование данных нескольких альтиметрических систем приводит к необходимости единовременного анализа нескольких тысяч измерений. В настоящий момент для расчета карт аномалий морской поверхности, в мире, чаще всего используется метод «объективного анализа (objective analysis)». Этот метод обработки данных, основанный на теореме Гаусса-Маркова [11], был изначально предложен Гандиным [12] для анализа метеорологических полей. Впоследствии для целей океанологии, он был модифицирован Бретертоном [11], а в настоящее время модернизировался Ле Трайоном [1] для применения непосредственно к альтиметрической информации.
Главным достоинством, а также и недостатком метода «объективного анализа» является то, что для его работы требуется иметь априорное статистическое представление об обрабатываемых данных, причем в случае если по каким то причинам эта оценка является неточной, метод теряет свои свойства «оптимального интерполятора» [11]. И соответственно использование этого метода становится нецелесообразным.
Главной целью настоящей работы является предложение альтернативного метода и разработка методики обработки альтиметрической информации. Предлагаемый метод основан на интерполяции исходных данных сплайном «тонкой пластины», основным достоинством которого является возможность его использования на хаотически распределенных узлах исходной информации. С основной теорией Dm - сплайнов можно ознакомиться в работах [13 - 16]. В исследовании [17] приводятся результаты восстановления некоторых математических функций, а также карты аномалий морской поверхности, восстановленные данным методом. Однако, в работе [17] как для моделирования, так и для построения реальных карт высот поверхности океана использовались данные только одной спутниковой системы - Topex/Poseidon. В настоящей работе проведена реконструкция различных функций, повторяющих поведение мезомасштабных структур в океане, причем исходное распределение восстанавливаемых данных повторяет реальные шаблоны спутниковых альтиметрических систем Topex/Poseidon, Geosat-FO и ERS-2. Рассчитаны невязки и ошибки восстановления как для случая использования двух, так и трех систем. Приведены ошибки восстановления функций в случае использования одной спутниковой системы Topex/Poseidon. Осуществлен расчет карт реальных аномалий морской поверхности для случая двух спутников (Topex/Poseidon и ERS-2), а так же трех спутников (Jason-1, Envisat и Geosat - FO), а так же осуществлено сравнение карт построенных при помощи сплайнов и методом «оптимальной интерполяции» по комбинации данных двух или трех спутниковых систем.
Данная работа состоит из двух основных частей.
В первой части представлено описание методики расчета карт аномалий морской поверхности.
Вторая часть работы содержит результаты моделирования восстановления различных поверхностей, анализ ошибок и невязок, а также графики их зависимости от некоторых параметров использованной методики. Для наглядного сравнения результатов работы метода сплайн интерполяции представлены карты аномалий морской поверхности, построенные методом «оптимальной интерполяции» и методом Dm - сплайнов.
В заключение представлены основные выводы данной работы.
Методика расчета карт аномалий морской поверхности
Расчет карт аномалии морской поверхности по данным спутниковой альтиметрии состоит из нескольких этапов.
Для начала необходимо провести обработку альтиметрических данных вдоль трека спутника. Провести учет поправок, осуществить отбраковку некорректных данных, т. е провести
все необходимые процедуры для получения значений аномалий морской поверхности. Подробнее со всеми вышеперечисленными процедурами можно ознакомиться в работе [18].
Далее необходимо провести фильтрацию данных вдоль трека спутника. Так как вдольтрековые данные, очевидно, содержат не только информацию о мезомасштабной изменчивости морской поверхности, но и более крупномасштабных явлениях, необходимо принять меры для исключения ненужной информации из дальнейших расчетов. Для этого можно, например, воспользоваться фильтрами с определенными значениями частот отсечки или медианной фильтрацией с определенным окном. После того, как данные отфильтрованы, необходимо провести прореживание информации, поскольку слишком близкие узлы интерполяции (расстояние между соседними подспутниковыми измерениями составляет порядка семи километров) могут привезти к ухудшению числа обусловленности матрицы интерполяции, а, следовательно, и к потере устойчивости и большим ошибкам при ее обращении.
Однако, даже после проведенной процедуры прореживания данных, количество альтиметрических измерений, например, для географического региона размером 30х20 градусов может достигать от нескольких тысяч до десятков тысяч точек в зависимости от количества спутниковых систем, используемых в расчете. Данный факт приводит нас к необходимости решения системы из тысяч линейных уравнений для расчета тысяч неизвестных коэффициентов Бт - сплайна. Естественно, решение системы подобного размера представляется крайне затруднительной задачей и поэтому предлагается использовать не все исходные данные в регионе одновременно, а для каждого узла к регулярной сетки, которую мы намеревались построить с использованием всех альтиметрических данных в регионе, проводить дополнительную выборку, чтобы уменьшить количество входной информации для расчета конкретного узла сетки. Подход с использованием подвыборок к основному набору используется достаточно часто [1], так как возможности персональных компьютеров пока ниже требований по быстродействию процессоров и количеству оперативной памяти, предъявляемых при обработке геофизических массивов информации. Методика расчета карт состоит в следующем:
1. Производим выборку исходных альтиметрических данных для выбранного региона. В соответствии с работой [17] видно, что минимальная ошибка при интерполяции с учетом времени получалась при восстановлении функции на центральный день изомаршрутного цикла спутника Topex/Poseidon. Соответственно выбираем альтиметрические данные по району симметрично по времени относительно выбранной даты Т, восстановления функции (Т -5, Т +4 для данных спутника Topex/Poseidon, Т -8, Т +8 для Geosat-FO и Т-17, Т +17 для БЯ8-2). Используя данные полных изомаршрутных циклов спутника, мы получаем максимальное покрытие альтиметрическими данными исследуемого региона;
2. Выбираем узел к из регулярной сетки в котором ищется значение сплайна;
3. Выбираем исходные альтиметрические измерения вокруг этого узла по следующему алгоритму. Сначала берем все исходные альтиметрические точки, которые попадают в квадрат со стороной Я градусов и нашим исследуемым узлом к посередине. Причем, центры квадратов выбираются таким образом, чтобы обеспечить «нахлест» узлов сетки внутри квадрата в одну строку и столбец с соседними. Соответственно, ранее рассчитанные граничные узлы на ребрах квадратов, так же отбираются для участия в процедуре интерполяции с целью обеспечения непрерывности значений восстанавливаемой функции в граничных узлах квадратов сетки. Затем проводим дополнительный поиск в круге с радиусом Яь градусов и центром в исследуемом узле к. В данном случае ограничиваемся каждой четвертой точкой, удовлетворяющей географическому условию. Ограничиваем общее количество точек, выбранных по приведенной выше схеме, числом N. К сожалению, введение ограничений на количество точек - это вынужденная мера, поскольку бесконтрольное увеличение точек интерполяции ведет к существенному увеличению времени расчета. Конкретные значения для параметров методики Я?, Яь и N выявит дальнейшее моделирование.
Точки, попавшие в рассмотрение, должны удовлетворять условию по времени \Т- Т\ < 20, где Т - время исходной точки, а Т - время узла.
4. Проводим расчет коэффициентов сплайна для рассматриваемой нами области. Так как сплайн описывает не только значение в узле к, но и во всей области, над которой мы провели процедуру интерполирования, рассчитаем значения восстанавливаемой функции не только в узле к, но и для узлов находящихся в квадратной области со стороной Я и центром в к.
5. Производим повтор процедуры для всех узлов, в которых необходимо восстановить исследуемую функцию.
Численный эксперимент
Основная цель численного эксперимента состоит в оценке параметров Яь и N значения которых необходимо знать для восстановления карт аномалий морской поверхности. Как и в исследовании [17] нас интересует акватория Охотского и части Японского морей, однако, в данной работе мы несколько увеличим регион исследования и определим его соответствующими вершинами четырехугольной области: 135 - 165 в.д. и 40 - 63 с.ш. Как и в работе [17], для создания «квазиреального» распределения моделируемых альтиметрических данных спутниковых систем Topex/Poseidon, ЕБ£-2 и Geosat-FO (спутники ^оп-1 и Envisat имеют исходные шаблоны данных, аналогичные шаблонам Topex/Poseidon и ЕЯБ-2, соответственно) использовалась модель SGP4/SDP4, которая часто применяется для предсказания местоположения в специальных программах слежения за космическими объектами.
В рамках данного эксперимента восстанавливались два вида функций. Первый - функция
вида:
/ (х, у, г) = 15 мп(2 х + 7) ^(2 у + 7), (1)
где х - широта, у - долгота и г - время.
Данная функция выбрана с целью проверки возможностей предложенного алгоритма к восстановлению «ковра» вихрей (рис. 1а), размерная константа «15» обеспечивает высоту изменение амплитуды вихрей от -15 до 15 см, а константа «2» для широты и долготы введена с целью получения диаметра структур порядка 1.5 градусов (~150 км). Временная константа подобрана таким образом, чтобы за 30 дней каждый вихрь перемещался на расстояние порядка 600 км, подобные параметры соответствуют ожидаемым характеристикам поведения реальных мезомасштабных структур в океане. Однако, следует отметить, что полная картина, которая описывается модельной функцией, конечно же, не является типической для поверхности океана и выбрана с целью максимального усложнения задачи восстановления для метода сплайн-интерполяции.
Второй вид функции представляет собой суперпозицию функции вида (1), но с амплитудой порядка одного сантиметра с более крупномасштабной структурой, которая к тому же имеет более сложный характер поведения (рис. 2а). Исследуемая структура перемещается в направлении северо-востока и одновременно расплывается (увеличивается ее радиус и уменьшается амплитуда). Точный вид функции представляет собой следующее:
,, ) 20 , (х- 4 -150)2 (У - 4 - 50)2,
/ (хуг)=тто+г^—т+1---7+т~);
0 < г < 17 а = 12; 17 < г < 25 а = 8;
25 < г < 30 а = 7; /(х, у, г) < 1 /(х, у, г) = /(х, у, г) + sin(х + 7) ^(у + 7).
Размерная константа «—20—» отвечает за изменение амплитуды вихря в зависимости от
г / а +1
времени, причем скорость ее изменения зависит от момента времени, за это отвечает переменная а. В остальном, смысл использованных констант аналогичен соответствующим константам функции (1).
а)
а)
160" 165*
в)
в)
135' 140' 145' 150
Рис. 1. Вид функции (1). а) - исходная функция, б) - восстановленная по двум альтиметрическим источникам и в) -восстановленная по трем альтиметрическим источникам.
Рис. 2. Вид функции (2). а) - исходная функция, б) - восстановленная по двум альтиметрическим источникам и в) -восстановленная по трем альтиметрическим источникам.
Для каждой из приведенных выше модельных функций проводился полный цикл моделирования с использованием различных значений параметров, при этом использовались данные одной, двух или трех альтиметрических систем. Также проводились оценки качества реконструкций, путем расчета невязок и вычислением ошибок, представленных следующим образом:
функции, а ~ - значения восстановленной функции в узлах регулярной сетки.
Суммирование подразумевается по всем значениям аргументов (по всем узлам) восстанавливаемой функции.
Результаты численного эксперимента и сравнение результатов различных методов
Для определения оптимальных параметров метода для каждой исследуемой функции был проведен цикл исследований, который включал в себя расчет ошибок и невязки от Я от 1 до 2 с шагом 0.1, Яь от 4 до 10 с шагом 1 и N от 100 до 350 с шагом 10. Проведенный эксперимент показал, что минимальные ошибки в совокупности соответствуют следующим параметрам Я = 1.4, Яь = 6 и N=150 при этом время расчета составляет 1 минуту независимо от количества исходных данных.
На всех рисунках 1 - 4, контуры, проведенные штрихованной линией, соответствуют отрицательным высотам, а проведенные сплошной - положительным.
В случае использования данных двух спутников (Topex/Poseidon и ЕЯ^-2) для восстановления функции (1) ошибки составили 5С = 1.41, 5Ь =0.46 и невязка 2.37. Для исходных
данных с трех альтиметров ошибки составили 5С = 1.12, 5Ь =0.25 и невязка 1.23, соответственно.
Видно, что добавление третьего спутника уменьшает 5С на четверть, 5Ь на половину, а невязка
сокращается почти вдвое. Причем следует отметить, что в процентном соотношении от всего диапазона высот невязка составляет 8% и 3% соответственно. Конечно, бросается в глаза довольно большая ошибка по метрике 5С , ошибка по метрике 5Ь составляет почти 50%, однако, как видно
на рисунках 1б и 1в восстановление, особенно если учесть сложность моделируемой функции, прошло неплохо.
При восстановлении функции (2) для шаблона альтиметрических данных двух спутников ошибки составили 5С = 0.28, 8Ь =0.18 и невязка 0.14, а для трех - 8С = 0.22, 8Ь =0.16 и невязка 0.11, соответственно. Ошибки и невязка существенно меньше, в первую очередь, в силу того, что сам по себе диапазон изменения высоты функции (1) составляет 30 см, а для функции (2) всего 8 см, так как она призвана имитировать относительно спокойную морскую поверхность. Улучшение за счет данных третьего альтиметра в случае низкоэнергетичного региона составляет - 30% для 5С (чуть больше, чем для функции (1)), 13% для дь и 25% для невязки. Визуально разницу в
восстановлении функции (2) с использованием двух и трех альтиметрических источников отследить сложнее, чем для функции (1) - рисунки 2б и 2в. В процентном соотношении невязка составляет 2% в случае двух спутников и 1.5% при использовании данных 3-х альтиметров для функции (2) от всего динамического диапазона высот. Время восстановления составило порядка 50 секунд независимо о количества источников информации.
Для сравнения приводим данные по моделированию восстановления функции (1) с использованием данных только системы Topex/Poseidon. 5С = 6.52, 8Ь =1.38 и невязка 7.21. Полученные ошибки восстановления чрезвычайно велики, а это говорит о том, что данных только одной спутниковой системы недостаточно для восстановления структуры поверхности океана с высокоэнергетичными мезомасштабными вихрями. Моделирование восстановления
, где 2 - значения в узлах исходной модельной
низкоэнергетичных структур с использованием данных только альтиметра с Topex/Poseidon проводилось в работе [17].
На рисунках 3 и 4 представлены результаты восстановления поверхностей комбинации реальных альтиметрических данных со спутников Jason-1 и ERS-2, а так же Jason-1, ERS-2 и Geosat-FO,
а)
а)
б)
135' 140' 145' 150 155' 160' 105'
Рис. 3. Аномалии морской поверхности по данным Jason-1 и БЯ8-2 на 2003/06/01. а) метод «объективного анализа», б) сплайн интерполяция_
Рис. 4. Аномалии морской поверхности по данным Jason-1, ERS-2 и Geosat-FO на 2003/06/01. а) метод «объективного анализа», б) сплайн интерполяция._
соответственно, методами интерполяции Бш-сплайнами и «объективного анализа». На двух типах карт (рисунки 3а и 3б) присутствуют вихри с координатами центров: 42 с.ш. ,147 в.д.; 42.5 с.ш., 150 в.д.; 47 с.ш., 154 в.д. и 49 с.ш. 158.5 в.д. Очень похожие структуры восстановились в районе Курильской гряды. Координаты центров соответствующих вихревых образований следующие: 47.5 с.ш., 147.5 в.д., 47 с.ш., 151 в.д. Однако, в данной зоне, на карте, построенной при помощи метода «оптимальной интерполяции» присутствует еще вихрь с координатами центра 46.5 с.ш., 149 в.д., что является отличием между двумя картами в этом регионе. Причина кроется, скорее всего, в том, что данные карты построены при помощи данных только двух альтиметрических источников, а как показало проведенное моделирование, разрешающая способность метода сплайн интерполяции просто не позволила бы четко выразить данное образование, что, по-видимому, и произошло. Подтверждение данной гипотезе можно наблюдать на рис. 4б (карта аномалии морской поверхности, восстановленная методом сплайн интерполяции, с использованием данных трех альтиметрических источников) где вышеописанная структура восстановилась полностью.
Некоторое расхождение структуры поверхности океана, на картах восстановленных различными способами, наблюдается в прибрежных зонах выбранного региона (например, северовосточная часть Охотского моря). Это может быть вызвано следующими причинами:
• к сожалению, в прибрежных зонах отсутствуют данные спутника ERS-2 (рис. 5), следовательно, количество альтиметрических данных резко уменьшается, что может приводить к потере разрешающей способности, что мы и видим на рисунках 3б, 4б для карт построенных по методу сплайн интерполяции. По определению, метод оптимальной интерполяции опирается на предварительные статистические исследования, следовательно, этот метод менее восприимчив к количеству спутниковых альтиметрических источников, что и подтверждают работы [4, 5] где переход от двух источников к трем гораздо менее заметен, чем для метода интерполяции Dm-сплайнами;
• исходные альтиметрические данные, которые использовались для восстановления карт на рисунках 3 и 4, представляют собой так называемую промежуточную (interim gdr) информацию, которая и отличается от финальной тем, что содержит аппроксимацию орбиты спутника худшего качества. Соответственно, подобные, ложные структуры, могут быть вызваны ошибками при определении орбиты. Сплайны четко отобразили структуру, описанную исходными данными, а «оптимальный интерполянт» мог их «не заметить» в силу своей природы, что, вообще говоря, можно отнести к заслугам метода, только в том случае, если исходные статистические исследования абсолютно верны [11].
135" 140" 145" 150" 155" 160" 165"
Рис. 5. Распределение данных альтиметрических спутников в регионе исследования. Круги - Jason -1, квадраты - ERS-2, треугольники - Geosat - FO.
В целом оба метода восстановления информации показывают сходный результат, причем для восстановления карты аномалий поверхности океана методом сплайн интерполяции требуется затратить порядка минуты на стандартном персональном компьютере, а при использовании метода «оптимальной интерполяции» необходимо около трех минут, на мощной рабочей станции, не считая времени, необходимого для проведения статистических исследований.
Электронный журнал «Исследовано в России»
Выводы
Найдены значения оптимальных параметров методики восстановления аномалий морской поверхности в выбранном регионе мирового океана: Ям = 1.4, Яб = 6 и N=150.
Результаты моделирования подтверждают, что невозможно восстанавливать высокоэнергетичные мезомасштабные структуры с хорошей точностью, основываясь на данных только одной альтиметрической миссии, что соответствует результатам работ [4, 5]. Ошибки восстановления при использовании только шаблона данных спутниковой системы Topex/Poseidon при восстановлении (1) составили: 5С = 6.52, 5Ь =1.38 и невязка 7.21.
В случае использования данных двух спутников Topex/Poseidon и БЯ8-2 метод продемонстрировал хорошие возможности по восстановлению функций. Ошибки восстановления составили 5С = 1.41, дь =0.46 и невязка 2.37. Для исходных данных трех спутников ошибки
составили 5С = 1.12, дь =0.25 и невязка 1.23 соответственно (при восстановлении (1)). Ошибка по метрике пространства Сх бесспорно велика, т. е существуют довольно серьезные ошибки при восстановлении, но их количество чрезвычайно незначительно, так как ошибка по метрике 5Ь не превышает 0.5, а средняя невязка составила всего 8% для случая двух спутников и порядка 3% для случая трех, что свидетельствует о достаточно успешной реконструкции заданной поверхности.
При исследовании восстановления (2) ошибки составили для случая двух спутников 5С =
0.28, 5Ь =0.18 и невязка 0.14, а для трех - 5С = 0.22, 5Ь =0.16 и невязка 0.11, соответственно.
Неплохое соответствие наблюдается при визуальном сравнении карт построенных по реальным спутниковым данным, причем добавление третьего источника альтиметрических данных оказывает серьезное влияние на качество восстановления структур на поверхности океана, что особенно заметно при использовании метода сплайн интерполяции.
Если использовать метод интерполяции Бт-сплайнами, то необходимо затратить примерно в три раза меньше времени на воссоздание аномалий морской поверхности в данном регионе, чем, при использовании метода «объективного анализа» для того же набора данных.
Электронный журнал «Исследовано в России» 463 http://zhurnal.ape.relarn.ru/articles/2004/042.html
Список литературы
1. LeTraon, P. Y., Nadal, F., Ducet, N. // J. Atm. Ocean. Tech. 1998. N. 15. P 522
2. Le Traon, P.Y. and Dibarboure G. // J. Atm. Ocean. Tech. 1996. N. 16. P 1208
3. Shen et. al // J. Phys. Oceangr. 1986. N. 16. P454
4. Greenslade et al. // J. Atm. Ocean. Techl. 1997. N. 14. P 849
5. Le Traon, P.Y., Faugere Y., Hernandez F., et. al // J. Atm. Ocean.Tech. 2002. N. 20. P 889
6. L-L Fue, Anny Cazenave Satellite Altimetry and Earth Science. A handbook of techniques and applications, Academic Press, 2001
7. Le Traon, P. Y. and Ogor, F. // J. Geophys. Res. 1998. V.95. N. C12. P 8045
8. Le Traon, P. Y., Gaspar P., Bouyssel F., et. al // J. Atm. Ocean. Tech. 1995 N. 12, 161-170, 1995.
9. Visser, P., ERS-1 Precise Orbit Determination Using TOPEX/ERS-1 Dual Satellite Crossover Differences, CSR Tech. Memo. CSR-TM-93-07, October 1993.
10. Nouel F., Berthias J. P., Deleuze M., et. al // J. Geophys. Res. 1994. V. 99. N. C12. P 24,405
11. Bretherton F.P., Davis R. E., Fandry C.B. // Deep Sea Res. 1976. N. 23. P 559
12. Gandin L. S. Objective analysis of meteorological fields, Israel Program for Scientific Translations, 1965
13. Завьялов Ю.С., Квасов Б.И., Мирошниченко В.Л. Методы сплайн-функций. М.: Наука, 1980
14. Василенко В.А. Сплайн функции: теории, алгоритмы, программы. Новосибирск: Наука, 1983
15. Ковалков А.В. Функции Грина и сплайн-аппроксимация в многомерных областях. Новосибирск, 1980
16. Ашкеназы В.О. Сплайн-поверхности: Методическая разработка. Тверь: Тверской гос. ун-т, 1993
17. Куницын В. Е., Романов А.А. // Радиотехника и Электроника. № 4. 2004. (в печати)
18. Benada R., PO.DAAC merged GDR (T/P) users handbook. Rep. JPL D-11007. Pasadena: Jet Propul. Lab., 1993
19. Le Traon P. - Y., Hernandez F. // J. Atm. Ocean. Tech. 1992. N. 9. P 687
20. Le Traon, P. - Y., Dibarboure G., Ducet N. // J. Atm. Ocean. Tech. 2001. N. 18. P 1277