Извлеките лесные метрики и отдельные древовидные атрибуты из воздушных данных о лидаре

В этом примере показано, как извлечь лесные метрики, и отдельное дерево приписывает из воздушных данных о лидаре.

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

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

Загрузите и визуализируйте данные

Разархивируйте forestData.zip к временной директории и загрузке данные об облаке точек из файла LAZ, forestData.laz, в рабочую область MATLAB®. Данные получены из Открытого Набора данных Топографии [1]. Облако точек, в основном, содержит точки растительности и земля. Загрузите данные об облаке точек в рабочую область с помощью readPointCloud функция lasFileReader объект. Визуализируйте облако точек с помощью pcshow функция.

dataFolder = fullfile(tempdir,"forestData",filesep);
dataFile = dataFolder + "forestData.laz";   
% Check whether the folder and data file already exist or not
folderExists = exist(dataFolder,'dir');
fileExists = exist(dataFile,'file');
% Create a new folder if it doesn't exist
if ~folderExists
    mkdir(dataFolder);
end
% Extract aerial data file if it doesn't exist
if ~fileExists
    unzip('forestData.zip',dataFolder);
end
% Read LAZ data from file
lasReader = lasFileReader(dataFile);
% Read point cloud along with corresponding scan angle information
[ptCloud,pointAttributes] = readPointCloud(lasReader,"Attributes","ScanAngle");
% Visualize the input point cloud
figure
pcshow(ptCloud.Location)
title("Input Point Cloud")

Земля сегмента

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

% Segment Ground and extract non-ground and ground points
groundPtsIdx = segmentGroundSMRF(ptCloud);
nonGroundPtCloud = select(ptCloud,~groundPtsIdx);
groundPtCloud = select(ptCloud,groundPtsIdx);
% Visualize non-ground and ground points in magenta and green, respectively
figure
pcshowpair(nonGroundPtCloud,groundPtCloud)
title("Segmented Non-Ground and Ground Points")

Нормируйте вертикальное изменение

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

Устраните дублирующиеся точки вдоль x-и осей Y, если таковые имеются, при помощи groupsummary функция.
Создайте interpolant использование scatteredInterpolant объект, чтобы оценить землю в каждой точке в данных об облаке точек.
Нормируйте вертикальное изменение каждой точки путем вычитания интерполированного наземного вертикального изменения из исходного вертикального изменения.

groundPoints = groundPtCloud.Location;
% Eliminate duplicate points along x- and y-axes
[uniqueZ,uniqueXY] = groupsummary(groundPoints(:,3),groundPoints(:,1:2),@mean);
uniqueXY = [uniqueXY{:}];
% Create interpolant and use it to estimate ground elevation at each point
F = scatteredInterpolant(double(uniqueXY),double(uniqueZ),"natural");
estElevation = F(double(ptCloud.Location(:,1)),double(ptCloud.Location(:,2)));
% Normalize elevation by ground
normalizedPoints = ptCloud.Location;
normalizedPoints(:,3) = normalizedPoints(:,3) - estElevation;
% Visualize normalized points
figure
pcshow(normalizedPoints)
title("Point Cloud with Normalized Elevation")

Извлеките лесные метрики

Извлеките лесные метрики из нормированных точек с помощью helperExtractForestMetrics функция помощника, присоединенная к этому примеру как вспомогательный файл. Функция помощника сначала делит облако точек на сетки на основе обеспеченного gridSize, и затем вычисляет лесные метрики. Функция помощника принимает что все точки с нормированной высотой ниже, чем cutoffHeight земля, и остающиеся точки являются растительностью. Вычислите эти лесные метрики.

Навес (CC) — Навес [2] является пропорцией леса, покрытого вертикальной проекцией древовидных корон. Вычислите его, когда отношение растительности возвращается относительно общего количества возвратов.
Часть разрыва (GF) — часть Разрыва [3] является вероятностью луча света, проходящего через навес, не сталкиваясь с листвой или другими элементами объекта. Вычислите его, когда отношение земли возвращается относительно общего количества возвратов.
Индекс [3] области Leaf area index (LAI) — Leaf является суммой односторонней листовой области на модуль земельного участка. Значение LAI вычисляется с помощью уравнения $LAI = - \frac{\cos (ang) * \ln (GF)}{k}$ , где ang средний угол сканирования, GF часть разрыва и k коэффициент исчезновения, который тесно связан с распределением листового угла.

% Set grid size to 10 meters per pixel and cutOffHeight to 2 meters
gridSize = 10;
cutOffHeight = 2;
leafAngDistribution = 0.5;
% Extract forest metrics
[canopyCover,gapFraction,leafAreaIndex] = helperExtractForestMetrics(normalizedPoints, ...
    pointAttributes.ScanAngle,gridSize,cutOffHeight,leafAngDistribution);
% Visualize forest metrics
hForestMetrics = figure;
axCC = subplot(2,2,1,Parent=hForestMetrics);
axCC.Position = [0.05 0.51 0.4 0.4];
imagesc(canopyCover,Parent=axCC)
title(axCC,"Canopy Cover")
axis off
colormap(gray)
axGF = subplot(2,2,2,Parent=hForestMetrics);
axGF.Position = [0.55 0.51 0.4 0.4];
imagesc(gapFraction,'Parent',axGF)
title(axGF,"Gap Fraction")
axis off
colormap(gray)
axLAI = subplot(2,2,[3 4],Parent=hForestMetrics);
axLAI.Position = [0.3 0.02 0.4 0.4];
imagesc(leafAreaIndex,Parent=axLAI)
title(axLAI,"Leaf Area Index")
axis off
colormap(gray)

Сгенерируйте Модель высоты навеса (CHM)

Модели высоты навеса (CHMs) являются растровыми представлениями высоты деревьев, созданий и других структур над землей топография. Используйте CHM в качестве входа для древовидного обнаружения и сегментации. Сгенерируйте CHM от своих нормированных значений вертикального изменения с помощью pc2dem функция.

% Set grid size to 0.5 meters per pixel 
gridRes = 0.5;
% Generate CHM
canopyModel = pc2dem(pointCloud(normalizedPoints),gridRes,CornerFillMethod="max");
% Clip invalid and negative CHM values to zero
canopyModel(isnan(canopyModel) | canopyModel<0) = 0;
% Perform gaussian smoothing to remove noise effects
H = fspecial("gaussian",[5 5],1);
canopyModel = imfilter(canopyModel,H,'replicate','same');
% Visualize CHM
figure
imagesc(canopyModel)
title('Canopy Height Model')
axis off
colormap(gray)

Обнаружьте вершины деревьев

Обнаружьте вершины деревьев с помощью helperDetectTreeTops функция помощника, присоединенная к этому примеру как вспомогательный файл. Функция помощника обнаруживает вершины деревьев путем нахождения локальных максимумов в переменных размерах окна [4] в CHM. Для обнаружения вершины дерева функция помощника рассматривает только вопросы с нормированной высотой, больше, чем minTreeHeight.

% Set minTreeHeight to 5 m 
minTreeHeight = 5;
% Detect tree tops
[treeTopRowId,treeTopColId] = helperDetectTreeTops(canopyModel,gridRes,minTreeHeight);
% Visualize treetops
figure
imagesc(canopyModel)
hold on
plot(treeTopColId,treeTopRowId,"rx",MarkerSize=3)
title("CHM with Detected Tree Tops")
axis off
colormap("gray")

Деревья индивидуума сегмента

Деревья индивидуума сегмента с помощью helperSegmentTrees функция помощника, присоединенная к этому примеру как вспомогательный файл. Функция помощника использует управляемую маркером сегментацию водораздела [5], чтобы сегментировать отдельные деревья. Во-первых, функция создает бинарное изображение маркера с местоположениями вершины дерева, обозначенными значением 1 . Затем функция фильтрует дополнение CHM наложением минимумов, чтобы удалить минимумы, которые не являются вершинами деревьев. Функция затем выполняет сегментацию водораздела на отфильтрованном дополнении CHM, чтобы сегментировать отдельные деревья. После сегментации визуализируйте отдельные древовидные сегменты.

% Segment individual trees
label2D = helperSegmentTrees(canopyModel,treeTopRowId,treeTopColId,minTreeHeight);
% Identify row and column id of each point in label2D and transfer labels
% to each point
rowId = ceil((ptCloud.Location(:,2) - ptCloud.YLimits(1))/gridRes) + 1;
colId = ceil((ptCloud.Location(:,1) - ptCloud.XLimits(1))/gridRes) + 1;
ind = sub2ind(size(label2D),rowId,colId);
label3D = label2D(ind);
% Extract valid labels and corresponding points
validSegIds = label3D ~= 0;
ptVeg = select(ptCloud,validSegIds);
veglabel3D = label3D(validSegIds);
% Assign color to each label
numColors = max(veglabel3D);
colorMap = randi([0 255],numColors,3)/255;
labelColors = label2rgb(veglabel3D,colorMap,OutputFormat="triplets");
% Visualize tree segments
figure
pcshow(ptVeg.Location,labelColors)
title("Individual Tree Segments")
view(2)

Извлеките древовидные атрибуты

Извлеките отдельные древовидные атрибуты с помощью helperExtractTreeMetrics функция помощника, присоединенная к этому примеру как вспомогательный файл. Во-первых, функция идентифицирует точки, принадлежащие отдельным деревьям от меток. Затем функция извлекает древовидные атрибуты, такие как древовидное местоположение вершины вдоль x-и осей Y, аппроксимированной древовидной высоты, древовидного диаметра короны и области. Функция помощника возвращает атрибуты как таблицу, где каждая строка представляет атрибуты отдельного дерева.

% Extract tree attributes
treeMetrics = helperExtractTreeMetrics(normalizedPoints,label3D);
% Display first 5 tree segments metrics
disp(head(treeMetrics,5));

    TreeId    NumPoints    TreeApexLocX    TreeApexLocY    TreeHeight    CrownDiameter    CrownArea
    ______    _________    ____________    ____________    __________    _____________    _________

      1          388        2.6867e+05      4.1719e+06       29.509          7.5325         44.562 
      2           22        2.6867e+05      4.1719e+06       21.464         0.99236        0.77344 
      3          243        2.6867e+05      4.1719e+06       24.201          5.7424         25.898 
      4          101        2.6867e+05      4.1719e+06       21.927          3.4571         9.3867 
      5           54        2.6867e+05      4.1719e+06       19.515          3.0407         7.2617

Ссылки

[1] Томпсон, Обзор Лидара Области с Ручьем. Иллилуетт, Йосемитская долина, CA 2018. Национальный Центр Бортового Лазера, Сопоставляющего (NCALM). Распределенный OpenTopography. https://doi.org/10.5069/G96M351N. Полученный доступ: 2021-05-14

[2] Мама, Цинь, Янджун Су и Цинхуа Го. “Сравнение Оценок Навеса От Бортового LiDAR, Воздушного Формирования изображений и Спутниковых снимков”. Журнал IEEE Выбранных Тем в Прикладных наблюдениях Земли и Дистанционном зондировании 10, № 9 (сентябрь 2017): 4225–36. https://doi.org/10.1109/JSTARS.2017.2711482.

[3] Ричардсон, Джеффри Дж., Л. Моника Москаль и Су-Хюнг Ким. "Моделирование Подходов, чтобы Оценить Эффективный Листовой индекс области от Дискретной Антенны - Возвращает ЛИДАР". Сельскохозяйственный и Лесная Метеорология 149, № 6-7 (июнь 2009): 1152–60. https://doi.org/10.1016/j.agrformet.2009.02.007.

[4] Pitkänen, J., М. Мэлтэмо, Дж. Хииппэ и Кс. Ю. "Методы адаптации для Отдельного Древовидного Обнаружения на Бортовой основанной на лазере Модели Высоты Навеса". Международные Архивы Фотограмметрии, Дистанционного зондирования и Пространственной Информатики 36, № 8 (январь 2004): 187–91.

[5] Чен, Ци, Деннис Болдоччи, Пэн Гун и Маджи Келли. “Изолируя Отдельные Деревья в Лесистой местности Саванны Используя Маленькие Данные о Лидаре Места”. Photogrammetric Engineering & Remote Sensing 72, № 8 (1 августа 2006): 923–32. https://doi.org/10.14358/PERS.72.8.923.

Документация