1. 서  론

지표면에서 발생하는 증발과 식물을 통한 증산을 포함한 증발산은 범지구적 수문순환을 파악하기 위한 중요한 변수이다. 증발산은 기상(기온, 상대습도, 풍속, 기압 등) 상태와 지표 환경(지표온도, 토양수분, 지표 반사도 등)에 영향을 받는다. 또한 증발산에 직접적인 관계가 있는 순복사량은 대기와 지표상태의 따른 단파, 장파 복사의 유출입 특성에 의해 결정된다. 이처럼 증발산은 대기와 주변 환경의 다양한 특성에 영향을 받기 때문에 지역의 에너지 수지를 알 수 있는 종합적인 지표이며 지역 환경을 이해하는 데에 기초가 된다(Sellers et al., 1996; Bonan, 2008; Im, 2013).

정량적인 증발산량을 측정하기 위해 Bowen Ratio Energy Balance Systems (BREBS), Eddy Correlation Systems (ECS), Lysimeter 등 물수지 방법과 에너지 수지 방법을 이용한 다양한 방법들이 개발되어 제시되었다. 물수지 방법은 강우량과 증발산량, 유출 및 침투 등 수문요소들의 이동에 따른 균형을 이용하여 증발산을 산정하며, 년 단위의 장기간의 증발산 산정에 적합하다. 반면 에너지 수지 방법은 유입 및 유출되는 단파, 장파 복사량의 균형을 이용해 증발산을 산정하는 방법으로 단기간의 증발산 산정에도 적합하다(Sur et al., 2012).

에너지 수지 방법을 이용해 증발산을 산정하는 대표적인 방법은 플럭스 타워(Flux tower)를 통한 관측과 인공위성 자료를 기반으로 한 원격탐사 산정 방법이 있다(Jeong et al., 2009). 플럭스 타워는 대기와 지표와의 에너지 및 물질 교환을 직접적으로 측정하여 증발산을 산정하며(Lee et al., 2004), 국내 관측지로는 경기도 포천시 광릉 국립수목원내 활엽수림/침엽수림, 경기도 설마천 혼효림, 경기도 청미천 논지 등이 있다. 플럭스 타워는 실시간 자료 수집, 전송을 통해 즉각적인 자료를 얻을 수 있는 장점이 있으나, 자료 처리를 위한 전문적인 지식이 필요하며 점(point) 추정 기법의 특성상 넓은 지역의 지형, 기상학적 특성을 반영할 수 없는 한계가 있다. 반면 원격탐사 산정 방법의 경우 관측 장비를 통해 넓은 지역의 지형적 특성 및 공간 분포를 반영할 수 있는 장점이 있다(Choi et al., 2009; Sur et al., 2012).

20세기 후반부터 원격탐사기법이 향상되면서 증발산량을 추정하기 위한 다양한 방법들이 연구되었다(Menenti and Choudhury, 1993; Norman et al., 1995; Anderson et al., 1997; Bastiaanssen et al., 1998; Allen et al., 2007). Bastiaanssen et al. (1998)이 제안한 Surface Energy Balance Algorithm for Land (SEBAL) 모형은 유역 내 지표상태가 불균일한 지역에서 최소한의 지상 관측 자료를 이용해 증발산량을 산정할 수 있어 유용성이 높고(Bastiaanssen, 2000), Lysimeter를 이용한 실측 증발산량과의 상관성이 높아 정확성 측면에서도 높은 평가를 받고 있다(Gowda et al., 2008). 또한 미국의 Idaho 주립대학과 Idaho Department of Water Resources에서 NASA와 기업의 지원을 받아 활발히 연구되고 있으며, 모형의 세부 알고리즘이 전부 공개되어 미국을 비롯한 터키, 프랑스, 스페인 등 30개국이 넘는 나라에서 서로 다른 시공간해상도에 대해 적용하는 연구가 활발히 진행되고 있다(Bastiaanssen et al., 2005; Singh et al., 2008; Teixeira et al., 2009). 본래 SEBAL 모형은 30 m 공간해상도와 16 일의 시간해상도를 가진 고해상도의 Landsat 위성영상을 이용해 지역적인 범위에 적용하기 위해 개발되었으나(Allen et al., 2002), 하루 혹은 이틀에 걸쳐 250, 500, 1000 m의 해상도로 전 지구를 관측하는 MODIS 위성과 같은 새로운 지구 관측 체계가 개발됨에 따라 중저해상도 위성을 SEBAL 모형에 적용하여 증발산량을 산정하는 연구가 다양하게 수행되고 있다(Nishida et al., 2003; Ahmad et al., 2006; Tasumi et al., 2008; Hong et al., 2009; Long et al., 2011; Ruhoff et al., 2012; Tang et al., 2013).

국내에서는 SEBAL 모형을 이용해 금강 상류의 보정천 유역(Yoo, 2003), 경기도의 경안천 유역(Ha et al., 2010), 왕숙천 유역(Im, 2013)에 대해 증발산량을 산정하고 플럭스 타워 관측 자료나 다른 모형과 비교하여 적용성을 평가한 연구가 진행된 바 있다. Lee et al. (2015)은 SEBAL 모형을 MODIS 위성영상에 맞게 알고리즘을 일부 수정하고 혼효림 지역인 용담댐 유역에 적용하여 시계열 공간증발산량을 산정하였다. 특히 위성 자료인 Albedo 및 NDVI에 보정계수를 도입하여 모의 증발산량의 보정을 시도하였으나, 에너지 수지 방정식의 각 항목의 시계열 분포가 연속적이지 않고 산발적으로 분포가 되는 문제가 발생하였다. 따라서 본 연구에서는 선행 연구(Lee et al., 2015)에서 적용한 월별 보정계수를 활용하지 않고, 월별 입력 자료를 통해 모형의 효율성을 개선한 뒤 국내의 플럭스 타워를 설치한 다른 유역에 적용하였다. 모의기간은 선행 연구와 동일하게 2012년부터 2013년까지로 하였으며, 산정된 모의 증발산량은 플럭스 타워 실측 자료를 통해 평가하고 이를 통해 모형의 효용성을 논해보고자 한다.

2. SEBAL 이론

SEBAL 모형은 에너지 수지 방정식을 기반으로 각 픽셀별로 계산하며, 에너지 수지 방정식의 각 항목인 순복사량(net radiation, ), 토양 열 플럭스(soil heat flux, ), 현열 플럭스(sensible heat flux, )를 개별적으로 계산한 뒤 총합하여 잠열 플럭스(latent heat flux, )를 계산한다. 증발산은 잠열 플럭스를 이용하여 계산하며, 위성영상 획득 순간의 증발산량을 계산한 뒤 일 증발산량으로 환산한다(Morse et al., 2000).

                          (1)

여기서, 은 순복사량(), 는 토양 열 플럭스(), 는 현열 플럭스()이며, 는 잠열 플럭스()이다. 은 지표상의 총 에너지 수지를 의미하며, 에너지의 총 유입량에서 유출량을 제하여 구할 수 있다(Eq. (2)).

  (2)

여기서, 은 지표로 유입되는 하향 단파 복사 에너지(), 는 지표에서 방출되는 상향 단파 복사 에너지(), 는 하향 장파 복사 에너지(), 는 상향 장파 복사 에너지(), 은 지표 반사율(albedo), 은 대기의 투과율(atmospheric transmissivity), 는 Stefan-Bolzmann 상수(), 는 기준점에서의 지표 온도(K), 는 지표 복사율(surface emissivity), 는 지표 온도(K)를 의미한다.

토양 열 플럭스는 전도현상(conduction)으로 인하여 토양 또는 식생에 저장되는 열의 양을 의미하며, NDVI와 지표 온도 및 반사율로부터 경험적으로 얻어진 Eq. (3)을 이용해 산출한다(Bastiaanssen, 2000).

 (3)

여기서, 는 지표 온도(°C), 는 지표 반사율이다.

현열 플럭스()는 지표면에서의 풍속, 지표 거칠기(surface roughness) 그리고 지표와 대기의 온도차를 이용해 계산한다.

                               (4)

여기서, 는 공기의 밀도(), 는 공기의 정압비열(), 는 지표와 대기의 온도차(K), 는 열 이동에 대한 공기역학적 저항()을 나타낸다. 지표와 대기의 온도차는 Anchor Pixels (Hot/Cold pixel)을 이용해 계산하며 이 기준점을 결정하는 것이 SEBAL 모형에서 증발산량을 산정하기 위한 가장 중요한 과정이다. Cold pixel은 작물로 덮인 잘 관개된 지역이며, 이 지점에서의 지표 온도는 공기 온도에 근사하다. Hot pixel은 잠열 플럭스가 0으로 추정되는 건조한 지역을 나타낸다. SEBAL 표준 모형에서는 사용자가 임의로 Anchor pixel을 산정하는데, 이는 연구자의 숙련도에 따라 기준점의 적정성이 결정되며 다수의 위성 자료를 활용할 경우 전체적인 결과물의 신뢰도가 낮아질 수 있다. 따라서 본 연구에서는 기준을 통해 Anchor pixel을 일관되게 선정할 수 있도록 자동화 기법을 도입하였다. Morse et al. (2000)에서 제시한 지표의 특성을 고려하기 위해 토지 이용도를 도입한 뒤, Cold pixel은 농업 지역(agricultural area) 중 NDVI의 상위 5%에 해당하는 pixel 중에서 지표 온도가 하위 15%이내에 속한 pixel을 선정하도록 하였다. Hot pixel은 나지나 시가지 중에서 NDVI 값이 하위 10%에 해당하고 지표 온도가 상위 15%에 속하는 지점을 선정하였다.  또한 지표의 가열로 인한 대기의 부양효과를 고려하기 위해 반복 순환 계산을 하여 오차를 최소화 한 뒤 현열을 산정하는데, 본 연구에서는 일별로 약 7 ~ 9회 순환 계산을 통해 일정한 현열 값으로 수렴하는 것으로 나타났다. 최종적으로 산정한 현열 플럭스와 순 복사량, 토양 열 플럭스를 에너지 수지 방정식(Eq. (1))에 대입하여 증발산량을 산정한다(Fig. 1).

3. 대상 유역

본 연구는 플럭스 타워가 설치된 설마천 유역과 청미천 유역을 대상유역으로 선정하였다(Fig. 2). 설마천 플럭스 타워(Seolmacheon mixed forest flux tower, SMK)는 경기도 파주시 적성면의 마지리와 설마리에 위치(북위 37°56′20″, 동경 126°57′17″)하고 있으며, 해발고도는 293 m로 설마천 유역의 중류부에 위치한다(Fig. 2a). 설마천은 북위 37°55′06″, 동경 126°57′25″에 위치하며, 임진강 하구에서 약 4.6 km 떨어진 임진강의 제 1지류로 유역면적은 8.5km² 유로연장은 5.8 km, 유로경사 2%, 해발고도 420 m의 전형적인 산지 하천이다(Korea Institute of Construction Technology, 2006). 이 유역은 90% 이상이 주로 20~40년 수령의 침엽수와 활엽수로 구성되어 있으며, 얕은 표토로 인해 수분 함양능력이 낮고 산지 사면에는 돌과 자갈이 산재하여 집중 호우시 토석류를 유발할 수 있다. 이 유역은 프런티어 사업과 관련하여 수문 관측이 활발하게 이루어지고 있다(Kwon et al., 2009).

청미천 플럭스 타워(Cheongmicheon farmland flux tower, CFK)는 경기도 여주군 점동면 농업기술원의 종자관리소 논지(북위 37°09′35″, 동경 127°39′10″)에 위치하며, 지점의 해발고도는 141 m 이다(Fig. 2b). 청미천은 북위 37°02′09″, 동경 127°39′26″에 위치하며, 여주군에서 남한강으로 유입하는 지방 2급하천으로 유역면적 569.6 km², 유로연장 60.8 km, 유역 평균경사는 17.01%이다. 유역의 식생 상태는 하천과 산악사이의 농경지로 이루어져 있으며, 산림이 48.5%, 농경지가 43.1%, 기타 지역(도심지, 도로 등)이 8.4%를 차지하고 있다(MLTM, 2008).

본 연구의 대상 유역은 모두 유역 내에 기상관측소가 위치하고 있지 않기 때문에 유역 인근의 관측소로부터 자료를 제공받아 모형의 기상입력 자료로 활용하였다. 설마천 유역은 인근에 위치한 문산(파주), 동두천 관측소를, 청미천 유역의 경우 수원, 이천, 원주, 천안, 청주 관측소의 자료를 제공받았으며, 유역 내 위치한 각 플럭스 타워 지점에서 일별 실측 증발산량 자료를 제공받아 이용하였다. 각 기상관측소 및 플럭스 타워에 위치는 Fig. 2에 나타내었다.

4. 입력 자료

본 연구에서는 미항공우주국(NASA)에서 제공하는 Terra MODIS 위성영상을 SEBAL 모형의 입력 자료로 사용하였다. Terra MODIS는 1999년 12월에 발사되어 전 지구를 하루 혹은 이틀에 걸쳐 관측하며, 36개의 분광 밴드로 0.4~15 ㎛의 파장영역을 관측한다. Terra 위성은 705 km 고도에서 일 2회 한반도 상공을 지나며, 북쪽에서 남쪽 방향으로 오전 10시 30분경에 적도를 통과한다. 관측범위는 남북방향으로 5000 km, 동서방향으로 2,300 km이다. MODIS의 공간해상도는 밴드 1과 2는 250 m × 250 m, 밴드 3부터 7은 500 m×500 m 및 밴드 8부터 36까지는 1 km×1 km등 세 가지 종류가 있다(Lee et al., 2015). 본 연구에서는 Land Processes Dirtributed Active Archive Center (LP DAAC)에서 제공하는 1 km 공간 해상도의 Land Surface Temperature/Emissivity (MOD11_L2)와 Normalized Difference Vegetation Index (MOD13A3), 500 m 공간 해상도의 Albedo (MCD43A3)를 구축하였다. 모든 영상은 2012년 1월부터 2013년 12월의 기간에 걸쳐 구축하였으며, MOD11_L2는 일별로, MCD43A3은 8일 간격으로 제공하는 영상으로 이를 월단위로 평균하여 월별 자료로 활용하였다. 또한 MOD13A3은 월별로 제공되는 영상으로 MCD43A3과 동일하게 월 자료로 활용하였다. Albedo 및 NDVI는 LST에 비해 시계열 변동성이 적고, 8일 혹은 16일 간격의 영상을 활용할 경우 영상 내에 이상치가 포함되어 있어 모형결과에 영향을 줄 수 있기 때문이다. 구축한 위성영상의 공간해상도는 각각 500 m (MCD43A3)와 1000 m (MOD11_ L2, MOD13A3)로 모의를 위해 공간해상도의 일치가 요구된다. 본 연구에서는 대상 유역의 특성을 고려하고 선행 연구와의 비교를 위해 500 m의 공간해상도로 모의를 진행하였으며, 이를 위해 1000 m 해상도의 위성영상은 ESRI ArcGIS 프로그램으로 resampling하여 이용하였다. Fig. 3은 본 연구에서 사용한 위성 영상을 각 유역별로 정리하여 나타낸 것이다.

MODIS 영상은 원격탐사장비의 특성상 구름과 대기의 교란 등으로 인한 결측값이 발생한다. 따라서 본 연구에서는 이러한 결측값 보정을 위해 대상 유역 인근의 기상관측소 지점(문산, 동두천, 수원, 이천, 원주, 천안, 청주)에서 제공하는 지면온도 자료를 이용하여 Gap-Filling을 실시하였다. Gap-Filling은 기상청에서 제공하는 지면온도 자료를 이용해 Alcamo et al(1997)이 제안한 Eq. (5)를 적용해 편이보정(Bias Correction) 한 뒤, 두 영상의 유효한 값을 각각 추출하여 통합하고 공간내삽법 중 하나인 역거리 가중법(Inverse Distance Weighted; IDW)을 적용해 공간분포화 하였다(Eq. (6)).

                   (5)

                               (6)

Eq. (5)에서 는 보정된 기상 관측소 지면온도, 는 기상관측소에서 관측된 지면온도, 은 MOD11_ L2 지면온도 영상 평균, 는 기상관측소의 지면온도 평균이며, Eq. (6)에서 는 추정치, 는 영향 범위 내에 해당하는 샘플 위치, 는 점 사이 거리()에 반비례()하는 가중치이다. MOD11_L2 지면온도와 기상청 지면온도의 편의보정 전후의 결정계수는 각각 0.79, 0.99로 편의 보정 후 결정계수(coefficient of determination, R²)가 향상됨을 보였다.

SEBAL 모형의 지형 자료로는 NGIS (National Geographic Information System)의 1:5,000 수치지도로부터 MODIS 위성영상의 크기와 동일한 격자크기 500 m의 DEM (Digital Elevation Model)을 제작하여 사용하였다(Fig. 4a). 토지이용도는 환경부에서 제공하는 2007년 중분류 토지이용도를 이용하였으며, 우리나라 위도를 맵으로 만들어 이용하였다(Fig. 4b, c). 모형의 기상자료로는 기상 관측소로부터 제공받은 풍속, 풍속측정 높이, 일사량을 활용하였다. SEBAL 모형은 현열플럭스 산정에 지표면으로부터 2 m 높이에서 측정한 풍속 자료를 이용하기 때문에 본 연구에서는 각 관측소의 풍속측정 높이를 이용해 2 m 높이의 풍속으로 환산하여 적용되도록 하였다. 또한 모형의 매개변수로 각 기상 관측소 지점의 위, 경도를 이용해 지점별 태양의 시간각(hour angle), 일출각(sunrise angle), 일몰각(sunset angle)을 일별로 계산하였으며, 태양의 시간각은 Terra 위성의 통과시각인 오전 10:30에 대해 계산하였다. 각 매개변수의 계산은 천문우주지식정보(http://astro.kasi.re.kr/)를 참고하였으며, 계산된 매개변수 및 풍속, 일사량 자료는 IDW를 통해 내삽하여 적용하였다. 설마천 유역의 경우 기상관측소가 유역에서 일정 거리 이상 떨어져 있어 유역 인근에 위치한 적성 AWS (Automatic weather station)에서 관측하는 풍속 자료를 이용해 입력 자료의 검증을 실시하였다. 2013년의 풍속 자료를 바탕으로 결정계수를 산정한 결과 0.65로 양호한 상관성을 보였다.

5. 결과 및 토의

본 연구에서는 SEBAL 모형의 검증을 위해 각 플럭스 타워 위경도 지점의 pixel로부터 모의 증발산량을 추출하여 비교하였다. Fig. 5는 SMK와 CFK에서 관측한 증발산량과 해당 지점에서 SEBAL 모형을 이용해 모의한 증발산량을 추출해 비교한 것이다. 각 플럭스 타워에서 2012년부터 2013년까지 실측증발산량 최대값은 CFK가 5.2 mm/day, SMK가 3.7 mm/day로 CFK의 최대값이 약 1.7배 높은 모습을 보이며 전반적인 증발산량 분포는 CFK에서 높은 모습을 보이고 있다. 이는 CFK의 경우 5월부터 9월까지의 논의 경작시기 동안 담수재배를 하며 수면으로부터 일어나는 증발의 양이 전체적인 증발산량의 증가의 원인이 되는 것으로 보인다. 같은 기간의 모의 증발산량의 최대값은 CFK에서 5.3 mm/day, SMK에서 3.4 mm/day로 실측 증발산량과 유사한 분포를 나타냈다. 몬순 기후에 속하는 우리나라의 특성상 여름에는 고온 다습의 남서 계절풍의 영향으로 장마기간이 뚜렷하게 나타난다. SMK는 2012년 6월과 2013년 7월에 장마 기간 동안 실측 증발산량과 모의 증발산량이 감소하는 모습을 보이며, CFK에서도 2012년 7월과 2013년 7월에 장마 기간으로 실측값과 모의값 모두 이 기간 동안 증발산량이 감소하고 장마 이후 증가하는 모습을 보였다. 한편, SMK에서 2013년 3, 4월에 SEBAL 증발산량이 실측증발산량에 비해 다소 높은 값으로 모의되는 모습을 보인다. 이는 2013년 3, 4월의 증발산량이 최대 0.7 mm/day까지 상승할 때, 2012년의 3, 4월 실측증발산량이 최대 1.4 mm/day까지 상승하는 경향을 나타내는 것으로 보아 2013년의 실측증발산량이 평소보다 낮게 측정된 것으로 판단된다.

Fig. 6은 SMK와 CFK에 대하여 SEBAL 모형의 주요 입력 자료인 NDVI 및 Albedo, 실측 증발산량과 모의 증발산량의 월별 값을 나타낸 것이다. CFK의 NDVI는 4 ~ 5월에 급격히 상승하고 8월 이후 감소하는 모습을 보이는 반면 SMK는 계절별로 NDVI가 일정한 분포를 나타낸다. 또한 전체적인 NDVI가 SMK가 CFK에 비해 높은 모습을 보이는데 이는 토지피복에 따른 산림과 논 지역의 특성이 반영된 것으로 보인다. 논 지역의 경작 시기에 따라 CFK의 NDVI가 급격히 상승하는 모습을 보이는 것이다. 2013년 6, 7월의 SMK에서 NDVI의 감소는 실제 식생의 변화가 아닌 장마의 영향에 따른 것으로 판단되며, 이에 따라 SEBAL 모형의 월 증발산량도 낮게 모의되는 것을 확인할 수 있다. 지점별 albedo의 경우 봄, 여름 및 가을에는 식생의 변화에 따라 증감하는 경향을 나타내고 겨울철에 급격히 상승하는 모습을 보인다. 이는 지표 위의 눈 및 얼음의 영향으로 인한 것으로 판단되며, 이에 따라 겨울철의 모의 증발산량은 두 지점 모두 매우 낮게 모의 되고 있는 것을 확인할 수 있다. 플럭스 타워 실측 증발산량 또한 겨울철은 거의 0에 가까운 모습을 보인다. 선행논문(Lee et al., 2015)의 경우 albedo 및 NDVI에 월별 보정 계수를 적용해 각 인자들의 월별 분포가 일정하지 않고 산발적인 모습을 보인 반면 본 연구에서는 월별 보정 계수를 적용하지 않고 입력 자료의 개선 및 보정, 현열 플럭스 산정 알고리즘의 개선을 통해 결과의 신뢰성을 높이고자 하였으므로, NDVI 및 albedo의 월별 분포는 시계열 변동성이 일정하게 나타난다.

Fig. 7은 각 플럭스 타워의 실측 증발산량과 모의 증발산량의 상관성을 분석하여 나타낸 것이다. 모의 증발산량과 실측 증발산량의 결정계수(R²)는 CFK가 0.78, SMK에서 0.54로 CFK의 상관성이 더 높은 것으로 나타났다.

Table 1은 모의기간동안의 플럭스타워 실측 증발산량과 SEBAL 모의 증발산량을 계절별로 비교하여 나타낸 것이다. 모의 결과에 따른 모형의 적합성과 상관성을 판단하기 위한 목적함수로는 결정계수(R²), 평균제곱근오차(average root mean square error, RMSE), IOA (index of agreement) 및 Nash and Sutcliffe (1970)가 제안한 모형효율성계수(Nash- Sutcliffe efficiency, NSE)를 사용하였다. 여기서 R², IOA 및 NSE는 값이 1에 가까울수록, RMSE는 0에 가까울수록 모의값과 실측값의 상관성이 높음을 의미한다(Moriasi et al., 2007).

2012년에 CFK에서 실측값과 모의값은 각각 496.1 mm와 388.4 mm, 2013년에 467.8 mm와 354.6 mm로 SEBAL 모형이 실측값에 비해 약 90 mm정도 과소 추정되는 모습을 보였다. 반면 2012년에 SMK에서 실측값과 모의값은 각각 311.7 mm와 192.4 mm, 2013년에 293.8 mm와 174.2 mm로 모의값이 실측값의 약 6/10 정도로 나타났다. 두 지점에서 SEBAL 모의 증발산량이 과소 추정된 주된 이유로는 일별 anchor pixels (hot pixel과 cold pixel)로부터 산정한 현열 플럭스의 과대추정으로 인한 것으로 판단된다. 본 연구에서는 선행연구(Morse et al., 2000)에서 제시한 토양조건을 고려하기 위해 토지 이용도를 이용하여 cold pixel은 관개지역(논)에서 hot pixel은 시가지와 나지에서 선정하도록 하였다. 그러나 선행연구와 달리 우리나라의 경우 논 지역은 5 ~ 9월 사이의 시기에서만 담수재배를 하고 그 외의시기에는 나지와 같아 anchor pixel의 선정에 개선이 필요하다. 최근에는 cold pixel을 관개 지역뿐만 아니라 수역에서 선정하는 방법들이 제안되고 있으며(Jiang and Islam, 2001; Bastiaanssen et al., 2005) 지표 에너지 수지 이론에 근거하여 지상에서 측정한 기상 자료 및 식생 지수를 통해 anchor pixel을 선정하는 등 (Moran et al., 1994; Long et al., 2012) anchor pixel 선정 방법의 개선을 위한 다양한 연구가 수행되고 있다(Tang et al., 2015). 따라서 본 연구의 활용한 anchor pixel 선정 알고리즘을 국내 실정에 맞게 더욱 개선한다면 보다 신뢰성 있는 결과를 얻을 수 있을 것으로 판단된다.

SMK의 계절별 증발산량은 실측값과 모의값 모두 여름에 가장 높고 봄과 가을에 근소한 차이를 보인 반면, CFK는 여름철에 가장 높고 봄에 두 번째로 높은 것으로 나타났다. 이는 각 지점의 토지피복 특성에 따라 혼효림 지역인 SMK는 봄/가을의 식생 변화가 크지 않기 때문에 각 계절의 증발산량도 큰 차이를 보이지 않으나, CFK는 작물의 경작시기에 따라 봄부터 여름까지 증발산량이 최대화 되었다가 작물의 수확으로 인한 식생의 급격한 변화로 인해 가을에는 증발산량이 감소하는 것으로 판단된다. 이로 인해 SMK의 가을철 증발산량은 CFK 지점에 비해 높게 나타난다. 또한 모내기를 위해 5월부터 논 지역은 물을 채워놓기 때문에 봄철의 증발산량이 전체적으로 높아지는 것으로 보인다.

통계집단의 산포도를 나타내는 표준편차는 0일 때 값이 모두 동일하고 표준편차가 클수록 평균에서 떨어진 값이 많음을 나타낸다. SMK에서 모의값의 표준편차는 2012, 2013년 모두 여름에 가장 높게 나타났고, 봄과 가을에는 0.01~0.02의 미세한 차이를 보였다. 반면 실측값의 경우 가을에 표준편차가 가장 높고 여름, 봄 순으로 나타났다. 2012년 여름의 경우 표준편차는 0.79로 가을(0.81)과 유사한 분포 경향을 보였으나, 2013년 여름에는 0.48로 가을(0.82)에 비해 낮은 모습을 보였다. 이는 2013년 여름에 강우와 그로 인한 NDVI 값의 저하로 인해 2012년에 비해 밀집되어 모의되었음을 나타내는 것이다. CFK 지점은 계절별로 실측값과 모의값의 표준편차가 유사하게 나타났다.

각 지점의 연별 R²는 CFK에서 2012년에 0.77, 2013년에 0.80으로 2013년이 2012년보다 상관성이 더 좋게 나타났으며, SMK의 R²는 2012년에 0.53, 2013년에 0.54로 CFK에 비해 약한 상관성을 보였다. 계절별 R²는 SMK에서는 2012년 봄에 0.50으로 가장 상관성이 좋게 나타났으며, 2013년에는 가을에 0.54로 가장 높았다. CFK에서는 2012년에는 SMK와 동일하게 봄에 0.66으로 가장 높았고, 2013년에는 여름에 0.73으로 상관성이 가장 좋은 것으로 나타났다. 대상 지역 인근의 기상관측소 중 운량 자료를 제공하는 4개 관측소(원주, 수원, 청주, 천안)로부터 계절별 운량 분석 결과, 여름철의 평균 운량은 2012년에 64.8%, 2013년에 63.6%로 가장 높았으며 봄, 가을에는 2012년에 각각 48.4%, 47.5%, 2013년에 43.1%, 43.7%로 50% 이하로 나타났다. 높은 운량은 MODIS LST 일 자료의 유효값 획득을 저하시키는 큰 원인이 되며, 여름철 상관계수가 낮게 나타나는 원인으로 작용했을 것으로 사료된다. 한편 겨울에는 상대적으로 높은 albedo와 낮은 NDVI로 인해 모의 값이 실측값에 비해 낮게 모의되고 있음을 확인하였다. SMK에서 겨울의 실측값은 2012년에 6.4 mm, 2013년에 11.4 mm인 반면 모의값은 각각 0.7 mm, 0.6 mm로 큰 차이를 보였으며, 이러한 양적 차이로 인해 겨울의 R²는 0으로 나타났다. CFK도 겨울의 모의값은 2012년에 12.4 mm, 2013년에 7.1 mm 낮게 모의되었으며 R²는 각각 0.1, 0.26으로 낮은 상관성을 나타냈다. 선행 논문(Lee et al., 2015)의 경우 동일한 연구기간을 대상으로 모의한 결과와 우리나라 혼효림 지역에 위치한 덕유산 플럭스 타워 실측 자료와의 R²는 0.45로 본 연구의 결과가 소폭 개선되었음을 보였으나, 겨울철의 R²는 0으로 겨울철의 증발산량에 대한 모형의 개선효과는 미미한 것으로 나타났다.

한편 R²를 제외한 다른 목적함수는 겨울의 모의 결과가 가장 좋은 것으로 나타났다. RMSE는 SMK와 CFK 모두 겨울철에 0에 가장 가깝게 계산되며 상관성이 가장 높은 것으로 계산되었고, IOA와 NSE 또한 겨울철에 1에 가깝게 나타나며 가장 높은 상관성을 보였다. 이는 표본 집단의 1:1 비례관계에 따른 일치 정도를 판별하는 R²와 달리 실측값의 평균을 기준으로 계산하는 IOA와 NSE 및 제곱근 오차를 계산하는 RMSE는 함수 계산상의 특성으로 인해 0 ~ 10 mm 이하로 매우 작게 모의되는 모의 값이 오히려 상관성을 더 좋게 나타내는 것으로 보인다(Moriasi et al., 2007).

Fig. 8과 9는 2012년의 공간 증발산량 월별 분포를 나타낸 것이다. Lee et al. (2015)에서 SEBAL 모형의 공간 증발산량 분포는 표고의 분포에 따라 표고가 높은 지역에서 증발산량이 낮고 표고가 낮은 지역에서 높게 나타내는 분포를 나타내는 모습을 보였는데, 청미천 유역의 경우 이러한 경향성이 동일하게 나타났다. 3월의 경우 유역의 중간 지점에서 특별히 높게 산정되는 모습을 보이고 있는데, 이는 MODIS 월별 알베도 영상에서 구름으로 인해 이상치(異常値)로 산정된 부분을 Kotoda (1989)의 토지피복별 알베도에 따라 보정을 하였는데, 보정된 부분이 토지피복에 따라 알베도 값이 일정하고 불연속적으로 산개되어 있기 때문인 것으로 보인다. 설마천 유역의 경우 500 m × 500 m의 해상도가 표고에 따른 증발산 분포 경향성이 표현하기에는 상대적으로 커서 다른 유역과 같은 분포 경향성이 잘 나타나지 않는 것으로 보인다.

6. 결  론

본 연구는 선행논문에서 적용했던 SEBAL 모형의 개선을 통해 우리나라의 혼효림, 논 지역에 증발산량을 산정하고 설치된 플럭스 타워 실측 증발산량을 이용해 모형의 효용성을 검토해보고자 하였다. 연구에 이용한 SEBAL 모형은 Lee et al. (2015)이 구축한 방법과 동일하게 MODIS 위성영상을 이용하여 증발산량을 산정할 수 있도록 일부 알고리즘을 수정하였으며 LST (MOD11_L2), NDVI (MOD13A3), Albedo (MCD43A3) 영상을 500m 해상도로 구축하여 활용하였다. albedo 및 NDVI 영상은 8일 혹은 16일 간격의 자료 활용 시 영상 내에 존재하는 이상치가 모형 결과에 영향을 줄 수 있기 때문에 보정된 월별 자료를 구축하여 활용하였으며, LST 영상은 구름이나 대기의 교란으로 인해 결측값이 발생하기 때문에 유역 인근의 기상관측소에서 제공하는 지면온도 자료를 활용해 공간 내삽법 중 하나인 IDW를 적용하여 gap-filling을 실시한 뒤 모의에 활용하였다.

2012년부터 2013년까지 증발산량을 산정한 결과 플럭스 타워 실측 증발산량과 시계열 분포는 유사하게 모의되었으나, 전체적인 모의 증발산량은 과소 추정되는 모습을 보였다. 이는 일별 anchor pixels (hot pixel과 cold pixel)로부터 산정한 현열 플럭스의 과대추정으로 인한 것으로 판단되며, 본 연구의 활용한 anchor pixel 선정 알고리즘을 국내 실정에 맞게 더욱 개선한다면 보다 신뢰성 있는 결과를 얻을 수 있을 것으로 판단된다. SMK의 계절별 증발산량은 실측값과 모의값 모두 여름에 가장 높고 봄과 가을에 근소한 차이를 보인 반면, CFK는 여름철에 가장 높고 봄에 두 번째로 높은 것으로 나타났다. 이는 지점별 토지 피복 특성에 따라 혼효림 지역인 SMK는 봄/가을의 식생 변화가 크지 않기 때문에 각 계절의 증발산량도 큰 차이를 보이지 않으나, CFK는 작물의 경작시기에 따라 봄부터 여름까지 증발산량이 최대화 되었다가 작물의 수확으로 인한 식생의 급격한 변화로 인해 가을에는 증발산량이 감소하는 것으로 판단된다. 특히, 여름철 장마기간에 나타나는 전체적인 증발산량의 감소와 장마기간 뒤에 반등하는 경향성을 따라가는 모습을 보여 위성영상을 통한 증발산량의 산정이 지면의 기상조건과 연관성을 가질 수 있음을 나타냈다. 각 지점에서의 albedo와 NDVI를 추출해 비교해본 결과, SMK에서 2013년 여름에 NDVI의 저하에 따라 모의 증발산량이 감소하는 모습을 보였고, 겨울에 albedo가 상승하고 NDVI가 감소함에 따라 증발산량이 적게 모의되는 모습을 나타냈다. 또한 지점별 NDVI의 값의 차이에 따라 모의 증발산량이 변화함으로써 위성 자료인 NDVI 및 albedo를 통해 지상의 토지 피복 특성을 고려하여 모의할 수 있음을 나타냈다. 각 지점에서 전체 모의 기간에 대한 R²는 CFK가 0.78, SMK에서 0.54로 선행연구에 비해 개선된 상관성을 보였다. SEBAL 증발산량의 공간분포는 선행연구(Lee et al., 2015)와 같이 표고의 분포에 따라 표고가 높은 지역에서 증발산량이 낮고 표고가 낮은 지역에서 높게 나타내는 분포를 나타내는 모습을 보였으나 설마천 유역과 같은 작은 유역의 경우 500 m 공간 해상도의 MODIS 영상은 표고에 따른 분포 경향성을 표현하기에 한계를 나타냈다.