Automatic Extraction of Tree Information in Forest Areas Using Local Maxima Based on Aerial LiDAR

항공 LiDAR 기반 Local Maxima를 이용한 산림지역 수목정보 추출 자동화

Abstract

Currently, the National Forest Inventory (NFI) collects tree information by human, so the range and time of the survey are limited. Research is actively being conducted to extract tree information from a large area using aerial Light Detection And Ranging (LiDAR) and aerial photographs, but it does not reflect the characteristics of forest areas in Korea because it is conducted in areas with wide tree spacing or evenly spaced trees. Therefore, this study proposed a methodology for generating Digital Surface Model (DSM), Digital Elevation Model (DEM), and Canopy Height Model (CHM) images using aerial LiDAR, extracting the tree height through the local Maxima, and calculating the Diameter at Breath Height (DBH) through the DBH-tree height formula. The detection accuracy of trees extracted through the proposed methodology was 88.46%, 86.14%, and 84.31%, respectively, and the Root Mean Squared Error (RMSE) of DBH calculated based on the tree height formula was around 5cm, confirming the possibility of using the proposed methodology. It is believed that if standardized research on various types of forests is conducted in the future, the scope of automation application of the manual national forest resource survey can be expanded.

현재 국가산림자원조사(National Forest Inventory, NFI)는 인력에 의한 수목정보를 수집하고 있어 조사 범위와 시간의 한계가 따른다. 항공 Light Detection And Ranging (LiDAR) 및 항공 사진 등을 이용하여 넓은 지역의 수목 정보를 추출하기 위한 연구가 활발하게 진행되고 있으나 수목의 간격이 넓은 지역이거나 수목의 간격이 일정하게 배치된 지역을 대상으로 이루어지고 있어 우리나라 산림지역 특성을 반영하지 못하고 있다. 이에 본 연구에서는 항공 LiDAR를 이용하여 수치표면모델(Digital Surface Model, DSM), 수치표고모델(Digital Elevation Model, DEM), 수목높이모델(Canopy Height Model, CHM) 영상을 생성한 후 local maxima 기법을 통해 수고를 추출하고 산정식을 통해 흉고직경(Diameter at Breast Height, DBH)을 산정하는 방법론을 제안하였다. 제안한 방법론을 통해 추출한 수목의 검출 정확도는 매목지구별 각 88.46%, 86.14%, 84.31%로 나타났으며, 수고 값을 기반으로 산정한 DBH의 평균제곱근오차(Root Mean Squared Error, RMSE)가 5 cm 내외로 나타나 제안한 방법론의 활용 가능성을 확인하였다. 향후 다양한 유형의 산림에 대한 표준화 연구를 진행한다면 수작업으로 이루어지는 국가산림자원조사의 자동화 적용 범위를 확대할 수 있을 것으로 사료된다.

1. 서론

수고, 흉고직경(Diameter at Breast Height, DBH), 임령 등은 산림 임목의 변화를 예측하는 임목 생장모델 개발에 필수적으로 요구되는 산림자원 정보이다(Kang et al., 2019). 산림자원 정보는 국가산림자원조사를 실시하여 전국 산림에 배치한 표본점을 대상으로 순환조사를 통해 수집된다. 표본조사 수행 시 수작업으로 가슴 높이에서 측정한 수목의 지름인 DBH를 측정하며 수고 정보의 경우 DBH-수고 산정식을 이용하여 산정한다(Korea Forest Service, 2016).

인력에 의해 수행되는 표본조사의 품질을 향상시키기 위해 Korea National Park Research Institute (2021)는 지상 LiDAR(Light Detection And Ranging) 데이터를 활용하여 수목의 구조정보를 분석하고자 하였다. 지상LiDAR로 취득한 수목의 구조정보는 기존 인력에 의해 측정된 수목의 구조정보에 비해 정확도가 향상되었으나 표본점 외의 수목 개체별 정보를 취득하기 어렵다는 한계는 해결되지 않고 있다.

항공 LiDAR 데이터를 활용한 수목 정보 추출 연구 사례를 살펴보면 도심지 가로수 또는 과수원 등 수목의 간격이 넓은 지역을 대상으로 연구가 활발히 수행되고 있다(Kim et al., 2009; Lee et al., 2014; Yang, J. et al., 2020; Liao et al., 2022). 또한, 수목 개체별 고품질 데이터 구축을 위해 무인항공사진측량 기법을 적용해 생성된 산출물과 기구축 공간정보를 활용해 수목 개체별 수고를 산정하는 연구가 진행되었다(Lee et al., 2017; Bang et al., 2018).

하지만 대다수 연구는 수목의 간격이 넓은 지역이거나 수목의 간격이 일정하게 배치된 지역을 대상으로 수목 정보 추출 연구가 진행되고 있어 우리나라 산림과 같이 수목의 밀집도가 높은 지역을 대상으로 수행된 연구는 미비한 실정이다. 이에 본 연구에서는 우리나라 산림지역을 대상으로 취득된 항공 LiDAR 데이터와 CHM을 기반으로 local maxima 기법을 적용하여 수목 개체별 수목 정보를 자동으로 추출하는 방법론을 제안하였다.

2. 대상지역 선정 및 자료 수집

2.1. 연구 대상 지역

본 연구에서는 Fig. 1과 같이 침엽수, 활엽수가 우거진 산림지역에서 연구를 수행하고자 치악산 국립공원 부곡지구 내 8.32 km² 면적에 해당하는 지역에 대하여 실험 데이터를 수집하였다.

Fig. 1. Study area.

산림 빅데이터 거래소에서 무료로 공개하고 있는 나무 밀도 지도를 이용하여 치악산 국립공원 부곡지구의 수목 밀도를 분석한 결과는 Fig. 2와 같다. 치악산 국립공원 부곡지고는 산림밀도 높음을 의미하는 C등급의 면적이 3,677,114.16 m², 보통을 의미하는 B등급의 면적이 17,211.30 m²로 나타났으며, 밀도가 낮음을 의미하는 A등급은 포함되지 않는 것으로 나타났다. 치악산 국립공원 부곡지구의 면적 대비 등급별 면적 비율을 계산한 결과 C등급이 66.76%, B등급이 0.31%로 나타났으며, 나머지 32.93%는 비산림지역으로 나타났다.

Fig. 2. Forest density map of study area.

본 연구 대상 지역은 산림밀도가 높은 C등급의 지역 중 30 × 30 m 범위의 국립공원관리 공단 매목조사 3개 조사구(침엽수림 2개, 활엽수림 1개)로 선정하였으며, 각 지구별 특징은 Table 1과 같다.

Table 1. Characteristics of study area

2.2. 항공 LiDAR 데이터 수집

연구에 활용한 항공 LiDAR 데이터는 Leica사의 City Mapper 장비(Table 2)를 활용하여 2020년 4월에 취득되었다. 항공 LiDAR 측량을 통해 수집된 데이터는 지상 해상도(Ground Sample Distance, GSD) 8 cm 급의 항공 정사영상(Fig. 3) 및 점밀도 40 pts/m²의 포인트 클라우드(Fig. 4)를 포함한다.

Fig. 3. Aerial orthophoto.

Table 2. Specification of CityMapper

Fig. 4. Aerial LiDAR point cloud.

2.3. 매목 조사 자료 수집

매목조사는 산림 내 입목의 재적(m³)의 양을 산출하기 위한 현장조사를 의미한다. 국립공원관리 공단에서는 매목조사 시 소형 LiDAR 장비를 이용해 지상에서 LiDAR 포인트 클라우드를 취득하고 수목의 위치(X, Y), DBH, 수고 정보를 산정한다. 본 연구를 위하여 수집된 매목 조사 자료는 국립공원관리 공단에서 지상 LiDAR 데이터 기반으로 산정한 수목의 위치, DBH 정보를 수집하였다. Fig. 5는 수집한 매목조사 데이터를 정사영상에 중첩하고 속성정보에 입력된 DBH 정보를 나타낸 것이다.

Fig. 5. Collected diameter measurement data: (a) Study area 1, (b) Study area 2, and (c) Study area 3.

3. 개체별 수목정보 추출

3.1. 항공 LiDAR 기반 산림지역 개체별 수목정보 추출 방법론

항공 LiDAR 기반 산림지역 개체별 수목정보 추출을 위한 절차는 Fig. 6과 같으며, 수목별 수고 정보 추출을 위한 전처리, local maxima 기반 수목별 수고 정보 추출, 그리고 DBH-수고 산정식 기반의 DBH 산정 세 단계로 구분할 수 있다. 제안한 방법론은 자동화를 위하여 파이썬(python) 기반의 스크립트로 작성하였으며 포인트 클라우드 처리를 위한 laspy 라이브러리와 다양한 영상 처리 기능을 제공하는 scikit-image 라이브러리를 활용하였다.

Fig. 6. Methodology of extract individual tree information in forested areas based on aerial LiDAR.

3.2. 개체별 수고 정보 추출을 위한 전처리

개체별 수고 정보를 추출하기 위해 지표면 표고에 대한 정보를 가진 수치표고모델(Digital Elevation Model, DEM), 지표면 표고와 자연·인공지물의 높이 정보를 가진 수치표면모델(Digital Surface Model, DSM), 그리고 수목의 높이 정보만을 가진 수목높이모델(Canopy Height Model, CHM)이 필요하다(Meng et al., 2010). 수관의 형상을 추출하기 위한 GSD를 고려하여 포인트 클라우드를 기반으로 GSD가 1 m인 DSM을 생성하였다. 또한, DEM 영상을 생성하기 위해 이상점이 제거된 항공 LiDAR 데이터에 Adaptive TIN (ATIN) 필터링 기법 적용하여 지면 점 분류를 수행한 후 GSD가 1 m인 DEM을 생성하였다. CHM은 지형과 수목의 높이를 지닌 DSM에서 지형의 높이를 지닌 DEM을 차분하여 생성하였으며 결과는 Fig. 7과 같다.

Fig. 7. Result of generating (a) DSM, (b) DEM, and (c) CHM.

생성한 DSM, DEM, CHM의 정량적인 분석을 위해 각 영상의 최대 높이와 최소 높이 값을 비교한 결과는 Table 3과 같다. 분석 결과 기존 DEM 및 DSM은 지형의 높이 즉, 표고에 대한 정보를 포함하고 있으나 CHM의 경우 순수한 수목에 대한 높이 값만을 포함하고 잇는 것을 확인할 수 있다.

Table 3. Comparison of height between DSM, DEM, CHM

수집한 LiDAR 데이터와 생성한 CHM 간 수고의 최대-최소 높이를 비교한 결과 3개의 매목지구 모두 높이 차가 10 cm 내로 나타났으며, 데이터 간 차이가 발생한 이유는 수관을 추출하기 위해 생성한 CHM의 GSD가 수집한 LiDAR 데이터의 해상도에 비해 낮아졌기 때문이다. 이를 통해 세 지역 모두 지면의 높이 영향을 제거한 CHM을 통해 순수한 수목의 높이 정보를 추출한 것을 확인할 수 있다.

지면 높이가 수고 추출 결과에 미치는 영향을 최소화하기 위해 CHM을 생성한 것과 동일하게 포인트 클라우드의 정규화를 수행하였다. 기존 포인트 클라우드(Fig. 8a)는 지형의 높이가 반영되어 지표면 표고에 따른 높이 차이가 존재하나, 정규화 된 포인트 클라우드(Fig. 8b)는 지면 높이가 평활화 되어 동일한 지표면 높이를 기준으로 수목의 형상이 나타났다.

Fig. 8. Result of normalization height of point cloud. (a) Before normalization. (b) After normalization.

높이 정규화에 대한 정량적 분석을 위해 Table 4와 같이 포인트 클라우드 높이 값의 최대, 최소, 평균, 표준편차를 계산하였다. 매목지구1 지역의 원본 포인트 클라우드의 최소-최대 높이는 각 562.34 m, 603.13 m로 나타났으나 정규화를 수행한 포인트 클라우드의 최소-최대 높이는 각 0 m, 27.38 m로 나타났다. 매목지구2 지역의 원본 포인트 클라우드의 최소-최대 높이는 각 547.20 m, 582.46 m로 나타났으나 정규화를 수행 이후 최소-최대 높이는 각 0 m, 26.83 m로 나타났으며, 매목지구3 지역의 원본 포인트 클라우드의 최소-최대 높이는 각 535.35 m, 556.63 m로 나타났으나 정규화를 수행 이후 최소-최대 높이는 각 0 m, 24.10 m로 나타났다. 이를 통해 포인트 클라우드를 높이 정규화함으로써 지형 높이의 영향을 제거한 수목의 높이만을 추출한 것을 확인할 수 있다.

Table 4. Quantitative analysis for height normalization

3.3. Local Maxima 기반 수목별 수고 정보 추출

수목별 수고 정보를 추출하기 위해 CHM을 기반으로 local maxima 기법을 적용하여 수목의 정점을 식별한다. Local maxima는 Fig. 9와 같이 탐색창(search window)을 이용하여 일정간격으로 이동하며 최소높이 값(minimum height)보다 크면서 국지 지역 내 가장 높은 극값을 탐색하는 알고리즘을 의미한다(Scholkmann and Boss, 2012;Lee et al., 2021).

Fig. 9. Concept of local maxima.

Local maxima를 수목 추출에 적용하는 경우 한 그루의 수목이 차지하는 영역 내에서 가장 높은 수고의 위치를 탐색할 수 있다는 장점이 있으나 이미지에 잡음이 포함되어 있을 경우 잘못된 위치에 극값이 추출될 수 있다는 단점이 있다. CHM 영상 내 극값을 효율적으로 탐색하기 위해 CHM에 median filter를 적용하여 이미지 내 잡음(noise)을 제거하였다. Median filter를 적용한 CHM영상을 기반으로 scikit-image 라이브러리에서 제공하는 local maxima 기법을 활용하기 위해 Table 5와 같이 산림 지역의 특성을 고려한 최소 수목 높이, 탐색창 크기 등 사용자 임계값(threshold)을 설정하였다.

Table 5. Threshold for local maxima

Local maxima를 적용하여 수목의 정점 위치를 추출하였으며(Fig. 10a), 추출한 수목의 정점 위치와 CHM을 활용하여 watershed 기반의 수관 형상을 검출하였다. 검출한 수관 형상은 높이 값이 정규화된 포인트 클라우드와 결합하여 개별 수목에 대한 포인트 클라우드를 추출하였으며 결과는 Fig. 10(b)와 같다.

Fig. 10. Result of detection of tree top and extraction of point cloud. (a) Detected tree top. (b) Extracted point cloud of individual tree.

3.4. DBH-수고 산정식 기반 DBH 산정

Local maxima 기법을 통해 추출한 수고 값을 기반으로 식(1)과 같은 DBH-수고 산정식을 이용하여 DBH를 역산하였다. 이때, H는 Local maxima 기법을 통해 추출한 수고 값을 의미하며, a와 b는 DBH-수고 간의 관계를 결정하는 계수로서 국가산림자원조사를 통해 실측된 데이터를 기반으로 비선형 회귀식을 통해 결정된 계수(Table 6)를 활용하였다(Jo et al., 2023).

DBH = a · b^H (1)

Table 6. Result of tree extraction accuracy test

4. 개체별 수목정보 추출 기법의 결과 분석

4.1. 수목 검출 정확도 검사

개체별 수목정보 추출 기법의 활용성을 확인하기 위해 수목의 검출 정확도 검사를 수행하였다. 매목 조사를 통해 수집된 수목의 위치는 수고의 위치가 아니므로 제안한 방법론으로 추출한 수고의 위치와는 위치적 차이가 발생한다. 따라서 매목 조사를 통해 수집된 수목의 개체 수와 제안한 방법론을 통해 검출한 수목의 개체 수를 비교하여 수목 검출 정확도 검사를 수행하고자 하였다. 제안한 방법론으로 수목의 개체를 검출한 결과는 Table 7과 같으며 매목지구별 수목 개체 검출률은 각 88.46%, 86.14%, 84.31%로 나타났다.

Table 7. Result of tree extraction accuracy test

4.2. DBH 산정 정확도 검사

DBH 산정 정확도 검사를 위해 매목 조사를 통해 수집된 DBH를 ground truth로 선정하고 제안한 방법론으로 산출한 수고에 DBH-수고 산정식을 적용하여 산출한 DBH를 비교함으로써 수목정보 정확도 검사를 수행하였다. 수목정보 정확도 검사를 위해 제안한 방법론을 통해 자동으로 추출한 수고의 위치와 매목 조사를 통해 수집된 수목의 위치 간 유클리디안 거리(euclidean distance)가 가장 근접한 쌍(pair)을 탐색하고 각 DBH 값을 비교하여 평균오차, 최대/최소오차, 그리고 표준편차를 비교하였다(Table 8).

Table 8. Result of DBH accuracy test

Table 8을 통해 매목지구별 DBH에 대한 평균오차는 각 4.842 cm, 3.827 cm, 1.955 cm로 나타났으며 3개 매목 지구에 대한 최소오차와 최대오차는 각 0.111 cm, 12.480 cm로 나타났다. 최대 오차의 경우 연구 대상 지역의 경계선 부분 즉, 외곽선 경계에서 추출된 수목에서 주로 발생한 것을 확인하였다. 매목조사를 통해 수집된 개체목의 매목지구별 DBH의 분포는 매목지구1의 경우 최소 7.05 cm에서 최대 33.47 cm, 매목지구2는 최소 6.77 cm에서 최대 36.08 cm 그리고 매목지구3은 최소 6.01 cm에서 최대 27.47 cm이며 매목지구별 DBH의 평균은 각 13.21 cm, 12.15 cm, 15.49 cm이다. 매목지구별 DBH에 대한 표준편차는 모두 5 cm 내외로 나타나 제안한 수목 정보 자동 추출 방법론의 활용 가능성을 확인하였다.

5. 결론

본 연구는 LiDAR를 기반으로 DEM, DSM, CHM 영상을 생성하여 local maxima 기법을 통해 산림지역의 수목별 수고 정보를 추출하고 DBH을 산정한 후 매목조사를 통해 수집한 수목별 DBH 정보 결과와 비교하여 다음과 같은 결론을 도출하였다.

첫째, 대규모 산림 지역의 수목별 수고 정보를 추출하기 지상 LiDAR를 활용하는 대신, 항공 LiDAR를 기반의 수목별 수고 정보 추출 방법론을 제시함으로써 좁은 규모에 국한된 현장조사에 비해 넓은 지역의 산림을 한 번에 조사할 수 있었다.

둘째, 항공 LiDAR를 기반으로 DEM, DSM, CHM 영상으로 변환하여 local maxima 기법을 이용하는 방법론을 제시하였으며, 제안한 방법론을 통해 추출한 수목의 검출 정확도는 매목지구별 각 88.46%, 86.14%, 84.31%로 나타났다. 또한, 추출한 수고 정보의 활용성을 확인하기 위해 DBH-수고 산정식을 통해 DBH를 역산하였다. 매목조사를 통해 산정한 DBH 정보와 비교하여 정확도 검사를 수행한 결과 DBH 값에 대한 평균오차는 각 4.842 cm, 3.827 cm, 1.955 cm로 나타났으며 표준편차는 3개 매목지구에서 모두 5 cm 내외로 나타나 제안한 수목 정보 자동 추출 방법론의 활용 가능성을 확인하였다.

셋째, 본 연구에서 사용한 항공 LiDAR 데이터는 국가 기본도를 생성하기 위해 활용하는 CityMapper 장비로 취득되었으며, 국가 기본도 갱신을 위해 전 국토를 대상으로 도시는 1년 주기, 비도시 지역은 2년 주기로 촬영을 수행한다. 따라서 2년 주기로 취득되는 비도시 지역의 항공 LiDAR 데이터를 수집할 경우 산림 지역의 광역적 모니터링이 가능할 것이라 생각된다.

본 연구에서는 3개의 매목지점을 선정하여 국한된 지역에 대해 연구를 수행하였으나, 향후 다양한 유형의 산림에 대한 항공 LiDAR를 수집하고 local maxima 기법에 대한 표준화 연구를 진행한다면 본 연구에서 제시한 방법론의 적용 범위를 확대할 수 있을 것으로 사료된다.