서 론
인간에 의한 종의 이동과 확산은 인류 역사와 함께 진행되지 만(Barnosky et al., 2017), 외래종 침입과 확산에 따른 생태 적·경제적 손실에 대한 부정적 인식이 점차 높아지고 있다 (Levine et al., 2003;Jones et al., 2010). 새로운 환경에 정착 한 종은 가용 자원을 두고 토착종과 직접적으로 경쟁하며, 교란 체제 및 생태 순환을 변형시킬 수 있다(Mack et al., 2000;Jones et al., 2010). 종의 생태적 지위 및 선호환경은 토착 및 외래 이입종 모두 서로 다를 수 있기 때문에, 관리 효율성을 높이기 위한 포괄적 생태 특성(부정 및 긍정 영향)을 규정하고, 가능한 수준에서 관리를 실시하는 것이 필요하다. 우선적으로, 침입에 따른 주요 생태 과정(예, 군집 생체량 및 토착종의 적합 도)의 변화에 대한 이해, 그리고 합리적인 관리가 필요하다. 외래종 관리는 이입을 막거나 이입 초기 발견과 제거가 가장 효율적이지만, 관리 범위 전체(예, 행정 단위 또는 국가)에 대한 정착 및 확산 관찰에 수반되는 비용 및 인력 운용을 고려해야 한다(Rew et al., 2006). 외래종의 정착 가능 수준을 예측한 결과에 기반하여 선정된 중요 위험 지역을 집중적으로 관리하 는 것은 그러한 문제들의 대안이 될 수 있다 (Benjamin and Hiebert, 2004).
외래식물은 주로 극한 기상 및 기후 변화, 화재, 병해충, 토지 이용의 변화, 지역 간의 무역 및 운송, 인간의 여행 등의 요인에 의해 교란이 발생한 지역에 침입한다(Dix et al., 2010). 선형 경관 요소들은 외래종의 주요 침입 경로 및 핵심 서식지로 기능 한다. 이것은 지속적 관리에 따른 교란, 이에 따른 충분한 광량, 토착 경쟁종의 상대적 낮은 역할에 따른 자원 확보의 용이성, 지속적인 새로운 흙 유입 및 바람길 형성에 따른 종자 산포의 용이성, 탐방객 및 차량에 의한 종자의 운반 등에 기인한다 (Säumel and Kowarik, 2010;Joly et al., 2011;Meunier and Lavoie, 2012). 이러한 측면에서 산림생태계 관리를 위한 필수 기반 시설로서 임도 역시 외래종의 정착 및 확산 측면에서 중요하지만(Mortensen et al., 2009), 외래종의 확산 및 그에 따른 위험성 평가에서 상대적으로 적은 관심을 받고 있다. 인간 활동 및 토지 이용 변형과 관련 높은 철로(Wołkowycki and Banaszuk, 2016;Benedetti and Morelli, 2017), 도로 (Meunier and Lavoie, 2012;Dar et al., 2015;Benedetti and Morelli, 2017) 및 하천(Richardson et al., 2007;Park et al., 2019) 같은 주요 경관 요소는 외래종의 주요 정착 및 확산의 핵심 경로이다.
특히, 한국은 국립공원 및 산림유전자원보호구역 같은 대부 분의 보호구역이 산을 중심으로 지정되어 있다. 따라서 보호구 역 같은 중요 산림생물다양성 지역에서 외래종의 분포 및 그것 의 확산 위험성 예측은 보전 관리의 기초 활동일 것이다. 여러 가지 예측모형들의 장점과 단점을 비교하여 산 지형, 그리고 숲, 벌채지, 그리고 임도 및 사방댐이 어우러진 고유 특성이 반영된 효과적인 모형 선택과 활용은 이러한 활동의 효율성과 효과성을 제고하기 위한 기반이 된다.
최근 지리정보시스템(Geographic Information System; GIS), 위성영상기술(Remote Sensing; RS) 및 공간통계모형 (Spatial Statistics Models; SSM)의 발달에 따라 이를 활용하 여 외래종 관리전략을 수립하고자 하는 시도가 점차 늘어나고 있다(Holcombe et al., 2010;Jones et al., 2010;Crall et al., 2013). 이 같은 시도는 주로 수집된 종의 위치 자료 및 지역의 환경특성을 종합하여 통계적인 절차에 따라 출현 확률 을 예측하는 종분포모형(Species Distribution Model; SDM) 에 기반하여 수행된다(Elith et al., 2010). 현재 다양한 종분포 모형이 개발되어 이용되고 있으며, 각 모형은 종의 출현·비출현 지역의 환경적 공간을 통계적으로 해석하고 출현확률을 예측하 기 위한 서로 다른 알고리즘에 기반하고 있다(Guisan and Thuiller, 2005;Holcombe et al., 2010).
이용빈도 및 예측력이 높은 것으로 평가되는 대표적인 종분포 모형으로는 bioclimatic model(Bioclim), generalized linear model(GLM), maximum entropy model(MaxEnt) 등이 있다. 각 모형은 서로 다른 특성을 나타내는데, Bioclim 모형은 전통 적인 기후 범위 모형으로서 profile 방법으로 분류되며 알고리즘 을 이해하기 쉽기 때문에 개발된 지 오랜 시간이 지났음에도 이용빈도가 높다(Hijmans and Graham, 2006). 일반적인 생물 조사 자료(장소 및 종 구성)에 기반한 예측에는 회귀 모형 (generalized linear or additive models, GLMs or GAMs; ensembles of regression trees; random forests or boosted regression trees)이 주로 이용되고 있다(Elith et al., 2011). MaxEnt 모형은 정보이론을 바탕으로 한 기계학습방법 (machine learning method)을 이용한다(Elith et al., 2010;Elith et al., 2011). 활용되고 있는 모델들은 서로 다른 특성을 나타내므로, 여러 모형을 통해 종의 출현확률을 예측한 후 이를 비교하여 가장 높은 정확도(accuracy) 및 예측력(performance) 를 나타내는 모형을 현장의 이해와 관리에 적용하는 것이 필요 하다.
본 연구에서 우리는 보호구역 내 외래식물 분포 자료를 활용 하여 예측력이 좋은 종분포모델(Jones et al., 2010)을 선발하 고자 하였다. 이를 위해 우리는 울진 소광리 유전자원보호구역 내 출현 외래종의 분포를 조사하여 만들어진 외래식물의 분포 점 자료를 대상으로범용적으로 활용되는 세가지 종분포모형을 적용하여 외래식물종의 출현 예측 분석을 실시하였다. 그리고 모형들의 수행 결과를 비교하여 최적 모형을 선정하였다.
연구방법
1. 연구대상지
본 연구는 한반도 남부 울진 지역의 산림유전자원 보호구역 (2,274 ha)(N 36° 59’~ 37︒ 02’, E 129° 14’~129° 09’)을 대상으로 수행되었다(Figure 1). 조사 지역은 목재로서의 우수성 이 높은 우점종인 소나무(Pinus densiflora)의 보전을 위해 산림 유전자원 보호구역으로 지정되었으며, 소나무의 생산을 촉진하 기 위해 다양한 기반 시설이 존재하며, 벌채 및 숲가꾸기 활동이 활발하다. 이는 연구지역이 보호구역으로서 외래종 확산 위험성 평가 및 확산 저감을 위한 이상적 장소가 되는 배경이 된다. 연구 지역의 산림 경관은 대부분 소나무림(Pinus densiflora for. Erecta forest)과 신갈나무림(Quercus mongolica forest) 으로 구성되었고, 조사지역의 평균 고도는 약 690m, 모암은 화강편마암으로 이루어졌다. 연평균 강수량은 1,059.1mm, 그리 고 연평균기온은 12.5℃이다(Bae et al., 2003).
2. 종 분포 좌표 수집 및 환경변수 생산
2014년 7월에 수행된 현장 조사에서 7과 22종의 외래식물이 확인되었다(Table 1). 현장 조사 결과를 바탕으로 138개 출현 지점의 좌표와 임도로부터의 인접한 벌채지 및 숲으로의 침입 거리를 고려하여 62개의 무작위 좌표(random point)를 생성하 였다. 또한 출현 지점을 제외한 지역에서 300개의 비출현 지점 좌표를 무작위로 추출하여 분석에 활용하였다(Figure 2).
외래식물 출현위치의 예측 모의에 사용된 독립변수로서 환 경변수는 12종으로, 지형변수 6종(해발고도, 경사도, 사면방향, 지형위치지수, 지형습윤지수 및 수계로부터의 거리), 기후변수 3종(태양복사량, 열부하지수 및 지표온도), 식생변수 2종(정규 식생지수 및 영급), 그리고 인위교란변수 1종(도로로부터의 거 리)을 고려하였다. 활용된 변수들의 상관성은 산점도 행렬을 통해 평가하였다(Figure 3).
수계로부터의 거리 변수를 제외한 5종의 지형변수(해발고 도, 경사도, 사면방향, 지형위치지수 및 지형습윤지수), 그리고 지표온도를 제외한 2종의 기후변수(태양복사량 및 열부하지 수)는 1:5,000 수치지형도를 활용한 수치고도모형(Digital Elevation Model; DEM)을 이용하여 생산하였다. 수계로부터 의 거리 변수는 1:5,000 수치지형도로부터 추출한 수계 선자료 를 대상으로 유클리드 거리를 분석하여 추출하였다. 지표온도 및 정규식생지수 변수는 2016년 6월 Landsat 8 영상으로부터 추출하였다. 영급 변수는 산림청 제작 5차 임상도로부터, 도로 부터의 거리 변수는 소광리 유전자원보호지역 임도 및 탐방로 분포 주제도로부터 구하였다. 모든 변수는 UTM 좌표계를 이 용하였고, WGS84 좌표계로 재투영하여 분석에 활용하였다. 환경변수 구축 및 좌표 투영은 ArcGIS 10.2 프로그램을 이용 하였다.
3. 모형 선정 및 구축
우리는 대표적인 종분포 예측 모형으로서 BIOCLIM model (Bioclim), Generalized Linear Model(GLM) 및 Maximum Entropy model(MaxEnt)을 예측 효과성 평가에 활용하였다. Bioclim은 출현빈도가 높은 지역에서 나타나는 환경 변수값을 이용하여 유사 지역을 추출하며, 0~1 사이의 예측확률을 나타 낸다(Hijmans and Elith 2013). 회귀모델로서 GLM은 서수 최소자승 회귀식(ordinary least squares regression)을 정규화 하여 예측하는 모형으로서 종속변수에 대한 독립변수의 연결함 수식으로 표현된다(Equation 1). GLM의 연결함수 계산에서 이항(binomial) 분포, gaussian 분포, 그리고 poisson 분포를 적용할 수 있는데, 본 연구에서는 0(비출현)과 1(출현)의 이항 분포를 적용하였다. GLM의 연결 방식은 logistic, identity, 그리고 log가 있으며 본 연구에서는 logistic 방식을 이용하였 다(Equation 1). MaxEent는 최대 엔트로피(최대 분포 가능 지역)를 추정하여 분포를 예측한다(Phillips et al., 2006). 활 용된 모형 중 GLM은 구축 시 출현지점 좌표와 비출현지점 좌표가 모두 필요하며(Presence/ Absence 방식), Bioclim과 MaxEnt는 출현지점 좌표(Presence only 방식)만을 활용한다.
모형 적용 및 검증, 그리고 일원분산 분석은 R Studio (RStudio Team, 2020) 및 dismo package를 이용하였다(Hijmans and Elith, 2013). 모형은 최적화(model fitting) 및 예측(prediction) 단계를 거쳐 생성하였다. 전체 종 표를 이용한 분석 시 출현 및 비출현 좌표를 5그룹으로 나눈 후(5-fold), 최적화 및 추정 을 다섯 차례 반복하여 도출된 결과를 생성하였다. GLM의 경우, 유의확률이 낮은 환경변수를 순차적으로 제거하면서 (stepwise removing) 재생성한 모형으로부터 산출된 AIC (Akaike information criterion)를 비교하여 가장 유의성이 높 은 모형을 선택하였다(Jones et al., 2010). AIC는 통계적 모형 의 질의 측정치로 활용되며 낮은 AIC를 나타낼수록 모형의 질 이 높다고 평가된다(Akaike, 1974). 종 별 GLM 모형 구축 시 차후의 모형 비교 단계를 위해 변수 제거 단계를 2단계로 한정하 였으며, 이 중 가장 낮은 AIC를 기록한 모형을 선택하였다.
예측된 확률을 출현지역과 비출현지역의 이진(binary) 형태 로 분류하기 위한 임계치(threshold value)로는 ROC(receiver operating characteristic) 곡선에서 특이성(specificity; true positive rate)와 민감성(sensitivity; true negative rate)의 합 이 최대를 기록하는 지점의 값을 선정하였다(Liu et al., 2013;Choe et al., 2016).
Equation 1. Generalized linear model의 확률추정식(p : 예측확률, var : 환경변수, a : 계수).
4. 모형 검증 및 비교
모형의 예측 적합도는 ROC 곡선 이하 지역의 넓이인 AUC(Area Under Curve)를 통해 평가하였다; AUC는 0~1의 분포를 나타내며 1에 가까울수록 적합함을 의미한다(Phillips et al., 2006). 모형 별 생성 결과의 비교는 통계적 방법과 현상학 적 방법을 이용하여 수행하였다. 통계적 비교는 일원배치 분산 분석(one-way ANOVA)을 활용하였고, 사후검증은 Tukey’s method를 이용하였다. 현상학적 비교를 위해 최종 산출된 예측 도의 현실성을 현장 조사 결과를 통해 검증하였다.
결과
1. 환경변수 별 출현 양상
현장 조사 결과, 외래식물은 해발고도는 442.2~1137.4m의 분포 범위를 보였으며, 상대적으로 해발고도가 낮은 지역에서 주로 출현하였다(Figure 4a). 전체 경사도 분포(0~69.5°) 중, 외래종은 20°이하의 비교적 평탄한 지역에서 출현빈도가 높았 다(Figure 4b). 외래식물이 출현하는 사면방향은 약 300°정도 의 서북사면에서 가장 높은 빈도로 출현하였다(Figure 4c).
연구지역의 지형위치지수는 –9.2~11.6의 분포를 보였으며, 외래종은 –1.5~1.5 사이의 비교적 평탄한 지역에서 출현빈도 가 높았다(Figure 4d). 연구지역의 지형습윤지수는 –1~24.5의 분포를 나타냈으며, 외래종은 2~6 사이에 정규형의 빈도 분포 를 나타내었다(Figure 4e). 연구지역은 수계로부터 최대 286m 떨어져 있었고, 수계로부터의 거리가 가까운 곳에서 외래종의 출현빈도가 높았다(Figure 4f).
연구지역은 임도 및 탐방로를 포함하는 도로에서 최대 1,540m 떨어져 있었으며, 외래종은 도로 인근에서 90% 이상 출현하였다(Figure 4g). 연구지역에서 태양복사량은 161,755~ 1,666,620WH m-2이었으며, 복사량이 많은 지역에서 외래종 의 출현빈도가 높았다(Figure 4h). 열부하지수는 0.16~1.11의 분포를 나타냈으며, 외래종 출현빈도는 열부하량과 특별한 상 관관계를 보이지 않았다(Figure 4i). 지표온도는 17.1~24.8℃ 의 분포를 나타냈으며, 20℃ 내외에서 가장 높은 출현빈도를 보였다(Figure 4j).
연구지역에서 정규식생지수는 0~0.76의 분포를 나타냈으 며, 외래식물은 중간 범위(0.4 내외)에서 가장 높은 출현빈도를 나타냈다(Figure 4k). 식물군락을 구성하는 우점종의 영급은 0(식물 없음)~6(51~60년)의 분포를 나타냈으며, 외래종은 0,2,3,4 영급 지역에서 출현하였다(Figure 4l).
2. 외래종 예측 모형
1) Bioclim
Bioclim 모형을 통해 예측한 외래종의 출현확률은 0~0.66 으로 나타났으며, 임계치는 약 0.017, AUC는 0.838을 기록하 였다(Figure 5-6).
2) GLM
GLM 모형 구축 시 사면방향(A), 해발고도(D), 지표온도 (Th), 경사도(Sl), 영급(Y) 변수를 제외한 모형에서 가장 낮은 AIC를 나타냈다(Table 2). 최종 모형에서 가장 높은 유의확률 을 나타낸 환경변수는 도로로부터의 거리와 정규식생지수였으 며(p<0.001), 수계로부터의 거리 및 태양복사량 역시 뚜렷한 상관성(p<0.01)을 나타냈다(Table 3). GLM에서 외래식물의 출현확률은 0~1의 값을 나타냈고, 임계치 및 AUC는 각각 0.52 및 0.976으로 계산되었다(Figure 5-6).
3) MaxEnt
MaxEnt에서 예측된 외래식물의 출현확률은 0~1의 분포를 나타냈고, 임계치 및 AUC는 0.14 및 0.979를 기록하였다 (Figure 5-6). 모형에 사용된 환경변수의 기여도는 도로로부터 의 거리, 정규식생지수, 경사도, 지표온도, 그리고 해발고도 순 으로 높았다.
3. 최적모형선정
각 모형으로부터 계산된 평균 AUC(각 모형 별 n=9)는 뚜렷 한 차이가 관찰되었고(F=14.96, p<0.001), 사후 검정한 결과 에서 Bioclim과 다른 두 모형은 차이를 보였으며, GLM과 MaxEnt 사이에는 유의한 차이가 없었다(Table 3). 세 모형 중 Bioclim로부터 도출된 AUC의 표준편차가 가장 컸으며, 이를 통해 표본 수에 따라 예측력의 차이가 가장 큰 모형은 Bioclim 인 것으로 나타났다. GLM의 표준편차는 0.02로 가 장 낮아 표본의 수에 따라 예측력의 차이가 크게 좌우되지 않는 일관성 높은 모형(robust model)으로 판단된다.
각 모형으로부터 도출된 평균 임계치(각 모형 별 n=9)를 유 의한 차이가 나타났다(F=7.44, p<0.01)(Table 3). 사후 검정 결과에서 Bioclim과 다른 두 모형은 차이를 보였지만, GLM과 MaxEnt 사이에는 뚜렷한 차이가 없었다(Table 3). 세 모형 중 MaxEnt 임계치의 표준편차가 가장 크게 나타나, 입력된 종에 따라 임계치가 가장 크게 달라지는 경향을 보였다.
모델에서 추출된 출현예측지역의 평균 면적(각 모형 별 n=9) 은 뚜렷한 차이가 있었다(F=14.14, p<0.001)(Table 3). 사후 검정 결과에서 GLM과 다른 두 모형은 차이를 보였지만 Bioclim과 MaxEnt에서는 유의한 차이가 없었다(Table 3). 세 모형 중 GLM으로부터 도출된 면적의 표준편차가 가장 컸으 며, 반면 MaxEnt의 출현예측지역 변동폭(SD=1.76)이 가장 작았다.
고찰
1. 보호구역 내 외래종 선호 환경
연구지역에서 확인된 모든 외래식물은 임도 및 탐방로, 그리 고 관리된 수계 같은 선 모양의 기반 시설에서 주로 관찰되었 다. 소광리 일대에서 출현한 외래종은 벌채지 일부를 제외하면, 숲 내부에서는 관찰되지 않았다. 보호구역의 숲 내부에서 외래 식물의 침입이 관찰되지 않은 것은 해당 종을 비롯한 대부분의 외래식물(77.3%, 17종)이 고립된 숲 내부의 환경보다 개발지 및 통로형 서식지를 선호하는 것에 기인한다(Table 1). 특별히, 연구지역에서 높은 빈도로 출현한 외래식물(개방초, 달맞이꽃, 서양민들레 등)이 상대적으로 긴 정착 시기를 가졌고, 모두 높 은 귀화도를 나타낸 것을 고려할 때(Table 1), 현 상태에서 숲 내부로의 외래식물 침입 가능성은 낮은 것으로 판단된다.
수관 닫힘 단계(closed canopy stage) 이전의 벌채 후 천이 초기 서식지들은 잔존하고 있는 지표면 교란 및 높은 광량 같은 교란 환경을 선호하는 식물(예를 들어, 천이초기 환경 선호 식 물)의 정착에 유리한 조건을 형성한다(Jung et al., 2017;Kim et al., 2017). 벌채에 따른 천이초기 환경은 토착종과 외래종 모두에게 정착 및 번성할 수 있는 기회가 되는데, 두 식물 기능 군 모두 자연 과정(매토종자 및 바람 등) 및 인간 도움(신발, 트랙터 같은 작업 장비 등)에 의해 유입될 수 있다. 따라서, 임도 및 사방댐 같은 산림생태계 관리의 기반 시설을 중심으로 자생종 위주의 초지 조성을 통한 외래식물 저감 활동이 필요하 다. 특히 숲가꾸기를 포함한 벌채지들은 외래식물이 이입하는 문으로 작용하는 임도와 접하는 구간을 최소화하는 활동이 필 요하다. 예를 들면, 벌채지와 임도가 접하는 면적을 최소화하 고, 벌채지 내부로 외래종이 유입되는 것을 감소시키기 위해 완충 역할의 숲 띠를 임도를 따라 잔존시키는 방법을 고려할 수 있다. 임도를 따라 잔존시키는 숲 띠는 국화과 같은 특히 풍부한 바람 산포 외래식물에 대한 여과지 역할이 기대된다. 더불어, 보호구역 내 온전한 산림 군반의 크기를 유지하고, 외 래식물의 정착시기와 귀화도를 고려한 외래식물 유형별 관리가 필요하다(Latombe et al., 2018).
2. 외래종 분포의 최적 예측 모형
통계 검정 결과, GLM과 MaxEnt는 Bioclim에 비해 상대적 으로 높은 예측력과 정확도를 보인 반면, Bioclim은 연구지역 같은 보호구역에 적용하기에 적합하지 않았다. 전체 종을 대상 으로 한 분석에서 예측된 출현지역은 도로를 중심으로 분포하 는 경향이었는데, Bioclim, GLM, 그리고 MaxEnt 순으로 출 현예상지역을 넓게 예측하였다(Figure 7).
통계적으로는 유의하지 않았지만, MaxEnt의 경우 모든 유 형의 분석에서 유사한 분포 예측 결과를 도출한 반면, GLM은 종에 따라 분포 예측 결과가 상이한 경향을 보였다. 특히, MaxEnt는 표본 수와 상관없이 가장 일관성 높은 모형이었다 (Wisz et al., 2008). Bioclim은 표본 수에 따라 예측력에서 큰 차이를 보였으며, 미국가막사리와 큰김의털의 경우 현장조 사로부터 종의 서식이 확인되었음에도 출현 예측 지역이 나타 나지 않았다. GLM의 결우, 통계적 검증에서는 높은 AUC 값 을 나타냈지만, 자료의 수가 적은 종에서는 다수 지역이 null data 형태로 추출되어 모형이 제대로 생산되지 않는 경향이었 다. 뿐만 아니라, GLM은 분포 좌표점 수가 적은 경우(n<100) 예측력이 10% 이상 감소하는 것이 보고된 바 있다(Wisz et al., 2008). 이러한 결과들을 종합하였을 때, 본 연구에 사용된 모형 중 연구지역의 특성이 반영된 외래종 분포 예측을 위한 최적 모형은 MaxEnt인 것으로 평가할 수 있다.
우리의 모형 효율성 검증은 보호구역 내 외래종 출현 자료를 활용한 예측 모델을 활용하였다. 모형 연구는 생태계 관리를 위한 실증적 검증(즉, 자연계에서 모형의 예측 효과가 실제 나 타나는 지에 대한 검증)이 거의 불가능한 것이 한계로 작용한 다. 보호구역 내 정밀 생물종 분포 자료를 이용한 우리의 모델 선발은 산림생태계로 보호구역의 보전 관리에 도움이 될 것이 다. 또한 지역 특성이 반영된 효율적이고 정교한 모델 발굴은 후속 생물종 확산 모형 연구에 도움이 될 것이다.