2.1 Data collection

본 연구에 이용된 자료는 한강물환경연구소의 “기후변화가 수생태계에 미치는 영향과 대응전략('10~'12)”, 환경부의 “환경생태유량 산정기준 연구 및 시범산정('14)”, 환경부⋅국립환경과학원의 “수생태계 건강성 조사 및 평가('12~'13)”, “하천 수생태계 현황 조사 및 건강성 평가('14~'21)”, “수생태계 참조하천 선정 및 활용방안 마련 연구('15~'18)”, “4대강 보 개방에 따른 수생태계 변화 조사('18~'20)”, 국립환경과학원의 “생물측정망 모니터링 및 평가기법 개발연구('17~'18)”의 지점 또는 정점 조사 결과였으며 총 조사지점 수는 3,695개, 총 표본단위는 25,370개였다. 본 연구의 조사지점들은 한강, 낙동강, 금강, 섬진강, 영산강, 제주 수계를 포함해 전국에 걸쳐 분포하고 있으며, 0~888 m의 다양한 고도를 가진 서식환경을 포함하고 있다. 모든 조사 결과는 줄날도래속 유충의 출현유무와 해당 조사지점의 고도에 대한 것으로, 줄날도래속 성충에 대한 조사결과는 포함되지 않은 것이다.

2.2 통계분석

2.2.2 연속적 확률분포 통계

고도구배에 따른 분류군의 출현도를 연속적으로 해석하기 위해서는 일차적으로 상대출현빈도를 이산적인 확률질량함수(Probability Mass Function, PMF)로 변환하고 이차적으로 연속적인 확률밀도함수(Probability Density Function, PDF)로 변환하는 과정이 필요하다. 이때 확률질량함수 또는 확률밀도함수는 그 자체로서 환경요인에 대한 분류군의 출현특성 또는 서식처 적합도(Habitat Suitability)로 이해될 수 있다. 통계학적으로 상대도수의 합은 1이고 확률질량함수는 상대도수를 급구간으로 나눈 값이다. 그러나 식 5에서 도출되는 급구간별 분류군의 상대출현빈도는 해당 급구간에서의 상대값이므로 이를 급구간의 총수 급구간까지의 전체 급구간에 대하여 합한 값은 1이 아니다. 따라서 급구간별 상대출현빈도와 급구간의 범위값()을 곱하고 이를 합한 값으로 해당 급구간의 상대출현빈도를 나눈 값을 확률질량함수 값으로 정하였다(식 13).

식 13

P M F i = R i ∑ i = 1 k R i Δ x i

식 13은 고정된 분모의 값으로 상대출현빈도를 나누는 것이므로 급구간별 상대출현빈도가 확률질량함수의 값과 비례하게 된다. 따라서 급구간의 크기가 달라진다고 하더라도 해당 구간에서 분류군의 출현빈도의 특성은 과대 또는 과소 평가되지 않고 유지된다. 식 13에 따라 급구간 까지의 누적질량함수(Cumulative Mass Function, CMF)는 식 14와 같다. 식 14에서 전체 급구간까지()의 누적질량함수 값은 1이 된다.

식 14

C M F m = ∑ i = 1 m P M F i Δ x i = ∑ i = 1 m R i Δ x i ∑ i = 1 k R i Δ x i

2.2.2.1 확률분포모형 선정

본 연구에서 적용한 확률분포모형은 1 인자 모형(지수분포)과 2 인자 모형(정규분포, 대수정규분포, 로지스틱분포, Weibull 분포, Gamma 분포, Beta 분포, Gumbel 분포)의 8개 유형이었으며, 이 중 4개의 모형(지수분포, 대수정규분포, Weibull 분포, Gamma 분포)에 대해서는 부적편포(Negatively Skewed Distribution)의 해석을 위한 역분포 모형(Inverted Distribution Type)을 추가하였다(Table 1). 이 외에도 2 인자의 지수분포를 비롯하여 유형별로 3 인자의 일반화 모형들이 있으며 이러한 모형들은 2 인자의 특수 모형에 비하여 유연성이 커서 실측치에 대한 적합도가 높지만 모수 추정에 어려움이 있기 때문에 적용하지 않았다.

Table 1. Probability density function (PDF), cumulative distribution function (CDF), mode, mean, and variance according to environmental variable ()

Exponential distribution (Type 1.1)			Type 1.1.1		Type 1.1.2
			, ,		, ,
	CDF		,
	PDF
	Mode		0
	Mean
	Variance
Inverted exponential distribution (Type 1.2)			, ,
	CDF		,
	PDF
	Mode
	Mean
	Variance
Normal distribution (Type 2)			,
	CDF
	PDF
	Mode
	Mean		,
	Variance
Lognormal distribution (Type 3.1)			Type 3.1.1		Type 3.1.2
			, ,		, ,
	CDF
	PDF
	Mode
	Mean		,
	Variance
Inverted lognormal distribution (Type 3.2)		, ,
	CDF
	PDF
	Mode
	Mean	,
	Variance
Logistic distribution (Type 4)		,
	CDF
	PDF
	Mode
	Mean
	Variance	,
Weibull distribution (Type 5.1)		Type 5.1.1		Type 5.1.2
		, ,		, ,
	CDF	,
	PDF
	Mode
	Mean	: upper incomplete gamma function
	Variance
Inverted weibull distribution (Type 5.2)		, ,
	CDF	,
	PDF
	Mode
	Mean	: upper incomplete gamma function
	Variance
Gamma distribution (Type 6.1)		Type 6.1.1		Type 6.1.2
		, ,		, ,
	CDF
	PDF
	Mode
	Mean	: upper incomplete gamma function
	Variance
Inverted gamma distribution (Type 6.2)		, ,
	CDF
	PDF
	Mode
	Mean	: upper incomplete gamma function
	Variance
Beta distribution (Type 7)		Type 7.1		Type 7.2
		, , ,		, , ,
	CDF
	PDF
	Mode	, , , ,		, , , ,
	Mean	: hypergeometric function
	Variance
Gumbel distribution (Type 8)		, ,
	CDF
	PDF
	Mode
	Mean	: Euler’s constant Ei: Exponential integral
	Variance

지수분포는 최빈값으로부터 분류군의 출현도가 급격하게 지수적으로 감소하는 경우에 적용할 수 있다. 정규분포와 로지스틱분포는 분류군의 출현도가 최빈값을 중심으로 대칭성을 가질 때 적합성이 높고 대수정규분포, Weibull 분포, Gamma 분포, Gumbel 분포는 최빈값을 중심으로 편포할 때 적합성이 높다. 반면 Beta 분포는 일반적으로 다양한 확률구조에 대하여 유연성이 높은 것으로 알려져 있다.

각 확률분포모형에서 확률변수의 정의역과 실제로 생물이 분포하는 범위는 다를 수 있다. 일례로 정규분포나 로지스틱 분포에서 확률변수의 정의역은 에서 의 범위를 갖지만 변환되지 않은 환경요인의 값은 음의 영역을 갖지 않기 때문에 확률분포모형으로부터의 평균이나 분산은 환경요인이 실제로 가질 수 있는 범위에서 산출되어야 한다. 이를 고려할 때 생물의 분포유형은 환경요인의 범위가 또는 인 조건에서 더미(dummy)분포를 가진 경우(Fig. 1a)와 더미분포를 갖지 않는 경우(Fig. 1b), 양의 영역에서 또는 (Fig. 1c), 음의 영역을 포함한 , , (Fig. 1d)의 경우로 구분될 수 있다. 환경요인이 고도, 유속, 수질의 경우에는 그 값이 0인 조건에서도 생물이 출현할 수 있기 때문에 Fig. 1a의 유형을 따를 수도 있고 역치 값을 가진 Fig. 1b의 유형을 따를 수도 있다. 유역면적, 하폭, 수심의 경우에는 그 값이 0인 조건에서는 생물이 출현할 수 없으므로 Fig. 1c 또는 역치 값을 가진 Fig. 1b의 유형을 따를 수 있다. 반면 [= : 입경(mm)] 값으로 환산된 하상의 평균입경은 음의 값을 가질 수 있으므로 Fig. 1d의 유형을 따를 수 있다.

Fig. 1. Distribution type of aquatic organisms according to environmental factors. a) with dummy, b) without dummy and non-shifted, c) without dummy and shifted, d) transformed, without dummy and shifted.

따라서 환경요인에 따른 각 분류군의 분포특성에 따라 그 유형이 세분화될 수 있는데, 고도의 경우에는 Fig. 1a와 Fig. 1b의 유형으로 나뉠 수 있다. 해안가에 있는 수체에서 수심이 깊은 곳은 해수면 아래에 위치하여 있을 수도 있으나 본 연구에서는 육상 기준의 고도에 따른 분포에 관심을 두고 있으므로 해수면 아래의 고도에 대한 확률분포는 더미(dummy)구간으로 평가하였다. 이때 Fig. 1b의 경우에는 확률변수의 정의역에서 통계량을 구할 수 있으나 확률분포가 Fig. 1a의 유형으로 나타나는 경우에는 더미구간을 제외한 범위에서 통계량을 구하여야 하며 본 연구에서 이에 대한 함수를 도출하여 Table 1에 수록하였다.

2.2.2.2 모수 추정

각 분류군의 분포에 적합한 확률분포모형을 알고 있는 경우에는 모멘트법(Moment method), 최대우도법(Maximum Likelihood Method), 베이즈 추정법(Bayesian Method) 등으로 확률밀도함수의 모수를 추정할 수 있다. 그러나 각 분류군의 분포에 맞는 모형을 사전에 알 수 없기 때문에 선정된 확률분포모형을 모두 적용하여 이산적인 확률질량함수 값에 대하여 최대의 적합도를 보이는 모수 값을 추정하였다. 이때 확률질량함수 값에 확률밀도함수 값을 적합시키는 방법(PMF-PDF)과 누적질량함수 값에 누적분포함수(Cumulative Distribution Function, CDF) 값을 적합시키는 방법(CMF-CDF)이 있을 수 있는데, PMF-PDF 방법은 적합한 확률모형이 확정되기 전 단계에서는 각 급구간의 계급치(class mark)를 결정할 수 없기 때문에 적용이 어렵다. 따라서 본 연구에서는 CMF와 CDF의 표준화 평균 제곱근 오차(Normalized Root Mean Squared Error, NRMSE) (식 15)가 최소가 되게 하는 모수를 Microsoft Excel 프로그램의 해 찾기 기능을 이용하여 도출하였다.

식 15

N R M S E   % = 1 k ∑ i = 1 k C M F i - C D F i 2 1 k ∑ i = 1 k C M F i × 100

3.1 상대출현빈도

줄날도래속 4종에 대한 출현고도의 범위는 줄날도래가 1~849 m, 동양줄날도래는 1~871 m, 흰점줄날도래는 3~849 m, 두줄날도래는 12~626 m였으며, 종별 분포는 Fig. 2와 같다. 총 표본단위에서 도출된 상대출현빈도는 줄날도래가 31.6%, 동양줄날도래가 10.4%, 흰점줄날도래가 11.5%, 두줄날도래가 0.5%였다(Table 2). 이러한 결과를 통해 줄날도래의 출현빈도가 줄날도래속 유충 4종의 분류군 중에 가장 높다고 볼 수 있으나 표본단위의 수가 고도구간에 따라 차이가 크기 때문에 절대적인 비교는 어렵다. 고도의 급구간을 100 m로 설정하였을 때 줄날도래는 200~600 m, 동양줄날도래는 600~700 m, 흰점줄날도래는 300~500 m, 두줄날도래는 100~300 m에서 상대출현빈도가 높았다.

Table 2. Number of sampling units, frequencies and relative frequencies of four species in genus Hydropsyche according to altitudinal gradient

Interval of altitude above sea level (m)	Sampling units ()	Frequency (relative frequency )
		H. kozhantschikovi	H. orientalis	H. valvata	H. formosana
0 – 100	16,916	3,984 (0.24)	713 (0.04)	1,547 (0.09)	33 (0.00)
> 100 – 200	4,498	2,053 (0.46)	646 (0.14)	700 (0.16)	52 (0.01)
> 200 – 300	2,147	1,102 (0.51)	487 (0.23)	367 (0.17)	26 (0.01)
> 300 – 400	965	505 (0.52)	348 (0.36)	190 (0.20)	3 (0.00)
> 400 – 500	347	173 (0.50)	185 (0.53)	62 (0.18)	0 (0.00)
> 500 – 600	256	135 (0.53)	122 (0.48)	39 (0.15)	0 (0.00)
> 600 – 700	143	50 (0.35)	100 (0.70)	4 (0.03)	1 (0.01)
> 700 – 800	82	26 (0.32)	37 (0.45)	8 (0.10)	0 (0.00)
> 800 – 900	16	1 (0.06)	9 (0.56)	1 (0.06)	0 (0.00)
Total	25,370	8,029 (0.316)	2,647 (0.104)	2,918 (0.115)	115 (0.005)

Fig. 2. Distribution of four species in genus Hydropsyche.

3.2 기술통계

Fig. 3의 상자수염그림은 각 분류군이 출현한 표본단위의 고도값만을 가지고 도출된 것이며, 개체수 및 상대출현빈도는 고려하지 않았다. 상자수염그림에서 보이는 각 분류군의 중위고도, 평균고도, 3 사분위고도는 각 분류군이 최대의 상대출현빈도를 보이는 고도구간(Table 2)에 비하여 낮은 수준으로 나타났으며 상위고도 값은 대부분 이상치(outlier)로 분석되었다. 이는 하위고도에 조사단위의 수가 집중되어 있기 때문으로 이러한 분석은 실제의 분류군 출현특성을 왜곡시킬 수 있음을 보여주는 것이다.

Fig. 3. Box-whisker plot for altitudinal gradient based on occurrence of four species in genusHydropsyche.

분류군이 출현한 조사단위의 평균고도와 중위고도는 물론 출현밀도 가중 평균고도는 상대출현빈도 가중 평균고도나 중위고도에 비하여 현저히 낮았다(Table 3). 이는 설명한 바와 같이 조사단위가 하위고도에 집중되었기 때문에 나타난 결과로 조사단위의 수가 고도구간별로 동일하지 않다면 식 1~5의 방법으로 도출된 통계량은 각 분류군의 고도에 따른 서식처 적합성을 대표할 수 없다.

3.3 연속 확률분포

연속확률분포 분석은 상대출현빈도를 기반으로 한 이산적인 분석에 비해서 환경요인의 전체 값에 대하여 줄날도래속 분류군의 출현특성을 연속적으로 해석할 수 있으며, 확률기반의 함수식으로 서식처 적합도를 제시할 수 있다. 본 연구에서 적용한 각 확률분포모형의 CDF와 실측자료를 기반으로 한 CMF 간의 NRMSE는 Table 4와 같다. 1 모수의 지수모형은 줄날도래속 유충 4종의 적합도가 전반적으로 매우 낮았다. 줄날도래, 동양줄날도래, 흰점줄날도래의 분포에는 I-Weibull 분포모형의 적합도가 가장 높았으며, Beta 분포모형의 적합도 역시 전반적으로 높았다. 두줄날도래의 분포에는 Gamma 분포모형의 적합도가 가장 높았으며, 비편포모형인 정규분포와 로지스틱분포의 적합도가 낮았다.

줄날도래속 분류군의 확률질량함수와 최적의 적합도를 보이는 확률분포모형의 확률밀도함수 및 누적분포함수는 Fig. 4와 같다. 고도의 급구간은 100 m로 동일하게 설정하였으나 다른 종에 비해 출현빈도가 현저히 낮은 두줄날도래는 300~700 m 구간에서 상대출현빈도의 값이 불연속적으로 나타나 병합구간으로 분석하였다.

Fig. 4. Cumulative mass function based on relative frequencies of four species in genus Hydropsyche according to altitudinal gradient, probability density function, and cumulative distribution function of the best-fit model.

줄날도래와 동양줄날도래, 흰점줄날도래의 분포는 모두 I-Weibull 분포를 따랐으나 줄날도래는 약한 부적편포(negatively skewed distribution), 동양줄날도래는 하위고도에 편향하는 강한 부적편포를 보이는 반면, 흰점줄날도래는 정적편포(positively skewed distribution)를 보였다. Gamma 분포를 따르는 두줄날도래는 상위고도에 편향하는 강한 정적편포를 보였다.

Fig. 5는 전체 분류군에 대한 확률밀도함수를 하나의 그래프로 나타낸 것이며, 각 분류군의 고도에 따른 출현구계의 차이를 더욱 뚜렷하게 나타낸다. 동양줄날도래와 두줄날도래는 고도 350~400 m를 경계로 상대출현빈도가 뚜렷하게 구분되었으며, 줄날도래와 흰점줄날도래는 고도 400~450 m를 경계로 상대출현빈도가 구분되었다. 줄날도래속 분류군의 고도에 대한 지위의 분할이 비생물적 환경요인에 의한 것인지 아니면 종간의 경쟁 또는 먹이망의 관계에 따른 것인지는 분명하지 않다.

Fig. 5. The best-fit probability density function of four species in genusHydropsyche according to altitudinal gradient.

Table 5는 최적의 적합도를 보이는 확률분포모형에서 도출된 각 분류군의 평균고도와 중위값, 표준편차이다. 줄날도래의 평균고도, 중위값, 표준편차는 각각 408 m, 406 m, 214 m, 흰점줄날도래는 각각 397 m, 373 m, 226 m로 두 종간의 큰 차이가 없어 유사한 결과를 보여주었으며, 동양줄날도래는 각각 565 m, 589 m, 210 m, 두줄날도래는 각각 245 m, 229 m, 109 m로 줄날도래와 흰점줄날도래의 결과와는 뚜렷하게 구분이 되는 차이를 보였다.

^{Kong et al. (2013)}은 고도와 위도에 따른 저서성 대형무척추동물의 분포를 바탕으로 한국온수생물지수(Korean Thermality Index)를 개발하면서 줄날도래의 온수치(thermal value)를 2.2, 동양줄날도래는 1.0, 흰점줄날도래는 2.1, 두줄날도래는 2.5로 적용한 바 있고, 온수계열에 따른 온수계급치의 변이를 기반으로 지표가중치(indicator value)를 줄날도래속에 속하는 각 분류군을 1로 적용한 바 있다. 본 연구의 결과와 비교할 때 ^{Kong et al. (2013)}의 지표치는 같은 경향을 보이고 있다.

줄날도래속에 속하는 각 분류군의 고도에 따른 분포는 다양한 요인(수질, 수온, 먹이원 등)에 영향을 받을 수 있으나 이들이 같은 섭식 및 서식 기능군으로 분류된다는 점을 감안할 때 그 주요 원인은 수온으로 추정된다. 고도에 따라 수온이 직접적인 관계를 가지며(^{Kong et al., 2013}; ^{Madsen et al., 2015}), 본 연구의 결과를 바탕으로 고도를 기준으로 각 분류군의 호온성(thermophility)을 유추할 수 있다. 본 연구에서 나타난 분포고도의 최빈값과 표준편차로 볼 때 두줄날도래는 호온성(thermophilic), 줄날도래와 흰점줄날도래는 중온성(mesophilic), 동양줄날도래는 혐온성(thermophobic)으로 추정될 수 있다.

3.4 서식처 적합도 평가

확률밀도함수 값을 나타낸 Fig. 4는 그 자체로 고도에 대한 각 종의 서식처 적합도 지수(Habitat Suitability Index, HSI)로 볼 수 있다. 그러나 각 종별로 확률밀도함수 값을 상대적으로 비교하기 위해서는 계량화가 필요한데 이는 최고값이 다르기 때문이다. 이에 대해서는 ^{Kong and Kim (2017)}이 가평천 저서성 대형무척추동물의 물리적 서식처 적합도 지수 산정 시 Weibull 분포의 최대치 확률밀도함수 값을 1로 계량화하는 방법을 제시한 바 있다. 마찬가지로 Fig. 4의 확률밀도함수의 최대값 함수를 구하고, 확률밀도함수를 최대값 함수로 나누면 고도에 대하여 최대값이 1이 되는 연속적인 HSI 함수식을 만들 수 있다.

복잡한 함수식으로 이루어진 HSI보다는 더 단순하고 가시적인 HSI가 더 유용할 수 있는데, 미국의 Instream Flow and Aquatic Systems Group (^{IFASG, 1986})에서는 환경요인에 대한 생물의 서식처 적합도 지수(Habitat Suitability Index)를 전체 생물분포의 50%, 75%, 90%, 95% 범위에 대하여 각각 1.0, 0.5, 0.1, 0.05의 값으로 부과하는 방법을 제안하고 있다. 확률밀도함수를 알면 ^{IFASG (1986)}의 기준에 따라 최빈값을 중심으로 고도적응 범위에 대한 임계값을 정량적으로 쉽게 산출할 수 있으며 그 결과는 Fig. 6과 같다. HSI 값이 1에 해당하는 범위를 최적 적응범위라 본다면 Fig. 5의 50% 범위의 하한값(low)과 상한값(high)은 최적 적응범위의 경계값이라 할 수 있다. 이상의 방법으로 도출된 결과로서 고도에 대한 줄날도래의 최적 적응범위는 242~574 m, 동양줄날도래는 500~793 m, 흰점줄날도래는 174~485 m, 두줄날도래는 137~272 m였다.

Fig. 6. Habitat suitability index of four species in genusHydropsychefor altitudinal gradient.