Использование данных SRTM в моделировании земной поверхности. Определение стоимости отрезка пути в рамках трехмерной модели рельефа.

Воловиков Олег Павлович

аспирант Кубанского государственного технологического университета

Одной из задач исследования Земли из космоса является подробное определение ее формы. Национальное управление США по аэронавтике и исследованию космического пространства (NASA) имеет в списке своих миссий исследование топографии с использованием космического радара. За время работы проекта спутникового радарного сканирования, имеющего в списке NASA название Shuttle Radar Topography Mission (сокращенно SRTM), в феврале 2000 г. были собраны данные о высотах в различных областях Земли. Подробнее о миссии можно узнать на сайте NASA по адресу http://www2.jpl.nasa.gov/srtm.

При помощи того же интернет-сайта можно получить обработанные данные результатов SRTM-сканирования. Данные представляют собой заархивированные файлы с именами NxxEyyy.hgt.zip, где xx – координата по широте, yyy – координата по долготе левого верхнего угла квадрата поверхности, представленного в файле. Буквы N и E в имени файла обозначают соответственно северное и восточное полушария. В архиве имеется одноименный файл с расширением .hgt, содержащий собственно результаты исследования.

Каждый ght-файл имеет размер 2884802 байт и представляет собой матрицу размером 1201x1201 из двухбайтовых значений. Каждому двухбайтовому элементу матрицы соответствует высота в метрах над уровнем моря узла решётки с шагом в три секунды дуги (1/1200 градуса) как по широте, так и по долготе. На цилиндрической равноугольной проекции Меркатора эта решётка прямоугольная, каждому файлу соответствует прямоугольник градус на градус. Отсчёт узлов в файле идёт с запада на восток и с севера на юг, строкам матрицы соответствуют параллели, столбцам – меридианы. Один файл покрывает квадрат геоида со сторонами в 1 градус, а последний 1201-й столбец и последняя 1201-я строка образуют перекрытие со смежными квадратами. Таким образом, в первую очередь выписаны данные вдоль самой северной широты (т.е. широта в названии файла плюс один градус), это 2402-байтный блок (строка матрицы). Затем идёт строка, которая на три секунды дуги ближе к югу. Затем следующая, и т.д. до 1201-й, самой южной, соответствующей названию файла. В пределах каждой строки долгота увеличивается.

Очерёдность байтов в файле обратная: сначала идёт старший байт числа, затем младший. Т.е., например, если весь файл прочитать в двухмерный массив 1200x1200 элементов типа word, а потом над каждым элементом массива для получения верных значений произвести операцию перестановки байтов, мы получим карту высот участка поверхности 1x1 градус с разрешением в 3 угловых секунды по обоим направлением (для территории России это означает шаг примерно в 60-70 м). Для быстрой перестановки байтов удобно использовать инструкции процессора x86 (функция на языке Pascal):

function sh(x: smallint): smallint;

begin

asm

mov ax,x

ror ax,8

mov result,ax

end;

end;

Полученный массив значений можно непосредственно использовать в качестве источника данных для схемы базы данных, описанной в «Журнале научных публикаций аспирантов и докторантов» (№10 за октябрь 2008 г.) в статье «Описание рельефа земной поверхности с помощью реляционной модели». При загрузке данных необходимо лишь учесть, что значение высоты, равное -32768, выражает неопределенную высоту (следует установить значение NULL в базе данных). Такое значение возникает в ситуациях, когда радар не смог точно определить высоту в данной точке.

Рассмотрим наполненную данными схему БД в качестве инструмента для определения стоимости кратчайшего пути между двумя произвольно выбранными точками на поверхности Земли. Для этого зададим две константы-параметра:

1.      dist_cost (стоимость перемещения на 1 м по горизонтали);

2.      alt_cost (стоимость перемещения на 1 м по высоте).

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

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

C = dist_cost · D + alt_cost · A,                                                                               (1)

где D – линейная дистанция между точками, A – разность высот.

Известно, что длина дуги в радианах, наикратчайшим образом соединяющая две точки на поверхности сферы, выражается зависимостью [1, 2]:

cos(d) = sin(φ_А)·sin(φ_B) + cos(φ_А)·cos(φ_B)·cos(λ_А – λ_B),                                            (2)

где φ_А и φ_B – широты, λ_А, λ_B – долготы данных пунктов, d – расстояние между пунктами, измеряемое в радианах длиной дуги большого круга земного шара. Для расчета расстояния между пунктами, расположенными в разных полушариях (северное-южное, восточное-западное), знаки (±) у соответствующих параметров (широт или долгот) должны быть разными. Расстояние между пунктами, измеряемое в метрах, определяется по формуле [2]:

L = d · R,                                                                                                                  (3)

где R – средний радиус земного шара (по современным представлениям около 6371023 м).

Тогда кратчайшее расстояние между двумя точками будет описываться формулой

L = R · acos(sin(φ_А)·sin(φ_B) + cos(φ_А)·cos(φ_B)·cos(λ_А – λ_B))                                     (4)

Разницу высот между смужными узлами сетки можно определить непосредственно из отношения Landscape по атрибуту height.

Пусть дистанция между точками выражается функцией distance(@point_a, @point_b), а высота точки функцией altitude(@point_id). Тогда суммарная стоимость C_summ пути P[i₁…i_n] , будет равна

C_summ = ∑_{i = 1..n} C(i, i + 1),                                                                                        (5)

в то время как C(i, i + 1) может быть определена по формуле

C(i, i + 1) = distance(i, i + 1) · dist_cost + ABS(altitude(i) – altitude(i + 1)) · alt_cost, (6)

где ABS – функция вычисления модуля числа. Для случая значений высоты NULL следует предусмотреть отдельный способ определения стоимости. Например, можно умножить значение (distance(i, i + 1) · dist_cost) на некий коэффициент сложности.

Сам кратчайший путь PATH может быть определен по следующему алгоритму:

1. Пусть P1 – начало, P2 – конец пути. Присвоим Pcur = P1, Pnext = NULL, Dmin = distance(P1, P2), N = 1, PATH[1] = P1.
2. Если Pcur ≠ P2 выполнить:
1. найти смежную с Pcur точку Pnext так, чтобы значение distance(Pnext, P2) было минимальным;
2. если distance(Pnext, P2) < Dmin то присвоить Dmin = distance(x, P2), иначе пути не существует – покидаем алгоритм с генерацией ошибки;
3. Присвоить N = N + 1; присвоить PATH[N] = Pnext;
4. присвоить Pcur = Pnext;
3. Конец алгоритма. Результат содержится в PATH (длина пути равна N).

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

Литература.

1. Лавский В.М., Кондратьев Н.Я., Жарков И.Н., Романов М.К., Буланов В.П., Лисодет В.Н., Исайчев В.В. Справочник авиационного штурмана. – М.: Военное издательство военного министерства СССР, 1952.

2. Определение расстояний на поверхности Земли. Код доступа: http://tvsh2004.narod.ru/geo_koor.htm.

Поступила в редакцию 02.02.2009 г.