Journal Search Engine

ISSN : 2671-9940(Print)
ISSN : 2671-9924(Online)

Journal of the Korean Society of Fisheries and Ocean Technology Vol.55 No.3 pp.206-216
DOI : https://doi.org/10.3796/KSFOT.2019.55.3.206

Comparative analysis of stock assessment models for analyzing potential yield of fishery resources in the West Sea, Korea

Min-Je CHOI

, Do-Hoon KIM^1*

, Ji-Hoon CHOI²

Student, Department of Marine & Fisheries Business and Economics, Graduate school, Pukyong National University, Busan 48513, Korea.
¹Professor, Department of Marine & Fisheries Business and Economics, Pukyong National University, Busan 48513, Korea
²Researcher, Coastal Water Fisheries Resources Research Division, National Institute of Fisheries Science, Busan 46083, Korea

Corresponding author: delaware310@pknu.ac.kr, Tel: +82-51-720-5954

Received 20190522 Review 20190624 Accepted 20190730

Abstract

This study is aimed to compare stock assessment models depending on how the models fit to observed data. Process-error model, Observation-error model, and Bayesian state-space model for the Korean Western coast fisheries were applied for comparison. Analytical results show that there is the least error between the estimated CPUE and the observed CPUE with the Bayesian state-space model; consequently, results of the Bayesian state-space model are the most reliable. According to the Bayesian State-space model, potential yield of fishery resources in the West Sea of Korea is estimated to be 231,949 tons per year. However, the results show that the fishery resources of West Sea have been decreasing since 1967. In addition, the amounts of stock in 2013 are assessed to be only 36% of the stock biomass at MSY level. Therefore, policy efforts are needed to recover the fishery resources of West Sea of Korea.

Key Words : State-space model , Bayesian inference , Observation-error model , Surplus production model , Maximum sustainable yield

서해 어획대상 잠재생산량 추정을 위한 자원평가모델의 비교 분석

최 민제

, 김 도훈^1*

, 최 지훈²

부경대학교 대학원 해양수산경영학과 학생
¹부경대학교 해양수산경영학과 교수
²국립수산과학원 연근해자원과 연구원

초록

키워드 :

This article has been cited by 0 article in crossref

Cited-By

Funding:

Ministry of Oceans and Fisheries

서 론

자원평가 모델 중 잉여생산량 모델(surplus production model)은 어획량과 어획노력량에 대한 시계열 자료만 존 재하면 최대 지속적 생산량(maximum sustainable yield, MSY)과 그에 상응하는 어획노력량을 추정할 수 있다는 장점이 있다. 그리고 계산절차의 간편함으로 수산자원의 평가에 널리 사용되고 있는 모델이다(Polacheck et al., 1993). 잉여생산량 모델은 자원동태 모델을 관측된 자료 를 이용하여 분석하는 방법에 따라 구분된다. 그 중 가장 널리 쓰이고 있는 모델은 과정오차(process-error) 모델 (Schnute, 1977;Clarke et al., 1992)과 관측오차(observationerror) 모델(Punt, 1990;Polacheck et al., 1993), 그리고 State-space 모델(Millar and Meyer, 2000;de Valpine and Hilborn, 2005;Winker et al., 2018) 등이다.

과정오차 모델은 자원동태 모델을 이용한 자원평가에 있어 관측치에서는 오차가 발생하지 않고 자원량의 변 화에서만 오차가 발생한다는 가정을 바탕으로 자원량을 추정하는 모델이다. 반대로 관측오차 모델은 자원동태 모델에서의 오차는 관측치인 단위노력당 어획량에서만 발생하며 자원량의 변화에서는 오차가 발생하지 않는다 고 가정하여 자원량을 추정하는 모델이다. 두 모델 중에 서는 통상적으로 관측오차 모델이 과정오차 모델에 비 해 보다 현실적인 추정치를 산출하는 것으로 평가되고 있다(Polacheck et al., 1993;Kwon et al., 2013;Kim et al., 2018). 하지만 두 모델 모두 관측오차와 과정오차 를 함께 고려하지 못한다는 한계가 있다.

이에 반해 State-space 모델은 과정오차와 관측오차를 동시에 고려할 수 있기 때문에 시간의 변화에 따른 자원 량의 변동을 추정하는 효과적인 방법으로 최근 널리 활 용되고 있다. State-space 모델은 마코브체인 몬테카를 로 기법을 이용한 베이지안 추론으로 수산자원을 평가 하는 방법으로 잉여생산량 모델 등을 활용한 자원평가 에 널리 사용되고 있다(McAllister et al., 1997;Meyer and Millar, 1999;de Valpine and Hastings, 2002;Chaloupka et al., 2007;Winker et al., 2018).

이들 모델 중 지금까지 국내에서 가장 널리 활용되고 있는 모델은 과정오차 모델(Zhang, 1992;Pyo, 2002)이 다. 비교적 최근에는 NOAA (National Oceanic and Atmospheric Administration)에서 제공하는 ASPIC (A Stock Production Model Incorporating Covariates) 프로 그램을 활용한 관측오차 모델(Kwon et al., 2013;Kim et al., 2018) 등도 사용되었다. 이에 반해 State-space 모 델을 활용한 자원평가 연구는 아주 제한적인 실정이다. 향후 보다 효과적인 자원평가를 위해서는 다양한 모델 들을 활용하고, 모델별 결과들을 비교해 가장 적합한 모델을 선택하는 것이 중요하다.

이러한 배경 하에서 본 연구에서는 Kim et al. (2018) 에서 서해 어획대상 잠재생산량의 추정에 사용된 어획 량과 표준화된 어획노력량 자료를 이용하여 과정오차 모델, 관측오차 모델, 그리고 State-space 모델의 결과를 비교ㆍ분석하고자 한다. 구체적으로 Kim et al. (2018) 의 연구에서는 Fox 모델, CYP 모델, ASPIC 모델, 그리 고 Maximum entropy 모델을 활용하여 서해 어획대상 잠재생산량을 추정하였다. 이에 본 연구에서는 관측오 차 모델과 State-space 모델을 추가하여 기존 모델 결과 들과 비교․분석함으로써 보다 유의한 자원평가 모델을 제시하고, 서해 어획대상 자원관리 및 회복 목표 설정에 있어 과학적인 의사결정에 기여하고자 한다.

재료 및 방법

분석 자료

관측오차 모델과 State-space 모델의 추정에는 Kim et al. (2018)의 서해 어획대상 잠재생산량의 추정에 사 용된 서해 어획대상 어종에 대한 어획량과 표준화된 어 획노력량 자료를 이용하였다(<Table-1> 참조). 본 연구 에서 사용된 자료는 1967~2013년까지 서해에서 조업 중인 어선들의 어획량과 총톤수(GT)이며, 연도별 어획 노력량(톤수)은 어군탐지기의 보급 또는 엔진출력의 증 가와 같이 톤수 자료에서 나타나지 않는 어획 능률의 변화를 고려하여 Fitzpatrick (1996)의 어업별 기술계수 를 활용하여 표준화된 자료이다(Kim et al., 2018).

자원동태 모델(잉여생산량 모델)

자원동태 모델은 현재 자원량의 수준을 과거의 자원 량과 해당 자원의 성장량 그리고 어획량의 관계로 설명 한다. 이 중 잉여생산량 모델은 자원의 가입량, 성장량, 그리고 자연사망률을 하나의 잉여생산 개념으로 합한 형태로 자원동태 함수는 다음과 같이 정의된다.

B_{y + 1} = B_{y} + g (B_{y}) - C_{y}

(1)

여기서, B_y 는 y년도 초의 자원량이고, g (B_y)는 자원 의 잉여생산 함수, 그리고 C_y 는 y년도의 어획량이다. 여기서, g (B_y)에 대한 가정에 따라 자원동태 모델의 형 태가 결정된다.

g (B_{y}) = r B_{y} (1 - \frac{B_{y}}{K})

(2)

g (B_{y}) = r B_{y} Ln (K) [1 - (\frac{ln (B_{y})}{ln (K)})]

(3)

g (B_{y}) = \frac{r}{p} B_{y} [1 - {(\frac{B_{y}}{K})}^{p}]

(4)

여기서, 식 (4)는 Pella-Tomlinson (1969)에서 고안된 보다 일반화된 형태의 성장함수로 r은 본원적 성장률, K는 환경수용력을 나타내며, 그리고 p는 성장함수의 형 태를 결정하는 매개변수이다. Pella-Tomlinson (1959) 함수는 p의 값의 변화에 따라 생산함수의 형태를 변화시 킬 수 있다. 예를 들어, p=1로 가정하면 식 (2)와 같은 형태가 되어 Schaefer (1954) 함수의 형태가 되고, p를 0에 근사시키면 식 (3)의 Fox (1970) 함수와 유사한 형태 가 된다. 위의 잉여생산 함수를 식 (1)에 대입하면 다음 과 같은 자원동태 함수를 도출할 수 있다.

B_{y + 1} = B_{y} + \frac{r}{p} B_{y} [1 - {(\frac{B_{y}}{K})}^{p}] - C_{y}

(5)

잉여생산량 모델에서는 어획능률과 자원량의 관계는 일정하다고 가정한다(Haddon, 2010). 이에 따라 자원의 관측치에 대한 함수는 다음과 같이 표현할 수 있다.

I_{y} = q B_{y}

(6)

여기서, I_y 는 자원의 상대적 풍요지수, q는 어획능률계 수를 나타낸다. 자원의 상대적 풍요지수는 단위 노력당 어획량(catch per unit effort, CPUE)으로 나타낼 수 있다. 식 (5)의 동태 모델은 과정오차 모델, 관측오차 모델, 그리 고 State-space 모델 등을 이용하여 추정할 수 있다.

관측오차 모델

관측오차 모델은 식 (1)의 함수에서는 오차가 발생하 지 않고 수산자원 평가에 있어서의 오차는 자원에 대한 관측치에서 발생한다고 가정한다(Polacheck et al., 1993;Haddon, 2010). 본 연구에서 사용된 관측오차 모 델과 ASPIC 모델의 주요 차이점은 본 연구에서는 관측 오차를 로그우도함수(log-likelihood function)를 사용하 여 확률적으로 고려했다는 점이다. 이를 함수식으로 나 타내면 다음과 같다.

{\hat{I}}_{y} = \frac{\hat{C_{y}}}{E_{y}} = q B_{y} e^{ε}, ε \sim ​ N (o; σ^{2})

(7)

여기서, e^ε 는 오차를 나타내고, 로그정규분포를 따른 다고 가정한다. 관측오차 모델 추정을 위한 우도함수는 다음과 같이 나타낼 수 있다.

L (d a t a | B_{0}, r, K, q) = \prod_{y} \frac{1}{I_{y} \sqrt{2 π} \hat{σ}} e^{\frac{- {(ln I_{y} - ln {\hat{I}}_{y})}^{2}}{2 {\hat{σ}}^{2}}}

(8)

{\hat{σ}}^{2} = \sum_{y} \frac{{(ln I_{y} - ln {\hat{I}}_{y})}^{2}}{n}

(9)

L L = - \frac{n}{2} (Ln (2 π) + 2 Ln (\hat{σ}) + 1)

(10)

여기서, $L (d a t a | B_{0}, r, K, q)$ 은 주어진 생물계수에서 나타나는 관측치의 우도(likelihood)이며, n은 관측된 단 위 노력당 어획량의 개수이다. 식 (8)의 함수식 양변에 로그를 씌워 정리하면 식 (10)과 같은 로그우도함수로 나타낼 수 있으며, 위 식을 최대화(maximum likelihood estimation)함으로써 관측오차 모델을 통한 생물계수를 추정할 수 있다.

관측오차 모델을 추정할 때 유의해야 할 점은 추정된 B₀ 이 K보다 높게 나타나는 비현실적인 경우가 발생한다 는 점이다. 이에 따라 B₀ 에 대한 추가적인 가정을 해야 한다. Punt (1990)에 따르면 관측오차모델을 추정할 때 B₀ 와 K를 동일하다고 가정하면 보다 좋은 결과치를 도출 할 수 있는 것으로 추정되었다. 이에 따라 본 연구에서도 관측오차 모델을 적용함에 있어 B₀ 를 K와 같다고 가정하 였다.

State-space 모델

State-space 모델은 베이지안 추론(Bayesian inference) 을 통해 식 (1)에서 발생하는 과정오차와 식 (6)에서 발 생하는 관측오차를 모두 고려하여 자원량을 추정할 수 있다. 베이지안 추론은 관측자료와 추정모수 모두에 확 률 분포를 가정하여 추론하는 방법으로 베이지안 추론 을 위해 식 (2)와 (6)을 재구성하여 자원량(B_y )을 자원량 의 비율(P_y = B_y/K)로 수정하여 추정하였으며, 추정에 있어 오차는 로그정규분포를 따르는 것으로 가정하였다 (Millar and Meyer, 2000).

\begin{array}{l} P_{1} | σ^{2} = e^{u_{0}} \\ P_{y} | P_{y - 1}, K, r, σ^{2} = (P_{y - 1} + r P_{y - 1} (1 - P_{y - 1}) - \frac{C_{y - 1}}{K}) e^{u_{y}} \end{array}

(11)

I_{y} | P_{y}, q, τ^{2} = q K P_{y} e^{υ_{y}}

(12)

베이지안 추론은 또한 관측된 자료와 사전정보를 결 합하여 사후분포를 추정하고 사후분포를 토대로 모수를 추정하는 방법이다. 우선 모수에 대한 사전분포는 다음 과 같이 나타낼 수 있다. 사전분포에 대한 가정은 Millar and Meyer (2000) 그리고 초기값에 대한 가정은 Kim et al. (2018)을 참고하였다.

\begin{array}{l} p (K, r, q, σ^{2}, τ^{2}, P_{1}, \dots, P_{N}) \\ = p (K) p (r) p (q) p (σ^{2}) p (τ^{2}) p (P_{1} | σ^{2}) \times \prod_{y = 1}^{N} p (P_{y + 1} | P_{y}, K, r, σ^{2}) \end{array}

(13)

그리고 생물계수가 특정한 값을 가질 때 관측된 단위 노력당 어획량이 발생할 확률은 다음과 같이 나타낼 수 있다.

p (I_{1}, \dots, I_{N} | K, r, q, σ^{2}, τ^{2}, P_{1}, \dots, P_{N}) = \prod_{y = 1}^{N} p (I_{y} | P_{y}, q, τ^{2})

(14)

베이즈 정리에 따르면, 사전확률 식 (13)과 우도함수 (14) 를 이용하여 식 (15)와 같은 사후분포를 추정할 수 있다.

\begin{array}{l} p (K, r, q, σ^{2}, τ^{2}, P_{1}, \dots, P_{N}, I_{1}, \dots, I_{N}) \\ = p (K) p (r) p (q) p (σ^{2}) p (τ^{2}) p (P_{1} | σ^{2}) \\ \times \prod_{y = 2}^{N} p (P_{y} | P_{y - 1}, K, r, σ^{2}) \prod_{y = 1}^{N} p (I_{y} | P_{y}, q, τ^{2}) \end{array}

(15)

추정된 사후분포로부터 각각의 모수를 산출하기 위해 서는 사후분포에 대한 적분계산이 필요하다. 그러나 위 와 같이 추정해야 하는 모수가 다차원적일 경우 이를 직접 계산하는 것은 거의 불가능하다. 이를 해결하기 위해 사용되는 방법 중 하나가 마코브체인 몬테카를로 (Markov Chain Monte Carlo, MCMC) 기법이다. 마코브 체인 몬테카를로는 변수들의 다차원적 사후분포를 통계 적으로 추정하는 방법으로, 추출된 변수들의 우도를 바 탕으로 사후분포를 확률론적으로 추정한다(Haddon, 2010;Kim, 2013).

마코브체인 몬테카를로는 매 단계에서 새로운 변수들 의 표본을 무작위로 생성한다. 그러나 다차원 변수에 대한 결합확률분포의 경우 랜덤표본을 생성하기 어려워 식 (16)과 같이 다른 변수들의 조건부확률분포로부터 변수 하나씩 순차적으로 표본을 생성하게 되는데 이를 깁스샘플링(Gibbs sampling)이라고 한다. 깁스샘플링을 이용한 마코브체인 몬테카를로 과정에는 사후분포를 추 정하는 다양한 알고리즘이 존재하는데, WinBUGS 프로 그램은 변수 각각의 분포형태에 따라 적절한 알고리즘 을 선정하여 사후분포를 도출한다(Gilks and Wild, 1992;Neal, 1997;Lunn et al., 2000).

\begin{array}{l} simulate K^{(1)} \sim f (K | r^{(0)}, \dots, τ^{2^{(0)}}, P_{1}^{(0)}, \dots, P_{n}^{(0)}, I_{1}^{(0)}, \dots, I_{n}^{(0)}), \\ simulate r^{(1)} \sim f (r | K^{(1)}, q^{(0)}, \dots, τ^{2^{(0)}}, P_{1}^{(0)}, \dots, P_{n}^{(0)}, I_{1}^{(0)}, \dots, I_{n}^{(0)}), \\ ⋮ \\ simulate I_{n}^{(1)} \sim f (I_{n} | K^{(1)}, r^{(1)}, \dots, τ^{2^{(1)}}, P_{1}^{(1)}, \dots, P_{n}^{(1)}, I_{1}^{(1)}, \dots, I_{n - 1}^{(1)}), \end{array}

(16)

위 과정을 m번 반복하면 각각의 변수에 대한 사후분 포를 식 (17)과 같이 추정할 수 있다.

\begin{array}{l} K (θ^{(m + 1)}, θ^{(m)}) = \\ \prod_{i = 1}^{n} f (K^{(m + 1)} | r^{(m + 1)}, \dots, τ^{2 (m + 1)}, P_{1}^{(m + 1)}, \dots, P_{n}^{(m + 1)}, I_{1}^{(m + 1)}, \dots, I_{n}^{(m + 1)}) \end{array}

(17)

깁스샘플링을 이용한 베이지안 추론을 위해서는 초기 값을 임의적으로 설정해 주어야 한다. 따라서 초기값이 정확하지 않을 경우, 최종적으로는 모델이 유의하게 도 출되더라도 추출과정의 초기 표본들은 사후분포에 수렴 하지 못하는 문제가 발생할 수 있다. 이러한 문제를 해결 하기 위해서 추출 표본 중 사후분포에 수렴하기 전에 추출된 표본들을 번인(burn-in) 과정을 통해 제거한 후 사후분포를 도출하면 초기값의 영향을 줄일 수 있다 (Bolker, 2008).

모델 검정

수산자원에 대한 모수를 파악하는 것은 거의 불가능 하기 때문에 분석 결과 어떠한 모델이 가장 정확한 추정 치를 반영하고 있는지를 평가하는 것은 상당히 어렵다. 이러한 문제를 해결하기 위한 대안으로 모델에서의 추 정치와 실제자료의 차이를 비교분석하는 방법이 널리 사용되고 있다(Kwon et al., 2013;Kim et al., 2018).

각 모델에 의해 추정된 연도별 단위 노력당 어획량 ( ${\hat{I}}_{y}$ )이 실제 관측된 단위 노력당 어획량(I_y )을 반영하고 있는 정도를 평가하기 위하여 식 (18)과 같이 오차의 제곱 평균 제곱근(root mean square error, RMSE)과 식 (19)의 결정계수(coefficient of determination, R²)를 계 산하였다.

RMSE = \sqrt{\frac{1}{n} \sum_{y} {(I_{y} - {\hat{I}}_{y})}^{2}}

(18)

R^{2} = 1 - \frac{\sum_{y} {(I_{y} - {\hat{I}}_{y})}^{2}}{\sum_{y} {(I_{y} - \bar{I})}^{2}}; \bar{I} = \frac{1}{n} \sum_{y} I_{y}

(19)

결과 및 고찰

관측오차 모델

관측오차 모델의 분석에 사용된 최대우도추정법의 경 우 로그우도 값이 클수록 추정치가 관측자료에 보다 잘 수렴되었음을 의미한다(Haddon, 2010). 관측오차 모델 추정 결과, Pella-Tomlinsom 함수를 가정한 경우 로그우 도값은 30.478로 나타났으며, Fox 함수의 경우는 30.479 로 나타났다. 이에 따른 생물계수는 Pella-Tomlinsom 함 수의 경우 r은 0.181, q는 1.31E-06, 그리고 K는 2,702,957톤으로 추정되었다. 그리고 Fox 함수의 경우 r은 0.180, q는 1.31E-06, 그리고 K는 2,708,372톤으로 추정되었다.

State-space 모델

State-space 모델의 추정은 깁스샘플링을 통한 베이지 안 추론에 주로 이용되는 프로그램인 WinBUGS (Bayesian inference Using Gibbs Sampling)를 활용하였 다. 분석에 있어서는 모델의 유의성을 확보하기 위해 1,250,000개의 표본을 추출하였고, 이 중 사후분포에 수 렴하지 못하는 초기표본들의 영향을 제거하기 위하여 초기 250,000개의 표본을 번인(burn-in)하였다. 그리고 추출 변수들의 자기상관을 제거하기 위해 매 50번째 표 본을 추출하여 최종적으로 20,000개의 표본을 이용하여 분석 결과를 도출하였다.

State-space 모델의 수렴 여부를 진단하는 방법 중 하 나는 변수들의 순차적인 추출결과를 나타낸 시도표 (trace plot)를 확인하는 것이다. 즉, 추정값들이 상승하 거나 하강하는 경향을 보이지 않고 백색잡음의 형태를 취하면 모델이 수렴한다고 판단한다. 서해 잠재생산량에 대한 State-space 모델의 추정 결과, Fig. 1에서 보는 바와 같이, 자원량(B)과 r, q, K 그리고 MSY 등 모든 변수들의 시도표가 특정한 패턴을 보이지 않는 것으로 나타났다.

또한 추정값들의 몬테카를로 오차가 작을수록 추정된 값이 보다 높은 정확도를 가지고 있다고 판단할 수 있다. 특히 추정값의 몬테카를로 오차가 사후표준오차의 5% 수준 미만일 경우 모델이 수렴한다고 판단할 수 있는데 (Spiegelhalter et al., 2003), State-space 모델의 추정 결 과, 모든 변수들의 몬테카를로 오차가 사후표준오차의 5% 수준보다 작은 것으로 나타났다(Table 2).

연도별 실제 단위 노력당 어획량과 사후분포에 의해 예측된 단위 노력당 어획량을 비교해 본 결과, 모든 단위 노력당 어획량의 실제 관측치가 사후 예측분포의 95% 구간 안에 포함되어 있는 것으로 나타났다(Fig. 2).

잉여생산 모델을 이용한 자원량 추정에 있어 유의해 야 할 점은 어획량과 단위 노력당 어획량의 자료만으로 다양한 생물계수들을 추정해야 되기 때문에 생물계수 간에 높은 상관관계가 발생한다는 점이다. 특히, 자원의 K와 r의 값이 높은 음의 상관관계를 갖게 되는데 K의 값이 과대 추정되면 r이 과소 추정되고, K가 과소 추정 되면 r이 과대 추정됨으로써 MSY는 상대적으로 정확한 값을 갖게 되지만 개별 생물계수의 추정에는 유의해야 한다(Hilborn and Walters, 1992).

State-space 모델의 경우 추정 모델의 자원량을 K에 대한 자원량의 비율로 재구성함으로써 생물계수 간의 상관관계를 감소시킬 수 있고(Millar and Meyer, 2000), 추정과정에서 설정한 생물계수에 대한 사전분포로 생물 계수 추정에 있어서의 불확실성을 감소시킬 수 있다 (McAllister et al., 2001). State-space 모델을 통해 추정 된 변수들의 사후분포는 Fig. 3과 같다. 95% 신뢰구간 범위는 r의 경우 0.212~0.570, K는 1,694,000~4,116,000톤, q는 9.40E-07~2.16E-06, 그리고 MSY는 194,400~269,500 톤의 범위를 갖는 것으로 나타났다.

분석 결과 비교

본 연구에서는 관측오차 모델과 State-space 모델로부 터 서해 어획대상 잠재생산량을 추정하였다(Table 3). 관측오차 모델에서 추정된 MSY는 Pella-Tomlinsom 함수를 가정한 경우 179,437톤 그리고 Fox 함수를 가 정한 경우는 179,291톤으로 나타났으며, Emsy는 각각 137,474 GT과 137,515 GT으로 추정되었다. 그리고 모 델 검정 결과, Pella-Tomlison 함수와 Fox 함수를 가정한 경우의 R²과 RMSE 값은 유사한 것으로 나타났다. 관측 오차 모델 분석 결과, ASPIC 모델보다는 관측 값의 변화 를 잘 반영하는 것으로 나타났으나 모델의 설명력은 최 대엔트로피 모델에 비해 떨어지는 것으로 나타났다.

State-space 모델의 R² 값은 0.999로 모형의 설명력이 가장 높은 것으로 나타났으며, RMSE 값 또한 0.038로 관측치의 변화를 가장 잘 반영하는 것으로 나타나 분석 모델 중 서해 어획대상 잠재생산량 분석을 위한 가장 적합한 모델로 평가되었다. State-space 모델에서 추정 된 MSY는 231,949톤 그리고 Emsy는 126,111 GT으로 나타났다.

결 론

본 연구에서는 서해 어획대상 자원의 지속적 생산량 추정을 위하여 어획량과 표준화된 어획노력량 자료(Kim et al., 2018)를 이용하여 관측오차 모델과 State-space 모델을 분석하였으며, 모델 분석 결과를 기존 연구결과 와 비교하여 보다 유의한 자원평가 방법을 제시하고자 하였다.

관측오차 모델은 수산자원의 평가에 있어 단위 노력 당 어획량과 같은 관측치에서만 오차가 발생한다고 가 정하여 수산자원의 지속적 생산량을 추정하는 방법으 로, 최우추정법(maximum likelihood estimation)을 통하 여 생물계수를 도출하고, 이를 바탕으로 MSY와 Emsy 를 산출하였다. 또한 보다 현실적인 추정치를 도출하기 위하여 추정 과정에서 B₀ 와 K가 동일하다는 가정과 관 측오차가 로그정규분포의 형태를 갖는다는 가정을 추가 하였다. 그리고 잉여생산 함수의 형태에 따른 영향을 추가적으로 고려하기 위해 Fox의 잉여생산 함수와 Pella-Tomlinsom의 잉여생산 함수의 두 가지 형태를 모 두 가정하여 추정하였다.

관측오차 모델 분석 결과, ASPIC 모델에 비해 관측치 변화에 대한 설명력이 더 높은 것으로 나타났다(Table 3). 따라서 B₀ 와 K가 동일하다는 가정과 관측오차에 대한 로그정규분포의 가정이 보다 효과적인 추정치를 도출하는데 유용한 것으로 나타났다. 반면에 잉여생산 함수에 대한 차이에도 불구하고 Pella-Tomlinson 함수를 이용한 관측오차 모델의 추정치와 Fox 함수를 이용한 관측오차 모델의 추정치는 크게 차이나지 않는 것으로 나타났다. 이는 Pella-Tomlinson 함수를 이용한 관측오 차 모델의 추정과정에서 Pella-Tomlinson 함수의 모형을 결정하는 p의 값이 0.0018로 추정되어 영(0)에 가까워짐 에 따라 상대적으로 Fox 함수와 비슷한 형태가 되었기 때문인 것으로 판단된다.

관측오차 모델이 ASPIC 모델에 비해 관측 자료를 잘 반영하는 것으로 평가되었지만, 최대엔트로피 모델과 비교했을 때는 어떤 모델이 보다 현실적인 결과치를 추 정했는지 판단하기 어려웠다. 반면, State-space 모델을 이용한 서해 어획대상의 잠재생산량 추정 결과는 Kim et al. (2018)에서 가장 적합한 것으로 평가된 최대엔트 로피 모델보다 높은 R² 값과 낮은 RMSE 값을 갖는 것으로 분석되었다(Table 3). 이에 따라 State-space 모델 이 서해 어획대상의 잠재생산량을 추정하는데 보다 유 용할 수 있는 모델로 평가되었다.

특히, 최대엔트로피 모델의 경우 MSY를 제외한 K, r, 그리고 q가 State-space 모델에서 추정한 각각의 생물 계수에 대한 95% 구간 안에 속하지 못하는 것으로 나타 났고, 다른 모델들과도 큰 차이를 보이는 것으로 나타났 다. 이에 따라 최대엔트로피 모델의 추정 결과 중 MSY 를 제외한 각각의 생물계수를 이용한 추정치 즉, 최대 지속적 자원량(Bmsy=459,762톤) 및 자원량 예측에 있 어서는 상당한 통계적 유의성 확보가 요구될 것으로 판 단된다.

State-space 모델을 이용한 서해 어획대상 잠재생산량 에 대한 MSY는 231,949톤으로 추정되었다. 그리고 서 해 어획대상의 연도별 자원량을 추정한 결과, 자원량이 1967년 이후 지속적으로 감소하여 2013년에는 최대 지 속적 생산을 위한 자원량(1,272,000톤)의 약 36% 수준 으로 추정되었다. 이를 95%의 구간으로 확대해도 여전 히 현재의 자원량이 MSY를 달성하기 위한 수준의 자원 량에 미치지 못하는 것으로 분석되었다(Fig. 4). 이에 따 라 State-space 모델에서 추정한 MSY를 현 시점에서 서 해 어획대상의 어획량 관리 수준으로 설정하면 현재의 자원량이 MSY를 달성하기 위한 자원량 수준으로 회복 하지 못하고 지속적으로 감소할 것으로 예상된다.

잉여생산량 모델을 활용한 수산자원의 평가에 있어서 생물계수에 대한 추정치의 해석에는 여전히 상당한 주 의가 필요하다. 특히 State-space 모델의 추정에는 쉐퍼 (Schaefer) 함수의 형태만 가정함으로써 자원의 성장함 수 형태의 보다 구체적인 변화는 고려하지 못하였다. 또한 모델 구조의 한계로 인해 어획량 자료의 오차까지 는 고려하지 못한 점이 본 연구의 한계로 남아 있다.

이러한 문제들을 해결하기 위해서는 향후 생물학적 조사가 보다 확대되어야 할 것이다. 이를 통해 보다 신뢰 가능한 사전분포의 자료를 확보하여 적용한다면 생물계 수들에 대한 보다 정확한 결과치를 도출할 수 있을 것이 다. 또한 관측오차 모델 간의 비교에서 Pella-Tomlinson 의 함수를 이용하면 모형의 변화까지 고려하여 보다 포 괄적인 형태의 자원평가가 가능한 것으로 나타나 차후 Pella-Tomlinson 함수를 활용하여 State-space 모델을 추 정한다면 보다 현실적인 자원평가가 가능할 것으로 기 대된다. 나아가 연령별 자료 등을 확보하여 연령구조 모델을 바탕으로 한 State-space 모델을 개발하여 평가 결과를 도출한다면 보다 적합한 결과치를 산출할 수 있 을 것이다. 마지막으로 어획량 자료에서 발생하는 오차 까지는 고려하지 못하는 자원평가 모델의 한계를 극복 하기 위한 노력이 더욱 필요할 것이다. 이를 위해 현재 사용하고 있는 어획량 자료 외에 직접적인 자원조사를 통한 자료를 이용하여 분석한다면 보다 정확한 수산자 원의 평가가 가능할 수 있을 것이다.

사 사

위 논문은 2019년도 정부(해양수산부) 재원으로 영남 씨그랜트센터의 지원을 받아 수행된 연구입니다. 그리 고 좋은 심사의견을 제시해 주신 익명의 심사위원님들 께 깊은 감사를 드립니다.

Figure

Fig. 1.

Trace plots of model parameters: Biomass (1967), Biomass (2013), Intrinsic growth rate (r), Catchability coefficient (q), Carrying capacity (K), and Maximum sustainable yield (MSY).

Fig. 2.

Observed CPUE(Kim et al., 2018) and the posterior predictive distribution of CPUE with the confidence intervals of 2.5% and 97.5% levels.

Fig. 3.

Posterior densities of model parameters: Biomass (1967), Biomass (2013), Maximum sustainable yield (MSY), Intrinsic growth rate (r), Catchability coefficient (q) and Carrying capacity (K).

Fig. 4.

Bmsy estimated and stock biomass trajectories with 95% confidence intervals (a), annual catch with MSY (b), and fishing effort with Emsy (c) for the Korean west coast fisheries by the State-space model from 1967 to 2013.

Table

Table 1.

Annual catch, standardized GT and CPUE for Korean west coast fisheries from 1967 to 2013 [data obtained from Kim et al. (2018)]

Table 2.

Standard error and MC error of model parameters: Biomass (1967), Biomass (2013), Intrinsic growth rate (r), Catchability coefficient (q), Carrying capacity (K), and Maximum sustainable yield (MSY)

Table 3.

Estimation results of potential yield of Korean Western sea. The estimation of Fox, CYP, ASPIC and maximum entropy models were quoted from Kim et al. (2018)

Reference

Bolker BM. 2008. Ecological models and data in R. Princeton University Press, 233-242.
Chaloupka M and Balazs G. 2007. Using Bayesian state-space modelling to assess the recovery and harvest potential of the Hawaiian green sea turtle stock. Ecological modelling 205, 93-109. (DOI:10.1016/j.ecolmodel. 2007.02.010).
Clarke RP , Yoshimoto SS and Pooley SG. 1992. A bioeconomic analysis of the Northwestern Hawaiian Islands lobster fishery. Marine Resource Economics 7, 115-140. (DOI:10.1086/mre.7.3.42629029).
de Valpine P and Hastings A. 2002. Fitting population models incorporating process noise and observation error. Ecological Monographs 72, 57-76. (DOI:10.1890/0012-9615(2002)072[0057:FPMIPN]2.0.CO;2).
de Valpine P and Hilborn R. 2005. State-space likelihoods for nonlinear fisheries time-series. Canadian Journal of Fisheries and Aquatic Sciences 62, 1937-1952. (DOI:10.1139/f05-116).
Fitzpatrick J. 1996. Technology and fisheries legislation. FAO Fisheries Technical Paper 350, 191-200.
Fox WW. 1970. An exponential surplus‐yield model for optimizing exploited fish populations. Transactions of the American Fisheries Society 99, 80-88. (DOI:10.1577/1548-8659(1970)99<80:AESMFO>2.0.CO;2).
Gilks WR and Wild P. 1992. Adaptive rejection sampling for Gibbs sampling. Applied Statistics, 337-348. (DOI:10.2307/2347565).
Haddon M. 2010. Modelling and quantitative methods in fisheries. CRC press, 285-333.
Hilborn R and Walters CJ. 1992. Quantitative fisheries stock assessment: choice, dynamics and uncertainty. Reviews in Fish Biology and Fisheries 2, 177-178. (DOI:10.1007/BF00042883).
Kim DH. 2013. Bayesian statistics using R and WinBUGS. Freedom Academy, 87-248.
Kim HA , Seo YI , Cha HK , Kang HJ and Zhang CI. 2018. A study on the estimation of potential yield for Korean west coast fisheries using the holistic production method (HPM). J Korean Soc Fish Ocean Technol 54, 38-53. (DOI:10.3796/KSFOT. 2018.54.1.038).
Kwon YJ , Zhang CI , Pyo HD and Seo YI. 2013. Comparison of models for estimating surplus productions and methods for estimating their parameters. J Korean Soc Fish Ocean Technol 49, 18-28. (DOI:10.3796/KSFT.2013.49.1.018).
Lunn DJ , Thomas A , Best N and Spiegelhalter D. 2000. WinBUGS-a Bayesian modelling framework: concepts, structure, and extensibility. Statistics and computing 10, 325-337. (DOI:10.1023/A:1008929526011).
McAllister MK , Pikitch EK and Babcock EA. 2001. Using demographic methods to construct Bayesian priors for the intrinsic rate of increase in the Schaefer model and implications for stock rebuilding. Canadian Journal of Fisheries and Aquatic Sciences 58, 1871-1890. (DOI:10.1139/f01-114).
McAllister MK and Ianelli JN. 1997. Bayesian stock assessment using catch-age data and the sampling-importance resampling algorithm. Canadian Journal of Fisheries and Aquatic Sciences 54, 284-300.
Meyer R and Millar RB. 1999. BUGS in Bayesian stock assessments. Canadian Journal of Fisheries and Aquatic Sciences 56, 1078-1087. (DOI:10.1139/f99-043).
Millar RB and Meyer R. 2000. Non‐linear state space modelling of fisheries biomass dynamics by using Metropolis‐Hastings within‐Gibbs sampling. Journal of the Royal Statistical Society: Series C (Applied Statistics) 49, 327-342. (DOI:10.1111/1467-9876.00195).
Neal RM. 1997. Markov chain Monte Carlo methods based onslicing' the density function. Department of Statistics, University of Toronto, Cananda, 1-27.
Pella JJ and Tomlinson PK. 1969. A generalized stock production model. Inter-American Tropical Tuna Commission Bulletin 13, 416-497.
Polacheck T, Hilborn R and Punt AE. 1993. Fitting surplus production models: comparing methods and measuring uncertainty. Canadian Journal of Fisheries and Aquatic Sciences 50, 2597-2607. (DOI:10.1139/f93-284).
Prager MH. 1994. A suite of extensions to a nonequilibrium surplus-production model. Fish Bull 92, 374-389.
Punt AE. 1990. Is B 1= K an appropriate assumption when applying an observation error production-model estimator to catch-effort data?. South African Journal of Marine Science 9, 249-259. (DOI:10.2989/025776190784378925).
Pyo HD. 2002. Determining appropriate bioeconomic models for stock assessment of aquatic resources. The Journal of Fisheries Business Administration 33(2), 75-98.
Robins CM, Wang YG and Die D. 1998. The impact of global positioning systems and plotters on fishing power in the northern prawn fishery, Australia. Canadian Journal of Fisheries and Aquatic Sciences 55, 1645-1651. (DOI: 10.1139/f98-037).
Schaefer MB. 1954. Some aspects of the dynamics of populations important to the management of the commercial marine fisheries. Inter-American Tropical Tuna Commission Bulletin 1, 23-56.
Schnute J. 1977. Improved estimates from the Schaefer production model: theoretical considerations. Journal of the Fisheries Board of Canada 34, 583-603. (DOI:10. 1139/f77-094).
Spiegelhalter D, Thomas A, Best N, and Lunn D. 2003. WinBUGS user manual. 1-60.
Winker H, Carvalho F and Kapur M. 2018. JABBA: Just Another Bayesian Biomass Assessment. Fisheries Research 204, 275-288. (DOI:10.1016/j.fishres.2018.03.010).
Zhang CI, KIM SA and YOON SB. 1992. Stock assessment and management implications of small yellow croaker in Korean waters. Korean J Fish Aquat Sci 25, 282-290.