В этом примере показано, как вычислить и визуализировать видимость линии видимости самолета с земли по местности. Во-первых, импортируйте данные местности для области и примените их к 3-D географическому глобусу. Затем выполните анализ видимости «точка-точка» от наземного местоположения до моделируемого угла тангажа и отобразите результаты на 3-D географическом глобусе. Наконец, выполните анализ видимости «точка-область» из наземного местоположения, соответствующего полету самолета на постоянной высоте, и отобразите результаты на 2-D географической оси.
Используйте анализ линии визирования для сценариев «земля-воздух», где важна беспрепятственная видимость, например, для радиолокационного наблюдения, связи и планирования пути БПЛА. Этот пример применяет анализ к радиолокационному наблюдению для аэропорта.
Укажите файл местности в формате 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 визуализацией.
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 картой.
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")
addCustomTerrain | campos | geoglobe | geoplot3 | los2 | readgeoraster