Импорт данных местности

Укажите файл местности в формате DTED для анализа данных и 3-D визуализации. Файл местности был загружен из набора данных «SRTM Void Fill», доступного из Геологической службы США (USGS).

dtedfile = "n39_w106_3arc_v2.dt1";
attribution = "SRTM 3 arc-second resolution. Data available from the U.S. Geological Survey.";

Импортируйте данные файла DTED в рабочую область как массив и географический объект растровой привязки, задавая тип возврата как double, чтобы данные работали со всеми функциями анализа.

[Zterrain,Rterrain] = readgeoraster(dtedfile,"OutputType","double");

Просмотрите географические пределы и разрешение выборки данных местности путем доступа к свойствам географического объекта растровой привязки. Пределы для файла соответствуют области вокруг Боулдера, Колорадо, США, и разрешение соответствует формату DTED уровня 1, который имеет разрешение выборки 3 угловых секунды, или около 90 метров.

latlim = Rterrain.LatitudeLimits;
lonlim = Rterrain.LongitudeLimits;
latspc = Rterrain.SampleSpacingInLatitude;
lonspc = Rterrain.SampleSpacingInLongitude;
disp("Latitude limits of terrain: " + mat2str(latlim) + newline + ...
    "Longitude limits of terrain: " + mat2str(lonlim) + newline + ...
    "Terrain resolution in latitude: " + latspc*3600 + " arc seconds" + newline + ...
    "Terrain resolution in longitude: " + lonspc*3600 + " arc seconds")

Latitude limits of terrain: [39 40]
Longitude limits of terrain: [-106 -105]
Terrain resolution in latitude: 3 arc seconds
Terrain resolution in longitude: 3 arc seconds

Визуализация линии зрения траектории самолета на карте 3-D

Создайте географический глобус с пользовательским рельефом

Добавьте пользовательские данные с помощью файла для использования с 3-D визуализацией.

addCustomTerrain("southboulder",dtedfile,"Attribution",attribution)

Задайте пользовательский рельеф с новым географическим глобусом. Сохраните пользовательский рельеф земного шара при добавлении данных путем установки hold on.

fig = uifigure;
g = geoglobe(fig,"Terrain","southboulder");
hold(g,"on")

Просмотр радиолокационного местоположения Земли

Определите радиолокационное местоположение земли в столичном аэропорту Рокки-Маунтин. Радар установлен на башне в 10 метрах над землей. Высота радара является суммой повышения земли и высоты радиолокационной башни, обозначаемой как средний уровень моря.

rdrlat = 39.913756;
rdrlon = -105.118062;
rdrtowerht = 10;
rdralt = 1717 + rdrtowerht;

Постройте график местоположения радара на географическом глобусе.

geoplot3(g,rdrlat,rdrlon,rdralt,"co", ...
    "LineWidth",6, ...
    "MarkerSize",1)

Симулируйте траекторию самолета

Симулируйте траекторию самолета, кружащего над горами.

Определите центральное положение кругового самолета.

tlat0 = 39.80384;
tlon0 = -105.49916;
tht0 = 3000;

Определить траекторные путевые точки для самолета с помощью Декартовых координат восток-север-вверх (ENU). Задайте кривую радиусом 5 км (5000 м) и вертикальным смещением 1 км (1000 м) на 1,5 оборота. Затем преобразуйте координаты ENU в геодезические координаты, которые привязываются к WGS84 эллипсоиду.

azs = 1:540;
r = 5000;
[X,Y] = pol2cart(deg2rad(azs),r);
Z = linspace(0,1000,numel(azs));
wgs84 = wgs84Ellipsoid;
[tlat,tlon,tht] = enu2geodetic(X,Y,Z,tlat0,tlon0,tht0,wgs84);

Просмотр траектории самолета по местности

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

traj = geoplot3(g,tlat,tlon,tht,"y", ...
    "HeightReference","ellipsoid", ...
    "LineWidth",3);

Просмотр 3-D местности и местоположения радара на расстоянии путем изменения положения камеры и углов поворота.

campos(g,39.77114,-105.62662,6670)
camheading(g,70)
campitch(g,-12)

Вычисление видимости линии визирования с помощью траектории самолета

Вычисление видимости линии визирования с помощью los2 функцию и данные DTED.

The los2 функция поддерживает ортометрическую высоту (высота над средним уровнем моря) или высоту над уровнем земли. Преобразуйте высоты траектории самолета с эллипсоидальной высоты на ортометрическую высоту. Затем вычислите линию зрения от местоположения радара аэропорта до каждой точки пути траектории самолета и преобразуйте результаты в логический массив.

numwaypts = numel(tlat);
isvis = zeros(1,numwaypts);
talt = tht - egm96geoid(tlat,tlon);
for k = 1:numwaypts
    isvis(k) = los2(Zterrain,Rterrain,rdrlat,rdrlon,tlat(k),tlon(k),rdralt,talt(k),"MSL","MSL");
end
isvis = logical(isvis);

Обратите внимание, что los2 вычисляет видимость линии визирования, принимая, что данные ссылаются на сферическую Землю, в то время как данные фактически ссылаются на WGS84 эллипсоид, и в результате могут быть незначительные расхождения. Расчет линии видимости также соответствует оптической линии видимости и не учитывает преломление через атмосферу.

Визуализация видимости линии зрения на местности

Постройте график видимости линии видимости. Используйте зеленые маркеры, где самолет виден из аэропорта, и пурпурные маркеры, где он не виден.

delete(traj)
geoplot3(g,tlat(isvis),tlon(isvis),tht(isvis),"og", ...
    "HeightReference","ellipsoid", ...
    "LineWidth",2, ...
    "MarkerSize",1)
geoplot3(g,tlat(~isvis),tlon(~isvis),tht(~isvis),"om", ...
    "HeightReference","ellipsoid", ...
    "LineWidth",2, ...
    "MarkerSize",1)

Просмотрите график видимости с точки зрения аэропорта. Получите геодезические координаты положения на 900 метров восточнее, 200 метров севернее и 100 метров вверх от радиолокационного местоположения. Затем установите положение камеры и углы поворота. Зеленые маркеры появляются на виду, но пурпурные маркеры либо полностью, либо частично препятствуют местности.

rdrht = rdralt + egm96geoid(rdrlat,rdrlon);
[camlat,camlon,camht] = enu2geodetic(900,200,100,rdrlat,rdrlon,rdrht,wgs84);
campos(g,camlat,camlon,camht)
camheading(g,-110)
campitch(g,0)

Визуализация контуров видимости линии видимости самолета на карте 2-D

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

Графическое изображение радиолокационного местоположения и пределов местности на карте 2-D

Постройте график местоположения радара на новом рисунке с топографической 2-D картой.

figure
geoplot(rdrlat,rdrlon,"co", ...
    "LineWidth",6, ...
    "MarkerSize",3, ...
    "DisplayName","Radar location")
hold on
geobasemap topographic
gx = gca;
gx.InnerPosition = gx.OuterPosition;

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

latmin = latlim(1);
latmax = latlim(2);
lonmin = lonlim(1);
lonmax = lonlim(2);
geoplot([latmin latmin latmax latmax latmin],[lonmin lonmax lonmax lonmin lonmin], ...
    "LineWidth",1, ...
    "Color","k", ...
    "DisplayName","Terrain limits")

Отображение легенды в северо-западном углу.

legend("Location","northwest")

График контуров видимости для самолета, летящего на постоянной высоте

Указание трех высот над средним уровнем моря для самолета. Для каждой высоты:

tgtalts = [3000 4000 5000];

minVertices = 10;
cfig = figure("Visible","off"); % Suppress contour plot using invisible figure
cax = axes("Parent",cfig);
for tgtalt = tgtalts
    vis = viewshed(Zterrain,Rterrain,rdrlat,rdrlon,rdralt,tgtalt,"MSL","MSL");
    
    C = contourm(vis,Rterrain,"LevelList",1,"Parent",cax);
    clat = C(2,:);
    clon = C(1,:);
    
    clats = [];
    clons = []; 
    k = 1;   
    while k < size(C,2)
        numVertices = clat(k);
        if numVertices > minVertices % Do not plot small segments 
            clats = [clats clat(k+1:k+numVertices) NaN]; %#ok<AGROW> 
            clons = [clons clon(k+1:k+numVertices) NaN]; %#ok<AGROW> 
        end
        k = k + numVertices + 1;
    end
    
    geoplot(gx,clats,clons,"LineWidth",2, ...
        "DisplayName", "Aircraft: " + string(tgtalt) + " m");
end

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

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

Как и в случае los2, а viewshed функция вычисляет видимость линии видимости, принимая, что данные ссылаются на сферическую Землю, в то время как данные фактически ссылаются на WGS84 эллипсоид, и в результате могут быть незначительные расхождения. Расчет линии видимости также соответствует оптической линии видимости и не учитывает преломление через атмосферу.

Очистка

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

if isvalid(fig)
    close(fig)
end
removeCustomTerrain("southboulder")