Obtencion y analisis de variables dendrometricas obtenidas mediante datos lidar en el ambito...

INGENIERÍA SUPERIOR EN GEODESIA Y

CARTOGRAFÍA

PROYECTO FINAL DE CARRERA

OBTENCIÓN Y ANÁLISIS DE VARIABLES

DENDROMÉTRICAS OBTENIDAS MEDIANTE

DATOS LIDAR EN ÁMBITO AGRÍCOLA

ESCUELA TÉCNICA SUPERIOR DE INGENIERÍA GEODÉSICA, CARTOGRÁFICA

Y TOPOGRÁFICA

Autor: Rubén Llorens Alonso

Tutores:

Luis Ángel Ruiz Fernández, Javier Estornell Cremades

Universidad Politécnica de Valencia

Septiembre, 2015

2

Resumen

Para la realización de este proyecto se seleccionó del término municipal de Viver, una zona

de olivos y se midieron en campo las siguientes variables para 25 árboles de estudio: diámetro

de tronco, diámetro de copa, altura total, altura de copa, peso de poda y producción.

El objetivo de este trabajo consiste en calcular modelos de regresión a partir de datos LIDAR,

que permitan estimar variables dendrométricas. Para ello, las variables tomadas en campo

fueron correlacionadas a partir de diferentes estadísticos de la distribución de puntos LIDAR

incluidos en cada árbol.

La metodología seguida en el procesado de los datos LIDAR consistió a) cálculo del modelo

digital del terreno (MDT) y su validación; b) extracción de árboles obteniendo la superficie de

cada copa; c) cálculo de los estadísticos de la distribución de las alturas de los puntos incluidos

en cada árbol; d) obtención de modelos estadísticos entre las variables medidas en campo y la

información anterior y análisis de resultados obtenidos.

Entre los modelos de regresión calculados destacar que la variable peso poda tuvo un buen

ajuste a partir de los estadísticos LiDAR (área copa automática, int máximum, kurtraster)

obteniendo un coeficiente de determinación R² de 0.91 y un RMSE de 2.35 kg. Del mismo modo

se obtuvieron buenos resultados en la variable diámetro de tronco, con un R² de 0,87 y RMSE

7,05 m.

Los resultados obtenidos para la variable peso de poda, podrían permitir estimar el potencial

de biomasa en plantaciones frutales, lo que puede ser utilizado como un elemento clave en la

cadena de gestión de biocombustibles.

3

Índice

1. Introducción .............................................................................................................................................. 7

2. Material y datos ........................................................................................................................................ 8

2.1 Zona de estudio ................................................................................................................................. 8

2.2 Datos ................................................................................................................................................. 9

2.2.1. Datos tomados en campo ...................................................................................................... 9

2.2.2. Clasificación de árboles de estudio......................................................................................... 9

2.3. Software .......................................................................................................................................... 12

3. Metodología ............................................................................................................................................ 13

3.1. Preprocesado de los datos LIDAR ..................................................................................................... 15

3.2. Clasificación puntos LIDAR ............................................................................................................... 17

3.2.1. Clasificación manual de puntos LIDAR .................................................................................. 17

3.2.2. Clasificación automática de puntos LIDAR ............................................................................ 19

3.3. Obtención del MDT, MDS y MDSN .................................................................................................... 21

3.3.1. Obtención de los modelos mediante FUSION ....................................................................... 23

3.3.2. Obtención de los modelos mediante lidex ............................................................................ 27

3.3.3. Obtención de los modelos mediante ArcGIS ......................................................................... 29

3.3.4. Obtención de los modelos mediante Lastools ...................................................................... 30

3.3.5. Comparación de métodos .................................................................................................... 31

3.4. Obtención de variables estadísticas .................................................................................................. 31

3.4.1. Cálculo de la variable área ................................................................................................... 31

3.4.2. Cálculo de las variables estadísticas relacionadas con las elevaciones e intensidades

y número de puntos LIDAR por copa de árbol ...................................................................... 33

3.4.3. Cálculo de la variable densidad ............................................................................................ 35

3.4.4. Cálculo de variable volumen obtenido mediante voxelización .............................................. 35

3.4.5. Variables tipo raster y frecuencia de cabida cubierta ........................................................... 36

3.5. Interpretación de los resultados de los modelos de regresión .......................................................... 36

3.5.1 Hipótesis sobre el término de perturbaciones ...................................................................... 37

4

3.5.2 Hipótesis sobre los parámetros del modelo ......................................................................... 39

3.5.3 Modelo matemático ............................................................................................................ 39

3.5.4 Significación estadística ....................................................................................................... 40

3.5.5 Inferencia estadística ........................................................................................................... 41

4. Resultados ............................................................................................................................................... 42

4.1. Clasificación LIDAR automática ......................................................................................................... 42

4.1.1. Bondad de la clasificación automática mediante la clasificación manual .............................. 45

4.2. MDT ................................................................................................................................................. 47

4.2.1. MDT mediante FUSION ........................................................................................................ 47

4.2.2. MDT mediante lidex ............................................................................................................. 49

4.2.3. MDT mediante Lastools y mediante ArcGIS .......................................................................... 50

4.3. Comparativa variable área manual vs área automática ..................................................................... 51

4.4. Variables estadísticas ....................................................................................................................... 53

5. Conclusiones ........................................................................................................................................... 63

5.1. Cálculo del MDT ............................................................................................................................... 63

5.2. Clasificación de puntos LIDAR manual vs automática ........................................................................ 63

5.3. Cálculo área de copa manual vs automática ..................................................................................... 63

5.4. Análisis de modelos obtenidos ......................................................................................................... 64

6. Referencias bibliográficas ........................................................................................................................ 65

5

Índice de figuras

Figura 1. Imagen de situación y detalle de la zona de estudio ............................................................................ 8

Figura 2. Esquema general con la metodología empleada ............................................................................... 14

Figura 3. Planta datos LIDAR con MDS…………………. ......................................................................................... 16

Figura 4. Vista cenital de los datos LIDAR…………….. ......................................................................................... 16

Figura 5. Perfil datos LIDAR con MDS .............................................................................................................. 16

Figura 6. Selección del perfil en 3D de los datos LIDAR con MDSN ................................................................... 17

Figura 7. Detalle de puntos LIDAR en uno de los árboles de estudio. Siendo el marrón puntos del suelo y

azul puntos no clasificados. ............................................................................................................................. 17

Figura 8. Vista del perfil seleccionado en figura 3.8 con árboles de estudio ..................................................... 18

Figura 9. Detalle árboles de estudio con puntos LAS reclasificados.

Código 5 en verde y código 2 en marrón. ........................................................................................................ 19

Figura 10. Detalle vista 3D de los puntos LIDAR reclasificados en suelo y copas árboles de estudio ................. 19

Figura 11. Extracto del fichero de cabecera de información datos LIDAR

tras la clasificación automática de puntos terreno .......................................................................................... 20

Figura 12. Extracto del fichero de cabecera de información datos LIDAR tras la normalización de alturas y

clasificación automática de puntos terreno y vegetación alta. Las unidades están en centímetros. ................. 20

Figura 13. Clasificación automática puntos LIDAR mediante Lastools. Los puntos de color marrón indican

puntos terreno, los de color verde puntos vegetación y los de color negro puntos no clasificados .................. 21

Figura 14. Esquema con la metodología de la obtención del MDT ................................................................... 22

Figura 15. Función de pesos. (McGaughey, 2010) ............................................................................................ 23

Figura 16. Metodología LIDEX para la obtención del MDT. (Adaptado de Estornell et al. 2011). ....................... 27

Figura 17. Metodología Lastools para la obtención del MDT. (Adaptado de Isenburg et al. 2006) .................... 30

Figura 18. MDSN y puntos LIDAR de arboles de estudio con zona de detalle previa a la clasificación.

Puntos LIDAR color verde son puntos de vegetación y puntos marrones, puntos terreno. ............................... 32

Figura 19. Esquema de la obtención de variables estadísticas con FUSION. ..................................................... 33

Figura 20. Modelo digital de superficies normalizado con shp de los árboles de estudio.................................. 34

Figura 21. Ecuación del coeficiente de determinación ajustado ....................................................................... 42

Figura 22. TIN de elevaciones de la clasificación automática de los puntos LIDAR mediante Lastools............... 43

Figura 23. Superposición del TIN de elevaciones de la figura 22 y los puntos LIDAR sin clasificar.

Y alzado de los puntos LIDAR sin clasificar ....................................................................................................... 44

Figura 24. Detalle zona de taludes junto a la carretera, ubicada a la parte derecha de la zona de estudio. ...... 45

Figura 25. Gráficas variación rmse en función de los parámetros de FUSION ................................................... 48

Figura 26. Resumen rmse de los mejores modelos con lidex ........................................................................... 50

Figura 27. MDS Lastools zona de estudio con mapa de sombras ...................................................................... 51

Figura 28. MDT Lastools zona de estudio......................................................................................................... 51

Figura 29. Superposición área automática y manual con puntos LIDAR clasificados.

Puntos LIDAR color marrón representan puntos terreno y color verde, vegetación. ........................................ 53

Figura 30. Histograma residuos estandarizados diámetro tronco .................................................................... 55

Figura 31. Gráfico de dispersión entre los residuos estandarizados y la variable predicha

para diámetro de tronco. ................................................................................................................................ 56

Figura 32. Gráfico de dispersión entre los residuos estandarizados y el número de fila

para diámetro de tronco. ................................................................................................................................ 57

6

Figura 33. Gráfica de valores observados frente a predichos para la variable Peso poda2013 Grupo 1 ............ 62

Índice de tablas

Tabla 1. Características de los datos LIDAR. ....................................................................................................... 9

Tabla 2. Estadillo de las variables tomadas en campo ...................................................................................... 10

Tabla 3. Estadísticos de las variables tomadas en campo para todos los árboles de estudio ............................ 11

Tabla 4. Estadísticos de las variables tomadas en campo con menor diámetro de fuste .................................. 11

Tabla 5. Estadísticos de las variables tomadas en campo con mayor diámetro de fuste ................................... 12

Tabla 6. Extracto del fichero de cabecera de información de los datos LIDAR .................................................. 15

Tabla 7. Clases LIDAR definida por (ASPRS), 2013 ............................................................................................ 18

Tabla 8. Variación de rangos de 0,5 para parámetro w .................................................................................... 25

Tabla 9. Variación de rangos de 0,5 para parámetro g ..................................................................................... 26

Tabla 10. Variación de rangos de 2 para parámetro iteraciones....................................................................... 26

Tabla 11. Variación de rangos específicos de 0,1 para parámetro w ................................................................ 26

Tabla 12. Variación de rangos específicos para parámetro g ............................................................................ 26

Tabla 13. Variación de rangos específicos para parámetro iteraciones ............................................................ 26

Tabla 14. Tamaño de ventanas aplicados para la búsqueda de puntos de cota mínima con LIDEX ................... 28

Tabla 15. Tolerancias aplicadas para la búsqueda de puntos de cota mínima con LIDEX .................................. 28

Tabla 16. Ampliación de valores de tamaño de ventanas aplicados para la búsqueda de puntos de cota

mínima con LIDEX ........................................................................................................................................... 29

Tabla 17. Tabla ANOVA ................................................................................................................................... 41

Tabla 18. Porcentaje de puntos LIDAR por árbol bien clasificados automáticamente mediante Lastools. ......... 46

Tabla 19. Porcentaje de puntos LIDAR total bien clasificados automáticamente mediante Lastools ................. 47

Tabla 20. Parámetros y rmse del mejor MDT ................................................................................................... 49

Tabla 21. Resumen con los mejores resultados rmse para cada grupo de ventanas ......................................... 49

Tabla 22. Rmse de los modelos generados en ArcGIS y Lastools ...................................................................... 50

Tabla 23. Comparativa área manual vs automática y estadísticos columna diferencia. Unidades en metros .... 52

Tabla 24. Matriz de correlación para variable diámetro de tronco en Grupo 1 ................................................. 53

Tabla 25. Resumen regresión múltiple para variable dependiente diámetro de tronco en Grupo 1 ................. 54

Tabla 26. Tabla ANOVA para variable diámetro de tronco en Grupo 1 ............................................................. 54

Tabla 27. Resumen modelo para producción ................................................................................................... 58

Tabla 28. Resumen modelo para peso poda .................................................................................................... 59

Tabla 29. Resumen modelo para altura total ................................................................................................... 60

Tabla 30. Resumen modelo para altura copa ................................................................................................... 61

Tabla 31. Resumen R2 ajustados de las variables dependientes de estudio .................................................... 62

7

1. Introducción

La tecnología LIDAR (Light Detection And Ranging) se basa en sistemas activos de

teledetección consistentes en representar la superficie terrestre en 3D mediante nubes de

puntos, obtenidos a partir de pulsos laser. En función del tiempo que tarda en retornar el pulso,

se puede saber la posición exacta de la superficie sobre la que impactan. El sistema mide la

distancia del sensor al terreno a partir del tiempo que tarda el rayo de luz en alcanzar el suelo y

regresar al sensor.

Entre los sistemas LIDAR se encuentra el LIDAR aerotransportado, consta de un

distanciómetro láser y de un espejo que desvía el haz perpendicularmente a la trayectoria del

avión. Este desplazamiento lateral combinado con la trayectoria del avión permite realizar un

barrido del terreno. Si conocemos las coordenadas, ángulos altitud del avión y el ángulo del

espejo, podemos calcular a partir de estos datos y de las distancias medidas las coordenadas de

los puntos. Para ello, el sistema LIDAR va equipado de GPS diferencial (DGPS) y un sistema

inercial de navegación (INS) (Ruiz, 2003).

Estudios precedentes (Holmgren, 2003) han demostrado el potencial de esta tecnología en el

ámbito forestal estimando variables forestales a nivel de árbol y a nivel de parcela: altura,

biomasa, volumen del árbol, área basimétrica o número de pies entre otros. (Zimble et al., 2003;

Andersen et al., 2005; Hall et al., 2005; Li et al., 2008; Estornell et al., 2011; Holmgren, 2003;

Lefsky et al., 1999; Naesset y Bjerknes, 2001; Riaño et al., 2004; Sociedad Española de Ciencias

Forestales y Teledetección, 2005). La información de estas variables es similar a la obtenida en

inventarios tradicionales. Por lo que este método supone un ahorro de tiempo y costes en la

toma de datos en áreas extensas, así como la posibilidad de medición de áreas inaccesibles

mediante métodos terrestres (Hall et al., 2005). Sin embargo pocos estudios se han centrado en

la aplicación del sistema LIDAR aéreo en agricultura pudiendo resultar una herramienta eficiente

en la gestión de las plantaciones frutales tanto en la caracterización morfométrica de los árboles

como en la estimación de biomasa residual (Estornell et al., 2013).

Para estimar las variables características de los árboles a partir de datos LIDAR es habitual

aplicar regresiones múltiples entre los datos medidos en campo (por parcelas o árboles) y los

estadísticos de la distribución de las alturas de los datos LiDAR en su interior (Holmgren, 2003).

Los resultados anteriores reflejan las ventajas que aportan el uso de esta tecnología en el

ámbito forestal siendo necesario analizar si los algoritmos y métodos utilizados en este ámbito

proporcionan resultados óptimos en el ámbito de la agricultura.

Dentro del estudio de los estadísticos mediante LIDAR se puede trabajar con LIDAR por

parcelas (variables dasométricas) o por individuos (variables dendrométricas). Para el caso de las

parcelas, se trabaja a partir de la selección de diferentes zonas de estudio, en las cuales se

establece un radio medio dentro del cual se extraen los estadísticos a partir del conjunto de

árboles que queden delimitados en ella.

8

En nuestro caso como ya se ha explicado se trabajará a nivel de árbol, donde se extrae una

nube de puntos para cada uno de ellos. Mediante las que se obtendrán los diferentes

estadísticos.

El objetivo de este trabajo consiste en calcular modelos de regresión a partir de datos LIDAR,

que permitan estimar variables dendrométricas. Para ello es necesario: a) obtener el MDT que

mejor se adapte a la zona de estudio. b) extraer el área de copa de cada uno de los árboles de

estudio. c) analizar los resultados de los modelos de regresión múltiple obtenidos y comprobar la

importancia de los datos LIDAR en la predicción de las variables dependientes.

2. Material y datos

2.1 Zona de estudio

La zona de estudio pertenece a la provincia de Castellón y se encuentra ubicada al suroeste

de la misma, en el municipio de Viver. La localización queda delimitada por las coordenadas

UTM de las esquinas superior izquierda (704103.36 m E, 4420764.34 m N) e inferior derecha

(704393.23 m E, 4420337.45 m N) en el sistema de referencia ETRS 89 huso 30 (Figura 1).

Figura 1. Imagen de situación y detalle de la zona de estudio

9

Se trata de una zona de cultivo de olivos, cuyas alturas van desde los 3 a los 4,5 m de altura y

los valores promedio para las siguientes variables son los que siguen: diámetro de tronco 29.78

cm, producción 28.89 kg, peso poda 9.26 kg y altura de copa 2.91 m. La morfología de las

parcelas agrícolas está caracterizada por la presencia de terrazas lo que se deberá tener en

cuenta al evaluar los modelos digitales del terreno que se van a calcular.

2.2 Datos

Los datos empleados para llevar a cabo este proyecto han sido datos LIDAR registrados en

julio de 2012, datos tomados en campo para el estudio estadístico, puntos GPS-rtk.

2.2.1. Datos tomados en campo

Respecto a los datos LIDAR, el vuelo fue realizado el 15 de septiembre de 2012. Las

principales características de estos datos se encuentran en la Tabla 1.

Especificación Áreas con densidad de 4 puntos/m²

Campo de visión transversal máximo (FOV) 18º

Frecuencia de escaneo 87,9 Hz

Frecuencia de pulso 150,6 Hz

Resolución espacial Densidad promedio 4ptos/m²

Tabla 1. Características de los datos LIDAR.

Los datos de campo de la zona de estudio fueron tomados el mismo año del vuelo LIDAR y las

variables medidas fueron las siguientes (altura total, diámetros de copa, peso de poda,

producción en kg, diámetro de tronco, alturas base por copa).

Para la evaluación de los diferentes MDT calculados se emplearon 62 puntos medidos en GPS

(Leica System 1200, Heerbrugg, Switzerland) basados en VRS (Virtual Reference Station) Internet

RTK.

2.2.2. Clasificación de árboles de estudio

Para este estudio se cuenta con 25 árboles, que presentan los siguientes datos tomados en

campo. (Tabla 2)

10

Identificador Diámetro

tronco(cm)

Altura

total

Diámetro

copa 1

Diámetro

copa 2

Producción

(kg)

Peso

poda (kg)

Altura de copa

(m)

1 15.51 3.28 2.98 3.45 17.05 3.56 2.54

2 17.2 3.49 3.22 4.09 22.1 4.66 2.64

3 18.2 3.81 3.96 4.28 30.9 3.64 3.12

4 14.4 3.48 3.46 3.59 18.95 2.39 2.8

5 15.3 3.82 3.76 3.92 26.6 3.92 3.1

6 18 3.48 3.92 4.05 28.9 7.54 2.86

7 17.1 3.91 3.32 3.39 31.65 2.39 3.03

8 16.3 3.84 2.86 3.87 29.25 5.34 3

88 16.6 4.23 3.47 3.79 32.5 3.94 3.38

17 54 3.77 3.74 4.39 21.32 7.95 2.96

18 80.7 4.08 5.06 6.09 24.75 21.2 3.1

19 52 3.64 4.73 5.34 27.3 26.64 2.58

20 46 3.87 5.13 5.16 38.9 19.86 2.96

21 66 4.28 5.89 6.86 35.55 31.4 3.21

22 47 3.96 4.44 4.59 31.25 18.5 2.93

23 52 4.48 4.89 5.38 54.65 19.4 3.5

24 53 3.49 4.63 4.9 22.75 12.4 2.89

25 15 3.72 2.89 3.23 45.2 2.06 2.99

26 17 3.56 3.47 3.76 30.1 4.84 2.75

27 20 3.7 3.78 4.13 34.8 8.4 2.93

28 18.6 3.7 3.32 3.6 36.3 5.28 2.89

29 14.5 3.12 2.92 3.16 29 3.91 2.38

30 17.6 3.87 3.21 3.42 13.25 3.18 3.19

31 25 3.27 3.34 3.35 19.05 2.57 2.47

32 17.4 3.45 2.71 3.41 20.1 6.42 2.64

Tabla 2. Estadillo de las variables tomadas en campo

11

Diámetro

tronco (cm)

Altura

total

Diámetro

copa 1

Diámetro

copa 2

Producción

(kg)

Peso

poda (kg)

Altura

de copa

Promedio 29.78 3.73 3.80 4.21 28.89 9.26 2.91

Desviación típica 19.72 0.33 0.84 0.95 9.17 8.44 0.27

Máximo 80.7 4.48 5.89 6.86 54.65 31.4 3.5

Mínimo 14.4 3.12 2.71 3.16 13.25 2.06 2.38

Tabla 3. Estadísticos de las variables tomadas en campo para todos los árboles de estudio

Los árboles se han clasificado en dos grupos atendiendo al parámetro diámetro de tronco.

Como se puede observar los valores de media, desviación, valor máximo y mínimo de los árboles

de cada grupo son diferentes y no habiendo superposición de clases. Además este parámetro no

está influenciando por las labores de los agricultores como por ejemplo la poda

En el primer grupo se han incluido aquellos árboles de estudio cuyo diámetro de fuste es

inferior a su respectivo promedio, que llamaremos árboles de menor fuste. Este grupo está

formado por 17 de los 25 árboles de estudio, siendo estos los citados a continuación; 1, 2, 3, 4,

5, 6, 7, 8, 88, 25, 26, 27, 28, 29, 30, 31 y 32.

Grupo 1 Diámetro

tronco (cm)

Altura

total

Diámetro

copa 1

Diámetro

copa 2

Producción

(kg)

Peso poda

(kg)

Altura

de copa

Promedio 17.39 3.65 3.35 3.69 28.04 4.41 2.89


Máximo 25 4.23 3.96 4.28 45.2 8.4 3.38

Mínimo 14.4 3.12 2.71 3.16 13.25 2.06 2.38

Tabla 4. Estadísticos de las variables tomadas en campo con menor diámetro de fuste

Y en el grupo dos estarán aquellos árboles de estudio que tengan un valor superior a su

respectivo promedio, que llamaremos árboles de mayor fuste. Este grupo está formado por 8 de

los 25 árboles de estudio, siendo estos los siguientes; 17, 18, 19, 20, 21, 22, 23 y 24.

12

Grupo 2 Diámetro

tronco (cm)

Altura

total

Diámetro

copa 1

Diámetro

copa 2

Producción

(kg)

Peso poda

(kg)

Altura de

copa

Promedio 56.34 3.95 4.81 5.34 32.06 19.67 3.02


Máximo 80.7 4.48 5.89 6.86 54.65 31.4 3.5

Mínimo 46 3.49 3.74 4.39 21.32 7.95 2.58

Tabla 5. Estadísticos de las variables tomadas en campo con mayor diámetro de fuste

Puesto que la muestra para el segundo grupo se considera pequeña, los análisis estadísticos

se realizarán considerando todos los árboles de estudio (25) y con el primer grupo (17). Estos

grupos se referenciarán como Grupo 1 y Grupo 2, respectivamente.

2.3. Software

Para llevar a cabo el siguiente proyecto, los programas informáticos empleados fueron los

siguientes:

ArcGIS Desktop:

Se empleó la versión de ArcGIS 9.3 para la consulta y edición de los datos LIDAR, así como

para obtener la precisión de los modelos digitales del terreno, y otros análisis espaciales como el

cálculo del área de las copas manual y automática.

FUSION/LDV:

Fue utilizado para la obtención del modelo digital del terreno y el modelo digital de

superficies (MDS), así como para preprocesar los datos LIDAR eliminando posibles puntos

anómalos. Y también en la obtención de estadísticos para cada árbol de estudio, mediante los

que se obtendrán los modelos de regresión múltiple.

Statistica:

Software estadístico empleado para calcular modelos de regresión múltiple entre las

variables medidas en campo y los estadísticos derivados de los datos LiDAR para la generación

de gráficas.

13

Envi + IDL:

Se ha utilizado un programa basado en procesos iterativos procesados en Interactive Data

Language (IDL, ITT Visual Systems). La extensión idl empleada “lidex” ha permitido generar

modelos digitales del terreno, de superficie y el normalizado. Así como la extracción de

características a partir del shp de las copas de los árboles de estudio.

Esta extensión es una aplicación experimental implementada por el grupo CGAT-UPV y

permitido obtener nuevas características a nivel raster, que no habían sido obtenidas mediante

el software FUSION.

Lastools:

Empleado para generación de variables estadísticas como la densidad de puntos por

intervalos. Así como para automatizar la clasificación de los puntos LIDAR.

Matlab:

Este software se empleó para obtener la variable volumen sólido a través de la

implementación de un código que recibía los archivos .txt de los puntos LIDAR como ficheros de

entrada y obtenía el volumen sólido de cada árbol de estudio.

3. Metodología

Inicialmente se comenzó por preprocesar los datos LIDAR y clasificarlos según el tipo.

Seguidamente se obtuvieron los MDT, MDS y MDSN y por último se obtuvieron las variables

estadísticas y se interpretaron los resultados obtenidos de los modelos de regresión.

El esquema de la figura 2 muestra la metodología general empleada en el proyecto:

Figura 2. Esquema general con la metodología empleada

Datos iniciales

Datos de campo regresión multiple y análisis de resultados

Datos LidarPreprocesado y

clasificación de datos Lidar

Generación de MDTs

Elección del dtm con mejor rmse y

generación de su dsm y ndsm

Digitalización del shp areas de copa

Normalización de puntos lidar y recorte de los mismos por el

área de las copas

3.1. Preprocesado de los datos LIDAR

Este paso es fundamental siempre que se trabaje con datos LIDAR, ya que se pueden generar

outliers durante el registro de los datos.

Para llevar a cabo esta tarea, se ha empleado el programa fusión con la extensión LDV y el

comando lasinfo de Lastools para comprobar si realmente existían puntos outliers.

Inicialmente se partía del archivo las original 000088_ort.las, el cual se recortó por la zona de

estudio mediante las coordenadas x, y de las esquinas: 704209 m; 4420359 m; 704365 m;

4420724.245 m. Se obtuvo un fichero de información de los datos LIDAR sin preprocesar, mediante

el comando lasinfo de Lastools (Tabla 6):

min x y z: 704.209 4.420.359 550,21

max x y z: 704.365 4.420.724 845,73

Reporting mínimum and máximum for all LAS point record entries …

X 0 156

Y 0 365

Z 550,21 845,72

Tabla 6. Extracto del fichero de cabecera de información de los datos LIDAR

Como se puede apreciar el intervalo de cotas oscila entre los 550 y 845 metros, lo cual ya induce

a pensar que existen puntos anómalos, ya que las alturas elipsoidales tomadas con el GPS oscilan

entre 555.176 y 572.57 y las alturas de los árboles de estudio van de lo 2 a los 5 metros

aproximadamente.

Para comprobar si existen puntos anómalos gráficamente, se ha empleado el software FUSION

con la extensión LDV.

Debido a la diferencia de altura respecto de la mayoría de puntos LIDAR con los outliers,

probablemente sean puntos generados mediante el retorno del pulso laser al interceptar algún ave

en vuelo.

Para corregir estos puntos outliers se ha realizado un filtrado de cota entre las cotas 550 y 600 y

se ha comprobado tanto numérica como gráficamente que ya no existen puntos outliers.

16

Tras llevar a cabo el preprocesado, como ya se aprecia en el fichero de información anterior, se

comprueba que dichos puntos outliers ya no están en los datos LIDAR preprocesados. Obteniendo

las siguientes figuras en el visualizador LDV en planta y perfil que lo corroboran, tal y como se refleja

en la leyenda del histograma. Por lo que ya no existirán esos valores de cota de los outliers que

podrían distorsionar los valores de cota del modelo al realizar la interpolación.

Figura 3. Planta datos LIDAR con MDS Figura 4. Vista cenital de los datos LIDAR

Figura 5. Perfil datos LIDAR con MDS

17

3.2. Clasificación puntos LIDAR

En este apartado se han analizado diferentes herramientas que se explican con mayor detalle en

los siguientes subapartados (3.2.1, 32.2.2 y 3.2.3). A partir de las cuales se ha llevado a cabo la

clasificación de los puntos LIDAR, siendo los puntos clasificados como vegetación los que se

utilizarán para la estimación de las diferentes variables dendrométricas de este trabajo.

3.2.1. Clasificación manual de puntos LIDAR

Este apartado se realizó para posteriormente poder obtener el área de copa de los árboles, lo

cual se explica en apartados posteriores. Y para ello, se realizó la clasificación de aquellos puntos que

estaban como “no clasificados”, empleando la herramienta llamada “las dataset profile view” de la

barra de herramientas “las dataset” de ArcGIS© versión 9.3. Mediante esta herramienta se

seleccionan los árboles de estudio y su contorno como se muestra en la figura 6. A partir de la cual

se obtiene la figura 7, en la que se ha cambiado el color de los puntos no clasificados de rojo a azul

para diferenciarlos mejor respecto a los puntos del suelo (color marrón).

Figura 6. Selección del perfil en 3D de los datos LIDAR con MDSN

Figura 7. Detalle de puntos LIDAR en uno de los árboles de estudio. Siendo el marrón puntos del suelo y azul puntos no clasificados.

18

Figura 8. Vista del perfil seleccionado en figura 3.8 con árboles de estudio

Una vez se ha seleccionado la ventana donde aparecen los árboles de estudio se procede a

asignarles su codificación. Para realizar esta función se han seguido los criterios definidos por la

ASPRS. La cual definió estos códigos para formatos de LAS 1.1 y 1.2 ó 1.3 y como nuestro archivo es

LAS 1.2, se podrán aplicar los siguientes códigos de clasificación para puntos LIDAR:

Código de clase Tipo de clasificación

0 Creado, no clasificado nunca

1 Sin clasificar

2 Suelo

3 Vegetación baja

4 Vegetación media

5 Vegetación alta

6 Edificio

7 Puntos bajos (ruido)

8 Clave de modelo

9 Agua

Tabla 7. Clases LIDAR definida por (ASPRS), 2013

Para ello se seleccionaron los puntos pertenecientes a las copas de los árboles de estudio y se les

asignó el código 5 “vegetación alta”. Se hizo lo propio para los puntos del suelo, asignándole el

código 2 y el resultado final se muestra en la figura 9:

http://www.asprs.org/

19

Figura 9. Detalle árboles de estudio con puntos LAS reclasificados. Código 5 en verde y código 2 en marrón.

Por último, para comprobar que la clasificación se hizo correctamente, se cambió el color de

fondo a negro para visualizar dicha zona en planta y con una cierta inclinación:

Figura 10. Detalle vista 3D de los puntos LIDAR reclasificados en suelo y copas árboles de estudio

El proceso explicado anteriormente se hizo para todos los árboles de estudio, ya que la muestra

explicada con figuras correspondía solo a la zona superior de estudio.

3.2.2. Clasificación automática de puntos LIDAR

Con el objetivo de automatizar la clasificación de los puntos LIDAR se empleó el software de

Lastools.

Este programa permite realizar la clasificación de los puntos LIDAR diferenciando entre árboles,

edificios, suelo y puntos sin clasificar. Para conseguir este proceso se han seguido los siguientes

pasos:

20

Primero se calculó un modelo digital del terreno que diferenció entre puntos del suelo y puntos

no clasificados mediante el comando “lasground”.

Figura 11. Extracto del fichero de cabecera de información datos LIDAR tras la clasificación automática de puntos terreno

Una vez ya se tienen diferenciados los puntos LIDAR pertenecientes al suelo de los que no, se

empleó el comando “lasheight” para normalizar las alturas y asignar a todos los puntos LIDAR que

estén entre 0.5 y 5 metros, la clase de vegetación alta. Se eligió este rango ya que ninguno de

nuestros árboles de estudio supera dicha altura.

Figura 12. Extracto del fichero de cabecera de información datos LIDAR tras la normalización de alturas y clasificación automática de puntos terreno y vegetación alta. Las unidades están en centímetros.

21

Por último se utilizó el comando “lasclassify” para reducir el número de puntos no clasificados,

diferenciando estos entre vegetación alta y edificios. Pero al no haber edificios en la zona de estudio,

no se redujo el número de puntos LIDAR no clasificados.

Figura 13. Clasificación automática puntos LIDAR mediante Lastools. Los puntos de color marrón indican puntos terreno, los de color verde puntos vegetación y los de color negro puntos no clasificados

3.3. Obtención del MDT, MDS y MDSN

Tras llevar a cabo el preprocesado de los datos LIDAR de la zona de estudio, en este apartado se

describirá cómo se calculó y evaluó la precisión de los diferentes modelos: MDT, MDS y el MDSN.

22

La metodología seguida para llevar a cabo el cálculo del mejor modelo digital del terreno ha sido

la siguiente:

Figura 14. Esquema con la metodología de la obtención del MDT

Tal y como se observa del esquema, se han obtenido los modelos digitales del terreno mediante

programas diferentes (LiDEX, FUSION, ArcGIS y Lastools) con la finalidad de conseguir el mejor MDT

que se adapte a las características de la zona de estudio. Tomando como premisa para todos los

softwares que el tamaño de la celda sea de 1 metro y el resto de parámetros se irán modificando en

función de lo que permita cada software. A continuación se procede a explicar el algoritmo de

cálculo utilizado por cada programa para el cálculo del MDT.

Datos lidar brutos

Preprocesado y clasificación de datos lidar zona de estudio

Batería de mdt con lidex

Batería de mdt con fusion

Obtención de rmse entre el mdt y puntos

rtk para cada mdt

Elección del mejor mdt

Obtención del mds y chm

Batería de mdt con arcgis

Batería de mdt con lastools

23

3.3.1. Obtención de los modelos mediante FUSION

Tras llevar a cabo el preprocesado de los puntos LIDAR, se emplea el nuevo archivo “.las“

generado para obtener el MDT mediante el software FUSION. Para ello se explica a continuación los

comandos empleados y que realiza cada uno de ellos:

-GroundFilter:

Se encarga de filtrar la nube de puntos LIDAR e identificar los retornos cuya probabilidad de

pertenecer a la superficie del terreno sea muy elevada. Pese a esta elevada probabilidad de filtrado,

este comando no es perfecto y deja puntos que no pertenecen al terreno sin filtrar. Como puedan

ser, zonas extensas en elevación que tienen valores de cota similares, como por ejemplo tejados de

edificios. Respecto a los retornos en vegetación cabe señalar que, la mayoría pueden ser eliminados

mediante una buena elección de coeficientes en la función de pesos y el número de iteraciones. En

general este comando produce el conjunto de puntos necesarios para obtener las alturas de

vegetación. Se debe tener presente que tanto el conjunto de puntos resultado como los resultados

de los modelos digitales pueden resultar inadecuados cuando se pretende obtener un modelo de

superficies de los puntos terreno, es por esto por lo que para obtener el modelo digital del terreno

se emplea el comando GridSurfaceCreate.

El algoritmo de filtrado (Kraus y Pfeifer, 1998) está basado en una predicción linear (Kraus y

Mikhail, 1972) con una precisión individual para cada medida. Está implementado como un proceso

iterativo en el que en un primer paso, una superficie es procesada con los mismos pesos para todos

los puntos LIDAR. Este resultado consiste en una superficie situada entre los puntos terreno y el

modelo de superficie. Donde los puntos del terreno son más propensos a estar por debajo de dicha

superficie y los puntos de vegetación por encima de la misma. La distancia y la dirección a la

superficie se utilizan para calcular los pesos para cada punto LIDAR a partir de la siguiente función de

pesos.

Figura 15. Función de pesos. (McGaughey, 2010).

Los parámetros a y b determinan los pasos o distancias de la función de pesos. Para la mayoría de

aplicaciones los valores de 1 y 4 para a y b respectivamente, han sido probados obteniendo buenos

resultados. El valor de g determina a qué puntos es asignado el valor del peso máximo, siendo estos,

aquellos puntos por debajo de la superficie generada y superiores al valor g y asignándoles un peso

igual a 1.

24

El parámetro de distancia de puntos situados sobre el terreno, w, se utiliza para establecer un

valor máximo de altura para los puntos que tengan una influencia sobre la superficie intermedia. Los

puntos sobre el nivel definido por (g + w) se les asigna un peso 0.

En el comando implementado, los valores de g y w no varían una vez se ha lanzado el comando.

Hay referencias en las que el parámetro g varía para cada iteración (Kraus y Pfeifer 1998)(13). Tras la

iteración final, los puntos terreno verdaderos son seleccionados utilizando la superficie intermedia.

Todos los puntos con elevaciones que satisfagan las dos primeras condiciones de la función de pesos

son considerados puntos terreno.

Los parámetros empleados para llevar a cabo este comando han sido los siguientes. Siendo estos

los que mejor se adaptan a la zona de estudio:

c:\fusion\GroundFilter /gparam:-1 /wparam:0.9 /aparam:5

C:\LIDAR_Fusion\groundfilter_dtmpts.lda 1 C:\LIDAR_Fusion\classify.lda

‐ GridSurfaceCreate:

El comando GridSurfaceCreate crea un modelo de superficie en formato imagen utilizando nubes

de puntos aleatorias. Dicho modelo de superficie se almacena en formato .dtm, cuyos valores de

elevación están en punto flotante. Por defecto, GridSurfaceCreate obtiene la elevación de cada celda

del grid a partir del valor medio de las elevaciones obtenidas a partir de los puntos que se

encuentran dentro de la propia celda. Este método parece que funciona bien con los datos LIDAR ya

filtrados, donde ya se tienen clasificados los puntos terreno.



c:\fusion\GridSurfaceCreate C:\LIDAR_Fusion\dtm1metro_fusion.dtm 1 M M 1 30 0 0

C:\LIDAR_Fusion\groundfilter_dtmpts-lda

La superficie de salida tendrá un tamaño de pixel de un metro, donde tanto las unidades

planimétricas (XY) como altimétricas Z están en metros (M). El sistema de coordenadas es el UTM en

huso 30 y debido a que FUSION no tiene programado el datum para el etrs89, se introduce datum

desconocido (0).

Gracias al uso de los comandos GorundFilter y en segundo lugar GridSurfaceCreate se consigue

obtener los modelos digitales del terreno. En cambio, si se quisiera obtener los modelos digitales de

superficie, bastaría con emplear el comando GridSurfaceCreate.

Una vez explicados estos comandos se procede a la explicación del comando CanopyModel, que

también permite obtener los modelos digitales de superficie, aunque la forma de obtenerlo difiere

de la explicada con el GridSurfaceCreate.

25

‐ CanopyModel:

CanopyModel crea el modelo de superficie a partir de las nubes de puntos LIDAR. Por defecto,

este comando asigna la elevación del retorno más alto de los puntos LIDAR que caen dentro de cada

celda al centro de la misma. Este comando también podrá crear un modelo digital de superficies

normalizado y un MDT, pero a diferencia del comando GridSurfaceCreate, lo obtiene a partir de los

mayores valores de elevación de la nube de puntos terreno LIDAR.



c:\fusion\CanopyModel C:\LIDAR_Fusion\dsm1metro_fusion.dtm 1 M M 1 30 0 0

C:\LIDAR_Fusion\classify.lda

La metodología seguida para obtener los MDT propios de FUSION ha sido la siguiente:

o Con los datos ya preprocesados, se comienza variando el parámetro w inicialmente

en rangos de 0,5 y dejando el resto de parámetros (g, a, b e iteraciones) constantes.

Esto se hace para acercarnos a nuestro mejor MDT.

o Una vez se ha obtenido el mejor MDT, se procede a variar el parámetro g con rangos

iniciales de 0,5 y dejando constantes el resto. Del mismo modo que en el paso

anterior y realizándose el mismo proceso para los parámetros a, b e iteraciones.

Siendo los rangos de 1 para a y b y de 2 para las iteraciones.

Cabe mencionar que tanto para los parámetros w como g, una vez se vieron entre cuales de sus

valores nos acercábamos al mejor MDT, se realizaron nuevos rangos hasta dar con el mejor MDT

posible. De la misma forma para el parámetro iteraciones, se realizó un nuevo rango de 1 con la

misma finalidad ya mencionada. Es por esto por lo que inicialmente se generaron 5 modelos para w

y g, 6 para a y b y 3 para iteraciones. Haciendo un total de 25 MDT. (Tablas 8, 9, 10, 11, 12 y 13).

g a b iteraciones

w

0,5 𝑔𝑐𝑡𝑒 𝑎𝑐𝑡𝑒 𝑏𝑐𝑡𝑒 𝑖𝑡𝑒𝑟𝑐𝑡𝑒

1 𝑔𝑐𝑡𝑒 𝑎𝑐𝑡𝑒 𝑏𝑐𝑡𝑒 𝑖𝑡𝑒𝑟𝑐𝑡𝑒


2 𝑔𝑐𝑡𝑒 𝑎𝑐𝑡𝑒 𝑏𝑐𝑡𝑒 𝑖𝑡𝑒𝑟𝑐𝑡𝑒


Tabla 8. Variación de rangos de 0,5 para parámetro w

26

w a b iteraciones

g

0,5 𝑤𝑐𝑡𝑒 𝑎𝑐𝑡𝑒 𝑏𝑐𝑡𝑒 𝑖𝑡𝑒𝑟𝑐𝑡𝑒

1 𝑤𝑐𝑡𝑒 𝑎𝑐𝑡𝑒 𝑏𝑐𝑡𝑒 𝑖𝑡𝑒𝑟𝑐𝑡𝑒


2 𝑤𝑐𝑡𝑒 𝑎𝑐𝑡𝑒 𝑏𝑐𝑡𝑒 𝑖𝑡𝑒𝑟𝑐𝑡𝑒


Tabla 9. Variación de rangos de 0,5 para parámetro g

w g a b

iteraciones

4 𝑤𝑐𝑡𝑒 𝑔𝑐𝑡𝑒 𝑎𝑐𝑡𝑒 𝑏𝑐𝑡𝑒



Tabla 10. Variación de rangos de 2 para parámetro iteraciones

Y tras observar en que rangos el rmse tenía una tendencia a aumentar la precisión del MDT, se

añadieron los siguientes valores, generándose 8 modelos más.

g a b iteraciones

w





Tabla 11. Variación de rangos específicos de 0,1 para parámetro w

w a b iteraciones

g

-0,75 𝑤𝑐𝑡𝑒 𝑎𝑐𝑡𝑒 𝑏𝑐𝑡𝑒 𝑖𝑡𝑒𝑟𝑐𝑡𝑒



Tabla 12. Variación de rangos específicos para parámetro g

w g a b

iteraciones 5 𝑤𝑐𝑡𝑒 𝑔𝑐𝑡𝑒 𝑎𝑐𝑡𝑒 𝑏𝑐𝑡𝑒

Tabla 13. Variación de rangos específicos para parámetro iteraciones

27

Finalmente se han generado 33 modelos, de los que cabe destacar que se han calculado los

RMSE para los modelos en los que ya se visualizaba el terreno sin puntos de copas de árboles.

También ha habido modelos donde al modificar los parámetros w y g se ha filtrado en exceso el

intervalo de cota del modelo obtenido, que no ha llegado a cubrir toda la zona de estudio.

Desestimando en este caso el cálculo del rmse también, debido a la falta de cota en puntos GPS de la

zona de estudio.

3.3.2. Obtención de los modelos mediante lidex

Los fundamentos de este algoritmo se basan en un proceso iterativo de búsqueda de puntos de

elevación mínima en ventanas decreciente (Estornell et al., 2011). Para entender mejor la

metodología que emplea el algoritmo nos apoyaremos de la figura 16, que muestra el proceso

mediante imágenes.

Figura 16. Metodología LIDEX para la obtención del MDT. (Adaptado de Estornell et al. 2011).

28

Como se aprecia en la figura 16, inicialmente LIDEX emplea ventanas para dividir el área de

estudio. A partir del primer tamaño de ventana (v1), obtiene los puntos de altura mínima de cada

ventana (v1) que va recorriendo el área de estudio. Obteniendo así el primer modelo digital del

terreno mediante el método de interpolación de la triangulación de Delaunay.

Tras la obtención del primer modelo digital del terreno, se escoge el tamaño de una segunda

ventana (v2) inferior a v1 y se realiza una nueva búsqueda de puntos mínimos a partir de los datos

originales. Tras la selección de todos los puntos mínimos empleando v2, se escogerán solo aquellos

que su diferencia con respecto al modelo digital del terreno anterior tengan un valor inferior al

umbral seleccionado (u1) y descartando aquellos que superen dicho umbral (u1). Finalmente se

obtiene el segundo modelo digital del terreno con aquellos puntos cuya diferencia era inferior a u1.

Por último, se vuelve a realizar el proceso anterior con un tamaño de ventana inferior a v2,

siendo v3. Donde se eliminarán aquellos puntos cuya diferencia con el segundo modelo digital del

terreno sean superiores a un segundo umbral (u2) establecido. Obteniendo el modelo digital del

terreno final con los puntos cuya diferencia sea inferior a u2.

Las tablas 3.8 y 3.9 contienen los valores de los diferentes parámetros analizados para el

cálculo de los MDTs. Se han probado 4 secuencias de ventanas (v1, v2 y v3). Para cada grupo de

ventanas se han analizado todas las posibles combinaciones de los parámetros u1 y u2 (25 posibles

combinaciones). Por tanto en total se han analizado 100 DTM

Tamaño de ventana (m)

v1 v2 v3

8 4 2 10 5 2,5 12 6 3 20 10 5

Tabla 14. Tamaño de ventanas aplicados para la búsqueda de puntos de cota mínima con LIDEX

Umbral 1 (m)

0,5 1 1,5 2 2,5

Umbral 2

(m)

0,5 𝑡11 𝑡12 𝑡13 𝑡14 𝑡15

1 𝑡21 𝑡22 𝑡23 𝑡24 𝑡25

1,5 𝑡31 𝑡32 𝑡33 𝑡34 𝑡35

2 𝑡41 𝑡42 𝑡43 𝑡44 𝑡45

2,5 𝑡51 𝑡52 𝑡53 𝑡54 𝑡55

Tabla 15. Tolerancias aplicadas para la búsqueda de puntos de cota mínima con LIDEX

29

Tras realizar el cálculo de estos 100 MDT, y teniendo en cuenta que la mayor variación del RMSE

variando las tolerancias, para el mismo tamaño de ventana fue igual a 0,0171 m. Se calcularon

nuevos MDT con los tamaños de ventana de la tabla 3.10.

Tamaño de ventana (m)

v1 v2 v3

4 2 1 6 3 1,5

10 5 2,5 12 6 3 20 4 2

Tabla 16. Ampliación de valores de tamaño de ventanas aplicados para la búsqueda de puntos de cota mínima con LIDEX

Respecto a la metodología seguida para la obtención de estos nuevos MDT. Se calculó un MDT

para cada grupo de ventanas de la tabla 16 con tolerancias de 0.5, puesto que se comprobó en los

cien MDT calculados de la tabla 14, que era con la que mejor RMSE se obtuvieron.

Y para aquellos MDT que obtuvieron un RMSE próximo al del mejor MDT obtenido con este

software hasta el momento, se varió el segundo umbral en rangos de 0.1 siendo estos siempre

inferiores de 0.5. Obteniéndose un total de 8 MDT, que sumados a los ya calculados hicieron un total

de 108 MDT.

Respecto al cálculo del MDS se seleccionaron aquellos puntos LIDAR cuyo valor de elevación fue

el mayor para cada celda. Estableciendo los mismo parámetros que para el MDT. Finalmente cabe

mencionar que el MDSN se obtuvo mediante la diferencia entre el MDS y el MDT.

3.3.3. Obtención de los modelos mediante ArcGIS

Para obtener el MDT y el MDS mediante este software la metodología a seguir consistió en pasar

el archivo “.las” preprocesado a vectorial, indicándole que se quede el primer o último retorno en

función de si lo que se quiere calcular es el modelo digital de superficies o el de terreno

respectivamente.

El tipo de interpolación que realiza este software es la interpolación por vecinos próximos, que

emplea ponderaciones basadas en áreas adyacentes (método de Voronoi).

Respecto al método de muestreo y la distancia utilizada para definir el tamaño de celda del

raster de salida, lo que controla la resolución horizontal del raster, se empleó 1 m. Se obtuvo 1

modelo, debido a los pocos parámetros que se permiten modificar.

30

3.3.4. Obtención de los modelos mediante Lastools

La metodología que lleva este software internamente para la interpolación y obtención del MDT

es la siguiente:

Figura 17. Metodología Lastools para la obtención del MDT. (Adaptado de Isenburg et al. 2006).

Primero la herramienta Spfinalize, se encarga de leer los puntos LIDAR del archivo *.las de

estudio ya preprocesado y genera un grid. A continuación la herramienta Spdelaunay2d, se encarga

de a partir del grid obtenido de spfinalize generar la triangulación mediante la técnica Delaunay. Y

por último la herramienta Tin2dem es la que genera un raster grid a partir de la triangulación de

Delaunay y para ello emplea una interpolación lineal de los valores de elevación.

Mediante este software se generó 1 modelo, ya que ocurre lo mismo que con el software de

ArcGIS.

31

3.3.5. Comparación de métodos

Tras explicar cómo se obtienen los MDTs mediante los diferentes programas citados, se procedió

a calcular el RMSE de cada uno de ellos con la finalidad de escoger el que mejor represente el

terreno y lo haga con un menor error.

Para ello se consideraron los puntos GPS medidos en la zona de estudio. Las elevaciones de

estos puntos se compararon respecto a las cotas de los MDTs y se calcularon los valores de RMSE.

RMSE= (𝑍𝑖−𝑧𝑖)

2𝑁𝑖=1

𝑁

Siendo N el número de puntos medidos, Zi el valor medido en campo con GPS y zi el valor

obtenido del MDE calculado.

3.4. Obtención de variables estadísticas

3.4.1. Cálculo de la variable área

Para obtener esta variable se hizo uso del software ArcGIS, siguiendo estas dos metodologías,

manual y automática. Las cuales se describen a continuación:

Para obtener el área de las copas manualmente, se digitalizaron aquellos árboles que se

pretenden estudiar, asignándole el identificador que le corresponde a cada olivo, según como se

explica en el apartado clasificación manual de puntos LIDAR.

A partir del archivo LIDAR y el MDSN obtenido del MDT de mejor RMSE de FUSION y su

correspondiente MDS, se procede a la digitalización del shp de las copas de los árboles de estudio a

partir del contorno que forman los puntos LIDAR clasificados como vegetación.

Para ayudar a la visualización y digitalización de ese shp se emplean las herramientas propias del

“las dataset”, que permiten filtrar los puntos LIDAR en función de su clase y del primer retorno. Con

ello se obtiene la figura 18, donde ya se aprecia que los puntos LIDAR clasificados como vegetación

caían dentro de las copas de los árboles del MDSN:

32

Figura 18. MDSN y puntos LIDAR de arboles de estudio con zona de detalle previa a la clasificación. Puntos LIDAR color verde son puntos de vegetación y puntos marrones, puntos terreno.

Para obtener el área de cada copa de árbol de forma automática, la metodología seguida ha sido

la siguiente:

1. Se han obtenido los puntos LIDAR normalizados mediante el comando lasheight de Lastools,

estableciéndose un umbral mínimo de 0.5 metros. Ya que aquellos puntos inferiores a dicha

altura serán vegetación del terreno y no forman parte de los puntos de copa de los árboles,

que es lo que se pretende filtrar.

2. Mediante la herramienta “las to multipoint” de ArcGIS se obtuvieron los puntos LIDAR de las

copas en un shp del “.las” obtenido en el anterior paso. Esto se consiguió especificando que

solamente tomase los puntos clasificados como vegetación (código ASPRS = 5), así se

consiguieron las copas delimitadas.

3. A partir del shp con multipuntos obtenido de las copas de los árboles, se empleó la

herramienta “multipart to singlepart” para pasar de entidades multipunto a puntuales. A

continuación ya se pudo trabajar a nivel puntual y obtener las coordenadas planimétricas de

cada uno de los puntos junto a su respectivo identificador.

4. Una vez se tuvieron las coordenadas x, y junto al campo idd tipo entero, se empleó la

herramienta “Grouping Analysis”. Esta herramienta se encargó de agrupar los puntos LIDAR

en función de la distancia euclidea a la que se encuentrasen unos de otros. Para ello se

apoyó de los campos x, y creados en el shp de singlepart.

33

Por último solo quedó especificar el numero de grupos que se querían obtener, que fue el

numero de copas de árboles de estudio (25).

5. Tras tener delimitadas las copas de los árboles de estudio a partir de los puntos LIDAR, se

empleó la herramienta “minimum bounding geometry”. La cual obtuvo la mínima

envolvente de cada grupo de puntos obtenido en el paso anterior.

6. Y para finalizar se calculó el área de cada poligono mediante la herramienta “calculate area”.

3.4.2. Cálculo de las variables estadísticas relacionadas con las

elevaciones e intensidades y número de puntos LIDAR por copa de

árbol

Para el cálculo de estas variables estadísticas se hizo uso del software FUSION. Estas variables

estadísticas relacionadas con las intensidades y elevaciones están compuestas por variables como

valores de máximo, mínimo, media, moda, desviación típica, varianza, coeficiente de asimetría y

kurtosis, así como los percentiles por intervalos. Tanto para variables de intensidad como de

elevación.

La metodología seguida para obtener estas variables fue la que se observa en la figura 19:

Figura 19. Esquema de la obtención de variables estadísticas con FUSION.

Para filtrar y normalizar los puntos LIDAR a partir del modelo digital del terreno escogido en el

anterior apartado, se empleó el comando “clipdata”, cuyo código empleado se muestra a

continuación:

clipdata /shape:0 /dtm:dtm1metro_fusion.dtm /zmin:0.5 /zmax:6 /height

filtradoZ_zonaEstudio.lda filtradoZ_zonaEstudio_clipdata.lda 704222 4420363 704359 4420712

Datos LIDAR preprocesados y clasificados

Filtrado y normalización de puntos LIDAR

Extracción de los puntos LIDAR de las copas de los árboles de estudio

34

Como lo que nos interesa es la copa de los

árboles, se seleccionó una cota mínima de 0,5

metros y una máxima de 6 metros. Para llegar

a la conclusión de estos valores, se tomaron

como referencia los datos de campo

referentes a la altura máxima de los árboles

de estudio y del MDSN. El cual se reclasificó

para contrastar que con esos valores

seleccionados, se seleccionaban aquellos

puntos LIDAR pertenecientes a las copas de

los árboles de estudio.

Figura 20. Modelo digital de superficies normalizado con shp de los árboles de estudio

Tras tener los puntos LIDAR normalizados y filtrados se empleó el comando “PolyClipData”, que

recorta el archivo .las obtenido mediante el clipdata a partir de un shp con las copas de los árboles a

estudiar. El código empleado dividió los puntos LIDAR en función de a qué copa de árbol de estudio

pertenecían, obteniendo un archivo de puntos LIDAR para cada copa de árbol.

PolyClipdata/multifile segmentacion_arboles_final.shp arboles_estudio.lda

filtradoZ_zonaEstudio_clipdata.lda

35

Y de cada copa de árbol se extrajeron las variables características con su correspondiente

identificador de árbol, mediante el comando “cloudmetrics”

Cloudmetrics/new /id arboles_estudio1.lda

C:\Lidar_Fusion\extraccion_caracteristicas_PolyclipData\estadisticas_arboles_estudio_1.csv

3.4.3. Cálculo de la variable densidad

Para obtener esta variable se empleó el software Lastools, mediante el cual, se obtuvieron los

puntos LIDAR que caían en cada uno de los intervalos fijados y sus correspondientes porcentajes de

densidades respecto a la totalidad de puntos LIDAR de cada copa de árbol de estudio.

Los siete intervalos fijados se han establecido con un rango de medio metro en medio metro,

empezando por 0.5 metros para no seleccionar puntos de vegetación baja, hasta los 4 metros, que

se corresponde con la altura más elevada de los árboles de estudio.

También se ha elegido el intervalo de 1.5 a 3.5 metros, el cual se introdujo debido a que en otros

estudios se ha comprobado que es en el intervalo donde aparecen la mayoría de puntos LIDAR de la

copa del árbol.

3.4.4. Cálculo de variable volumen obtenido mediante voxelización

Para obtener esta variable se utilizó el software de Matlab, a partir del cual, se obtuvo la variable

volumen sólido de las copas de los árboles de estudio, para emplearla como una variable

independiente en el análisis estadístico.

Esta variable se obtuvo a partir de los puntos LIDAR de las copas de los árboles ya normalizados.

Se hizo uso del software FUSION para realizar la conversión de archivo .las a archivo .xyz, para

obtener información en un txt del punto LIDAR y sus respectivas coordenadas x, y, z, separadas

mediante comas.

Una vez obtenidos los ficheros de coordenadas x, y, z para cada copa de árbol, mediante Matlab

se obtuvo el volumen sólido de cada copa mediante el método voxel. Este método consiste en

determinar el sumatorio de los volúmenes de todos los voxels existentes en las copas de los árboles.

Un voxel no es más que una superficie en tres dimensiones, con un valor de tamaño

determinado. En nuestro caso cada vez que haya un punto LiDAR con una altura superior a 0,5 m

dentro de un cubo con dimensiones por lado de 0,5 m se creará un voxel (Fernández-Sarría et al.,

2013).

36

3.4.5. Variables tipo raster y frecuencia de cabida cubierta

Se trata de una aplicación basada en IDL que permite procesar datos LIDAR. Esto se lleva a cabo

mediante el archivo con los puntos LIDAR filtrados, el fichero vectorial en formato shape con las

áreas de las copas y el modelo digital del terreno seleccionado anteriormente como el elegido para

el trabajo, siendo los resultados obtenidos un extenso listado de variables. De las cuales se elegirán

únicamente aquellas que no se hayan obtenido previamente mediante cualquier otro software de

los citados, siendo estas las siguientes:

Meanraster: media de los valores raster

Stdvraster: desviación típica de los valores raster

Rangeraster: rango de oscilación de los valores raster

Kurtraster: coeficiente de kurtosis de los valores raster

Skewraster: coeficiente de asimetría de los valores raster

Edgmeanraster: media de los valores raster de los bordes

Edgstdvraster: desviación típica de los valores raster de los bordes

Dens0: tanto por cien o proporción de puntos al nivel del terreno

Hdmax: Altura donde mayor número de puntos LIDAR se registran

Pdmax: Proporción de puntos LIDAR registrados en hdmax respecto al total

Hdmax2: Se trata de la misma variable que Hdmax, pero para un segundo lugar.

Pdmax2: Se trata de la misma variable que Pdmax, pero para un segundo lugar.

Fcc: frecuencia cabida cubierta, es la superficie que cubren los puntos LIDAR frente a la que

se queda sin cubrir.

3.5. Interpretación de los resultados de los modelos de regresión

En este apartado se explicó en qué consiste el tipo de análisis realizado, regresión lineal múltiple,

como se analizan los resultados y la metodología seguida.

37

Las variables dependientes tomadas en campo que se estimaron para todos los árboles de

estudio, fueron;

Altura de copa

Altura total

Peso poda 2013 (kg)

Producción 2012

Diámetro de tronco

Para el conjunto de 17 árboles de estudio del Grupo 2 fueron las mismas variables exceptuando

el diámetro de tronco.

Donde para llevar a cabo los modelos de regresión múltiple se empleo el método “paso hacia

delante”, poniendo un límite de tres variables independientes para explicar el modelo. Llevándose a

cabo la evaluación de los modelos por medio del coeficiente de determinación R2.

3.5.1 Hipótesis sobre el término de perturbaciones

Los supuestos de un modelo estadístico se refieren a una serie de condiciones que deben darse para garantizar la validez del modelo. Al efectuar aplicaciones prácticas del modelo de regresión, nos veremos en la necesidad de examinar muchos de estos supuestos.

1. Linealidad:

La ecuación de regresión adopta una forma particular. En concreto, la variable dependiente es la suma de un conjunto de elementos: el origen de la recta, una combinación lineal de variables independientes y los residuos. El incumplimiento del supuesto de linealidad suele denominarse error de especificación.

Se emplea el diagrama de dispersión de los residuos estandarizados frente a los valores de la

variable predichos. Y debe de ser una nube de puntos sin forma entre el intervalo [-2, 2], cuya media de los valores residuos sea cero.

2. Independencia:

Los residuos son independientes entre sí, es decir, los residuos constituyen una variable aleatoria (recordemos que los residuos son las diferencias entre los valores observados y los pronosticados). Es frecuente encontrarse con residuos autocorrelacionados cuando se trabaja con series temporales.

Se emplea el coeficiente de Durbin-Watson y se acepta la independencia de residuos cuando se

obtengan valores entre 1,5 y 2,5.

38

𝐷𝑊 = (𝑒𝑖−𝑒𝑖−1)2𝑛

𝑖=2

𝑒𝑖2𝑛

𝑖=1

También se emplea la gráfica de residuos estandarizados frente al número de fila, donde debe de

aparecer una nube de puntos sin forma, para que los residuos sean independientes.

3. Homocedasticidad: Para cada valor de la variable independiente (o combinación de valores de las variables

independientes), la varianza de los residuos es constante. Y se empleará la misma gráfica que para la linealidad, pero esta vez se estudiará que los puntos no tengan tendencia a aumentar tanto sus valores positivos como negativos a medida que nos alejamos del origen de coordenadas.

4. Normalidad:

Para cada valor de la variable independiente (o combinación de valores de las variables independientes), los residuos se distribuyen normalmente con media cero.

Se emplean los histogramas de los residuos para comprobar la normalidad de los residuos

visualmente y se realiza el test de Kolmogoroff-Smirnoff-Lilliefors y Shapiro-Wilk. 5. No-colinealidad:

No existe relación lineal exacta entre ninguna de las variables independientes. El incumplimiento

de este supuesto da origen a colinealidad o multicolinealidad. Para comprobar la colinealidad se emplea el coeficiente r de Pearson entre las variables

dependientes y la independiente. Debiendo de presentar alta correlación entre la variable independiente y las variables dependientes, pero baja correlación entre las variables dependientes del mismo modelo.

6. Residuos studentizados:

Estos residuos son estandarizados de una manera singular: cada residuo i-ésimo es dividido por la

desviación estándar de todos los residuos salvo el i-ésimo. Studentizacion es una vieja palabra que significa estandarización. Es una forma de estandarizar no muy lejana al cálculo de los valores z. En la estandarización z, a cada valor se le resta la media y se divide por una única desviación estándar del conjunto de residuos. En la studentización de los residuos no es necesario restar la media, ya que la media de los residuos es cero. Se divide por una desviación estándar distinta para cada elemento. La desviación utilizada, se calcula utilizando todos los residuos, salvo el que está siendo considerado.

En conclusión, los residuos studentizados permiten obtener los outliers de la relación.

39

3.5.2 Hipótesis sobre los parámetros del modelo

Al realizar el contraste de hipótesis estadístico, se siguieron los siguientes pasos:

1. Se especificó un nivel de significación , lo que cuantifica la tolerancia para un error de tipo I.

2. Se empleó un valor crítico, obtenido mediante tablas estadísticas en función del valor de , para contrastar las hipótesis formuladas.

3. Y por último se definió una regla de decisión para determinar si la hipótesis nula debe ser,

o no, rechazada en función del valor que tome el estadístico construido frente al valor crítico.

Dicho contraste de hipótesis pueden derivar en estos dos tipos de errores:

1. Error tipo I: es la probabilidad de rechazar la H₀ cuando es cierta. Para el estudio, se

estableció un =0.05, lo que quiere decir que se acepta rechazar H₀ falsamente en un 5%

de los casos.

2. Error tipo II: β es la probabilidad de rechazar H₀ cuando es falsa.

La hipótesis de permanencia estructural, fue la establecida para los parámetros del modelo. La

cual cita que los parámetros poblacionales, βj, se mantienen constantes a lo largo de toda la

muestra. Siendo estas las siguientes en función del estadístico empleado:

3.5.3 Modelo matemático

Mediante el modelo de regresión lineal múltiple se trata de explicar el comportamiento de una

determinada variable, llamada variable dependiente y representada mediante la letra Y. Esta variable está en función de un conjunto de k variables explicativas 𝑥0, 𝑥1, 𝑥2, ..., 𝑥𝑘 mediante una relación de dependencia lineal (suponiendo 𝑥0= 1):

𝑌 = 𝛽0 + 𝛽1*𝑥1 + 𝛽2*𝑥2 +…+ 𝛽𝑘*𝑥𝑘 + 𝜀

H₀: 𝛽𝑗 = 0

H₁: 𝛽𝑗 ≠ 0

Para H₁ con dos colas en el

estadístico T-Student

H₀: 𝛽1 = 𝛽2 = ⋯ = 𝛽𝑘 = 0

H₁: 𝛽𝑗 ≠ 0 para algún j = 1, … , k

Para el estadístico F-Snedeccor

40

Donde 𝜀 es la variable que representa los residuos, siendo estos la diferencia entre el valor

estimado por el modelo y el verdadero valor de Y. Y los parámetros 𝛽𝑗 , j = 1, … , k se llaman

coeficientes de regresión y han sido estimados mediante el método de mínimos cuadrados.

Para n observaciones se tiene:

𝑦1 = 𝛽0 + 𝛽1*𝑥11 + 𝛽2*𝑥12 +…+ 𝛽𝑘*𝑥1𝑘 + 𝜀1 . . . . . . . . . . . . . . . . . . 𝑦𝑛 = 𝛽0 + 𝛽𝑛*𝑥11 + 𝛽2*𝑥𝑛2 +…+ 𝛽𝑘*𝑥𝑛𝑘 + 𝜀𝑛

Lo que matricialmente queda como:

𝑌 = 𝑋 ∗ 𝛽 + 𝜀

𝑦1

𝑦2

⋮𝑦𝑛

=

1 𝑥11 𝑥12 ⋯ 𝑥1𝑘

1 𝑥21 𝑥22 ⋯ 𝑥2𝑘

⋮ ⋮ ⋮ ⋮ ⋮1 𝑥𝑛1 𝑥𝑛2 ⋯ 𝑥𝑛𝑘

*

𝛽0

𝛽1

⋮𝛽𝑘

+

𝜀0

𝜀1

⋮𝜀𝑛

Donde el vector 𝛽 estimado es:

𝛽 = 𝑋𝑇 ∗ 𝑋 −1 ∗ (𝑋𝑇 ∗ 𝑌) Quedando así la ecuación ajustada:

𝑌 = 𝑋 ∗ 𝛽

3.5.4 Significación estadística

En este apartado se diferenció entre la significación estadística individual y la global. Empleando las siguientes reglas de decisión para contrastar las hipótesis:

T-Student: Para analizar la significación estadística individual de cada variable independiente del modelo.

41

o Si 𝑡𝑗 ≥ 𝑡𝑛−𝑘 ,∝/2 , el estadístico cae fuera de la región de aceptación, por lo que se

rechaza la H₀.

o Si 𝑡𝑗 < 𝑡𝑛−𝑘 ,∝/2 , el estadístico cae dentro de la región de aceptación, por lo que no

se puede rechazar la H₀. Por ello, el parámetro no es individualmente significativo.

Siendo n el numero de observables (25 árboles de estudio) y k el número de parámetros del modelo ( 𝛽 ). Donde:

𝑡𝑗 =𝛽𝑗 −𝛽𝑗

𝑜

𝑒𝑒 (𝛽𝑗 )

, siendo ee(βj ) el error estándar del coeficiente 𝛽𝑗

F-Snedeccor: Para analizar la significación estadística global del modelo.

o 𝐹 ≥ 𝐹𝑞 ,𝑛−𝑘∝ , se rechaza H₀

o 𝐹 < 𝐹𝑞 ,𝑛−𝑘

∝ , no se rechaza H₀

Siendo q los grados de libertad del numerador y n-k los grados de libertad del denominador.

Donde el estadístico F, se obtuvo a partir de la tabla 17:

Variación Suma de cuadrados g.d.l Media cuadrática F

Explicada 𝑆𝐶𝐸 = (𝑦 𝑖 − 𝑦 )2

𝑛

𝑖=1

k SCE/k (𝑆𝐶𝐸/𝑘)

𝑆𝐶𝑅/(𝑛 − 𝑘 − 1)

Residual 𝑆𝐶𝑅 = (𝑦𝑖 − 𝑦 𝑖)2

𝑛

𝑖=1

n-k-1 SCR/(n-k-1)

Total 𝑆𝐶𝑇 = (𝑦𝑖 − 𝑦 )2

𝑛

𝑖=1

n-1

Tabla 17. Tabla ANOVA

3.5.5 Inferencia estadística

Tras calcular los modelos, el siguiente paso fue analizar si la regresión era válida para predecir las variables dependientes a partir de la influencia de las variables independientes. Para ello se contrastó que la correlación entre ambas variables fuese próxima a 1 o -1, lo que parece indicar una existencia de correlación lineal entre los valores de la muestra. Pero este valor del coeficiente de correlación no garantiza que también estén correlacionadas en la población.

42

Para poder contrastar este caso, se partió de la estimación de la recta de regresión y de los parámetros del modelo

𝑌 = 𝛽 ₀ + 𝛽 1*𝑥1

Y se comprobó que el valor de los parámetros (𝛽) fuesen significativamente distintos de cero.

Existiendo en tal caso una correlación lineal entre ambas variables poblacionales.

Para representar la bondad del ajuste se empleó el coeficiente de determinación 𝑅2 ajustado, ya que solo aumentará su valor al introducir variables independientes que reduzcan el cuadrado medio de los residuales.

𝑅𝑎𝑑𝑗2 = 1 −

𝜀𝑖 2

𝑖

𝑦𝑖2

𝑖 *

(𝑛−1)

(𝑛−𝑘)

Figura 21. Ecuación del coeficiente de determinación ajustado

Siendo;

𝑦𝑖 = 𝑌𝑖 − 𝑌 𝑖 , donde su sumatorio será la suma de los cuadrados totales

𝑦 𝑖 = 𝑌 𝑖 − 𝑌 𝑖 , donde su sumatorio será la suma de los cuadrados residuales

ε i = Yi − Y i , donde su sumatorio será la suma de los cuadrados de la regresión

Este coeficiente puede interpretarse como un porcentaje de reducción de la incertidumbre

cuando las variables independientes son conocidas, teniendo mayor poder explicativo cuanto más

próximo a 1 esté. Es por esto, por lo que se analizaron los valores de este coeficiente con cada

variable independiente introducida en el modelo. Para así ver la cantidad de la muestra que se

representaba con cada variable independiente y estudiar si merecía la pena introducir más variables,

o la cantidad de muestra que mejora, no era significativa para introducir una variable independiente

más.

4. Resultados

4.1. Clasificación LIDAR automática

A continuación se procede a mostrar un TIN de elevaciones obtenido a partir de los puntos

LIDAR para demostrar que los puntos que tienen una cota normalizada por debajo de 0 metros son

tan escasos que no se aprecian visualmente.

43

Figura 22. TIN de elevaciones de la clasificación automática de los puntos LIDAR mediante Lastools

Respecto a la cantidad de puntos LIDAR que se quedan sin clasificar mediante este método, cabe

señalar que serán puntos que sean inferiores a 0.5 y no se han clasificado como suelo. Como se

puede apreciar en la figura 23 dichos puntos corresponden con las alineaciones que forman las

terrazas de la zona de estudio y los taludes que existen junto a la carretera en la parte derecha de la

figura 22. La cual se muestra en la figura 24.

44

Figura 23. Superposición del TIN de elevaciones de la figura 22 y los puntos LIDAR sin clasificar. Y alzado de los puntos LIDAR sin clasificar

45

Figura 24. Detalle zona de taludes junto a la carretera, ubicada a la parte derecha de la zona de estudio.

También habrá otro bloque de puntos sin clasificar pertenecientes a la zona donde la vegetación

alcanza sus máximas alturas, sobrepasando los 5 metros que se habían marcado como máximo en el

intervalo de vegetación para los árboles de estudio. Lo cual se aprecia perfectamente en la figura 23.

Y el resto de puntos aislados podrían deberse a hierbas bajas inferiores a 0.5 metros.

4.1.1. Bondad de la clasificación automática mediante la clasificación

manual

Una vez realizada la clasificación automática se procede a contrastar los puntos LIDAR

pertenecientes a las copas de los árboles de estudio con los obtenidos en la clasificación manual

para comprobar la bondad del ajuste automático. Para ello se tomará el shp de las copas de los

árboles de estudio como recorte, obteniendo los siguientes resultados:

46

Id Nºpuntos LIDAR_manual Nºpuntos LIDAR_automaticos % bien

clasificado

1 107 107 100

2 116 116 100

3 124 124 100

4 107 107 100

5 92 92 100

6 100 100 100

7 94 94 100

8 97 97 100

88 93 93 100

17 38 38 100

18 58 58 100

19 61 61 100

20 51 51 100

21 68 68 100

22 54 54 100

23 84 84 100

24 56 56 100

25 23 23 100

26 29 29 100

27 42 42 100

28 86 86 100

29 124 124 100

30 105 105 100

31 86 86 100

32 88 88 100

Tabla 18. Porcentaje de puntos LIDAR por árbol bien clasificados automáticamente mediante Lastools.

47

El porcentaje total será igual a:

Total Nºpuntos LIDAR_manual Total Nºpuntos LIDAR_automaticos % bien clasificado

1983 1983 100

Tabla 19. Porcentaje de puntos LIDAR total bien clasificados automáticamente mediante Lastools

4.2. MDT

Este ha sido un apartado fundamental para obtener unos estadísticos de la distribución de las

alturas sin errores asociadas a la superficie topográfica El MDT se ha utilizado para normalizar los

puntos LIDAR, de los que se obtendrán las variables estadísticas empleadas en el cálculo de los

modelos de regresión como variables independientes.

Es por todo esto por lo que se han obtenido un total de 143 MDTs entre los cuatro programas

citados.

4.2.1. MDT mediante FUSION

A continuación se va a mostrar un resumen de los resultados estadísticos de los modelos más

significativos, que han hecho que lleguemos a obtener el óptimo:

Figura 25. Gráficas variación rmse en función de los parámetros de FUSION

Donde finalmente, el mejor MDT obtenido mediante este software tiene un rmse de 0,2486 m con los siguientes parámetros y se puede visualizar en la tabla 20:

rmse = 0,2486 m

Parámetros empleados:

tamaño celda = 1 metro

w = 0.9

g = -1

a = 5

b = 4

iteraciones = 5

Tabla 20. Parámetros y rmse del mejor MDT

4.2.2. MDT mediante lidex

En este apartado se va a representar un resumen con la gráfica y rmse más precisos obtenidos

para cada una de las variaciones de las tolerancias explicadas en el apartado “Obtención de los

modelos mediante lidex”.

Tabla 21. Resumen con los mejores resultados rmse para cada grupo de ventanas

Mejores resultados de cada grupo de ventanas

identificador tamaño de celda ventana1 ventana2 ventana3 tolerancia1 tolerancia2 rmse

1 1 20 10 5 0,5 0,3 0,2892 2 1 12 6 3 0,5 0,5 0,303 3 1 10 5 2,5 0,5 0,5 0,3145 4 1 8 4 2 1,5 0,5 0,3169 5 1 25 12,5 6,25 0,5 0,5 0,3808 6 1 20 4 2 0,5 0,5 0,314 7 1 4 2 1 0,5 0,5 0,3234 8 1 6 3 1,5 0,5 0,5 0,3194

50

Figura 26. Resumen rmse de los mejores modelos con lidex

4.2.3. MDT mediante Lastools y mediante ArcGIS

Este apartado se va a redactar junto, ya que solo se ha generado 1 modelo con cada software

debido a la escasa opción de modificar los parámetros para su generación.

Los resultados obtenidos respecto al rmse de cada MDT se muestran en la siguiente tabla:

MDT Lastools MDT ArcGIS

Tamaño celda 1 1

Rmse (m) 0,3406 0,3462

Tabla 22. Rmse de los modelos generados en ArcGIS y Lastools

Finalmente se muestra el MDT y el MDS obtenido mediante Lastools, para que se pueda apreciar

que ya se han conseguido dos modelos donde se distinguen claramente las copas de los árboles en

el MDS y las terrazas en el MDT.

Cabe mencionar que al MDS se le ha generado un mapa de sombras para conseguir una óptima

visualización y que los puntos que aparecen en ambas imágenes hacen referencia a los árboles de

estudio.

51

Figura 27. MDS Lastools zona de estudio con mapa de sombras Figura 28. MDT Lastools zona de estudio

4.3. Comparativa variable área manual vs área automática

En este apartado se procede a mostrar la comparativa entre ambas áreas, mediante la tabla 23

que refleja los valores de diferencia obtenidos para cada árbol de estudio, así como el RMSE

obtenido.

52

Comparativa áreas

id área automática área manual diferencia

1 7.85665000063 8.64207948859 0.78542948796

2 11.83765000050 11.49167056300 -0.34597943750

3 12.89915000010 13.78311759330 0.88396759320

4 12.54949999910 11.94715352570 -0.60234647340

5 10.28310000000 10.69436590640 0.41126590640

6 13.02934999990 11.87629198630 -1.15305801360

7 11.41530000030 11.98577281520 0.57047281490

8 11.04585000040 11.88190631050 0.83605631010

17 11.89544999980 14.38496657570 2.48951657590

18 20.84935000000 22.24137706960 1.39202706960

19 20.64310000030 24.48807120880 3.84497120850

20 17.11099999980 19.97114595940 2.86014595960

21 27.21644999990 28.38198624090 1.16553624100

22 19.70234999980 21.83396346770 2.13161346790

23 24.17945000000 26.36820820330 2.18875820330

24 13.04580000030 15.26477988380 2.21897988350

25 6.36825000005 8.96392500554 2.59567500549

26 8.95650000060 12.65694206430 3.70044206370

27 13.89364999990 14.24018290550 0.34653290560

28 13.01060000070 14.84979901870 1.83919901800

29 7.74835000017 10.19909849790 2.45074849773

30 9.20310000030 10.11028130520 0.90718130490

31 8.54249999987 9.94937699119 1.40687699132

32 10.08140000010 11.04577504950 0.96437504940

88 12.33235000120 11.93481483910 -0.39753516210

Tabla 23. Comparativa área manual vs automática y estadísticos columna diferencia. Unidades en metros

desv.tipica 1.292

promedio 1.340

máximo 3.845

mínimo -1.153

53

Como muestran los resultados el

área manual de color rosa en la figura

29, presenta valores mayores que el

área automática.

Figura 29. Superposición área automática y manual con puntos LIDAR clasificados. Puntos LIDAR color marrón representan puntos terreno y color verde, vegetación.

4.4. Variables estadísticas

Tras tener agrupados los árboles de estudio según el apartado “Clasificación de árboles de

estudio”, se procede a obtener las regresiones múltiples para cada uno de los grupos mencionados.

Seguidamente se muestran los resultados estadísticos obtenidos para cada variable del Grupo 1:

Diámetro de tronco:

Al contrastar las correlaciones de las variables del modelo se comprueba que existe una

correlación baja entre las variables independientes del modelo, dándose la mayor correlación entre

la variable área copa automática e int mínimum con un valor de 0.35. Lo cual demuestra que no

existe colinearidad entre las variables independientes del modelo. Sin embargo se obtienen

correlaciones elevadas entre las variables independientes respecto a la variable dependiente,

dándose la mayor correlación con la variable área copa automática con un valor de 0.80.

Variables

Diámetro tronco

Área copa automática

Int mínimum

Int L swekness

Diámetro tronco

1 0.796 0.660 -0.152


0.796 1 0.350 -0.027

Int minimum 0.660 0.350 1 0.265

Int L skewness -0.152 -0.027 0.265 1

Tabla 24. Matriz de correlación para variable diámetro de tronco en Grupo 1

54

Al introducir la primera variable del modelo área copa automática, se obtiene un valor de R²

ajustado de 0.618. Introduciendo la segunda variable Int mínimum se mejora hasta 0.781 e

introduciendo la última variable Int L skewness se llegó a 0.848.

N=25

R= 0.93118699 R²= 0.86710922 Adjusted R²= 0.84812482

F(3,21)=45.675 p<0.00000 Std.Error of estimate: 7.6860

Beta Std.Err. of Beta

B Std.Err. of B

t(21) p-level

Intercept

-14.0174 4.41419 -3.17554 0.004555

Área copa automática 0.606532 0.085692 2.2485 0.31767 7.07804 0.000001

Int minimum 0.519555 0.088836 5.6926 0.97335 5.84850 0.000008

Int L skewness -0.273419 0.083236 -72.6979 22.13122 -3.28486 0.003532

Tabla 25. Resumen regresión múltiple para variable dependiente diámetro de tronco en Grupo 1

Tras obtener los modelos se calcula la tabla ANOVA para contrastar el modelo de regresión lineal

global y se cumple que el p-valor obtenido de la tabla 26 es inferior al valor de significancia

introducido para obtener el modelo de 0.05. Lo que implica que este modelo sirve para explicar la

variable dependiente.

Variación Suma de

cuadrados Grados de

libertad Media

cuadrática F P-valor

Explicada

8094.569 3 2698.190 45.675 0.000000

Residual 1240.551 21 59.074

Total 9335.120

Tabla 26. Tabla ANOVA para variable diámetro de tronco en Grupo 1

Seguidamente se realizan los contrastes de significatividad individual de los parámetros mediante

el estadístico t-student y se obtiene que los p-valor de cada una de las variables independientes son

menores que el valor de significación introducido inicialmente para crear el modelo (0.05), como se

observa en la tabla 25. Por lo que finalmente se rechazan las Ho, lo que significa que la variable

dependiente depende linealmente de cada una de las variables independientes del modelo.

55

Tras pasar los contrastes de significatividad se analizan los residuos para comprobar si siguen

una distribución normal. Como se observa a través del histograma de la figura 30 y a partir de los

contrastes de normalidad de Kolmogorov-Smirnov con corrección de Lilliefors y Shapiro-Wilk, los

residuos estandarizados siguen una distribución normal, ya que su p-valor es mayor al valor de

significatividad introducido en el modelo, para ambos contrastes de normalidad (0.57771 > 0.05 y

0.20 > 0.05). Por lo que se rechaza la Ho.

Figura 30. Histograma residuos estandarizados diámetro tronco

Seguidamente se comprueba la linealidad y la homocedasticidad de los residuos a partir de la

gráfica de dispersión de los residuos estandarizados frente a la variable predicha. Y como se observa

en la figura 31, los puntos se encuentran distribuidos indistintamente arriba y abajo del 0 en el

intervalo [-2, 2] y no presentan una tendencia a aumentar los valores de los residuos a medida que

aumentan los valores de la variable predicha, lo cual indica una varianza constante en los residuos.

Por lo que se cumple la homocedasticidad.

56

Figura 31. Gráfico de dispersión entre los residuos estandarizados y la variable predicha para diámetro de tronco.

A continuación se va a estudiar la independencia de los residuos a partir del coeficiente Durbin

Watson, el cual debe de estar próximo a 2 para que esto se cumpla. Siendo para los residuos de la

variable diámetro de tronco igual a 2.11, con un coeficiente de correlación próximo a cero (-0.06), lo

cual indica que son independientes. También es de utilidad la figura 32, la cual indica la presencia de

residuos independientes si la nube de puntos no tiene forma.

57

Figura 32. Gráfico de dispersión entre los residuos estandarizados y el número de fila para diámetro de tronco.

Por último se comprueban los valores anómalos mediante los residuos estudentizados y se

obtiene que solamente 3 residuos presentan valores superiores a 2, que son los árboles 17, 26 y 27

cuyos valores respectivos son (2.22, 2.11 y 2.14). Al eliminar estos árboles con valores de residuos

estudentizados anómalos, se generan árboles cuyos valores son mayores a los ya obtenidos. Por lo

que se decide no eliminarlos ya se parte de una muestra de datos pequeña y los valores están muy

próximos al límite, considerando que no afectan en exceso al conjunto.

Producción:

Tabla 27. Resumen modelo para producción

Cabe mencionar que se ha realizado el mismo estudio para la variable producción con el grupo 2 y empleando las variables total return, elev

mínimum e int p75 se ha obtenido un 𝑅2 ajustado de 0.67.

Variables

modelo

Matriz de correlación

Método

paso hacia

delante

Contraste

global

Contraste

individual

Contraste de

normalidad

residuos

Contraste

independencia

residuos

Gráfica residuos

estandarizados vs

valor predicho

Residuos

anómalos

Máxima correlación

entre variables

independientes

Máxima

correlación entre

variable

dependiente-

independiente

R cuadrado

ajustado

P-valor

tabla

ANOVA

p-valor

t-student

p-

valor

K-S

p-valor

Shapiro

Wilk

Durbin

watson Coef.corr

Homocedasticidad,

linealidad

ID

árbol valor

Área copa

automática Área copa

automática e int

mode (0.17)

Área copa

automática-

producción (0.42)

0.143

0.008

0.009

0.20 0.93 1.53 0.22

Se obtiene un

modelo lineal y

homocedástico

4

25

2.09

2.72 Int mode 0.299

0.02

Constante 0.000048

Peso poda:

Tabla 28. Resumen modelo para peso poda

Variables

modelo


Método

paso hacia

delante

Contraste

global

Contraste

individual

Contraste de

normalidad

residuos

Contraste

independencia

residuos

Gráfica residuos

estandarizados vs

valor predicho

Residuos

anómalos


entre variables

independientes

Máxima

correlación entre

variable

dependiente-

independiente

R cuadrado

ajustado

P-valor

tabla

ANOVA

p-valor

t-student

p-

valor

K-S

p-valor

Shapiro

Wilk

Durbin

watson Coef.corr

Homocedasticidad,

linealidad

ID

árbol valor


Área copa

automática-

Kurtraster (0.25)

Área copa

automática-peso

poda (0.92)

0.835

0

0

0.20 0.6 2.14 -0.08

Se obtiene un

modelo lineal y

homocedástico

20

3.26

Int maximum

0.88

0.001861

Kurtraster

0.901

0.027

Constante 0.24

Altura total:

Tabla 29. Resumen modelo para altura total

Cabe mencionar que se ha realizado el mismo estudio para la variable altura total con el grupo 2 y empleando las variables elev máximum, elev P01 y

elev P05 se ha obtenido un 𝑅2 ajustado de 0.799.

Variables

modelo


Método

paso hacia

delante

Contraste

global

Contraste

individual

Contraste de

normalidad

residuos

Contraste

independencia

residuos

Gráfica residuos

estandarizados vs

valor predicho

Residuos

anómalos


entre variables

independientes

Máxima

correlación entre

variable

dependiente-

independiente

R cuadrado

ajustado

P-valor

tabla

ANOVA

p-valor

t-student

p-

valor

K-S

p-valor

Shapiro

Wilk

Durbin

watson Coef.corr

Homocedasticidad,

linealidad

ID

árbol valor

Volumen

chglobal

Kurtraster-volumen

chglobal (0.19)

Volumen

chglobal-altura

total (0.73)

0.514

0.00002

0.000083

0.20 0.13 2.48 -0.27

Se obtiene un

modelo lineal y

homocedástico

19

29

2.57

2.50

Kurtraster No converge

0.022

Dens0

0.630

0.042

Constante 0

61

Altura copa:

Tabla 30. Resumen modelo para altura copa

En esta variable también se ha realizado el mismo estudio con el grupo 2 y empleando las variables elev máximum, elev P01 y elev P05 se ha

obtenido un 𝑅2 ajustado de 0.773.

Variables

modelo


Método

paso hacia

delante

Contraste

global

Contraste

individual

Contraste de

normalidad

residuos

Contraste

independencia

residuos

Gráfica residuos

estandarizados vs

valor predicho

Residuos

anómalos

Máxima correlación entre variable

dependiente-independiente

R cuadrado

ajustado

P-valor

tabla

ANOVA

p-valor

t-student

p-

valor

K-S

p-valor

Shapiro

Wilk

Durbin

watson Coef.corr

Homocedasticidad,

linealidad

ID

árbol valor

Elev maximum

Elev maximum-altura copa (0.60)

0.334

0.0015

0.0015

0.20 0.47 2.87 -0.44

Se obtiene un

modelo lineal y

homocedástico

23 2.99

19 2.58

Constante 0.24 29 2.01

Seguidamente se muestra una tabla resumen con los resultados más significativos obtenidos de

las rectas de regresión múltiple. En ella también se van a mostrar los resultados obtenidos para el

grupo de 8 árboles (grupo 3), ya que será un punto a citar en las conclusiones finales:

Variable Grupo 1 Grupo 2 Grupo 3

Producción (kg) 0.299 0.667 0.986

Peso poda 2013(kg) 0 .901 0.592 0.974

Altura total (m) 0.630 0.799 0.953

Altura copa (m) 0.334 0.773 0.997

Diámetro tronco (cm) 0.848 ---------- ----------

Tabla 31. Resumen 𝐑𝟐 ajustados de las variables dependientes de estudio

Tras obtener los resultados estadísticos se procede a comprobar las rectas de regresión múltiple

obtenidas para los dos grupos de árboles de estudio, Grupo 1 y Grupo 2. Con esto se pretende

comprobar que los valores obtenidos de campo frente a los predichos por ambos ajustes se

obtienen con una confianza del 95%. Donde únicamente se visualizará la gráfica del mejor modelo

obtenido para la mayor muestra (n=25), que se corresponde con la variable Peso de poda.

Figura 33. Gráfica de valores observados frente a predichos para la variable Peso poda2013 Grupo 1

63

Tal y como se aprecia en la figura 33, tan solo existe un árbol de estudio que se aleja de la

obtención de los resultados a un 95% de probabilidad. Este árbol con un valor observado

próximo a 10 y predicho próximo a 15, coincide con el único árbol de estudio que presentó un

residuo estudentizado para la variable peso de Poda. Es por esto por lo que podría de haberse

dado algo error accidental en la toma de este dato en campo.

5. Conclusiones

5.1. Cálculo del MDT

Para la obtención del mdt que mejor se adaptase a las condiciones del terreno, se emplearon

cuatro programas diferentes, sumando un total de 143 MDTs. De entre los diferentes programas se

concluye que aquellos que permite modificar mayor número de parámetros a la hora de generar los

MDTs son LIDEX y FUSION. Obteniendo el mejor MDT mediante el software FUSION con un RMSE de

0.2486 metros.

5.2. Clasificación de puntos LIDAR manual vs automática

Tras realizar la clasificación de forma manual y de forma automática mediante el software

lastools, se concluye que mediante el intervalo fijado de 0.5 a 5 metros en el proceso automático, se

consigue representar el 100% de los puntos LIDAR clasificados de forma manual con código 5 del

ASPRS de vegetación. Lo cual representaría un ahorro de tiempo para estudios con mayor número

de árboles tomados en campo.

5.3. Cálculo área de copa manual vs automática

Como se ha mostrado en el apartado de resultados, los valores de las áreas manuales eran

mayores que los de las automáticas. Esto es debido a que el método automático traza la linea del

poligono sobre los puntos de vegetación que delimitan el contorno y en el manual se tomaba mayor

área. Por lo que se concluye que la obtención del área automatica además de ser un proceso más

rápido y cómodo a la hora de obtener las áreas, también es mas preciso. Obteniendo una desviación

tipica para el campo de diferencias de áreas de 1,29 m y un RMSE de 1.84 m.

64

5.4. Análisis de modelos obtenidos

Tras la obtención de los cinco modelos obtenidos para las variables dependientes se extraen

diferentes conclusiones, citadas a continuación:

1. La variable más representativa en los modelos ha resultado ser el área de copa. De ahí la

importancia que ha tenido el poder obtener esta variable de forma automatizada.

2. El modelo que ha permitido predecir la variable dependiente con mayor probabilidad, ha

sido el obtenido para la variable Peso de poda, con un R² de 0.91. Lo cual resulta muy

interesante, ya que puede permitir estimar el potencial de biomasa en plantaciones

frutales, lo que puede ser utilizado como un elemento clave en la cadena de gestión de

biocombustibles.

3. En lo referente a la variable producción, que presenta un R² de 0.667, se concluye que para

que tenga mayor representatividad, dicha variable no podrá obtenerse exclusivamente

mediante variables obtenidas a partir de los datos LIDAR. También habrá que tener

presentes factores externos que afectan a la producción de los árboles, como puedan ser

enfermedades de los mismos, el tipo de abono o la meteorología que se dé en dicho año.

La tesis de Álamo Romero, S. (2004)(18), llega a obtener unas relaciones de variables

muy interesantes para el análisis de las zonas que presentan mayor producción. Obteniendo

menor cosecha en las zonas de mayor cota y deficiencia de nitrógeno y mayor cosecha en

zonas de menor cota y pendientes suaves inferiores al 10%.

4. Pese a que el grupo 3, formado por los árboles de estudio de mayor fuste fue descartado

debido a un tamaño de muestra insuficiente. Cabe destacar que se han predicho las

variables dependientes de los modelos alcanzando valores de R² ajustado superiores al 0.95.

Es por esto por lo que sería interesante la posibilidad de realizar una nueva toma de datos

en campo de árboles que entren en este grupo para aumentar la muestra y contrastar si

continua obteniendo valores tan elevados.

5. Por todo esto, se concluye que la tecnología LIDAR resulta de gran utilidad a la hora de

predecir variables dendrométricas para árboles frutales a raíz de los resultados obtenidos en

este proyecto.

65

6. Referencias bibliográficas

Álamo Romero S. Efectos económicos de la aplicación de la agricultura de precisión en una explotación de olivar en la provincia de Jaén. 2004.

Andersen H, McGaughey RJ, Reutebuch SE. Estimating forest canopy fuel parameters using LIDAR data. Remote Sens Environ 2005;94(4):441-449.

Estornell J, Ruiz L, Velázquez-Martí B, Fernández-Sarría A. Estimation of shrub biomass by airborne LiDAR data in small forest stands. For Ecol Manage 2011;262(9):1697-1703.

Estornell J, Velázquez-Martí B, López-Cortés I, Salazar D& Fernández-Sarría A. Estimation of wood volume and height of olive tree plantations using airborne discrete-return LiDAR data. 2013. p. 17-29.

Fernández-Sarría A, Velázquez-Martí B, Sajdak M, Martínez L, Estornell J. Residual biomass calculation from individual tree architecture using terrestrial laser scanner and ground-level measurements. Comput Electron Agric 2013;93:90-97.

Hall S, Burke I, Box D, Kaufmann M, Stoker J. Estimating stand structure using discrete-return lidar: an example from low density, fire prone ponderosa pine forests. For Ecol Manage 2005;208(1):189-209.

Holmgren J. Estimation of forest variables using airborne laser scanning. ; 2003.

Isenburg M, Liu Y, Shewchuk J, Snoeyink J, Thirion T. Generating raster DEM from mass points via TIN streaming. Geographic Information Science: Springer; 2006. p. 186-198.

Kraus K, Mikhail E. Linear least-squares interpolation 1972.

Kraus K, Pfeifer N. Determination of terrain models in wooded areas with airborne laser scanner data. ISPRS Journal of Photogrammetry and remote Sensing 1998;53(4):193-203.

LAS Specification Version 1.4-R13. : American Society for Photogrammetry and Remote Sensing Bethesda, Maryland; 2013.

Lefsky MA, Cohen W, Acker S, Parker GG, Spies T, Harding D. Lidar remote sensing of the canopy structure and biophysical properties of Douglas-fir western hemlock forests. Remote Sens Environ 1999;70(3):339-361.

Li Y, Andersen H, McGaughey R. A comparison of statistical methods for estimating forest biomass from light detection and ranging data. Western Journal of Applied Forestry 2008;23(4):223-231.

McGaughey RJ. Fusion. LDV: Software for Li-DAR Data Analysis and Visualization 2010.

Næsset E, Bjerknes K. Estimating tree heights and number of stems in young forest stands using airborne laser scanner data. Remote Sens Environ 2001;78(3):328-340.

66

Riaño D, Valladares F, Condés S, Chuvieco E. Estimation of leaf area index and covered ground from airborne laser scanner (Lidar) in two contrasting forests. Agric For Meteorol 2004;124(3):269-275.

Ruiz A, Kornus W. Experiencias y aplicaciones del Lidar. V Semana Geomática de Barcelona 2003;11(03).

Sociedad Española, y Teledetección de Ciencias Forestales, Grupo de Trabajo Inventario, Pando FJS, de Toda, Santiago Saura Martínez. Actas de la I Reunión de Inventario y Teledetección Forestal INVITEL 2004: Lleida, 25-26 marzo 2004. : Sociedad Española de Ciencias Forestales; 2005.

Zimble DA, Evans DL, Carlson GC, Parker RC, Grado SC, Gerard PD. Characterizing vertical forest structure using small-footprint airborne LiDAR. Remote Sens Environ 2003;87(2):171-182.

Obtencion y analisis de variables dendrometricas obtenidas mediante datos lidar en el ambito...

Technology

Transcript of Obtencion y analisis de variables dendrometricas obtenidas mediante datos lidar en el ambito...