2.1 자료수집

본 연구에 이용된 자료는 한강물환경연구소의 “기후변화가 수생태계에 미치는 영향과 대응전략('10~'12)”, 환경부의 “환경생태유량 산정기준 연구 및 시범산정('14)”, 환경부⋅국립환경과학원의 “수생태계 건강성 조사 및 평가('12~'13)”, “하천 수생태계 현황 조사 및 건강성 평가('14~'21)”, “수생태계 참조하천 선정 및 활용방안 마련 연구('18)”, “4대강 보 개방에 따른 수생태계 변화 조사('18~'19)”, 국립환경과학원의 “생물측정망 모니터링 및 평가기법 개발연구('17)”의 지점 또는 정점 조사 결과였으며 총 조사지점 수는 6,787개, 총 표본단위는 23,846개였다. 대상 표본단위는 남북으로는 민통선 이남부터 제주도, 동서로는 서해 접경부터 동해 접경, 고도는 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)는 식 14와 같다. 식 14에서 전체 급구간까지()의 누적질량함수 값은 1이 된다.

(14)

C M F l = ∑ i = 1 l P M F i Δ x i = ∑ i = 1 l 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)의 해석을 위한 역분포 모형(reversed 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
Reversed 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
Reversed 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
Reversed 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
Reversed 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) 값을 적합시키는 방법(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 상대출현빈도

분류군별 출현고도의 범위는 동양하루살이-사할린하루살이 군이 1~700 m, 무늬하루살이는 5~871 m, 가는무늬하루살이는 7~846 m였다. 총 표본단위에서 도출된 상대출현빈도는 동양하루살이-사할린하루살이 군이 약 30%, 무늬하루살이는 약 7%, 가는무늬하루살이는 약 1%였다(Table 2). 이를 통해 볼 때 동양하루살이-사할린하루살이 군의 출현빈도가 무늬하루살이에 비해 높고 가는무늬하루살이에 비해서는 더욱 높다고 볼 수 있으나 표본단위의 수가 고도구간에 따라 차이가 크기 때문에 절대적인 비교는 어렵다. 고도의 급구간을 100 m로 설정하였을 때 동양하루살이-사할린하루살이 군은 100~200 m, 무늬하루살이는 300~500 m, 가는무늬하루살이는 600~900 m에서 상대출현빈도가 높았다.

Table 2. Number of sampling units, frequency and relative frequency of four species in genusEphemeraaccording to altitudinal gradient

Interval of altitude above sea level (m)	Sampling units ()	Frequency (relative frequency )
		E. orientalis-sachalinensis	E. strigata	E. separigata
0 – 100	15,521	4,597 (0.30)	431 (0.03)	53 (0.00)
> 100 – 200	4,585	1,681 (0.37)	530 (0.12)	52 (0.01)
> 200 – 300	2,102	652 (0.31)	302 (0.14)	37 (0.02)
> 300 – 400	938	224 (0.24)	226 (0.24)	21 (0.02)
> 400 – 500	319	47 (0.15)	73 (0.23)	10 (0.03)
> 500 – 600	256	25 (0.10)	53 (0.21)	11 (0.04)
> 600 – 700	125	5 (0.04)	25 (0.20)	13 (0.10)
> 700 – 800	83	-	14 (0.17)	5 (0.06)
> 800 – 900	19	-	1 (0.05)	3 (0.16)
Total	23,948	7,231 (0.30)	1,655 (0.07)	205 (0.01)

3.2 기술통계

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

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

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

3.3 연속 확률분포

상대출현빈도를 기반으로 한 이산적인 분석에 비하여 연속 확률분포 분석은 환경요인의 모든 값에 대하여 분류군의 출현특성을 연속적으로 해석할 수 있고 서식처 적합도를 확률기반의 함수식으로 제시할 수 있다는 점에서 이점이 있다. 본 연구에서 적용한 각 확률분포모형의 CDF와 실측자료를 기반으로 한 CMF 간의 NRMSE는 Table 4와 같다. 1 모수의 지수모형은 전반적으로 적합도가 낮았다. 동양하루살이-사할린하루살이 군과 무늬하루살이의 분포에는 Beta 분포모형의 적합도가 가장 높았으며, Weibull 분포모형의 적합도 역시 전반적으로 높았다. 가는무늬하루살이의 분포에는 역대수정규분포모형(reversed lognormal distribution model)의 적합도가 가장 높았으며 정적편포(positively skewed distribution) 모형들과, 비편포모형(정규분포, 로지스틱분포)은 적합도가 낮았다.

Fig. 3은 각 분류군의 확률질량함수와 최적의 적합도를 보이는 확률분포모형의 확률밀도함수 및 누적분포함수를 도시한 것이다. 고도의 급구간은 100 m로 동일하게 설정하였으나 출현빈도가 낮은 가는무늬하루살이는 700~900 m 구간에서 상대출현빈도의 값이 불연속적으로 나타나 병합구간으로 분석하였다.

Fig. 3. Cumulative mass function based on relative frequency of four species in genusEphemeraaccording to altitudinal gradient, probability density function and cumulative distribution function of the best-fit model.

동양하루살이-사할린하루살이 군과 가는무늬하루살이의 분포는 모두 Beta 분포를 따랐으나 동양하루살이-사할린하루살이 군은 상위고도로 치우치는 정적편포(positively skewed distribution)를 보이는 반면 무늬하루살이는 하위고도로 치우치는 약한 부적편포(negatively skewed distribution)를 보였다. 역대수정규분포를 따르는 가는무늬하루살이는 강한 부적편포를 보였다. 서론에서 언급된 바와 같이 본 연구에서 사용된 동양하루살이-사할린하루살이 군의 자료는 동양하루살이와 사할린하루살이를 포함한 결과인데 그 분포는 단봉형(unimodal)으로 나타나고 있어 두 종류가 동소서식하는 경향이 있음(^{Hwang et al., 2013})을 결과론적으로 입증하고 있다.

모든 분류군에 대한 확률밀도함수를 함께 도시하였을 때 각 분류군의 고도에 따른 출현구계는 더욱 뚜렷하게 대비되었다(Fig. 4). 동양하루살이-사할린하루살이 군과 가는무늬하루살이는 고도 450~500 m를 경계로 상대출현빈도가 뚜렷하게 구분되었으며, 무늬하루살이는 이 구간에서 최대의 출현도를 보였다. 하루살이속 분류군들의 고도에 대한 지위의 분할이 비생물적 환경요인에 의한 것인지 아니면 종간의 경쟁 또는 먹이망의 관계에 따른 것인지는 분명하지 않다.

Fig. 4. The best-fit probability density function of four species in genusEphemeraaccording to altitudinal gradient.

최적의 적합도를 보이는 확률분포모형에서 도출된 각 분류군의 평균과 중위값 및 표준편차는 Table 5와 같다. ^{Lee et al. (1995)}은 남양주군의 수동천에서 고도 40~350 m에 위치한 12개 지점을 대상으로 조사한 결과 동양하루살이-사할린하루살이 군이 출현하는 고도의 범위는 40~200 m, 평균고도와 최빈고도는 160 m였으며, 무늬하루살이가 출현하는 고도의 범위는 140~350 m, 평균고도와 최빈고도는 각각 279 m와 250 m였다고 보고한 바 있다. ^{Lee et al. (1995)}의 결과는 동양하루살이-사할린하루살이 군의 최빈고도를 제외하고 모두 본 연구의 결과에 비해 출현범위와 평균 및 최빈값이 크게 낮은 것이다. 이와 반면 ^{Bae (1995)}와 ^{Lee et al. (1999)}은 가는무늬하루살이가 대략 고도 500 m 이상에서 제한적으로 출현하는 것으로 보고하였는데, 이는 본 연구의 결과와 비교할 때 그 경향은 일치하지만 그 분포구계가 지나치게 상위고도로 치우쳐 해석된 것으로 판단된다. 하루살이속 분류군의 고도에 따른 분포구계에 대한 기존의 연구결과는 조사지점이 일부 지역에 국한되었다는 한계로 인해 각 분류군의 고도에 대한 고유의 서식지위를 파악하기에는 충분하지 않았던 것으로 보인다.

그러나 본 연구에서 사용된 자료에는 집중강우 시 빠른 유속에 의해 하상이 교란되어 상류로부터 부동(drift)된 부분이 포함되었을 가능성을 배제할 수 없으므로 가는무늬하루살이의 하위고도에 대한 분포에 대해서는 추후 연구가 필요하다. 부동의 정도가 어떠한지, 부동된 약충들이 하류 수체에서 성충으로 우화하여 생활사를 이어가는지에 대한 확인이 필요하다.

^{Kong et al. (2013)}은 고도와 위도에 따른 저서성 대형무척추동물의 분포를 바탕으로 한국온수생물지수(Korean thermality index)를 개발하면서 동양하루살이의 온수치(thermal value)를 2.9, 무늬하루살이는 1.0, 가는무늬하루살이는 0.4로 적용한 바 있고, 온수계열에 따른 온수계급치의 변이를 기반으로 지표가중치(indicator value)를 동양하루살이와 무늬하루살이는 1, 가는무늬하루살이는 4로 적용한 바 있다. 본 연구의 결과와 비교할 때 ^{Kong et al. (2013)}의 지표치는 같은 경향을 보이고 있으나 동양하루살이가 무늬하루살이에 비해 고도에 따른 분포의 표준편차가 낮은 결과로 볼 때 ^{Kong et al. (2013)}의 동양하루살이의 지표가중치가 상대적으로 과소평가된 것으로 보인다.

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

3.4 서식처 적합도 평가

Fig. 3의 확률밀도함수 값은 그 자체로 고도에 대한 각 종의 서식처 적합도 지수(habitat suitability index, HSI)로 볼 수 있다. 그러나 각 종별로 확률밀도함수 값은 최고 값이 다르기 때문에 상대적인 비교를 위해서는 계량화가 필요하다. 이에 대해서는 ^{Kong and Kim (2017)}이 가평천 저서성 대형무척추동물의 물리적 서식처 적합도 지수 산정 시 Weibull 분포의 최대치 확률밀도함수 값을 1로 계량화하는 방법을 제시한 바 있다. 마찬가지로 Fig. 3의 확률밀도함수의 최대값 함수를 구하고, 확률밀도함수를 최대값 함수로 나누면 고도에 대하여 최대값이 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. 5와 같다. HSI 값이 1에 해당하는 범위를 최적 적응범위라 본다면 Fig. 4의 50% 범위의 하한값(low)과 상한값(high)은 최적 적응범위의 경계값이라 할 수 있다. 이상의 방법으로 도출된 결과로서 고도에 대한 동양하루살이-사할린하루살이 군의 최적 적응범위는 42~257 m, 무늬하루살이는 335~644 m, 가는무늬하루살이는 641~824 m였다.

Fig. 5. Habitat suitability index of four species in genusEphemerafor altitudinal gradient.