(Chun Gun Park)
박천건1
(Chang-Soo Rim)
임창수2†
(Kwang-Suop Lim)
임광섭3
(Hyo-Sok Chae)
채효석4
-
경기대학교 수학과 조교수
(Kyonggi University)
-
경기대학교 토목공학과 부교수
(Kyonggi University)
-
K-water 연구원 책임연구원
(K-water)
-
K-water 연구원 수자원연구소 수석연구원
(K-water)
Key words (Korean)
잠열플럭스, 이상치, 에디공분산, 통계적 분석, IQR 규칙
Key words
Latent heat flux, Outliers, Eddy covariance, Statistical analysis, IQR rule
-
1. 서 론
-
2. 연구자료 및 방법
-
2.1 연구지역 및 자료
-
2.2 연구방법
-
2.2.1 잠열플럭스 측정
-
2.2.2 IQR 규칙 및 대체방법
-
3. 연구결과
-
3.1 이상치 탐색 및 대체
-
3.2 실제 증발산량 비교
-
4. 요약 및 결론
1. 서 론
증발산은 수문순환과정에서 중요한 요소 중의 하나로서 실제증발산량을 이해하는 것은 지구상의 물순환 과정 중에서 지구상에 존재하는 수자원의 손실량을 이해하고
수자원을 효율적으로 관리하는데 필요하다. 실제증발산량 관측을 위해서 가장 많이 적용되는 방법은 에디 공분산(eddy covariance)방법을 이용하는
것으로 전 세계적으로 많은 실험적 관측이 수행되고 있다(Kustas et al., 1994; Humes et al., 1994; Oak Ridge
National Laboratory Distributed Active Archive Center for Biochemical Dynamics, 2014;
FLUXNET Oak Ridge National Laboratory, 2014; US-China Carbon Consortium, 2014; AsiaFlux,
2014; Max Planck Institute for Biogeochemistry, 2014; Global Energy and Water Exchange
Project-World Climate Research Programme, 2014; CARBO-AFRICA, 2014; CarboEurope-IP,
2014; OzFlux, 2014).
국내의 경우 광릉 국립수목원에서는 2000년 4월부터 에디 공분산법(eddy covariance method)을 이용하여 낙엽활엽수림과 대기간의 탄소,
에너지 및 물의 교환과정을 이해하기 위한 연구를 수행하였다. 이를 위하여 30m 관측 타워가 설치되었고, 현열, 잠열 및 운동량 플럭스의 관측이 시작되었다(Park,
2001). 설마천 유역은 유역면적이 8.50km2이고, 1995년부터 한국건설기술연구원에서 시험유역으로 운영하고 있으며, 증발산은 에디 공분산법을 사용하여 관측되고 있다. 기상연구소에서는 2002년
7월부터 해남의 기상 관측소 부근의 농경지에 20m 플럭스 관측 타워를 설치하여 플럭스 관측을 수행하고 있다(Kim et al., 2003; Lee
et al., 2003). 하리 농경지는 강화도 삼산면 석모도에 위치하고 있으며, 이 지역은 넓고 균질한 벼논으로 이루어져 있고, 에디공분산 방법으로
현열과 잠열 플럭스를 측정하고, 순복사, 기온, 수증기압, 지중열 등도 측정하고 있다(National Research Foundation of Korea,
2004). Kang et al. (2013)은 낙동강 강정 고령보 근처에서 관측된 플럭스 자료를 이용하여 농경지와 도시의 플럭스 차이를 분석한 바
있다.
에디 공분산방법을 이용하여 잠열 플럭스를 측정하는 경우 플럭스 원시자료의 품질관리가 필요하며(Vickers and Mahrt, 1997), 특히 우리나라와
같이 산지가 많은 지역에서는 평지에서 관측하는 방법과 달리 산지지형을 고려하여 잠열 플럭스를 측정하여야 한다. Wilczak et al. (2001)은
이러한 산지지형을 고려하여 잠열플럭스를 측정하는 방법을 제안한 바 있다. 또한 잠열 플러스의 측정에 필요한 요소들 중에 풍속의 측정이 매우 중요하다.
만약 측정된 풍속이 측정기기의 내부적인 오류 및 비, 바람 등과 같은 외부적인 충격에 의해서 자료가 오염되었다면 올바른 잠열 플러스의 측정이 어려울
것이다. 오염된 자료라는 것은 자료에 이상치가 포함되어 있다는 것을 의미하며 과학적인 분석을 위해서는 필히 이상치에 대한 처리가 필요하다.
과거에 수행되었던 플럭스 관측 시 발생되는 연직 풍속의 이상치의 탐색 및 대체방법과 관련된 연구사례를 분석한 결과 대부분의 경우 관측된 자료의 10Hz
관측값과 30분 평균값과의 표준편차를 비교하여 이상치를 탐색하고 제거한 바 있다(Lim, 2007; Korean Meterological Administration,
2013; Kang et al., 2013). 한편 Papale et al. (2006)은 연속된 3개 관측값의 차이를 이용하여 이상치를 탐색하는
방법을 제안하였으며, Hong et al. (2009)은 광릉 수목원 플럭스타워에서 관측된 연직 풍속자료의 이상치를 탐색하기 위해 Papale et
al. (2006)이 제안한 연속된 3개 연직 풍속 관측값의 차이를 이용한 방법을 보완한 방법을 제안한 바 있다. 또한 관측 기기에 따른 연직 풍속오차를
분석하기 위해 Frank et al. (2013)은 직각 음파 풍속계(orthogonal sonic anemometer)와 비직각 음파 풍속계(non-orthogonal
anemometer) 사용에 따른 연직 풍속 관측값의 차이를 연구한 바 있다. 하지만 과거 연구사례를 분석한 바에 의하면 연직 풍속 이상치에 대한
대체방법에 대한 연구가 미미하며, 앞으로 플럭스 관측 시 연직 풍속의 이상치 탐색과 대체 방법에 따른 증발산량 변화에 대한 더 많은 연구가 수행되어야
할 것으로 판단된다. 한편 증발산량 관측 이외의 다른 많은 응용분야에서는 이상치가 분석 및 해석의 방향을 반대로 이끌 수 있기 때문에 이상치의 처리에
많은 관심을 두어왔다. 일반적으로 이상치는 대부분의 자료가 생성된 체계가 아닌 것으로부터 온 자료라고 정의할 수 있다(Hawkins, 1980).
이상치의 처리과정은 탐색부터 시작되어 대체까지를 말하다. 탐색기법은 크게 통계적 모형을 기반으로 한 공간 근접도 그리고 다변수 접근 등으로 구성된다.
특히 단일 확률변수(또는 자료)인 경우에는 모수적인 방법과 비모수적인 방법으로 이상치 탐색방법이 나누어진다. 본 연구에서는 플럭스의 분포가 알려져
있지 않기 때문에 비모수적인 방법 중 interquantile range (IQR) 규칙을 이용하여 이상치를 탐색하였다(Moore and McCabe,
1999). 또한 플럭스 이상치 탐색 결과를 이용한 이상치를 보완하기 위해 이상치로 확인된 값들을 실시간으로 대체할 수 있는 통계적 방법들을 제시하고,
이상치 대체 전후의 플럭스 변화량을 비교 검토하였다.
2. 연구자료 및 방법
2.1 연구지역 및 자료
본 연구의 지역은 전라북도 안성면의 덕곡제 유역(동경 127° 42′23″~127° 44′53″, 북위 35 °50′47″∼35° 52′50″)이며,
유역 내 토지이용분포를 보면 대부분 산림으로 이루어져 있다. 임상분포를 분석한 바에 의하면 유역 내 대부분의 경우 침엽수와 활엽수의 혼합림과 낙엽송림이
분포하고 있다. 유역의 평균표고는 1007.58m이고, 유역의 평균경사는 26.89°이다. 플럭스 관측을 위해 군락의 높이를 고려하여 플럭스 타워
높이를 25m로 하여 타워본체 및 케이블을 설치하였고, 본 연구자료의 경우 19m 지점에 설치된 에디공분산 시스템인 3차원 풍속계(CAST3)와 /가스분석기(EC 155) 그리고 자동기상시스템(Automatic Weather System: AWS, 모델 WXT510-NC) 1기로부터 측정된 자료를
이용하였다.
2.2 연구방법
2.2.1 잠열플럭스 측정
본 연구에서는 잠열 플럭스를 측정하기 위해 에디 공분산방법을 적용하였으며, 연구지역이 평지가 아니고 경사가 있는 산지라는 점을 고려하여 평판회전(planar-fit
rotation)을 실시하였다. 관측지점의 지형경사를 하나의 평판으로 가정하여 360° 전체를 하나의 평판으로 가정하는 경우(Wilczak et al.,
2001) 3차원풍속계로부터 얻어진 0.1초 간격 z, x, y 방향 풍속을 이용하여 30분 평균 z, x, y방향 평균 풍속 , , 을 계산하였다. 또한 과거 연구사례에 의하면 30분 평균값을 적용하는 경우 타당한 것으로 나타났다(Choi et al., 2003; Göckede
et al., 2004; Yuan et al., 2007; Kwon et al., 2009; Frank et al., 2013; Korean Meteorological
Administration, 2013). 30분 평균 풍속자료 336개(7일 자료)를 이용하여 Eq. (1)과 같은 다변량 회귀분석을 실시하여 계수
, , 를 산정하였으며, 산정된 계수 , , 를 이용하여 다음 Eq. (2)로부터 연직방향 실제 풍속()을 측정하였다.
(1)
(2)
여기서, , , 은 평판회전을 위한 좌표변환 변수이다(Wilczak et al., 2001).
바람의 난류적 특성에 관련되어 있는 성분들은 수증기 밀도() 및 연직풍속()이며 각 성분들은 어느 일정기간의 평균값과 변동하는 부분으로 나타내질 수 있으므로 대기 중의 바람에 의한 잠열 플럭스는 다음 Eq. (3)과 같이
표현되어질 수 있다(Oke, 1987).
(3)
만일 30분 동안 대기성분의 수직변동 부분을 평균화시키면 Eq. (3)은 다시 Eq. (4)와 같이 간략하게 나타낼 수 있다. 따라서 Eq. (4)을
이용하여 잠열플럭스()를 계산하였다.
(4)
여기서 =잠열 플럭스(latent heat flux, W/m2), =잠재증기화열(J/kg), = 0.1초 간격 연직방향 풍속(m/sec) , = 30분 평균 연직방향 풍속(m/sec), = 0.1초 간격 공기의 밀도(kg/m3) , = 30분 평균 공기의 밀도(kg/m3)
위에서 언급되어진 잠열 플럭스를 관측하기 위해 지표면으로부터 일정한 높이에서 수직 풍속과 수증기 변화를 매우 빨리 감지할 수 있는 관측장비가 필요하며,
지표면에 가까울수록 와류(eddy)의 크기가 작아지므로 이를 감지할 수 있도록 보다 민감하게 작용하는 관측장비를 필요로 한다.
2.2.2 IQR 규칙 및 대체방법
IQR 기법은 상자그림(box plot)에서 1사분위수 및 3사분위수에서 거리보다 더 멀리 떨어진 자료를 이상치라고 한다(Frigge et al., 1989). 여기서 는 통상적으로 1.5 또는 3으로 설정되며, Fig. 1에서 1.5와 3 사이에 위치한 관측값들을 온순한(mild) 이상치 또는 의심되는 이상치이고
3 이상에 위치한 값들을 심한(severe) 이상치라고 구분한다.
|
|
|
Fig. 1. Outlier Locations in Box Plot
|
Fig. 2는 자료의 분포가 알려져 있지 않은 경우에 이상치를 탐색하고 대체하여 증발량을 계산하는 방법까지를 순서도로 표현했다. 우선 0.1초 간격의
x, y, z 방향 풍속자료를 입력하고, 30분 간격으로 x, y, z 방향 평균풍속자료를 구분한 후 z 방향의 풍속자료를 IQR규칙에 의해서 이상치를
탐색하였다. z 방향뿐만 아니라 동일한 위치인 x, y 방향에 대해서 삭제 또는 30분 평균으로 대체하고, 회귀계수산정과 평판회전을 수행한 후 Eq.
(4)로부터 잠열플럭스를 계산하였다.
3. 연구결과
|
|
|
Fig. 2. Flow Chart for Outlier Detection and Replacement
|
3.1 이상치 탐색 및 대체
풍속자료의 이상치 탐색은 x, y, z 축 모두에 IQR 규칙이 적용되어 수행되었다. 따라서 어느 한 방향이라도 이상치라고 판별이 되면 다른 방향의
풍속값들도 동일한 이상치로 간주되고, 탐색된 이상치들은 해당되는 x, y, z축 모두 30분 평균으로 대체하거나 삭제하였다. Fig. 3은 2011년
6월 3일(비강우일) 및 6월 25일(강우일)의 0시부터 24시까지 0.1초당 연직방향(z 방향) 풍속 자료이다. 패널 (a), (d)는 원자료의
산점도로 나타낸 것이며, (b), (e)는 IQR 규칙을 로 적용하여 탐색된 이상치만 표시했고 그리고 (c), (f)는 로 적용하여 탐색된 심한 이상치만을 표시한 것이다. 강우일(2011년 6월 25일)인 경우 0.1초 단위로 측정된 풍속의 산점도에 의하면 새벽과 오전까지는
풍속이 잔잔했으나 오후이후부터는 점차 심한 것으로 나타났다. 특히 대략 밤 20시 이후에는 풍속이 더욱 강하게 나타났다. 따라서 이상치도 풍속이 심한
밤 시간대에 많이 나올 것으로 추측될 수 있다. 반면 비강우일(2011년 6월 3일)인 경우 대체적으로 풍속의 변동이 심하지 않은 것으로 나타났다.
|
|
|
|
(a)
|
(b)
|
|
|
|
|
(c)
|
(d)
|
|
|
|
|
(e)
|
(f)
|
|
Fig. 3. Outlier Detection and Replacement (June 3, 2011, Sunny and June 25, 2011,
Rainy) (a) & (d) Raw Data, (b) & (e) Outlier Detected by IQR(), (c) & (f) Outlier Detected by IQR()
|
IQR 규칙()에 의해서 이상치를 탐색하고 30분 간격에 따른 이상치의 비율은 강우일이 시간에 관계없이 대략 3%보다 높았으며 최대 11.5%까지도 나온 반면
비강우일에 대한 이상치 비율이 3.4%미만으로 기록됐다. 더욱이 강우일에 대한 IQR 규칙()에는 시간에 관계없이 이상치의 비율이 1.5% 미만의 결과가 도출되었고 비강우일은 0.22%로 상당히 작은 이상치 비율로 도출되었다.
3.2 실제 증발산량 비교
이상치 판별 방법과 이상치 대체 방법에 따른 증발산량의 차이를 비교하기 위해 비강우일 혹은 적은 강우일과 강우일로 구분하여 측정된 일별 증발산량을
비교분석하였다(Table 1 and Fig. 4). Fig. 4의 경우 첫 번째는 원자료와 IQR 규칙 을 적용한 이상치 탐색과 평균 대체, 두 번째는 원자료와 IQR 규칙 을 적용한 이상치 탐색과 평균 대체, 세 번째는 원자료와 IQR 규칙 을 적용한 이상치 탐색과 삭제, 네 번째는 원자료와 IQR 규칙 을 적용한 이상치 탐색과 삭제된 경우를 보여준다. 비교분석을 위해 비강우일 혹은 아주 적은 강우일로 선정된 날은 2011년 5월 28일, 6월 3일,
6월 28일, 그리고 7월 18일이고, 강우일로 선정된 날은 2011년 6월 25일, 6월 26일 그리고 7월 12일이다. 일별 실제증발산량 측정결과
비강우일의 경우 2011년 5월 28일은 2.40mm/day, 6월 3일은 2.43mm/day, 6월 28일은 2.35mm/day, 7월 18일은
3.09mm/day를 보였다. 강우일의 경우 2011년 6월 25일은 0.71mm/day, 6월 26일은 1.84mm/day, 7월 12일은 0.98
mm/day를 보였다.
분석결과에 의하면 IQR규칙을 적용하여 이상치를 판단하고 이상치를 대체(평균, 삭제)하여 증발산량을 관측하는 경우 증발산 원자료와 비교하여 비강우시보다는
강우 시에 보다 큰 차이를 보였다(Table 1). 이러한 결과는 Fig. 4에서도 나타나고 있다. 그림에서 보는 바와 같이 강우 시(Figs. 4(a)~(d))의
경우 비 강우 시(Figs. 4(d)~(f))와 비교하여 시간별로 증발산변화가 크게 나타나고 있다. 이는 강우 시 혹은 적은 강우 시에 증발산량 관측과정에서
비강우시보다 측정 상 불확실성이 크다는 것을 보여준다. 또한 IQR (k=1.5)를 적용하는 경우 강우 시(June 26, 2011)를 제외하고 평균을
사용하여 이상치를 대체하는 경우가 이상치를 삭제하는 경우보다 원자료와 비교하여 보다 큰 차이를 보였다. 반면에 IQR (k=3.0)을 적용하는 경우
강우 시(June 25, 2011; July 12, 2011)에 이상치를 삭제하는 경우보다 평균치로 이상치를 대체하는 경우에 원자료와 비교하여 보다
큰 차이를 보였다. 또한 같은 평균치를 적용하거나 이상치를 배제하고 증발산량을 관측하는 경우 모두에서 IQR (k=1.5)를 적용하는 경우가 IQR
(k=3.0)을 적용하는 경우보다 증발산 원자료와 비교하여 보다 큰 차이를 보였다.
|
Table 1. Comparisons of Outlier Detection and Replacement Methods for Evapotranspiration
Measurement (mm/day)
|
|
Measurement day
(weather condition)
|
Raw data
|
IQR(k=1.5)
|
IQR(k=3.0)
|
|
Mean
|
% change
|
Deleting
|
% change
|
Mean
|
% change
|
Deleting
|
% change
|
|
May 28, 2011 (no rain)
|
2.40
|
2.07
|
-13.73
|
2.19
|
-9.02
|
2.38
|
-1.07
|
2.43
|
1.04
|
|
June 3, 2011 (no rain)
|
2.43
|
2.13
|
-12.38
|
2.18
|
-10.40
|
2.42
|
-0.49
|
2.42
|
-0.31
|
|
June 28, 2011 (light drizzle)
|
2.35
|
2.11
|
-10.49
|
2.14
|
-9.09
|
2.34
|
-0.64
|
2.33
|
-0.77
|
|
July 18, 2011 (no rain)
|
3.09
|
2.47
|
-20.11
|
2.67
|
-13.42
|
3.01
|
-2.42
|
3.00
|
-2.82
|
|
June 25, 2011 (rain)
|
0.71
|
0.40
|
-43.01
|
0.53
|
-25.67
|
0.63
|
-10.86
|
0.73
|
2.81
|
|
June 26, 2011 (rain)
|
1.84
|
1.44
|
-21.72
|
1.11
|
-39.41
|
1.78
|
-2.85
|
1.67
|
-9.15
|
|
July 12, 2011 (rain)
|
0.98
|
0.70
|
-28.33
|
1.04
|
5.56
|
0.93
|
-5.30
|
0.99
|
0.85
|
|
|
|
|
|
Fig. 4. Comparisons of Hourly Latent Heat Flux Density According to Outlier Detection
and Replacement Methods (a) May 28, 2011 (b) June 3, 2011 (c) June 28, 2011 (d) July
18, 2011 (e) June 25, 2011 (f) June 26, 2011 (g) July 12, 2011
|
이러한 분석의 차이는 이상치 탐색에 적용된 IQR규칙이 이상치의 크기뿐만 아니라 개수와 위치에 기인한다. 즉 k=1.5가 적용되어 탐색된 이상치의
비율은 대략 30분 구간 자료에서 5%∼6%차지하고 있지만 k=3은 1% 내외인 것으로 나타났으며 Eq. (1)과 Eq. (2)에서 추정된 계수가
이상치의 탐색 및 대체에 의존적이다. 더욱이 Fig. 1로부터 IQR 규칙에 의한 이상치의 탐색에서 k=1.5와 k=3 사이에서 발생된 이상치는 의심의
여지가 있지만 정상적인 자료라고 확신할 수 없다. 따라서 이상치로 간주하는 것이 올바른 분석의 정도라고 할 수 있다. 따라서 분석결과에 의하면 증발산량
관측에 있어서 비강우시보다 강우 시에 이상치를 고려한 적절한 관측 및 분석 방법에 대한 세밀한 연구가 필요하다.
4. 요약 및 결론
본 연구에서는 에디 공분산방법을 이용하여 실제 증발산량을 측정하는 경우 연직 풍속의 이상치 탐색 결과에 따른 이상치를 보완하기 위해 이상치로 확인된
값들을 실시간으로 대체할 수 있는 통계적 방법들을 적용하여 이상치 대체 전후의 플럭스 변화량을 비교 검토하였다. 분석결과에 의하면 IQR규칙을 적용하여
이상치를 판단하고 이상치를 대체하여 증발산량을 관측하는 경우 증발산 원자료와 비교하여 비강우시나 혹은 적은 강우 시 보다는 강우 시에 보다 큰 차이를
보였다. 이는 강우 시에 증발산량 관측과정에서 비강우시보다 측정하는 과정에서 불확실성이 크다는 것을 보여준다. 따라서 증발산량 관측에 있어서 비강우시보다
강우 시에 이상치를 고려한 적절한 관측방법이 필요하다.
IQR (k=1.5)를 적용하는 경우 강우 시(June 26, 2011)를 제외하고 평균을 사용하여 이상치를 대체하는 경우가 이상치를 삭제하는 경우보다
원자료와 비교하여 보다 큰 차이를 보였다. 반면에 IQR (k=3.0)을 적용하는 경우 강우 시(June 25, 2011; July 12, 2011)에
이상치를 삭제하는 경우보다 평균치로 이상치를 대체하는 경우에 원자료와 비교하여 보다 큰 차이를 보였다. 또한 같은 평균치를 적용하거나 이상치를 배제하고
증발산량을 관측하는 경우 모두에서 IQR (k=1.5)를 적용하는 경우가 IQR (k=3.0)을 적용하는 경우보다 증발산 원자료와 비교하여 보다 큰
차이를 보였다. 따라서 적용된 이상치 탐색과 대체방법에 따라서 증발산량 측정값에 차이를 보이고, 이상치 대체 시 증발산량 측정값이 원자료 증발산량과도
차이를 보여서, 증발산 관측 시 연직풍속 이상치를 고려한 증발산관측이 이루어져야 할 것으로 판단된다. 또한 연직 풍속 이상치의 문제점을 개선하고 이를
증발산량 관측에 적용하기 위해서는 보다 많은 이상치 탐색 및 대체방법들을 적용하여 이들 방법들이 증발산 측정에 적합한지를 검증하는 추가 연구가 필요하다.