вторник, 7 августа 2018 г.

Определение нормали ЦМР и уклона в заданном направлении

Для задач ГИС, использующих уклоны, их величины и направления, необходимо по растровым данным ЦМР (цифровая модель рельефа) определять условную плоскость и ее нормаль в заданной точке.

Поверхность ЦМР задается в виде матрицы высот, с заданным размером пиксела по горизонтальным осям $X$ и $Y$ и содержащей значения высоты $Z$.

Методы оценивания в точке используют окрестность точки размером $3\times 3$ с центром в заданной точке. Обозначим значения высот через $z_i$: $$ \begin{array}{ccc} z_1 & z_2 & z_3 \\ z_4 & z_5 & z_6 \\ z_7 & z_8 & z_9 \end{array} $$ Здесь $z_5$ обозначает высоту центральной точки, в которой определяется уклон и направление нормали. Будем полагать что ось X направлена слева направо, а ось Y снизу вверх как это принято в математике. Нужно отметить что при применении в реальных системах оси и направление отсчетов нужно учитывать согласно правилам этих систем. В машинной графике традиционно ось X направлена слева направо, а ось Y сверзу вниз, а в географических и геодезических координатах в зависимости от используемой системы координат, например в системе координат Гаусса-Крюгера ось X направлена снизу вверх (по меридиану от экватора к полюсу) и ось Y направлена слева направо (по параллели с запада на восток).

Охват центральной точки бОльшим количеством пикселов, например, $5 \times 5$, не используют по причине того, что ЦМР задаются пикселами, имеющими относительно большой размер на Земле, и пикселы далее, чем непосредственно окружающие, не считаются характеризующими уклон в центральной точке.

По методике, предложенной в работе Зевенбергена-Торна в 1987, используется воспроизведение частных производных поверхности в точке на основе креста прилегающих пикселов $$ \begin{array}{ccc} & z_2 & \\ z_4 & z_5 & z_6 \\ & z_8 & \end{array} $$ Уклону в направлении оси X соответствует величина $$ \frac{\partial z}{\partial x} = \frac{1}{2}\left( \frac{z_6 - z_5}{l} + \frac{z_5 - z_4}{l} \right) $$ и уклону в направлении оси Y величина $$ \frac{\partial z}{\partial y} = \frac{1}{2}\left( \frac{z_2 - z_5}{l} + \frac{z_5 - z_8}{l} \right) $$ Здесь l -- размер пиксела по осям X и Y в тех же единицах что и отсчет высот по оси Z.

В формуле уклонов значения в центральном пикселе сокращаются: $$ \frac{\partial z}{\partial x} = \frac{1}{2}\left( \frac{z_6 - z_4}{l} \right) $$ $$ \frac{\partial z}{\partial y} = \frac{1}{2}\left( \frac{z_2 - z_8}{l} \right) $$ В программе ArcGIS уклоны вычисляются с учетом также угловых пикселов $z_1$, $z_3$, $z_7$ и $z9$. Значения в угловых пикселах считаются не в равных пропорциях к крестовым (горизонтальным и вертикальным соседям), и используется формула: $$ \frac{\partial z}{\partial x} = \frac{(z_3 + 2z_6 + z_9) - (z_1 + 2z_4 + z_7)}{8l} $$ $$ \frac{\partial z}{\partial y} = \frac{(z_1 + 2z_2 + z_3) - (z_7 + 2z_8 + z_9)}{8l} $$ Здесь в формуле уже участвуют угловые пикселы, но взвешиваются неравномерно по отношению к горизонтальным и вертикальным.

Средневзвешенный расчет уклонов предполагает что пикселы справа, слева, сверху и снизу образуют линию, и в качестве значения пиксела использует усреднение из трех, воспроизводя билинейную интерполяцию. При средневзвешенной оценке уклоны имеют вид: $$ \frac{\partial z}{\partial x} = \frac{1}{2}\left( \frac{\frac{z_3 + z_6 + z_9}{3} - \frac{z_1 + z_4 + z_7}{3}}{l} \right) $$ $$ \frac{\partial z}{\partial y} = \frac{1}{2}\left( \frac{\frac{z_1 + z_2 + z_3}{3} - \frac{z_7 + z_8 + z_9}{3}}{l} \right) $$ Разные системы воспроизведения плоскости по ЦМР также могут иметь разные трактовки как вычислять уклоны для краевых пикселей. Обозначим через NaN специальное значение Not-A-Number, обозначающее пиксел выходящий за пределы предоставленной матрицы высот. Могут быть два случая выхода охватывающей матрицы на край матрицы ЦМР -- выход на линейную границу и выход на угловую границу.

Выход на линейную границу, например на верхнюю границу матрицы высот, определяется: $$ \begin{array}{ccc} NaN & NaN & NaN \\ z_4 & z_5 & z_6 \\ z_7 & z_8 & z_9 \end{array} $$ Для этого случая расчет средневзвешенного уклона по горизонтали использует лишь горизонтальные пикселы и по вертикали лишь одну примыкающую (нижнюю) линейку: $$ \frac{\partial z}{\partial x} = \frac{1}{2}\left( \frac{z_6}{l} - \frac{z_4}{3l} \right) $$ $$ \frac{\partial z}{\partial y} = \left( \frac{z_5 - \frac{z_7 + z_8 + z_9}{3}}{l} \right) $$ Аналогичным образом должны строиться уклоны по осям X и Y при выходе центрального пиксела на другие линейные границы.

Выход на угловую границу, например на верхний левый, определяется: $$ \begin{array}{ccc} NaN & NaN & NaN \\ NaN & z_5 & z_6 \\ NaN & z_8 & z_9 \end{array} $$ В этом случае для определения уклонов используются лишь горизонтальные и вертикальные примыкающие пикселы: $$ \frac{\partial z}{\partial x} = \frac{z_6 - z_5}{l} $$ $$ \frac{\partial z}{\partial y} = \frac{z_5 - z_8}{l} $$ Аналогичным образом должны строиться уклоны по осям X и Y при выходе на другие углы. Нужно отметить, что при выходе на угол не учитываются диагональные пикселы и при выходе на линейную границу диагональные пикселы учитываются лишь частично.

В программе ArcGIS случай выхода центрального пиксела на линейную и угловую границу рассматривается иначе. Производится неявное дополнение отсутствующих пикселов центральным значением и рассчет уклонов производится как для полной матрицы окружения. Выход на верхнюю линейную границу эквивалентен дополнению: $$ \begin{array}{ccc} z_5 & z_5 & z_5 \\ z_4 & z_5 & z_6 \\ z_7 & z_8 & z_9 \end{array} $$ И выход на верхний левый угол эквивалентен дополнению: $$ \begin{array}{ccc} z_5 & z_5 & z_5 \\ z_5 & z_5 & z_6 \\ z_5 & z_8 & z_9 \end{array} $$ По этой причине при работе в ArcGIS возникают гарантированные погрешности для оценки уклонов на краевых пикселах. При средневзвешенной оценке такая погрешность не является гарантированной.

Далее будем полагать, что сделан выбор метода оценивания уклонов по осям X и Y и способ определения уклонов для краевых пикселов. Будем полагать вычисленными значения $$ \frac{\partial z}{\partial x} = \lambda_x $$ $$ \frac{\partial z}{\partial y} = \lambda_y $$ По этим значениям далее надо воспроизвести плоскость. Метод воспроизведения будет полагать, что если в пространстве 3D есть две пересекающиеся прямые перпендикулярные друг другу, то через них всегда можно провести плоскость. Систематическая погрешность определения плоскости через уклоны как частные производные по осям состоит в том, что возможны совершенно разные матрицы высот ЦМР, дающие одинаковые оценки уклонов по осям, например отличающиеся расположением угловых пикселов. В частности, матрицы $$ \begin{array}{ccc} 1 & 3 & 5 \\ 3 & 3 & 3 \\ 5 & 3 & 1 \end{array} $$ $$ \begin{array}{ccc} 5 & 3 & 1 \\ 3 & 3 & 3 \\ 1 & 3 & 5 \end{array} $$ дают одинаковые значения уклонов по осям и воспроизведение плоскости по уклонам будет одинаковым, в то время как характер наклонов по диагоналям совершенно противоположный.

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

Получив значения $\lambda_x$ и $\lambda_y$, воспроизведем по ним плоскость в общем виде, чтобы впоследствии получить уклон в заданном направлении.

Двумерная плоскость в трехмерном пространстве имеет каноническое уравнение $$ ax + by + cz + d = 0 $$ Здесь коэффициенты a, b и c образуют коэффициенты вектора нормали к плоскости и значение d характеризует насколько плоскость отстоит от начала координат в направлении нормали.

Поместим воспроизводимую плоскость так чтобы в начало координат попал центральный пиксел $z_5$, для этого проведем коррекцию: $$ \begin{array}{ccc} z_1 - z_5 & z_2 - z_5 & z_3 - z_5 \\ z_4 - z_5 & 0 & z_6 - z_5 \\ z_7 - z_5 & z_8 - z_5 & z_9 - z_5 \end{array} $$ Если плоскость проходит через начало координат, а в начале координат значения x, y и z равны нулю, то в этом случае $$ a0 + b0 + c0 + d = 0 $$ И для этого случая плоскость не отстоит ни на какой расстояние от начала координат: $$ d = 0 $$ $$ ax + by + cz = 0 $$ Все три коэффициента a, b и c имеют возможность быть умноженными на любое ненулевое число. Эти три коэффициента задают направление вектора нормали, и обычно в качестве вектора нормали используют единичный вектор, получаемый нормированием: $$ \frac{ai + bj + ck}{\sqrt{a^2 + b^2 + c^2}} $$ Если нам известно сечение по оси X (вычисленное ранее значение $\lambda_x$), то оно считается проходящим в плоскости через начало отсчета, следовательно значение Y равно 0 и уравнение плоскости принимает вид: $$ ax + b0 + cz = 0 $$ Следовательно, коэффициенты a и c связаны с вычисленным ранее уклоном в направлении оси X: $$ \frac{a}{c} = \lambda_x $$ Аналогично для уклона в направлении оси Y значение X равно 0: $$ a0 + by + cz = 0 $$ Следовательно, коэффициенты b и c связаны с вычисленным ранее уклоном в направлении оси Y: $$ \frac{b}{c} = \lambda_y $$ Поскольку значения a, b и c всегда могут быть умножены и разделены на одно и то же ненулевое число, мы можем выбрать для решения системы уравнений $$ \left\{ \begin{aligned} & \frac{a}{c} = \lambda_x \\ & \frac{b}{c} = \lambda_y \end{aligned} \right. $$ значение c = 1, и значения коэффициентов нормали плоскости определять как $$ \lambda_xx + \lambda_yy + c = 0 $$ При этом нормированная нормаль имеет значение $$ n_xi + n_yj + n_zk = \frac{\lambda_xi + \lambda_yj + k}{\sqrt{\lambda_x^2 + \lambda_y^2 + 1}} $$ Далее, имея каноническую нормаль, надо определить уклон в направлении, задаваемом горизонтальной составляющей вектора направления. Если уклон задается матрицей ЦМР, то направление уклона задается вектором ветра. У ветра не используют составляющую по высоте и считают что вектор ветра задается лишь в плоскости осей X и Y: $$ v = v_xi + v_yj $$ Поскольку рассматриваемая нами воспроизводимая плоскость проходит через начало отсчета, от этого же начала отсчета отложим вектор v. Проекция его концевой точки на плоскость должна удовлетворять уравнению плоскости. Поэтому высота точки проекции на плоскость удовлетворяет уравнению $$ \lambda_xv_x + \lambda_yv_y + z = 0 $$ Вектор скорости и точка поднятия на плоскости образуют треугольник с основанием (горизонтальный катет) $$ v_xi + v_yj $$ имеющий величину $\sqrt{v_x^2 + v_y^2}$

Второй, горизонтальный, катет имеет высоту $$ z = - \lambda_xv_x - \lambda_yv_y $$ Поэтому угол наклона плоскости в этом направлении определяется через арктангенс отношения катетов: $$ \alpha = arctg\left( \frac{- \lambda_xv_x - \lambda_yv_y}{ \sqrt{v_x^2 + v_y^2} } \right) $$

PS. Тег "Космонавтика" ставлю потому что работа выполнялась в рамках работ ДЗЗ (дистанционное зондирование Земли) для задачи перколяционного прогнозирования распространения очагов лесных пожаров. Это часть задачи анализа информации со спутниковых снимков. В предыдущей методике МЧС рельев и направление уклонов по ветру не учитывались, только сила и направление ветра (в дополнение к классам горимости, пожароопасности погоды и вида пожара - верховой или низовой). В нынешней уже учитываются. Также стали учитываться данные по влажности. Но, на мой взгляд, дополнительно нужно учитывать также и точку росы как результат оценивания температуры, давления и влажности (на ночные часы прогнозируемого периода). Также другими средствами прогнозируется класс пожароопасности погоды по предыдущим осадкам и температуре, эти материалы легко найти в Интернете.

Комментариев нет:

Отправить комментарий