LiDAR의 DBH 추정식, 제대로 논문을 찾아보니 크게 달랐다【임업×AI×GIS #6】

서론

지난 회차(#5)에서는 LiDAR의 수고(tree height) 데이터로부터 탄소량을 계산했다.

계산식의 핵심은 다음과 같다.

수고(LiDAR) → DBH(흉고직경) → 재적 → 탄소량

이때 사용한 DBH 추정식은 DBH = 0.7 × H^1.1이라는 수고만을 이용한 1변수 식이었다.

이 기사에서는 지난번의 식을 개량한 이야기를 쓰고자 한다.

제대로 논문을 찾아 계수를 구현했더니, 탄소 추정치가 4배나 변했다.

문제: 수고만으로 DBH를 추정할 수 있는가?

지난 모델을 되돌아본다.

DBH = 0.7 * H ** 1.1 # 단위: cm, m

이는 "수고가 높은 나무는 굵다"라는 직관을 식으로 만든 것이다.

솔직히 말하자면, 근거가 되는 논문을 확인하지 않고 구현했었다.

실측값과 비교한 연구에 따르면, 수고만을 이용한 모델의 오차는 **±40~50%**에 달하는 경우가 많다.

이 정도로는 탄소 크레딧(carbon credit) 산정에 사용할 수 있는 수준이 아니다.

Watershed로 무엇이 바뀌었나

08_watershed.ipynb의 구현을 통해 각 나무의 **수관 반지름(crown_radius_m)**을 취득할 수 있었다.

944본 × {height_m, crown_area_m2, crown_radius_m}

수고에 더해 수관 크기를 사용할 수 있게 되었다.

같은 수고라도 가지가 크게 뻗어 있는 나무는 줄기도 굵다. 이것을 활용하지 않을 이유가 없다.

논문 찾기: 글로벌 모델 (Jucker et al., 2017)

"tree allometry height crown diameter"로 검색하니 자주 인용되는 논문을 발견했다.

Jucker, T. et al. (2017)

Allometric equations for integrating remote sensing imagery into forest monitoring programmes

Global Change Biology, 23(1), 177–190.

전 세계 10만 8,753본의 데이터로부터 글로벌 추정식을 도출하고 있다.

PubMed Central에서 전문을 무료로 공개하고 있었다.

식의 내용

논문에 게재된 식:

DBH = 0.557 × (H × CD)^0.809 [cm]

H：수고(m)

CD：수관 직경(m) = 2 × 수관 반지름

• RMSE = 9.7cm (전 바이옴 평균)

자체 데이터는 "수관 반지름(Cr)"이므로, CD = 2 × Cr로 치환하여 전개한다.

DBH = 0.557 × (H × 2Cr)^0.809
= 0.557 × 2^0.809 × H^0.809 × Cr^0.809
≈ 0.976 × H^0.809 × Cr^0.809

Python 구현:

JUCKER_A = 0.557 * (2 ** 0.809) # ≈ 0.976
DBH_v2 = JUCKER_A * (H ** 0.809) * (Cr ** 0.809)

지역 특화 모델 찾기: 시마다(2011)

Jucker의 글로벌 모델은 모든 수종을 대상으로 한다.

다마(多摩)의 인공림은 삼나무(Sugi)와 편백나무(Hinoki)가 중심이므로, 일본의 삼나무에 특화된 연구가 없는지 조사했다.

시마다 히로미츠 (島田博匡, 2011)

三重県の高齢人工林における胸高直径・樹高・樹冠幅の関係

三重県林業研報 (3), 19–26.

미에현의 삼나무 고령림(68~157년생, 41 임분)에서의 다중 회귀식(multiple regression equation)을 발견했다.

식의 내용

논문에 기재된 다중 회귀식:

D = 1.3866×Ht + 4.8520×W - 22.6459 [cm]

Ht：수고(m)

W：수관 폭(m) = 2 × 수관 반지름

• R² = 0.808

W = 2 × Cr로 치환:

DBH = 1.3866×H + 4.8520×(2×Cr) - 22.6459
= 1.3866×H + 9.7040×Cr - 22.6459

Python 구현:

DBH_v3_raw = 1.3866 * H + 9.7040 * Cr - 22.6459
DBH_v3 = np.maximum(DBH_v3_raw, 0) # 유령목에서 음수값이 나오므로 클리핑(clip)

이 식은 임분 평균값으로부터 도출된 다중 회귀식이므로, 개별 나무에 적용할 경우 오차가 커질 수 있다.

또한 적용 권장 범위는 80~150년생의 고령 삼나무(Sugi) 숲으로 한정된다.

3개 모델의 비교

실제로 데이터에 적용한 결과는 다음과 같다.

==============================================================
v1(수고만) v2(Jucker) v3(미에현 삼나무)
==============================================================
...

동일한 수고(H) 및 수관 반지름(Cr) 데이터를 넣어도 모델에 따라 크게 다르다.

탄소량에 미치는 영향

나아가 탄소량으로 환산하면 다음과 같이 되었다.

==============================================================
v1(05) v2(Jucker) v3(미에현)
==============================================================
...

v1의 4~5배 차이가 발생했다.

왜 이렇게 차이가 나는가

직관적으로 이해하기 위해 수고별 DBH 추정치를 플롯(plot)해 보았다 (수관 반지름을 중앙값 5.3m로 고정).

h_range = np.linspace(H.min(), H.max(), 200)
cr_fixed = np.median(Cr) # = 5.3m
dbh_v1_line = 0.7 * h_range ** 1.1
...

v1은 "수고 30m에서 약 30cm"라고 추정하고 있다.

v2·v3는 "수고 30m, 수관 반지름 5m라면 50~70cm 이상"이라고 추정하고 있다.

다마(Tama)의 오타케산(Mitake) 지역은 고령 삼나무·편백나무 숲이 많으며,

교목(High tree)이 될수록 줄기도 굵어지는 것이 실제 현상에 가깝다.

수관 반지름을 더함으로써 이 관계가 올바르게 반영되고 있다.

민감도 분석: 지수의 영향을 확인하기

Jucker의 글로벌 계수(지수 0.809)가 얼마나 중요한지 확인했다.

for exp in [0.6, 0.7, 0.809, 0.9, 1.0, 1.1]:
a = 0.557 * (2 ** exp)
dbh_e = a * (H ** exp) * (Cr ** exp)
...

결과:

exponent 평균DBH(cm) 총탄소량(t) vs v1
--------------------------------------------------
0.600 17.5 391.0 -67.7%
...

지수가 0.1 변하면 탄소량이 한 자릿수 차이로 변한다.

계수 선정이 추정 정밀도에 극히 큰 영향을 미친다는 것을 알 수 있다.

향후 과제

이번 2변수 모델로 정밀도는 높아졌겠지만, 과제도 남아 있다.

• 실측 데이터와의 비교 필요

다마 지역의 실측 DBH가 있다면 v2·v3·실측의 세 가지를 비교하여 최적 모델을 선택할 수 있다.
scipy.optimize.curve_fit을 사용하여 다마 고유의 계수를 최적화할 수도 있다.

• v3의 적용 범위

시마다(Shimada, 2011)의 식은 미에현 삼나무 고령림에 특화되어 있다.

유령림(Young forest)이나 침엽수 이외의 수종에 적용하려면 별도의 식이 필요하다.

• 복수 타일로의 전개

이번에는 오타케산의 1개 영역뿐이다.

다마 전역으로 확대하면 더 큰 규모의 탄소량 추정이 가능하다 (다음 회차 예정).

요약

모델	식	근거	탄소량
v1(기존)	0.7 × H^1.1	—	1,212 t
v2(Jucker)	0.976 × H^0.809 × Cr^0.809	Jucker et al. (2017) 전 수종	4,933 t
v3(미에 삼나무)	1.3866×H + 9.704×Cr - 22.65	시마다 (2011) 미에현 삼나무	6,586 t

계수를 "그럴듯한 값"으로 넣었을 때와 논문의 실측치를 사용했을 때, 탄소량이 4배 차이 났다.

DBH는 재적(Volume)·탄소량 계산에서 2제곱·3제곱으로 작용하므로, 추정 정밀도가 최종 결과에 매우 큰 영향을 미친다.

"일단 돌아가는 식"이 아니라 "근거가 있는 식"을 사용하는 것의 중요성을 실감했다.

참고문헌

• Jucker, T. et al. (2017). Allometric equations for integrating remote sensing imagery into forest monitoring programmes. Global Change Biology, 23(1), 177–190. https://doi.org/10.1111/gcb.13388 - 島田博匡 (2011). 三重県の高齢人工林における胸高直径・樹高・樹冠幅の関係. 三重県林業研報(3), 19–26.

데이터 라이선스

본 기사에서 사용한 데이터:

도쿄도 다마·도서 지역 3차원 점군 (Point Cloud) 데이터

출처: 도쿄도 오픈 데이터 카탈로그

라이선스: CC BY 4.0

저자의 다른 실험

LiDAR 점군 (Point Cloud)을 이용한 산림 해석·블록체인×AI 등, 다른 프로젝트도 모아두었습니다.