최근 12년간 아까시나무 만개일의 변화와 과정기반모형을 활용한 지역별 만개일 예측

아까시나무 개화시기 예측 모형의 모수 추정 및 검증을 위해서 2006년부터 2017년까지 전국 각지에서 수집된 아까시나무의 만개일 자료를 활용하였다. 개화 관측은 매년 최남단 지역의 아까시나무 꽃이 만개한 시점에 해당 지역을 시작으로 북부 방향으로 이동하며 진행되었다. 조사지는 수평적 산지구분도에 근거하여 기후대별로 선정되었고, 그 지역을 대표할 수 있는 아까시나무 임분이 조사구로 선택되었다(가로수 및 공원수는 관측 대상에서 제외). 조사 지점마다 흉고직경 10~15 cm 사이의 건전한 평균목 5본이 관측대상목으로 지정되었으며, 남쪽으로 뻗은 가지 중에서 수관을 3등분 하였을 때 하부에 위치한 1차 가지를 대상으로 개화 상태가 측정되었다. 관측 시점에 측정된 개화 상태로부터 만개일을 추정하기 위하여 아까시나무의 전체 개화기를 7단계로 구분하였다(1단계: 꽃망울이 맺힌 상태, 2단계: 꽃망울이 벌어지기 직전 상태, 3단계: 꽃잎이 반 정도 열린 상태, 4단계: 꽃잎이 완전히 열린 상태, 5단계: 꽃잎이 변색하는 상태, 6단계: 꽃잎이 떨어지고 있는 상태, 7단계: 꽃잎이 완전히 떨어진 상태). 평년기온에서 한 단계는 약 2~3일간 지속되며, 전체 단계 중 만개기는 4단계에 해당된다. 관측 시 조사목마다 개화 단계를 측정한 후 평균 개화 단계를 계산하였으며, 이로부터 만개일을 역으로 추정하였다.

전체 조사지점 중, 결측년이 1년 이하인 남부 지역 8개소(목포, 땅끝, 해남, 완도, 남해, 진주, 밀양, 청도), 중북부지역 7개소(춘천, 홍천, 충주, 횡성, 원주, 양구, 제천), 서울 도심지역 11개소(우장산공원, 서오릉, 강변북로, 관악산, 현충원, 북악산, 남산, 정릉, 홍릉, 태릉역앞, 어린이대공원)에서 관측된 자료를 바탕으로, 지난 12년 간의 아까시나무 만개일의 변화 경향을 분석하고 만개일 예측을 위한 모형을 개발하였다. Table 1에 각 조사 지점의 경위도, 측정 고도 및 관측 기간 동안의 평균만개일을 수록하였으며, Figure 1에 해당 조사 지점들의 지리적 위치를 나타내었다.

Figure 1. Geographical locations of 26 site, where the flowering of Robinia pseudoacacia was observed.

Download Original Figure

아까시나무 개화 관측 지점의 일기온정보를 추정하기 위해 1 km 해상도로 상세화된 과거 일기온 복원자료를 활용하였다. 해당 자료는 종관기상관측시스템(Automated Surface Observing System, ASOS), 방재기상관측시스템(Automatic Weather System, AWS), 산악기상관측시스템으로(Automatic Mountain Meteorology Observation System, AMOS)에서 관측된 기온정보를 융합하여 공간 분석을 수행한 후, 일별 기온 3요소(일최고기온, 일평균기온, 일최저기온)를 래스터 정보로서 산출한 자료이다. 공간 분석은 역거리 가중치법(Inverse Distance Weighting Method, IDW Method)을 바탕으로 하였으며, 고도편차에 의한 일별 기온감률을 적용하였다(Yun et al., 2001). 개화 관측 지점이 위치한 격자에서 일별 최고기온, 평균기온, 최저기온 정보를 추출하였다. 땅끝과 해남 지점은 산출된 래스터 자료의 격자 범위 외부에 위치하므로, 두 지점의 경우 실제 관측 지점와 가까우면서도 고도가 비슷한 범위 내에 있는 격자를 선택, 해당 격자에서 일기온 정보를 추출하였다. 전체 관측 기간 동안의 지점별 겨울철(12~2월) 평균 기온 및 봄철(3~5월) 평균기온을 Table 1에 나타내었다. 남부지역의 조사 지점 중 목포와 땅끝, 완도, 남해는 해안으로부터의 거리가 1 km 이내에 불과하며, 겨울철 평균기온은 2℃ 이상, 봄철 평균온도는 12℃ 이상으로 다른 지점보다 상대적으로 기온이 높았다. 중북부 지역의 경우, 모든 관측 지점들은 해발고도가 100 m 이상으로 다른 지역보다 상대적으로 높은 지대에 위치하고 있으며 낮은 기온 분포를 보였다(Table 1).

본 연구에서 활용한 식물계절모형은 Cesaraccio et al. (2004)이 제안한 저온요구량모형(Chill days model)이다. 이 모형은 온대 낙엽수의 겨울철 휴면생리에 근거한 과정기반모형으로 내생휴면과 환경휴면 두 단계의 휴면 과정을 모의한다. 모형에 의하면, 겨울철 냉각온도의 적산량이 저온요구량(Chilling requiremnet)에 도달하면 내생휴면이 해제되고 그 이후, 생육촉진온도의 적산량이 같은 크기의 임계값에 도달하면 환경휴면이 타파되고 눈(bud)은 생장을 재개한다(2단계 순차휴면모형, Sequential model). 냉각온도 및 생육촉진온도의 적산률은 단일삼각도일 계산법(Single triangle degree day calculation, 단위: °D(degree days))으로 구하였다(Table 2).

아까시나무 만개일 모형의 모수는 내생휴면 시작일과 저온요구량, 기준온도, 고온요구량(Heating requirement)이다. 고온요구량은 내생휴면이 해제된 후 꽃이 만개하기까지 생육촉진온도의 총 적산량으로 정의하였다. 온대 낙엽수목이 휴면을 시작하는 시점은 직접 관측이 불가능하므로, 대부분의 저온요구량모형에서 휴면 시작일은 가을철 식물계절현상이 시작되는 날짜(Cesaraccio et al., 2004), 임의로 가정한 특정한 날짜 등으로 대체된다(Jung et al., 2005; Matzneller et al., 2014; Chmielewski and Götz, 2016). 본 연구에서는 10월 1일을 휴면의 시작점으로 가정하였으며, 결과적으로 휴면 시작일을 제외한 저온요구량, 기준온도, 고온요구량을 모형의 모수로서 추정하였다.

12년 동안의 만개일 자료는 년도별로 7:5 또는 6:5의 비율로 나누어 각각 모수 추정을 위한 자료(훈련자료, calibration dataset)와, 개발 모형의 성능을 검증하는 자료(검증자료, validation dataset)로서 활용하였다. 모형의 과적합을 최소화하기 위하여 자료의 분리는 무작위적으로 이루어졌으며, Wilcoxon 순위합 검정을 시행하여 최종적으로 유사한 분포를 갖는 두 집단으로 분리하였다(Gehan, 1965; Bogawski et al., 2019; Liu et al., 2019). 결과적으로 2006, 2008, 2009, 2010, 2015, 2016, 2017년의 만개일 자료를 모수 추정에 사용하였고, 그 외의 자료는 모형 검증에 사용하였다.

모수 추정에 R 패키지 GenSA를 활용(Xiang et al., 2013), 최적화 과정에 담금질 기법(generalized simulated annealing)을 적용하였다(Chuine et al., 1998; Kalvans et al., 2015; Roberts et al., 2015). 알고리즘 적용에는 모수의 최적화 범위를 설정하는 단계가 선행되어야 한다. 최적화 범위가 주어지면, 예측 만개일과 관측 만개일의 평균제곱근오차(Root mean squared error, RMSE)가 최소가 되는 모수 추정치가 해당 범위 내에서 결정된다. 저온요구량, 기준온도, 고온요구량의 최적화 범위는 각각 -200 이상 -40 이하(단위: °D), 2 이상 15 이하(단위: ℃), 40 이상 600 이하(단위: °D)로 설정하였다. 모형마다 난수의 씨앗(seed) 값을 1부터 30까지 지정하여 30번 담금질 과정을 시행하였으며, 가장 낮은 RMSE를 갖는 모수 조합을 최적 모수로 선택하였다.

또한, 전 지역에 대해 단일 모수를 추정하는 전국 통합 모형(Single Model, SM)과 전국 통합 모형에 지점별 보정계수를 도입한 수정된 통합 모형(Modified Single Model, MSM), 권역 내에서 단일 모수를 추정하는 권역별 통합 모형(Group Model, GM), 관측 지점마다 개별적으로 모수를 추정하는 개별 지역 모형(Local Model, LM)까지 총 4가지 방식으로 모수를 추정, 모형을 개발하였다. 수정 통합 모형은 전국 통합 모형의 기준온도와 저온요구량, 고온요구량의 추정값과 동일한 값을 이용하여 각 지점의 만개일을 계산한 후, 지점별 보정계수 값을 합하였다. 수정 통합 모형의 보정계수는 전국 통합 모형에서 훈련자료를 대상으로 계산된 만개일의 지점별 편의(Bias) 값에 음의 부호를 취하여 각 지점별 보정계수로 정의하였다.

모형의 성능을 평가하기 위하여 평균제곱근오차(Root mean square error, RMSE) 편의(Bias), 피어슨 상관계수(Pearson correlation coefficient, r)의 지표를 사용하였다. 피어슨 상관계수는 프로그래밍 언어 R의 기본 패키지에 내장되어 있는 cor 함수를 활용하였다. 이를 바탕으로 훈련자료와 검증자료에 대한 전국 통합 모형, 수정된 통합 모형, 권역별 통합 모형, 개별 지역 모형의 성능을 비교하였다.

N：분석기간 동안의 관측년도 수

d _pi：i 번째 년도의 예측 만개일

d _oi：i 번 째 년도의 관측 만개일

d _oi：분석기간 동안의 평균 관측 만개일

만개시기에 영향을 크게 미치는 기간을 탐색하기 위해 26개 관측 지점을 대상으로 12년 간의 만개일과 월별 평균기온의 상관관계를 조사하였다. 그 결과, 12개 지점에서는 4월 평균기온이, 그 외 11개 지점에서는 3월 평균기온이 만개일과 가장 높은 상관계수를 보였다(p-value<0.05). 이에 따라, 26개 관측 지점의 과거 12년 간(2006년~2017년)의 만개일의 변화와 3, 4월 평균기온의 변화 경향을 권역별로 구분하여 분석하였다. 권역별 구분 기준은 지점별 일평균기온의 상관관계를 바탕으로 계층적 군집화를 통해 4개 권역으로 분류하였다(남서부(목포, 땅끝, 해남, 완도), 남동부(남해, 진주, 밀양, 청도), 중북부(춘천, 홍천, 충주, 횡성, 원주, 양구, 제천), 서울 도심(우장산공원, 서오릉, 강변북로, 관악산, 현충원, 북악산, 남산, 정릉, 홍릉, 태릉역앞, 어린이대공원)). 또한, 시간의 흐름에 따라 권역 간 만개일의 차이가 어떻게 변화하는지 조사하였다. 통계 분석법으로 맨-켄달 경향 분석법(Mann-Kendall trend test)을 이용하였으며(Xiang et al., 2013), R 패키지 zyp와 wql를 활용하였다.

모형 선택을 위해 각 모형별로 전체 지역에 대한 모형 성능, 권역별 모형의 성능, 26개 관측 지점의 오차(RMSE, 편의) 분포를 비교하였다. 오차의 분포 형태가 모형 간에 차이를 보이는지 검정하기 위하여 일원배치 분산분석(One-way ANOVA)을 수행하였으며, 모형 간 차이가 통계적으로 유의한 경우 사후분석으로 Tukey 검정(Tukey's test)을 시행하였다. 또한, Student의 T 검정법(Student’s T test)을 이용, 오차의 분포에 있어 훈련자료와 검증자료 사이에 통계적으로 유의미한 차이가 존재하는지 확인하였다. 비교해야하는 집단 간의 분포가 등분산성을 만족하지 못하는 경우, 비모수 검정법을 활용하였다(Kruskal-Wallis 검정, Wilcoxon 순위합 검정)(Chuine et al., 1998; Kalvans et al., 2015; Roberts et al., 2015). 모든 통계 분석은 프로그래밍 언어 R로 수행하였다.

과거 12년(2006년~2017년)동안의 아까시나무 만개일의 변화 및 봄철(3~4월) 평균 기온의 변화 경향을 4개 권역별로 분석하였다(Table 3, Figure 2). 분석 기간 동안 봄철 평균기온은 모든 지역에서 상승하는 경향을 보였다[Table 3(a), Figure 2(a)]. 서울 도심지역의 봄철 기온 변화의 기울기(Sen’s slope)는 0.2123℃/yr로 전체 지역 중 가장 컸으며, 분석 기간 중 기온이 가장 높았던 2014년 3, 4월 평균기온(9.72℃)은 기온이 가장 낮았던 2011년 3, 4월 평균기온(5.11℃)에 비해 4℃ 이상 높은 값을 나타냈다. 중북부 내륙 지역의 경우, 기온 변화의 기울기(Sen’s slope)는 0.2033℃/yr로 전체 지역 중 서울 도심 다음으로 두번째로 컸으며, 통계적으로 유의한 기온 증가 추세를 보였다(p-value=0.0327). 한편, 남서부 및 남동부 지역은 기온 변화의 기울기가 다른 두 지역의 50% 이하 수준으로 낮게 나타났으며, 특히 남서부 해안지역은 네 지역 중 기온 변화가 가장 적었다(Sen’s slope=0.0886℃/yr).

Figure 2. Changes of mean temperature in spring (March and April) and (b) the bloom dates of Robinia pseudoacacia by region over the past 12 years (2006~2017).

Download Original Figure

만개일의 변화 경향은 봄철 평균기온의 변화 경향과 달리 지역에 따라 상반된 양상을 보였다[Table 3(a), Figure 2(b)]. 만개일은 중북부 내륙 지역과 서울 도심 지역에서 점차 앞당겨지는 경향을 보인 반면(Sen’s slope＜0), 남부 지역에서는 관측 기간 동안 만개일 변화의 방향성이 뚜렷하게 나타나지 않았고, 남동부 지역에서는 오히려 변화의 기울기가 양의 값을 갖는 것으로 나타났다. 하지만 이러한 지역별 만개일 변화의 경향성은 관측 기간 내 만개일의 큰 연별 변화(annual variation)로 인해 통계적으로 유의하지는 않았다.

결과적으로, 지난 12년 동안 지역 간의 만개일 차이는 점차 감소하였다[Table 3(b), Figure 2(b)]. 중북부 내륙 지역과 남서부 해안 및 남동부 지역의 평균 만개일 차이는 분석 기간 동안 매년 0.7일 이상의 속도로 줄어들었는데, 이는 중북부 내륙 지역의 급격한 봄철 기온 상승에 의한 것으로 판단된다. 특히 남서부 해안 지역과 중북부 내륙 지역의 평균 만개일 차이는 12년의 비교적 짧은 기간에도 통계적으로 유의미한 감소를 보였다(p-value=0.0417). 2007년, 중북부 지역과 남서부 해안 지역의 평균 만개일은 각각 144.29 DOY(Day Of Year), 132.25 DOY로 두 지역간의 만개일은 평균적으로 12일 정도의 차이를 보였으나, 2014년 이후부터 두 지역 간의 평균 만개일 차이는 4일 이하로 지속되었다. 중북부 내륙 지역과 남동부 지역의 경우도, 2011년 까지 평균 만개일 차이가 10일 내외로 유지되었으나, 최근 들어 급격히 줄어드는 양상을 보였다. 2007년, 15.78일까지 차이를 보였던 두 지역의 평균 만개일은 2014년에 불과 3.21일만의 차이를 나타냈고, 2014년부터 2017년까지 3~5일 정도의 차이를 유지하였다.

여러 지역을 대상으로, 보다 정확도 높은 개화시기 예측을 위해 서로 다른 4가지 방식으로 모수를 추정하였다. 전 지역에 대해 단일 모수를 추정하는 SM에서 최적의 기준온도와 저온요구량, 고온요구량은 각각 8.88℃, -138.95 °D, 247.73°D로 결정되었다. 그리고 SM에 지점별 보정계수를 도입한 MSM, 권역 내에서 단일 모수를 추정하는 GM, 관측 지점마다 개별적으로 모수를 추정하는 LM에서 추정된 모수를 Table S1에 나타내었다.

MSM에서 지역별 보정계수는 대략 -11에서 7까지의 분포를 나타냈다. 양의 부호를 갖는 보정계수 중 절대값이 가장 컸던 지점은 남해였으며(보정계수=7.1429), 음의 부호를 갖는 보정계수 중 절대 값이 가장 컸던 지점은 관악산이었다(보정계수=-11.1429). 이는 SM이 남해 지점에서 만개일을 실제보다 7일 정도 빠르게, 관악산 지점에서 만개일을 실제보다 11일 정도 느리게 예측한다는 것을 의미한다. 남부지역의 경우 보정계수는 모두 0보다 컸으며, 남해를 포함하여 진주와 밀양은 5 이상의 보정계수 값을 가졌다.

GM에서, 서울 도심지역 및 남서부 지역은 기준온도가 5℃ 가까이 크게 감소하였고, 그와 동시에 고온요구량의 크기는 2배 이상 증가하였다. 남동부 지역은 모든 지역 중 추정 모수의 변화가 가장 적었다. 한편 LM 방식으로 추정된 모수는 지점별로 다양한 분포를 보였다.

본 연구에서는 지난 12년(2006~2017년)동안 관측된 아까시나무 만개일 자료를 활용하여 최근 들어 더욱 심화되고 있는 온난화로 인해 아까시나무의 개화 시기가 어떻게 변화하고 있는지 살펴보았다. Figure 2에 나타난 바와 같이, 시간에 따른 권역별 만개일의 변화 경향은 봄철(3, 4월) 기온의 변화 경향과 상응한다. 봄철 기온이 증가하면 만개일은 빨라지는 음의 상관관계를 보이며(Figure S2), 이러한 변화 경향은 분석 기간 동안 4개 권역에서 일치한다.

아까시나무 만개일의 변화 경향을 권역별로 조사한 결과에 따르면(Table 3), 봄철(3, 4월) 기온 상승률이 낮았던 남부 지역에서는 아까시나무의 만개일 변화 경향성이 뚜렷하게 나타나지 않은 반면, 남부 지역보다 2배 이상 높은 봄철 기온 상승률을 보였던 중북부 내륙 지역은 분석 기간 동안 네 지역 중 가장 빠른 속도로 만개일이 앞당겨졌다(Sen’s slope=-0.6508 days/yr). 그 결과 남부 지역과 중북부 지역에서 아까시나무가 만개하는 시기의 간격이 크게 단축 되었으며, 특히 남서부 해안 지역과 중북부 내륙 지역의 평균 만개일 차이는 통계적으로 유의하게 감소하였다(p-value=0.0417). 이와 유사하게, 봄철 개화종 7종을 대상으로 국내 14개 지점에서 관측된 개화일의 장기(1945~ 2007년) 시계열 변화를 분석한 선행 연구 결과에서, 해안 지역보다 내륙 지역의 봄철 식물계절 시기가 더 빠른 속도로 앞당겨져 해안, 내륙 지역간 개화일 차이가 줄어들었음이 밝혀졌다(Lee et al., 2009). 국외에서도 지역간 불균형한 기온 상승으로 인해 봄철 식물계절 시기의 지리적 변동이 감소한 사례가 보고된 바 있다. Ma et al.(2018)은 유럽에 분포하고 있는 온대 수종 16종의 1980~2013년 개화, 개엽 시기와 겨울 및 봄철 기온을 분석하여, 기온이 낮은 지역에서 온난화가 더욱 빠르게 진행되어 따뜻한 지역보다 개화, 개엽 시기가 더 앞당겨졌다는 사실을 확인하였다. 또한 Vitasse et al.(2018)은 고도가 다른 알프스 128개 지역에서 1960~2016년에 걸쳐 4수종의 개엽 시기가 어떻게 변화하는지 조사하였는데, 그 결과 고도가 높을수록 개엽 시기가 더 크게 앞당겨졌음이 드러났고, 이는 겨울과 늦은 봄의 높은 기온 상승으로 인한 것이라고 주장하였다.

한편, 서울 도심 지역의 경우 시간에 따른 봄철 평균 기온 및 만개일의 맨-켄달 경향 분석결과, 기온의 변화율은 0.2123℃/yr로 네 지역중 가장 높았던 반면, 만개일의 변화율은 -0.3028 days/yr로 중북부지역에 비해 50% 이상 낮았다(Table 3). 즉 서울의 3, 4월 기온 상승률은 중북부지역보다 높았지만 만개일의 변화율은 비교적 낮았고, 그 결과 남부 지역과의 만개일 차이 또한 크게 감소하지 않았다. 이는 중북부 내륙 지역이 다른 지역에 비해 아까시나무가 개화하는 시기가 늦어[Figure 2(b)], 개화 직전 및 개화 후 만개에 이르기까지의 기간인 5월 중의 기온 또한 중북부 지역의 만개 시기에 큰 영향을 미치기 때문인 것으로 판단된다. 실제로 분석 기간 중 중북부 지역의 평균 만개일은 140.71일(±6.96)로 서울 도심(137.09±5.93) 지역이나 남서부(134.63±3.61) 및 남동부(132.25±4.33) 지역에 비해 만개 시기가 평균적으로 4~8일 가량 늦었다. 특히 양구에서는 2009년, 2017년을 제외한 모든 년도에서 5월 25일 이후 만개한 것으로 관측되었으며, 2010년과 2011년에는 중북부 내륙 지역에 속하는 7개 지점 모두 5월 20일 이후에 아까시나무가 만개한 것으로 관측되었다. 또한, 맨-켄달 경향 분석을 통해 2006~2017년 간 5월 평균기온의 변화 경향을 권역별로 살펴보았을 때, 네 지역 모두에서 통계적으로 유의한 온도 상승이 나타났는데, 그 중 중북부 내륙 지역의 기온 상승률이 가장 높았다(중북부: +0.2027℃/yr (p-value=0.0019), 서울 도심: +0.1732℃/yr (p-value=0.0240), 남동부: +0.1030℃/yr (0.0240), 남서부: +0.1255℃/yr (p-value= 0.0087)).

우리나라의 기후변화 특성을 분석한 결과에 따르면(NIMS, 2018), 과거 106년(1912~2017년) 동안의 계절별 기온 상승률은 겨울(+0.25℃/10년), 봄(+0.24℃/10년), 가을(+0.16℃/10년), 여름(+0.08℃/10년) 순으로 높았으나, 최근 10년(2008~2017년) 동안 겨울철 기온 상승이 둔화되어 겨울 이외의 계절의 기온 상승이 더 뚜렷하게 나타났다. 또한 기온의 월별 변화를 보았을 때, 최근 30년(1988~2017년)에 비해 최근 10년 간의 기온은 5~8월에 상승, 12~1월에 하강, 2~4월에는 큰 차이가 없는 것으로 드러났다. 본 연구 결과에 의하면, 지역에 따라 5월 초에서 6월에 걸쳐 개화하는 아까시나무는 최근 12년(2006~2017년) 간의 짧은 관측 기간에도 불구하고 일부 지역 간의 만개일 변이가 크게 감소하였으며, 2014년 이후부터 이러한 양상이 더욱 심화되었다[Figure S3(a)]. 이러한 연구 결과는 앞서 언급했던 국내 기온 변화의 최근 동향이 반영된 것으로 보인다. Lee et al.(2014)는 1951~1980년, 1981~2010년 두 기간 동안 국내 6개 지역에서 관측된 개나리와 벚꽃의 개화시기 변화를 조사하였는데, 그 결과 수종 및 지역 별로 개화시기 변화의 폭이 달리 나타났다. 개화시기가 상대적으로 늦은 벚꽃의 경우, 개나리보다 개화시기가 크게 단축되었으며 관측 지점 간의 개화일 변이도 크게 감소하였는데, 이 또한 최근 한반도에서 관측되는 늦봄의 급격한 기온 상승으로 인한 것으로 해석된다(Kim et al., 2015; Lee et al., 2014).

전 지역에 대해 단일 모수를 추정하는 SM과 함께, SM에 지점별 보정계수를 도입한 MSM, 권역 내에서 단일 모수를 추정하는 GM, 관측 지점마다 개별적으로 모수를 추정하는 LM을 각각 개발하고 성능을 비교하였다. 훈련 결과, 전체 모형 중 LM이 가장 우수한 성능을 보였고, 검증 결과, 네 모형의 성능 차이가 거의 나타나지 않았던 남서부 지역을 제외하면 MSM과 LM이 다른 두 모형에 비해 양호한 성능을 보였다. 또한, Figure S3에 지역 평균 만개일의 표준편차와 모형의 성능을 비교한 결과를 연도별로 나타내었다. 지역 평균 만개일의 표준편차, 즉 만개일의 지역 편차[Figure S3(a)]는 2015년(1.58)에 가장 작았고, 2007년(6.97)과 2013년(6.97)에 가장 컸다. 2007~2017년의 기간 중 가장 최근 4년에 해당하는 2014, 2015, 2016, 2017년의 만개일 지역 편차는 나머지 7개 년도보다 작았다. 만개일의 지역 편차가 작았던 최근 4년 중 2017년을 제외한 3개 년도에서 LM의 성능이 가장 뛰어난 것으로 나타났다[Figure S3(b)]. 만개일의 지역 편차가 가장 컸던 2007, 2013년과 두번째로 컸던 2011년에는 MSM이 LM보다 더 좋은 성능을 보였다. 하지만 전체적으로 보았을 때, 특정 년도를 제외하면 LM과 MSM이 GM 및 SM에 비해 더 좋은 성능을 보였다. 특히, LM은 2012년과 2013년 두 해를 제외하면, GM과 SM보다 항상 더 좋은 성능을 보였다.

SM, MSM, GM, LM 중 비교적 좋은 성능을 보인 MSM과 LM을 선택하여 해당 모형의 예측 알고리즘을 전국 범위로 확대, 적용하였다. 알고리즘 적용 시, 위도 또는 고도가 높은 곳에 위치하여 다른 지역에 비해 아까시나무의 개화 시기가 상대적으로 느린 지역(5월 말~6월 개화 지역)까지 포괄해야 하므로, 봄철 평균기온으로써 3~5월 평균기온을 활용하였다.

MSM의 예측 알고리즘을 전국 범위로 확대하기 위하여, 모형을 구축한 26개 지점을 대상으로 분석 기간 중의 봄철 평균기온과 지점별 보정계수의 상관관계를 조사하였다. 그 결과, 두 변수는 양의 상관관계를 나타냈으며, 봄철 기온이 높은 지역일수록 보정계수가 증가하는 경향을 확인할 수 있었다[Figure 4(a)]. 이러한 상관관계를 바탕으로, 봄철 평균기온을 설명변수 x, 보정계수를 반응변수 y로 하는 선형회귀모형을 도출하였다. 모형의 적합도를 향상시키기 위해 R 패키지 car 을 활용하여 이상치 검정을 수행하였으며, 최종적으로 회귀직선의 방정식은 y=3.6034x‒40.8797로 추정되었다(결정계수 R²=0.7114)[Figure 4(a)].본 연구에서 모형의 입력자료로서 활용했던 과거 일기온 복원 래스터자료와 앞서 도출된 회귀식을 이용하여 MSM 보정계수의 전국 분포도를 1 km 해상도로 작성하였다. 작성된 보정계수 분포도를 활용하면 MSM을 전국적으로 적용하여, 아까시나무의 만개일을 l km 격자 단위로 예측할 수 있다.

한편, 26개 지점에서 독립적으로 구축된 LM을 전국 범위의 예측에 적용하기 위하여 공간보간(spatial interpolation)을 통해 모형 개발 이외의 지역에서도 예측 만개일이 산출되도록 하였다. 이를 위해 26개 지점에서 각각 추정된 모수를 이용하여 지점별 만개일을 계산한 후, 코크리깅(co-kriging)을 수행하여 만개일의 전국 분포를 추정하였다(Figure 5). 만개 시기에 영향을 미치는 해당 년도의 봄철(3~5월) 평균기온을 코크리깅의 이차변수로 적용하였으며, 과거 일기온 복원 래스터 자료를 활용하여 1 km의 공간 해상도로 예측을 상세화 하였다.

Figure 5. Spatial distribution of bloom dates of Robinia pseudoacacia in 4 years (2014~2017) estimated using co-kriging interpolation of Local Model (LM). (a) 2014, (b) 2015, (c) 2016, and (d) 2017.

Download Original Figure

MSM, LM의 예측 알고리즘을 전국 범위로 확대 적용 후, 최종적으로 모형이 개발된 이외 지역을 대상으로 모형의 예측력을 평가하기 위해 4년 간(2014~2017년) 16개의 추가 관측 지점에서 동일한 방법으로 관측된 만개일 자료를 활용하여, 외부 검증을 수행하였다. 그 결과, LM에 코크리깅 기법을 적용한 방법이 지점별로 보정계수를 추정하는 MSM보다 더 높은 정확도를 보였다(Table 5). 16개 지점 중 북강릉, 이천, 평택, 부안을 제외한 모든 곳에서, LM은 MSM보다 더 낮은 RMSE를 보였고 편의(Bias)는 0에 더 가까운 것으로 나타났다. 추가 관측 지점의 평균 RMSE는 LM이 MSM보다 1.78일 낮았으며(MSM：5.08, LM：3.30), 편의의 크기는 MSM이 LM에 비해 평균적으로 1.59일 더 컸다(MSM：-4.50, LM：-2.91). 또한 T 검정 결과, 두 모형의 RMSE, 편의의 분포는 통계적으로 유의한 차이를 보였다(RMSE: p-value=0.0118, Bias: p-value=0.0471). 즉, 개별 지역 모형은 수정 통합 모형보다 모형 개발 이외의 지역을 대상으로 더 나은 예측력을 보였다.