DOI QR코드

DOI QR Code

A Study for Estimation of High Resolution Temperature Using Satellite Imagery and Machine Learning Models during Heat Waves

위성영상과 머신러닝 모델을 이용한 폭염기간 고해상도 기온 추정 연구

  • Lee, Dalgeun (Researcher Officer, National Disaster Management Research Institute, MOIS) ;
  • Lee, Mi Hee (Senior Researcher, National Disaster Management Research Institute, MOIS) ;
  • Kim, Boeun (Senior Researcher, National Disaster Management Research Institute, MOIS) ;
  • Yu, Jeonghum (Principal Researcher, National Disaster Management Research Institute, MOIS) ;
  • Oh, Yeongju (Researcher, National Disaster Management Research Institute, MOIS) ;
  • Park, Jinyi (Principal Researcher, National Disaster Management Research Institute, MOIS)
  • 이달근 (행정안전부 국립재난안전연구원 공업연구사) ;
  • 이미희 (행정안전부 국립재난안전연구원 선임연구원) ;
  • 김보은 (행정안전부 국립재난안전연구원 선임연구원) ;
  • 유정흠 (행정안전부 국립재난안전연구원 책임연구원) ;
  • 오영주 (행정안전부 국립재난안전연구원 연구원) ;
  • 박진이 (행정안전부 국립재난안전연구원 책임연구원)
  • Received : 2020.10.05
  • Accepted : 2020.10.22
  • Published : 2020.10.31

Abstract

This study investigates the feasibility of three algorithms, K-Nearest Neighbors (K-NN), Random Forest (RF) and Neural Network (NN), for estimating the air temperature of an unobserved area where the weather station is not installed. The satellite image were obtained from Landsat-8 and MODIS Aqua/Terra acquired in 2019, and the meteorological ground weather data were from AWS/ASOS data of Korea Meteorological Administration and Korea Forest Service. In addition, in order to improve the estimation accuracy, a digital surface model, solar radiation, aspect and slope were used. The accuracy assessment of machine learning methods was performed by calculating the statistics of R2 (determination coefficient) and Root Mean Square Error (RMSE) through 10-fold cross-validation and the estimated values were compared for each target area. As a result, the neural network algorithm showed the most stable result among the three algorithms with R2 = 0.805 and RMSE = 0.508. The neural network algorithm was applied to each data set on Landsat imagery scene. It was possible to generate an mean air temperature map from June to September 2019 and confirmed that detailed air temperature information could be estimated. The result is expected to be utilized for national disaster safety management such as heat wave response policies and heat island mitigation research.

본 연구에서는 지상기상센서가 설치되지 않은 미 관측지점의 기온정보를 추정하기 위하여 K-최근접 이웃, 랜덤 포레스트, 신경망 알고리즘을 대상으로 위성영상을 이용하여 기온자료를 산출하고 그 정확성을 평가·분석하고자 하였다. 위성영상자료는 2019년에 취득된 Landsat-8과 MODIS Aqua/Terra을 이용하였으며, 기상자료는 기상청과 산림청의 AWS/ASOS 자료를 이용하였다. 또한 추정 정확도를 향상시키기 위하여 수치표면 모델, 일사량, 경사방향, 경사도를 생성하여 이용하였다. 머신러닝 알고리즘 정확도 비교는 10-fold 교차검증을 통하여 R2(결정계수) 및 RMSE(평균제곱근오차)의 통계량을 계산하여 대상지역별 추정결과를 비교하였다. 그 결과 신경망 알고리즘이 R2=0.805, RMSE=0.508로 세 알고리즘 중 가장 안정적인 결과를 나타내었다. 신경망 알고리즘을 구축된 위성영상 데이터셋에 적용하여 2019년 6월부터 9월까지의 평균기온 지도를 생성할 수 있었으며 세밀한 기온 정보를 관측할 수 있음을 확인하였다. 연구 성과는 폭염 대응 정책, 열섬완화 연구 등 국가재난안전 관리에 활용 될 수 있을 것으로 기대된다.

Keywords

요약

본 연구에서는 지상기상센서가 설치되지 않은 미 관측지점의 기온정보를 추정하기 위하여 K-최근접 이웃, 랜덤 포레스트, 신경망 알고리즘을 대상으로 위성영상을 이용하여 기온자료를 산출하고 그 정확성을 평가·분석하고자 하였다. 위성영상자료는 2019년에 취득된 Landsat-8과 MODIS Aqua/Terra을 이용하였으며, 기상자료는 기상청과 산림청의 AWS/ASOS 자료를 이용하였다. 또한 추정 정확도를 향상시키기 위하여 수치표면 모델, 일사량, 경사방향, 경사도를 생성하여 이용하였다. 머신러닝 알고리즘 정확도 비교는 10-fold 교차검증을 통하여 R2(결정계수) 및 RMSE(평균제곱근오차)의 통계량을 계산하여 대상지역별 추정결과를 비교하였다. 그 결과 신경망 알고리즘이 R2=0.805, RMSE=0.508로 세 알고리즘 중 가장 안정적인 결과를 나타내었다. 신경망알고리즘을 구축된 위성영상 데이터셋에 적용하여 2019년 6월부터 9월까지의 평균기온 지도를 생성할 수 있었으며 세밀한 기온 정보를 관측할 수 있음을 확인하였다. 연구 성과는 폭염 대응 정책, 열섬완화 연구 등 국가 재난안전 관리에 활용 될 수 있을 것으로 기대된다.

1. 서론

전 세계적으로 2019~2020년 호주 산불, 2018년 동아시아 폭염, 2017년 미국 허리케인 하비, 2015년 네팔 지진, 2011년 동일본 대지진 등과 같이 각종 대형재난이 빈번 하게 발생하고 이에 따른 인명 피해와 경제적 손실의 규모가 급격하게 증대되었다. 현대의 재난은 과거와 비교 하여 발생 주기가 단축되었고, 양상 및 피해규모가 복잡화·대형화 되고 있음을 알 수 있다.

우리나라의 경우도 2017년 포항지진, 2018년의 기록적인 폭염, 2019년 이례적인 태풍에 이어 2020년은 유례없는 장마가 지속되는 등 재난의 빈도와 강도가 증가하고 있다(Kim, 2020). 기상청은 최근 한반도의 기온 및 강수 변동성이 전 지구적인 온난화 현상 및 장기적 기후 변동성의 직접적인 영향을 받고 있는 근거로 전 지구 평균 지표온도가 1880~2012년 동안 0.85°C 상승한 반면, 우리나라는 1912~2017년 동안 약 1.8°C 상승하였다고 분석하였다(KMA, 2020). 관련된 사례로 폭염이 심했던 2018년 서울시의 최고기온은 111년의 기상관측 역사상 최고의 기록인 39.6°C이고 열대야 지속일수도 1973년 이래 최장기간인 26일을 기록하였다(Cho and Lee, 2018). 우리나라 대부분의 지역에서 열대야 발생빈도가 증가 현상이 나타나고 있으며, 1990년대 후반 이후에 발생빈도, 강도, 지속 기간이 증가하고 있다(Park and Suh, 2011). 2010년대 이후 가속화 되고 있는 전 지구 평균기온의 상승은 세계적으로 폭염 발생 빈도와 강도를 증가시키고 있으며 이에 따라 우리나라에서도 보건, 환경, 기상 등 다양한 측면에서 폭염에 대한 연구가 활발히 진행되고 있다(KMA, 2020).

폭염과 열대야 변화 경향을 나타내는 중요한 지표인 여름철 평균기온의 변화는 지상에서 관측된 기상자료를 이용하여 위험수준을 판단하고 있다. 하지만 공간적인 불확실성을 포함한 5~11 km 간격의 점적인 기온정보는 인구, 토지이용, 시설자산 등 공간적으로 분포하고 있는 사회 취약성(Vaulerable)을 고려하는데 어려움이 있으며 결과적으로 종합적인 위험성(Risk)을 고려한 대응전략 수립에 한계가 있다(NDMI, 2017). 이에 대한 방법적 대안으로 폭염에 대응하기 위하여 도시 내 사회적 취약성들을 분석하기 위하여 광역적·주기적 관측이 가능한 위성영상을 활용한 연구들이 수행되고 있다. Keramitsoglou et al. (2013)은 NOAA/AVHRR와 MODIS Terra/Aqua의 열적외 영상을 이용하여 대도시 내 공간 분포와 사회 경제적 변수를 이용하여 폭염 위험성을 제한하였으며, Jedlovec et al. (2017)은 MODIS와 VIIRS의 1 km 해상도의 열적외 영상에서 지표온도를 산출하고 해상도를 30 m로 다운스케일링하여 폭염 위험성을 분석 하였다. Yoo et al. (2018)은 MODIS Terra/Aqua의 해상도 1 km 지표온도 자료와 랜덤 포레스트 기법만을 이용하여 도시 내 기온을 추정하는 연구를 수행하였으며, Noi et al. (2017)은 MODIS Terra/Aqua의 지표온도 자료와 여러 개의 선형회귀, 큐비스트회귀, 랜덤포레스트 기법을 이용하여 기온자료를 추정하고 비교하였다. 하지만 위성영상을 이용한 연구들은 1 km의 낮은 해상도와 학술적인 연구를 목적으로 특정지역 및 기간을 대상으로 수행되고 있어 국가적 차원에서의 지속적인 활용에는 많은 제약요소들을 가지고 있다. 폭염과 같이 전국적으로 발생하면서 국소적인 다양한 발현양상의 특징을 분석하기 위해서는 전국을 대상으로 행정단위별로 정밀한열 분포 특성을 파악하여야 한다(NDMI, 2017).

본 연구에서는 기존 연구에서 MODIS Terra/Aqua 위성만을 사용하였을 경우 1 km의 낮은 해상도를 갖는 제약사항을 30 m 해상도를 갖는 Landsat 위성영상을 적용하여 해상도를 향상시켰으며, 기존 도시 지역의 연구에 국한되어 국가 단위의 데이터 산출에 갖는 한계성을 전국 대상으로 관측된 위성영상과 다양한 머신러닝 기법을 이용하여 전국규모 데이터 생산 기법을 제안하였습니다. 폭염기간인 6월부터 9월까지의 전국단위의 30 m 고해상도의 일정한 분포를 갖는 평균기온 정보를 산출하기 위하여 위성영상과 머신러닝 기법을 이용하여 적합한 평균기온 추정 모델을 개발하고자 하였다. 이를 위해 MODIS AQUA/TERRA 위성에 Landsat 위성영상에서 관측된 지표온도를 함께 사용하여 해상도를 향상시켰으며, 기상자료로는 기상청과 산림청에서 관측된 AWS와 ASOS의 기상자료를 이용하였다. K-최근접 이웃, 랜덤 포레스트, 신경망의 3가지 기법을 이용하여 평균기온을 추정한 이후 정확도 검증을 통하여 최적의 알고리즘을 선정하고 공간적으로 분석 할 수 있도록 평균기온 지도를 생성하였다.

2. 연구방법

본 연구는 위성에서 산출할 수 있는 다양한 변수들과 머신러닝 기법을 이용하여 여름철인 6월부터 9월까지의 평균기온을 추정하는 모델을 구축하고자 한다. 연구 방법은 Fig. 1과 같이 크게 3가지로 나눌 수 있다. 첫 번째는 연구대상지역의 위성영상 전처리를 포함한 데이터세트 구축이다. 데이터세트별로 모든 픽셀에 대하여 동일한 조건으로 구성이 되어야하기 때문에 데이터세트는 입력자료 중 가장 국소적인 관측영역을 가지고 있는 Landsat 위성의 영상 크기를 기준으로 하였다. 두 번째는 K-최근접 이웃, 랜덤 포레스트, 신경망 기법을 이용하여 평균기온 추정 모델을 구축하는 것이다. 파이썬의 싸이킷런(Scikit-learn) 라이브러리를 이용하여 머신러닝 모델을 구축하고 결정계수(R2)와 평균제곱근오차(RMSE)를 이용하여 모델별 정확성을 분석하였다. 세 번째는 최적의 모델을 선정하여 평균기온 지도를 생성 하는 것이다. 기상관측지점이 없어 평균기온 측정이 불가능했던 지역을 관측한 위성영상기반의 데이터세트에 머신러닝 모델을 적용하여 전 지역의 평균기온을 산출하고 영상 간 모자이크를 통하여 전국단위의 평균기온 지도를 생성하였다.

OGCSBN_2020_v36n5_4_1179_f0001.png 이미지

Fig. 1. Flowchart of this study.

1) 연구대상지역

Landsat 시리즈 위성영상들은 세계표기시스템(WRS, Worldwide Reference System)을 사용하여 각 영상에 Path와 Row 번호를 통해 촬영 영역을 나타낸다. WRS는 WRS-1, WRS-2로 두 가지가 있으며 Landsat 8호는 WRS-2 표기법을 사용한다. Path는 Landsat 위성이 지나 가는 방향에서 행(동쪽에서 서쪽)을 나타내며 Row는 열(북쪽에서 남쪽)을 의미한다. 본 연구지역은 제주도를 포함한 남한 전체이며, 연구지역을 포함하는 Landsat 8호 영상은 총 10장이며, Path/Row는 Fig. 2와 같다.

OGCSBN_2020_v36n5_4_1179_f0002.png 이미지

Fig. 2. Study areas and the location of weather stations.

2) 연구자료

(1) 지표온도 자료

본 연구에서는 폭염기간인 6월부터 9월까지의 전국 단위의 상세한 평균기온을 산출하기 위하여 위성영상과 머신러닝 기법을 이용하여 적합한 평균기온 추정 모델을 개발하고자 하였다. 이를 위해 Landsat 위성과 MODIS AQUA/TERRA 위성에서 관측된 지표온도를 사용하였다. Landsat 영상의 지표온도 자료는 Landsat-8 Collection 1 Level-2 자료 중 열 밴드를 밝기온도(BT, Brightness Temperature)자료로 변환한 자료를 이용하였다. BT자료에 NDVI범위에 따라 계산된 지표면 방출률을 이용하여 식 (1)과 같이 지표온도로 산출하여 사용하였다(Zhang et al., 2006).

Ts = є1/4T       (1)

여기서, Ts는 지표온도, є은 지표면 방출률, 그리고 T는 밝기온도(BT)를 나타낸다.

Landsat-8 Collection 1 Level-2 자료는 기존 DN자료로 제공되는 Level-1자료에서 사용자의 분석 편의를 위해 대기효과를 보정해준 자료이며(Teixeira Pinto et al., 2020), United States Geological Survey EROS Science Processing Architecture(USGS ESPA)에서 사용자 요청하에 제공 받을 수 있다(https://espa.cr.usgs.gov/). Landsat-8자료는 2019년 6월부터 9월동안 취득된 자료 중 구름의 영향을 적게 받은 자료를 각 Path-row당 1장씩 사용하였으며, 6월부터 9월까지 촬영된 영상 중에 구름의 영향이 많은 경우에는 5월 영상을 사용하였다. 본 연구에 사용된 Landsat-8 자료는 Table 1과 같다. Landsat-8 위성영상은 광학 위성영상 특성상 구름 등 대기의 영향을 많이 받기 때문에 구름 아래의 지표정보가 차단되기 때문에 영상에서 구름의 영향이 있을 경우, Landsat-8 QA밴드를 이용하여 구름 및 구름그림자 영역을 제거하고 제거된 영역에 대해 변형된 Spectral Similarity Group(SSG) 알고리즘을 통해 영상복원을 수행하였다. 변형된 SSG 알고리즘은 Jin et al. (2013)가 제안한 SSG 알고리즘에서 Kim et al. (2014)이 변형한 방법이다. 기존 영상복원에 주로 사용하는 방법은 다른 날 취득한 영상 또는 타 위성 영상으로부터 화소값을 참조하는 것이 일반적이었지만 변형된 SSG 알고리즘은 타 영상자료의 화소값을 사용 하는 것이 아닌 참조영상의 화소값 위치만 참조할 뿐 화소값은 복원하고자 하는 영상에서 참조하는 방법이다. MODIS 자료는 지표온도 산출물(Product)인 MOD11A1, MYD11A1(Terra, Aqua)를 사용하였다. MOD11A1, MYD 11A1 자료는 하루에 2번씩(주간, 야간)으로 1 km 공간 해상도로 지표온도를 산출한 자료이다(Wan, 2006). MODIS 자료는 1 km 해상도의 영상으로서 학습시 30 m 해상도를 갖는 Landsat 위성영상과 동일한 해상도의 자료를 사용하기 위하여 30 m 해상도로 리샘플링 처리를 수행하였다. 1일 주기로 제공되는 지표온도 자료는 2019년 6월부터 9월까지의 평균을 산출하여 대표 평균 기온자료로 사용하였다.

Table 1. Description of variables in this study

OGCSBN_2020_v36n5_4_1179_t0001.png 이미지

(2) 입력변수

① 정규식생지수(NDVI, Normalized Difference Vegetation Index)

정규식생지수는 적색 파장 대역과 근적외선 파장 대역에서의 식물의 반사율 차이를 이용하여 영상 내 식생의 상대적인 활력도를 알 수 있는 지수로 식 (2)와 같다(Rouse et al., 1974). 정규식생지수는 가장 일반적으로 사용되는 식생지수로 -1부터 1까지의 값을 나타낸다. 구름, 눈, 물 등과 같이 수분을 포함하는 경우에는 음수를 나타내며 일반적으로 0.1~ 0.6 사이의 값을 나타내고 식생이 우거질수록 1에 가까운 값을 나타낸다(Jung et al., 2012)

\(N D V I=\frac{\rho_{N I R}-\rho_{R e d}}{\rho_{N I R}+\rho_{R e d}}\)       (2)

여기서, ρNIR은 근적외선 파장 밴드, 적색 파장 밴드를 나타낸다.

본 연구에서는 Landsat-8 Collection 1 Level-2자료에서 대기보정된 표면 반사율(Surface Reflectance)로 변환된 근적외선 파장 밴드와 적색 파장 밴드를 연산하여 NDVI를 산출하였으며, 구름의 영향이 있을 경우 지표 온도자료와 동일하게 변형된 SSG알고리즘을 이용하여 구름 영역을 복원하여 사용하였다.

② 수치표면모델(DSM, Digital Surface Model)

수치표면모델은 인공지물 및 식생 등을 포함한 지구 표면에서부터의 고도를 나타낸 자료이다(Lee and Son, 2016). 수치표면모델은 일본 Japan Aerospace Exploration Agency(JAXA)에서 제공하는 ALOS World 3D(AW3D30)를 이용하였다. AW3D30은 2006년부터 2011년까지 일본 Advanced Land Observing Satellite(ALOS) 위성으로부터 관측된 자료를 통해 제작된 고도 자료이며 해상도는 1 arc(약 30 m)로 무상으로 제공되고 있다(JAXA EORC, 2020).

③ 수치표고모델(DEM, Digital Elevation Model)

대상지역에 대한 총일사량(Solar Area Radiacne), 경사방향(Aspect), 기울기(Slope)를 추출하기 위해 수치표고모델을 사용하였다. 수치표고모델은 Shuttle Radar Topography Mission(SRTM) DEM 자료를 사용하였다. STRM은 미국 국립지리정보국(NGA, National Geospatial)과 미국 항공우주국(NASA, National Aeronautics and Space Administraion)의 협력으로 추진된 전 세계의 수치표고 모델을 제작하는 프로젝트로, 2000년 2월 11일부터 22일까지 11일동안 취득한 자료를 통해 DEM을 제작하여 1 acr(약 30 m)로 제공되고 있다.

④ 총 일사량(Solar Radiance)

총 일사량은 특정 기간 동안 태양으로부터 지표면에 도달하는 태양복사량을 산출한 자료이며 지리적 특성에 따라 얼마나 많은 태양복사량을 받는지에 대한 정보를 나타낸다. 지표면에 도달하는 태양복사에너지인 직사광(Direct Radiation), 반사광(Reflected Radiation), 산란광(Diffuse Radiation) 중 직사광과 산란광을 이용하여 일사량을 산출하며 단위를 가진다(Fu and Rich, 2000; Fu and Rich, 2002). 총 일사량 자료는 DEM자료를 입력자료로, ArcGIS의 Solar Radiance Tool을 이용하여 2019년 폭염대책기간(6월~9월) 동안의 태양고도각에서 받는 일사량을 산출하였다. 태양복사량 자료는 작물 재배, 자원관리, 토목공학 등 다양한 분야에서 최적의 부지를 결정하는데 사용되는 중요한 입력변수가 된다.

⑤ 경사방향(Aspect)

경사방향은 경사면이 어디를 향하고 있는지 나타내는 경사면의 방향을 의미한다. 경사면이 향하고 있는 방향이 북쪽이면 0°이며 북쪽을 기준으로 시계방향으로 동쪽은 90°, 남쪽은 180°, 서쪽은 270°로 0°에서 360°로 표현되며 경사면의 방향이 없는 평평한 지역에서는 -1값을 나타낸다. 경사방향은 ArcGIS Aspect Tool을 이용 하여 산출하였으며, 산출한 Aspect 자료는 분석에 용이 하게 사용하기 위해 Beers et al. (1966)에 의해 제안된 Transformed Aspect로 재산출하여 사용하였다(Yoo et al., 2018). Transformed Aspect는 식 (3)과 같다. Transformed Aspect는 0°부터 360° 값의 범위를 가지던 Aspect를 0부터 2로 나타낸 자료이다.

Transformed Aspect = cos(45° – Aspect) + 1       (3)

⑥ 경사도(Slope)

경사도는 수평을 기준으로 기울기를 나타내며, 일반적으로 고저차를 수평거리로 나누어 백분율(%) 또는 각도(°)로 표현한다. 본 논문에서는 경사도를 각도로 나타냈으며 식 (4)와 같다. 경사 값이 높을수록 지형이 가파르고, 낮은 경사값은 비교적 평평한 지형을 나타내며 0°에서 90° 범위를 가진다(Burrough et al., 2015). 경사도는 ArcGIS의 Slope Tool을 이용하여 산출하였다.

\(\text { Slope }=\tan ^{-1}\left(\frac{v d}{h d}\right)\)      (4)

여기서, vd는 고저차, hd는 수평거리를 나타낸다.

(3) 지상 관측자료

본 연구에서 위성영상을 이용하여 평균기온을 추정하기 위하여 목표값(참값)으로 지상에서 관측된 기상자료를 사용하였다. 기상자료는 기상청과 산림청에서 운영 중인 종관기상관측장비(ASOS, Automated Synoptic Observing System)와 자동기상관측장비(AWS, Automatic Weather System)의 지점별 자동으로 측정된 기온을 연구대상기간 동안 평균 기온 값을 산출하여 사용하였다. 위성영상에서 관측된 영역 내 위치한 지상기상관측지점을 하나의 데이터셋으로 구성하였다(Table 2).

Table 2. Acquisition date and weather stations of satellite scene in study areas

OGCSBN_2020_v36n5_4_1179_t0002.png 이미지

3) 기온 추정방법

(1) K-최근접 이웃(KNN, K-Nearest Neighbor)

K-NN 알고리즘은 값을 구하고자 하는 목표 지점과 가장 근접한 k개의 표본 지점을 그룹화하여 계산을 통해 결과 값을 추정하는 방법이다. 식 (5)과 식 (6)와 같이 목표 지점과 표본지점과의 거리 계산을 통하여 비례한 가중치를 부여하여 계산한다. 가중치는 거리역산에 의해 가중치를 부여하는 방법을 사용한다(Gu et al., 2006; Nilson et al., 2005)

\(\omega\left(p_{i}\right) p=\frac{1}{d_{t, r}} / \sum_{j=1}^{k}\)        (5)

\(\hat{y}_{t}=\sum_{i=1}^{k}\left(w_{t, r} \times y_{r}\right)\)        (6)

식 (5)에서 dt,r은 구하고자 하는 값(t)과 관측값(r)간의 거리를 나타내며, k는 사용되는 데이터 영상 개수를 나타낸다. 각 데이터별 값의 차이는 모두 합산되어 최종값 차이를 계산한다. 두 값이 근접 할수록 보다 큰 가중치를 부여한다. 거리 역산에 의해 가중치(wi,t)의 공식은 식 (5)와 같다(Jung et al., 2010; Gu et al., 2006; Lee, 2005). 여기서 k는 표본 지점의 지정할 개수를 나타낸다. k-NN 알고리즘을 적용하는데 있어서 k 값에 따라 최종 결과 값에 큰 영향을 미친다. 일반적으로 지정한 k 값이 적으면 민감하지만 에러 발생 확률이 높지만, k 값이 늘어날수록 계산량이 증가되고 추정 값이 특성을 따르지 않고 데이터의 평균에 수렴하게 되는 단점이 있기 때문에 최적의 k 값 선정은 연구대상에 맞게 변동 가능하다. 본 연구에서는 k 값을 추정하기 위하여 1에서 50까지 변화시켜 가며 정확도를 비교하였고, 본 연구에 적합한 최적 K 값을 선정하였다. 최대값을 30까지 제한시킨 것은 영상별로 참조점의 수가 30~233개로 많지 않기 때문이다. 영상별로 비교해본 결과, k 값이 증가 할수록 R2 및 RMSE가 정확성이 증가하다가 k=5 이후에는 결과값이 전체 평균값으로 수렴하며 정확성이 감소하는 변화를 보였다. k 값이 증가되면 계산량이 증가되고 결과값이 평균으로 수렴하는 단점이 있으므로 최종적으로 k=5로 선정하여 모델을 구축하였다(Yim et al., 2007; Reese et al., 2002; Katila and Tomppo, 2001)

(2) 랜덤포레스트(RF, Random Forest)

랜덤포레스트는 다수의 의사결정나무를 결합하여 정확도를 높이는 앙상블 기법이다(Lee et al., 2019; Han et al., 2018). 단일 의사결정 나무는 규칙 기반 나무로써 결과를 해석하기 쉽다는 장점이 있지만 학습 데이터에 따라 생성되는 매우 모델이 다르기 때문에 결과값이 과적합되는 경향이 있다(Yoo et al., 2018). 랜덤포레스트는 과적합 문제를 해결하기 위하여 무작위로 다수의 의사결정 나무를 생성하여 포레스트를 구성하는 방법을 사용한다(Breiman, 2001). 각 의사결정 나무는 노드(node)들과 에지(edge)들로 구성되어 있으며, 모든 의사결정 나무들은 독립적으로 훈련 단계를 거친다. 랜덤 포레스트의 최종 결과는 모든 의사결정 나무의 다수결(분류)이나 평균(회귀)을 통하여 생성된다. 이러한 앙상블 과정을 거쳐 단일 의사결정 나무의 과적합을 줄인 결과를 산출한다(Noi et al., 2017; Yoo et al., 2018).

(3) 신경망(NN, Neural Network)

신경망은 인간의 정보처리 패턴과 유사한 방식의 시스템으로 비선형적인 특징을 갖는 문제를 해결할 수 있을 뿐만 아니라 대량의 데이터를 처리 할 수 있는 장점이 있어 보다 일반적이고 보편화된 모형으로 연구되고 있다(Lee et al., 2008). 신경망은 하나의 입력층(input layer)과 하나 이상의 은닉층(hidden layer), 마지막 층인 (output layer)로 구성되어 있다. 출력층을 제외한 모든 층은 뉴런을 포함하며 다음 층과 연결되어 있다. 신경망 학습의 목적은 기댓값(참값)과 모델에서 계산된 출력 값의 차이를 최소화시키는 가중치를 계산하는데 있다. 기댓값과 출력 값을 비교하면서 역전파 알고리즘을 이용한 반복적인 훈련과정을 통해 최적의 가중치와 편차로 갱신한다(Koo et al., 2019; Bishop, 1995). 인공신경망 기법은 이론을 기반으로 한 물리식을 세우기 어려운 문제에 적용하면 기존의 통계적인 방법보다 정확하고 효율적으로 문제를 해결할 수 있다(Gardner and Dorling, 1998).

신경망은 목적에 따라 다양한 종류로 나눌 수 있으며, 본 논문에서는 다층퍼셉트론으로 훈련하는 기법을 이용하여 1개의 층을 이용하였으며 적합한 뉴런의 개수를 선정하기 위하여 위성영상별로 1에서 100까지 뉴런 개수를 변화시켜 정확도를 조사하였습니다. 위성영상별 뉴런 개수가 30개 이후로는 R2 및 RMSE가 일정수준의 정확성에 수렴하는 경향을 나타내었다. 또한, 뉴런 개수가 지나치게 증가하면 계산량을 증가시키고 추정 값의 정확성 향상에 큰 영향을 미치지 않는 단점을 나타내기 때문에 최종적으로 뉴런 개수를 30개로 선택하였다. 또한, 최종 값은 활성함수에 따라 다양한 값으로 추정하게 되는데 입력 값이 증가할수록 출력 값도 증가 하는 ReLU 활성함수를 선택하여 기온 추정 모델을 구축하였다.

3. 실험결과 및 분석

본 연구는 위성영상별 머신러닝기반 평균 기온 추정 모델은 기온 관측이 이루어지고 있지 않은 영역을 위성 영상과 기상자료를 이용하여 평균기온을 추정하는 것이다. 따라서 평균기온과 위성영상 변수들과의 상관관계를 분석하고 연구모델이 기상 관측소가 없는 지점의 평균기온을 추정할 수 있는 머신러닝 모델을 구축하고 정확도를 평가하기 위하여 검증을 수행하였다. Table 2 에 나타난 위성영상별 자료를 기반으로 머신러닝 모델을 구축하였으며, 학습자료를 10개의 그룹으로 나누어 학습 및 검증데이터의 일정비율(90:10)로 10번 반복하여 오류를 평가하는 교차검증을 수행하였다. 검증을 통하여 위성영상별 머신러닝 기법 정확도를 비교·분석하였으며 정확도가 가장 높은 머신러닝 기법을 이용하여 평균기온 지도를 생성하였다.

1) 입력변수 상관관계

본 연구에서는 평균기온에 영향을 미치는 입력변수와의 상관관계를 분석하기 위하여 피어슨 상관관계 분석법을 사용하였다. 평균기온은 기상청의 AWS/ASOS와 산림청의 AWS에서 관측된 자료를 이용하였으며, 위성영상자료는 위성영상별 지표온도, 식생지수, 수치표면모델, 총 일사량, 경사방향, 경사도를 사용하였다. 피어슨 상관계수는 다음 방정식으로 계산 할 수 있으며, 상관계수의 값이 -1과 +1 사이의 값을 가지는데 +1에 가까울수록 양의 선형 상관관계를 나타내고, -1에 가까울수록 음의 선형 상관관계를 나타낸다(Choi et al., 2019).

\(r=\frac{\sum_{i}^{n}\left(x_{i}-\bar{x}\right)\left(y_{i}-\bar{y}\right)}{\sqrt{\sum_{i}^{n}\left(x_{i}-\bar{x}\right)^{2}} \sqrt{\sum_{i}^{n}\left(y_{i}-\bar{y}\right)^{2}}}\)         (7)

식 (7)에서 r은 평균기온에 대한 피어슨 상관계수를 나타낸다. xi는 위성영상 변수이며, yi는 평균기온 데이터이다. n은 데이터 수를 나타낸다.

Table 3은 총 10개의 위성영상별 10개의 변수와 평균 기온 값의 상관계수(r)을 나타낸 표이다. Table 3에서 LSTA-D가 10개 영상 중 5개에서 가장 강한 양의 상관계수 0.67~0.85를 나타내었으며, DSM은 모든 7개 영상에서 가장 강한 음의 상관계수 -0.73~-0.96을 나타내었다. 또한, 각 위성영상별 변수와 평균기온의 상관관계가 다르게 나타난다. 지표온도는 모두 양의 상관계수를 보이고 있는 반면 총 일사량은 위성영상별로 양과 음의 상관계수를 보이는 특징을 보여주고 있다. Fig. 3는 변수와 평균기온과의 상관계수를 바 그래프 형태로 표현한 것이다. 본 연구에서는 다양한 변수의 상관성이 최종 결과 값에 미치는 영향성 분석을 고려하여 선행연구로 선정된 변수를 포함하여 모델을 구축하였습니다.

Table 3. Correlation coefficients between the ten variables and Tmean in satellite images

OGCSBN_2020_v36n5_4_1179_t0003.png 이미지

OGCSBN_2020_v36n5_4_1179_f0003.png 이미지

OGCSBN_2020_v36n5_4_1179_f0004.png 이미지

Fig. 3. Bar plots between the ten variables and Tmean in satellite images.

2) 머신러닝 정확도 평가

평균기온 추정 모델의 성능을 분석하기 위해서 본 연구에서는 10-fold 교차 검증을 수행하였다(Table 4, Fig. 4). 10-fold 교차검증은 평가하고자 하는 모델의 성능 측정을 위하여 전체 샘플을 10등분하여 9개의 부분 샘플이 학습 자료인 모델을 구축하고, 나머지 1개의 부분 샘플을 검증 자료로 이용하는 방법이다(Han et al., 2018). 반복적으로 10개의 부분 샘플을 검증 자료로 교체하면서 모델을 평가하기 때문에 모델의 정확도를 추정하는데 다양한 분야에 이용되고 있으며 입력변수는 표준화(Z-score)된 표준점수를 사용하였다(Cao et al., 2016). 본 연구에서는 머신러닝 기법의 적합성을 판단하기 위해서 다음과 같이 두 변수의 선형상관관계에 대한 척도를 측정하기 위한 R2(결정계수, 식 (8))과 추정오차를 계산하기 위한 RMSE(평균제곱근오차, 식 (9))를 계산하였다(Wu et al., 2016).

\(R^{2}=\frac{\sum_{j=1}^{n}\left(\hat{y}_{i}-\bar{y}\right)^{2}}{\sum_{j=1}^{n}\left(y_{i}-\bar{y}\right)^{2}}\)        (8)

\(R M S E=\sqrt{\frac{1}{n} \sum_{j=1}^{n}\left(\hat{y}_{i}-y_{i}\right)^{2}}\)          (9)

Table 4. 10 (Ten)-fold cross-validation results (R2 and RMSE) for satellite images

OGCSBN_2020_v36n5_4_1179_t0004.png 이미지

OGCSBN_2020_v36n5_4_1179_f0005.png 이미지

OGCSBN_2020_v36n5_4_1179_f0006.png 이미지

OGCSBN_2020_v36n5_4_1179_f0007.png 이미지

Fig. 4. Scatter plots between estimated and observed air temperatures from 10(Ten)-fold cross-validation results (R2 and RMSE) for satellite images.

여기서 \(\hat y\)은 i번째 추정 값, yi은 실제 관측 값, \(\bar y\)는 실제 관측 값의 평균, n은 전체 표본 개수를 나타낸다.

Table 4는 Landsat 위성영상 관측 크기 기준으로 총 10개 지역을 머신러닝기법을 적용하여 평가된 결과 값이다. 결정계수 R2은 신경망 기법이 10개 지역 중 7지역에서 R2(0.733 ~ 0.938)이 가장 높게 나타났으며, 랜덤 포레스트 기법은 1개 지역(R2 =0.879), K-최근접 이웃 기법도 1개 지역(R=0.712)에서 가장 높게 나타났다. 또한 각 기법의 평균을 기준으로는 정확도가 신경망(R2mean=0.805) > 랜덤 포레스트(R2mean=0.789) > K-최근접 이웃(R2mean=0.745) 기법 순으로 나타났다. 평균제곱근오차RMSE는 신경망기법이 10개 지역 중 8개 지역에서 RMSE(0.444~0.586)가 가장 낮게 나타났으며, K-최근접 이웃은 1개 지역 RMSE(0.472)로 가장 낮게 나타났다. 그리고 머신러닝 기법별 평균을 보면 오차가 신경망(RMSE=0.508) < 랜덤 포레스트(RMSE=0.551) < K-최근접 이웃(RMSE=0.630) 기법 순으로 나타났다. R2과 RMSE를 확인해보면 모두 신경망이 높은 정확도를 보이는 것을 확인할 수 있다.

Fig. 4는 위성영상 별 머신러닝 기법을 적용하여 산출 된 결과물의 산포도로서 파란점은 평균기온 값, 회색선은 기준선, 빨간선은 회귀선을 나타낸 그림이다. 평균기온 추정 값과 관측 값이 모두 선형적인 특징을 나타내는 것을 확인 할 수 있다. 그러나 10개의 영상 중 9개 영상에서 모든 알고리즘에서 R2이 0.7 이상을 나타나지만 1개(116-036)의 영상에서 0.3~0.5 정도의 값을 나타내어 현저히 떨어진 성능을 확인할 수 있었다. 이는 해당 영상 내의 기상자료 관측점의 위치, 입력변수에서 나타내는 이상 특징 등을 면밀히 파악하여 원인을 분석할 필요가 있다고 판단된다.

3) 평균기온 지도 생성

본 연구에서 구축된 머신러닝 모델 중 R2, RMSE의 정확도 평가에서 적합성이 가장 높게 나온 신경망 기법을 이용하여 연구대상지역에 적용한 후 영상 간 모자이크를 통하여 평균기온 지도를 생성하였다. 생성된 지도는 가시성을 확보하기 위하여 0 ∼ 29°C 간 색상을 반영 하였다. 연구대상지역 토지피복도(Fig. 2)와 평균기온 지도(Fig. 5)를 비교하면 도시지역, 농림지역에서 주로 상위 온도 분포를 가지고 있는 것으로 보이며, 산림지역에서 하위의 온도 분포가 나타나는 것을 확인할 수 있다. 이는 토지피복이나 고도가 평균기온 값을 결정하는데 중요한 변수라고 보여 질 수 있지만 보다 정확한 정량적인 분석을 위해서는 머신러닝 모델에 입력되는 변수와 결과의 상관성 분석 등을 통한 추가적인 연구가 필요하다고 판단된다.

OGCSBN_2020_v36n5_4_1179_f0008.png 이미지

Fig. 5. The mosaic map of mean air temperature images calculated through machine learning.

4. 결론

본 연구에서는 동일 입력 변수를 이용하여 평균기온을 추정할 수 있는 머신러닝 모델을 구축하여 정확도를 비교·분석하였다. K-최근접 이웃, 랜덤 포레스트, 신경망 3개의 기법을 이용하여 모델을 개발하고 평가한 결과 신경망 기법이 K-최근접 이웃이나 랜덤 포레스트 기법의 모델보다 모델 정확도(R2)가 높았으며, 오차를 나타내는 RMSE가 낮은 것으로 나타나 모델의 예측력이 우수한 것으로 분석되었다. 이는 10개의 영상에서 1개 지역을 제외한 9개 지역에서 동일하게 나타났으며 신경망을 최적의 모델로 선정하고 대상지역에 해당되는 위성영상에 적용을 통하여 전국단위의 평균기온 지도를 제작하였다. 이러한 연구방법 및 결과를 이용하여 목적에 맞게 입력변수를 조정하여 모델을 구축한다면 폭염의 중요한 요소인 최고기온, 최저기온 등 다양한 시기 및 기온 정보를 산출할 수 있을 것으로 판단된다.

하지만 본 연구에서는 훈련 표본 수로 사용된 기상관측자료는 주로 시가지, 초지, 산림 지역에 위치한 자료의 제한된 결과로 보다 다양한 토지피복을 가지고 있는 지표에 적용하기 위한 머신러닝 기법의 효용성을 결정하기에는 다소 무리가 있을 것으로 판단된다. 따라서 향후 연구에서는 다양한 토지피복에서 관측된 기상자료를 대상으로 한 추가 검증이 필요한 것으로 판단된다. 본 연구의 경우 3개의 머신러닝 기법의 정확도 차이를 분석하는데 초점을 두고 진행하였기 때문에 이번 연구 결과물을 기반으로 향후 추정된 평균기온의 분포나 도시지역의 미세한 변화 패턴을 공간적인 관점에서 분석 한다면 국가나 지자체의 폭염 대응 정책을 수립하는데 과학적인 근거로서 활용 할 수 있을 것으로 기대한다.

사사

본 연구는 행정안전부 국립재난안전연구원의 주요 사업(위성자료 활용 현업지원 기술개발(NDMI-주요-2020-03-01))의 지원으로 수행되었습니다.

References

  1. Beers, T. W., P. E. Dress, and L. C. Wensel, 1966. Notes and observations : aspect transformation in site productivity research, Journal of Forest, 64(10): 691-692.
  2. Bishop, C. M., 1995. Neural Networks for Pattern Recognition, Oxford University Press, New York, p. 482.
  3. Breiman, L., 2001. Random Forests, Machine Learning, 45(1): 5-32. https://doi.org/10.1023/A:1010933404324
  4. Burrough, P. A., R. McDonnell, R. A. McDonnell, and L. D. Lloyd, 2015. Principles of Geographical Information Systems, Oxford University Press, NY, USA, p. 190.
  5. Cao, X. H., I. Stojkovic, and Z. Obradovic, 2016. A robust data scaling algorithm to improve classification accuracies in biomedical data, BMC Bioinformatics, 17(1): 359. https://doi.org/10.1186/s12859-016-1236-x
  6. Cho, H. M. and Y. H. Lee, 2018. The Improvement plan for Seoul heat response, Policy Report, 257: 4.
  7. Choi, M. H., N. J. Jung, K. C. Lee, J. S. Jeong, and I. Y. Seo, 2019. Development of artificial neural network algorithm for the prediction of power failures by natural disaster, Transactions of the Korean Institute of Electrical Engineers, 68(9): 1085-1093. https://doi.org/10.5370/KIEE.2019.68.9.1085
  8. Fu, P. and P. M. Rich, 2000. The solar analyst 1.0 user manual, Helios Environmental Modeling Institute (HEMI), USA, p. 1616.
  9. Fu, P. and P. M. Rich, 2002. A geometric solar radiation model with applications in agriculture and forestry, Computers and Electronics in Agriculture, 37(1-3): 25-35. https://doi.org/10.1016/S0168-1699(02)00115-1
  10. Gardner, M. W. and S. R. Dorling, 1998. Artificial neural networks (the multilayer perceptron) - A review of applications in the atmospheric sciences, Atmospheric Environments, 32(14-15): 2627-2636. https://doi.org/10.1016/S1352-2310(97)00447-0
  11. Gu, H., L. Dai, G. Wu, D. Xu, S. Wang, and H. Wang, 2006. Estimation of forest volumes by integrating Landsat TM imagery and forest inventory data, Science in China Series E: Technological Sciences, 49(1): 54-62. https://doi.org/10.1007/s11431-006-8107-z
  12. Han, D., Y. J. Kim, J. Im, S. Lee, Y. Lee, and H. Kim, 2018. The estimation of arctic air temperature in summer based on Machine Learning approaches using IABP Buoy and AMSR2 satellite data, Korean Journal of Remote Sensing, 34(6-2): 1261-1272. https://doi.org/10.7780/KJRS.2018.34.6.2.10
  13. JAXA EORC (Japan Aerospace Exploration Agency Earth Observation Research Center), May 2020. ALOS Global Digital Surface Model (DSM) ALOS World 3D-30 m (AW3D30), Product Description Ver.3.1.
  14. Jedlovec, G., D. Crane, and D. Quattrochi, 2017. Urban heat wave hazard and risk assessment, Results in Physics, 7: 4249-4295
  15. Jin, S., C. Homer, L. Yang, G. Xian, J. Fry, P. Danielson, and P. A. Townsend, 2013. Automated cloud and shadow detection and filling using two-date Landsat imagery in the USA, International Journal of Remote Sensing, 34(5): 1540-1560. https://doi.org/10.1080/01431161.2012.720045
  16. Jung, J. H., J. Heo, S. H. Yoo, K. M. Kim, and J. B. Lee, 2010. Estimation of aboveground biomass carbon stock in Danyang area using kNN algorithm and Landsat TM seasonal satellite images, Journal of Korean Society for Geospatial Information Science, 18(4): 119-129.
  17. Jung, M. H., S. H. Lee, E, M. Chang, S. W. Hong, 2012, Method of monitoring forest vegetation change based on change of MODIS NDVI time series pattern, Journal of Korea Spatial Information Society, 20(4): 47-55. https://doi.org/10.12672/ksis.2012.20.4.047
  18. Katila, M. and E. Tomppo. 2001. Selecting estimation paramters for the finish multisource national forest inventory, Remote Sensing of Environment, 767: 16-32. https://doi.org/10.1016/S0034-4257(00)00188-7
  19. Keramitsoglou, I., C. T. Kiranoudis, B. Maiheu, K. Ridder, I. A. Daglis, P. Manunta, and M. Paganimi, 2013. Heat wave hazard classification and risk assessment using artificial intelligence fuzzy logic, Environmental Monitoring and Assessment, 185: 8239-8258. https://doi.org/10.1007/s10661-013-3170-y
  20. Kim, B. H., Y. H. Kim, Y. K. Han, W. S. Choi, and Y. I. Kim, 2014. Fully automated generation of cloud-free imagery using Landsat-8, Journal of the Korean Society of Surveying, Geodesy, Photogrammetry and Cartography, 32(2): 133-142 (in Korean with English abstract). https://doi.org/10.7848/ksgpc.2014.32.2.133
  21. Kim, M., 2020. The hazard Viz-platform for the establishment of heatwave response strategies, Journal of Korea Multimedia Society, 23(5): 683-699. https://doi.org/10.9717/KMMS.2020.23.5.683
  22. KMA (Korea Meteorological Administration), 2020. Korean Climate Change Assessment Report 2020, Korea Meteorological Administration, KR, pp. 307-328.
  23. Koo, Y. H., S. M. Kim, M. Oh, and H. D. Park, 2019. Estimation of solar irradiance at weather stations in Korea using regionally trained artificial neural network models, Journal of the Korean Society of Mineral and Energy Resources Engineers, 56(2): 155-171. https://doi.org/10.32390/ksmer.2019.56.2.155
  24. Lee, C. H., 2005. Calculating attribute weights in k-nearest neighbor algorithms using information theory, The Korean Institute of Information Scientists and Engineers, 32(9): 920-926.
  25. Lee, E. J., C. H. Min, and T. S. Kim, 2008. Development of the KOSPI (Korea Composite Stock Price Index) forecast model using neural network and statistical methods, The Institute of Electronics and Information Engineers, 45(5): 95-101.
  26. Lee, G. W. and S. W. Son, 2016. Geo-Spatial Information System, Seoul, Goomibook.
  27. Lee, M., Y. Kim, Y. Jun, and Y. Shin, 2019. Random forest based prediction of road surface condigion using spatio-temporal features, Journal of Korean Society of Transportation, 37(4): 338-349 (in Korean with English abstract). . https://doi.org/10.7470/jkst.2019.37.4.338
  28. NDMI (National Disaster Management Research Institute), 2017. Development of Analysis Technique Using Land Surface Temperature Based on Satellite Data, Research Report, National Disaster Management Research Institute, KR.
  29. Nilson, M., J. Bohlin, H. Olsson, S. A. Svensson, and M. Haapaniemi, 2005. Operational use of remote sensing for regional level assessment of forest estate values, New Strategies for European Remote Sensing, 24: 263-268.
  30. Noi, P. T., J. Degener, and M. Kappas, 2017. Comparison of multiple linear regression, cubist regression, and random forest algorithms to estimate daily air surface temperature from dynamic combinations of MODIS LST data, Remote Sensing, 9(5): 398. https://doi.org/10.3390/rs9050398
  31. Park, W. S. and M. S. Suh, 2011. Characteristics and trends of tropical night occurrence in South Korea for recent 50 years (1958-2007), Atmosphere, 21(4): 361-371. https://doi.org/10.14191/ATMOS.2011.21.4.361
  32. Reese, H., M. Nilson, P. Sandstrm, and H. Olsson. 2002. Applications using estimates of forest paramters derived from satellite and forest inventory data, Computers and Electronics in Agriculture, 37(1): 37-55. https://doi.org/10.1016/S0168-1699(02)00118-7
  33. Rouse, J. W., R. H. Haas, J. A. Schell and D. W. Deeringm, 1974. Monitoring vegetation systems in the Great Plains with ERTS, Proc. of 1974 3rd Earth Resource Techonology Satellite (ERTS) Symposium, Washington D.C, USA, Dec. 10-14, pp. 309-317.
  34. Teixeira Pinto, C., X. Jing, and L. Leigh, 2020. Evaluation analysis of Landsat level-1 and level-2 data products using In situ measurements, Remote Sensing, 12(16): 2597. https://doi.org/10.3390/rs12162597
  35. Wan, Z., 2006. MODIS land surface temperature products users' guide, Institute for Computational Earth System Science, University of California: Santa Barbara, CA, USA.
  36. Wu, C., H. Shen, A. Shen, J. Deng, M. Gan, J. Zhu, H. Xu, and K. Wang, 2016. Comparison of machine-learning methods for above-ground biomass estimation based on Landsat imagery, Journal of Applied Remote Sensing, 10(3): 03510.
  37. Yim, J. S., G. S. Kong, S. H. Kim, and M. Y. Shin, 2007. Forest thematic maps and forest statistics using the k-nearest neighbor technique for Pyeongchang-gun, Gangwon-do, Journal of Korean Society of Forest Science, 96(3): 259-268 (in Korean with English abstract).
  38. Yoo, C., J. Im, S. Park, and L. J. Quackenbush, 2018. Estimation of daily maximum and minimum air temperatures in urban landscapes using MODIS time series satellite data, ISPRS Journal of Photogrammetry and Remote Sensing, 137: 149-162. https://doi.org/10.1016/j.isprsjprs.2018.01.018
  39. Zhang, J., Y. Wang, and Y. Li, 2006. A C++ program for retrieving land surface temperature from the data of Landsat TM/ETM+ band6, Computers & Geosciences, 32(10): 1796-1805. https://doi.org/10.1016/j.cageo.2006.05.001