A State-space Production Assessment Model with a Joint Prior Based on Population Resilience: Illustration with the Common Squid Todarodes pacificus Stock

자원복원력 개념을 적용한 사전확률분포 및 상태공간 잉여생산 평가모델: 살오징어(Todarodes pacificus) 개체군 자원평가

Abstract

It is a difficult task to estimate parameters in even a simple stock assessment model such as a surplus production model, using only data about temporal catch-per-unit-effort (CPUE) (or survey index) and fishery yields. Such difficulty is exacerbated when time-varying parameters are treated as random effects (aka state variables). To overcome the difficulty, previous studies incorporated somewhat subjective assumptions (e.g., B₁=K) or informative priors of parameters. A key is how to build an objective joint prior of parameters, reducing subjectivity. Given the limited data on temporal CPUEs and fishery yields from 1999-2020 for common squid Todarodes pacificus, we built a joint prior of only two parameters, intrinsic growth rate (r) and carrying capacity (K), based on the resilience level of the population (Froese et al., 2017), and used a Bayesian state-space production assessment model. We used template model builder (TMB), a R package for implementing the assessment model, and estimating all parameters in the model. The predicted annual biomass was in the range of 0.76×10⁶ to 4.06×10⁶ MT, the estimated MSY was 0.13×10⁶ MT, the estimated r was 0.24, and the estimated K was 2.10×10⁶ MT.

서 론

한국 연근해 어장에서 주요 개체군의 경우, 정부는 해당 개체군의 TAC (total allowable catch)를 산정하여 어업관리를 하고 있다(Sim et al., 2020). 살오징어(Todarodes pacificus) 어업관리 또한 TAC를 계산하여 이루어진다. 그러나 TAC의 산정에 이용가능한 자료는 연도별 어획량과 단위노력당어획량 (catch per unit effort, CPUE) 에만 한정되어 있어서 잉여생산모델을 적용하여 자원평가를 해야하는 상황이다. 고전적인 잉여생산모델에 수반된 모수의 개수는 5개(b, σ_o, r, K, q; Table 1) 인데, 연도별 자원크기[생체량(biomass)]를 random effects로 다룬다면[즉, 상태공간(state-space) 잉여생산모델을 적용할 경우] 추정해야 할 모수가 하나 더 추가(σ_p; Table 1) 되고 더욱이 연도별 생체량도 예보해야만 한다. 연도별 어획량과 CPUE 자료만 이용할 경우, 과모수화(over-parameterization) 문제가 발생하여 고전적인 잉여생산모델에 수반한 모수들 조차도 추정이 어렵다. 이런 어려움을 극복하고자 과거 연구자들은 다소 임의적인 가정을 적용하여 모수들을 추정하였다. 종종 이용되는 대표적인 가정은 초기 시작시기의 생체량과 환경수용력은 같다고(B₁=K) 가정하여 추정할 모수의 개수를 줄이는 것이다 (Polacheck et al., 1993; Millar and Meyer, 2000; Punt, 2003; Zeller et al., 2008; Carruthers et al., 2011; Hilborn and Wal- ters, 2013). 어획자료가 수집되기 전에 어획활동이 전무하였다면 이 가정(B₁=K)이 객관적인 논리일 수 있으나 그러한 상황이 아니라면 이 가정은 옳지 않다.

Table 1. Points estimates of parameters and their standard errors from the state-space surplus production model for common squid Todarodes pacificus stock in Korean waters

본 연구에서는 1999–2020년 시기의 살오징어 어획량과 CPUE를 사용하였다. 1999년 이전에도 한국 연근해 해역에서 살오징어 어업은 활발히 이루어졌기에 우리는 흔한 가정 (B₁=K)을 배제하고 대신에 모수의 베이지언 사전확률분포 (prior probability distribution)를 적용하여 모수들을 추정하고 연도별 생체량을 예보하고자 하였다. 이용되는 자료(어획량과 CPUE)와 무관한 객관적인 사전확률분포를 적용하는 것이 베이지언 방법론의 핵심인데, 산정 결과치인 베이지언 사후확률분포(posterior probability distribution)에 미치는 사전확률분포의 영향을 줄이기위해서 가능한 정보사전확률분포(informa- tive priors)를 덜 적용해야 한다. 그런 이유로 두 개의 모수(r과 K) 에만 정보사전확률분포를 적용하고 나머지 모수들에겐 무정 보 사전확률분포(non-informative priors)를 적용하였다. 특히 r과 K에 적용하는 정보사전확률분포에 객관성을 높히기 위해서 자원복원력(resilience) 개념(Froese et al., 2017)을 적용하는 방법을 사용하였다. 본 연구에서는 r과 K추정치 사이에 상관관계가 존재하므로 정보사전확률분포를 적용할 때 이들 간의 상관관계도 고려하였다. FishBase (https://www.sealifebase. ca/search.php)에는 Musick (1999)과 Froese et al. (2000)이 규정한 각 수산생물종에 대한 자원복원력 기준이 등재되어 있다. 이들은 자원복원력을 ‘high’, ‘medium’, ‘low’, ‘very low’의 네 가지 수준으로 구분하였으며, 각각의 항목에 상응하는 r의 범위가 제시되었다(Table 2 in Froese et al., 2017). 자원복원력이 ‘high’ 에 가까울수록 r은 높은 경향이 있다.

본 연구의 목적은 과거 잉여생산모델의 자원평가에 흔하게 사용되었던 임의적 가정을 배제하고 객관적인 모수의 사전확률분포를 적용하여 한국 살오징어 개체군의 자원평가를 수행하는 것이다. 특히 적합도(goodness of fit)를 향상시키기 위해서 상태공간 잉여생산모델을 적용하였다. 최근에 개발된 TMB (template model builder; Kristensen et al., 2015) 언어를 이용하여 본 연구의 상태공간 잉여생산량모델을 실행하였다.

재료 및 방법

어획량 및 CPUE 자료

한국 살오징어 자원의 최대지속가능어획량(maximum sustainable yield, MSY)를 추정하기 위해 한국 살오징어의 연도별 총 어획량 자료(1999–2020) (KOSTAT, 2022)와 오징어를 주로 어획하는 근해채낚기 어업의 연도별 단위노력당 어획량자료(1999–2020)를 사용하였다(Fig. 1 and Fig. 2). 근해채낚기 어업의 살오징어 어획량은 2006년부터 2020년까지 전체 살오징어 어획량의 약 28%를 차지해 왔다. 2020년에 근해채낚기어업은 전체 살오징어 어획량(56, 621 MT)의 32%를 차지하는 18, 194 MT를 어획함으로서 살오징어를 어획하는 어법 가운데 가장 많은 어획량을 기록하였다. 따라서 본 연구에서는 근해채낚기 어업의 연도별 CPUE 자료를 한국 살오징어 자원량의 상대적인 크기를 나타내는 지수로 간주하였다.

Fig. 1. Observed yield data of the common squid Todarodes pacific us from 1999 to 2020 and the estimated MSY. The circles and broken line denote the data and MSY, respectively.

베이지언 상태공간 잉여생산량 모델

한국 살오징어 자원의 MSY를 추정하기 위해 베이지언 상태 공간 Schaefer 잉여생산량 모델을 사용하였다. 상태 공간 모델은 자료에 수반되는 관측오차(observation error)와 상태변수 (random effects, state variable, hidden variable)에 수반되는 과정 오차(process error)를 함께 고려하는 모델이다. 상태공간 Schaefer 잉여생산량 모델의 구조를 수식으로 나타내면 다음과 같다.

\(B_{t+1}=\left[B_{t}+r \cdot B_{t}\left(1-\frac{B_{t}}{K}\right)-Y_{t}\right] \cdot \exp \left(\varepsilon_{t}^{p}\right), \text { where } \varepsilon_{t}^{p} \sim N\left(0, \sigma_{p}^{2}\right)\)

Fig. 2. Predicted and observed CPUE from the jigger fisheries in Korean waters. The circles and solid line denote the data and the predicted values, respectively. CPUE, catch per unit effort.

\(I_{t}=q \cdot B_{t} \cdot \exp \left(\varepsilon_{t}^{o}\right), \text { where } \varepsilon_{t}^{o} \sim N\left(0, \sigma_{o}^{2}\right)\)

여기서 B_t는 연도별 자원량, Y_t는 연도별 어획량 자료, I_t는 연도별 CPUE 자료이며, r은 내재성장률, K는 환경수용력, q는 어획능률계수를 의미한다. 과정오차 (ε_t ^p)와 관측오차(ε_t ^o)는 평균이 0이고 분산이 각각 σ_p ², σ_o ²인 정규분포를 따른다고 가정하였다. 여기서 B_t는 모수(fixed parameters)가 아닌 상태변수로 취급되었다.

그런데 본 연구에서 사용한 방법은 위 수식과 몇 가지 차이가 있다. 먼저, 수치최적화를 안정시키기 위하여 과정오차의 분산 (σ_p ²)과 관측오차의 분산(σ_o ²)을 서로 동일한 값으로 가정하였다(σ_p ²=σ_o ²=σ²) (Auger-Methe et al., 2016). 둘째, 상태변수 B_t를 K로 나눈 값, P_t로 재모수화하여(P_t =B_t /K), 숫자의 스케일을 낮추었으며, 결과적으로 P_t에 로그를 취한 값을 상태변수로 취급하였다. 마지막으로 자원량의 초기값(B₁₉₉₉)을 추정하기 위해, 자원량의 초기값을 K로 나눈 값을 b로 정의하여(b=B₁₉₉₉/K), 이값을 추정하였다. 위 내용들을 정리하면 다음과 같다.

\(\log P_{1} \sim N\left(\log (b), \sigma^{2}\right)\)

\(\log P_{t+1} \sim N\left(\log \left[P_{t}+r \cdot P_{t}\left(1-P_{t}\right)-\frac{Y_{t}}{K}\right], \sigma^{2}\right)\)

\(\log I_{t} \sim N\left(\log \left(q \cdot K \cdot P_{t}\right), \sigma^{2}\right)\)

사전확률분포

본 연구는 베이지언 방법론을 바탕으로 모델에서 추정하고자 하는 모수 r과 K에 사전확률분포를 적용했다. Froese et al. (2017)은 자원복원력(resilience) 개념을 이용하여 어떤 수산생물 개체군의 모수 r과 K에 대한 사전확률분포의 정보를 결정하는 방법을 제시하였다. 전세계의 여러 어종에 대한 정보를 보유하고 있는 FishBase (https://www.sealifebase.ca/summary/ todarodes-pacificus.html)에 따르면, 살오징어의 자원복원력은 ‘high’이다. 이 경우, Froese et al. (2017)에 따른 r의 범위는 0.6 에서 1.5이다. 한편 K는 다음 식에 의해 계산된다(equation 3 in Froese et al., 2017).

\(K_{\text {low }}=\frac{\max \left(Y_{t}\right)}{r_{\text {high }}}, K_{\text {high }}=\frac{4 \cdot \max \left(Y_{t}\right)}{r_{\text {low }}}\)

여기서 K_low와 K_high는 각각 K의 하한값과 상한값이며, max (Y_t)는 주어진 연도별 어획량 자료의 최대값, r_high와 r_low는 앞서 살 오징어의 자원복원력에 근거한 r의 상한값과 하한값으로서 r_high=1.5, r_low=0.6 이다.

본 연구에서는 r과 K의 사전확률분포로서 로그 정규분포를 가정하였다(Meyer and Millar, 1999; Brodziak and Ishimura, 2011). 사전확률분포 각각의 최빈값(mode)은 앞서 계산된 r 과 K의 상한값과 하한값의 평균값을 사용하였다(r_mode=1.05, K_mode=0.92×10⁶ MT). 한편 사전확률분포의 변동계수(CV)는 가능한 한 큰 값을 취함으로써 사전확률분포의 주관성을 약화시키고자 하였다. 이를 위해 사전확률분포 각각의 변동계수를 수치적 최적화가 실패하기 직전까지 증가시켰다. 결과적으로 r과 K의 사전확률분포의 CV는 각각 110%로 설정되었다.

목적함수 및 모수의 추정

모수 추정을 위한 모델의 목적함수는 다음과 같다. L₁ 과 L₂는과정오차 방정식에 대한 가능도 함수이며, L₃는 관측오차 방정식에 대한 가능도 함수이다.

\(L_{1}=\frac{1}{\sqrt{2 \pi \sigma^{2}}} \exp \left[-\frac{\left(\log P_{1999}-\log (b)\right)^{2}}{2 \sigma^{2}}\right]\)

\(L_{2}=\prod_{1999}^{2020} \frac{1}{\sqrt{2 \pi \sigma^{2}}} \exp \left[-\frac{\left(\log P_{t+1}-\log \left[P_{t}+r \cdot P_{t}\left(1-P_{t}\right)-\left(Y_{t} / K\right)\right]\right)^{2}}{2 \sigma^{2}}\right]\)

\(L_{3}=\prod_{1999}^{2020} \frac{1}{\sqrt{2 \pi \sigma^{2}}} \exp \left[-\frac{\left(\log I_{t}-\log \left(q \cdot K \cdot P_{t}\right)\right)^{2}}{2 \sigma^{2}}\right]\)

위의 세 가능도 함수는 상호독립을 가정하여 아래와 같이 결합 가능도 함수의 형태로 표현할 수 있다.

L(K, r, q, b, ^σ2, P_~|I_~, Y_~)=L₁·L₂·L₃

r과 K에 대한 사전확률분포도 마찬가지로 서로 독립이라고 가정하였다. 따라서 다음과 같은 결합사전확률분포를 얻는다.

Pr(r,K)=Pr(r)·Pr(K)

모수 추정

모델의 상태변수 P_~ 와 모수들 (K, r, q, b, σ²)은 R 소프트웨어의 패키지인 TMB를 사용하여 추정하였다(Kristensen et al., 2015). TMB에 의해 모델의 모수들은 최대가능도법으로 추정 되며, 상태변수 P~ 의 예측치는 사후분포(posterior probability distribution)의 최빈값으로 예측된다. R 패키지 tmbstan을 이용하여 MCMC (Markov chain Monte Carlo) 샘플링을 실시하였으며, 총 1, 000, 000개의 표본을 추출하였다. 그리고 최종적으로 총 5, 000개의 샘플을 이용하여 사후분포를 형성하였다.

각각 모수의 사후분포가 형성된 후, 사후분포의 수렴 여부 및 자기 상관(auto-correlation)의 정도를 진단하는 과정을 실시하였다. 본 연구에서는 4가지 방법을 바탕으로 사후분포의 수렴 여부 및 자기상관성의 정도를 진단하였다: (i) Raftery and Lewis statistic, (ii) lag 1 auto-correlation, (iii) the ratio between the time series standard error and the naïve standard error, (iv) 사후분포의 형태. 진단에는 R 소프트웨어 CODA 패키지를 이용하였다. 일반적으로 Raftery and Lewis statistic의 의존성 인자 (dependence factor) 값이 5보다 작거나, lag 1 auto-correlation 의 결과값이 0에 가까울 때, 그리고 the ratio between the time series standard error and the naïve standard error의 결과값이 1 에 가까우면 추정된 사후분포들이 수렴한다고 진단할 수 있다.

결 과

모수의 점추정치와 추정치의 불확실성

모델의 모든 모수들은 성공적으로 추정되었으며, 모수의 점 추정치와 추정치의 표준오차는 Table 1에 정리하였다. 한국 살오징어 자원의 MSY 추정치는 0.1310⁶ MT로서, 1999년부터 2015년 까지의 한국 살오징어 어획량은 MSY 값을 상회하였으나, 최근 5년(2016–2020) 동안의 살오징어 어획량은 MSY 를 하회하였다(Fig. 1)

CPUE의 적합값, MSY의 추정치, 자원량의 예측치

CPUE 자료에 관측오차가 고려되어 예측된 적합값은 자료 값의 추세를 잘 반영하였다(Fig. 2). 모델에서 개체군의 자원량은 CPUE 자료값과 비례한다고 가정하였기 때문에 예측된 자원량의 추세는 CPUE 자료의 시계열과 유사하다. 예측된 살오징어 자원량은 1999년에는 4.0×10⁶ MT을 상회하였으나 2008년에 1.0×10⁶ MT 수준으로 감소한 뒤, 2021년까지 0.8×10⁶ MT 내외의 연간 자원량을 보였다(Fig. 3). 한편, 베이지언 MCMC 샘플링을 통해 생성된 MSY의 사후분포는 MSY의 추정치를 중심으로 하나의 봉우리를 이루는 형태로 나타났다(Fig. 4).

Fig. 3. Predicted annual biomass of the common squid stock in Korean waters. The circles denote the predicted values. The error bars denote plus/minus one standard error.

Fig. 4. Posterior probability distribution of MSY obtained from Markov chain Monte Carlo sample sets. The vertical line denotes the estimate of the MSY.

r과 K에 대한 사전확률분포 및 사후확률분포

본 모델에서는 모수 r과 K의 수치적 최적화를 돕기 위해 Froese et al. (2017)가 제시하는 방법에 근거하여 두 모수에 각각 사전분포를 가정하였다. 결과적으로 생성된 각 모수의 사후분포는 사전분포와 형태는 다르지만 모두 단봉 그래프의 형태를 보였다(Fig. 5). 모델 추정 결과는 사후분포의 수렴 여부 및 자기 상관(auto-correlation)의 정도의 진단기준을 모두 만족하였다. r과 K에 대한 사전확률분포 및 사후확률분포의 최빈값은 r은 1.05와 0.24, K는 0.92×10⁶ MT과 2.10×10⁶MT로 나타났다.

Fig. 5. Prior probability distributions (solid line) and posterior probability distribution (histogram) for parameters r and K.

고 찰

본 연구에서는 잉여생산량 모델에 입력되는 자료 중 어획량은 한국 살오징어의 총 어획량, CPUE는 근해채낚기 어업으로부터 얻어진 자료를 사용하였다. 근해채낚기 어업은 살오징어를 주로 어획하는 어업이며, 2005년 이후 살오징어의 총 어획량 중 1/4 이상을 어획해왔다. 따라서 근해채낚기 어업의 CPUE는잉여생산량 모델에서 평가하고자 하는 살오징어 자원량의 상대적인 수치로서 사용되기에 타당하다고 판단하였다. 그러나 명백히 근해채낚기 어업이 살오징어 전체 어획량 중 과반을 차지하지는 못하기에, 근해채낚기 CPUE 자료의 대표성에 대한 문제가 제기된다. 살오징어의 어업별 어획 비중은 상황에 따라 변할 수 있으므로 자료의 대표성 문제는 차후에도 발생할 수 있다. 만약 어떤 어종의 어획량을 몇 개의 주요 어업이 상당량을 차지하고 각 어업의 CPUE 자료가 존재한다면, 잉여생산량 모델을 이용하여 복수어업을 고려한 자원평가를 수행할 수 있다. 이러한 방법으로 본 연구에서 제기되는 CPUE 자료의 대표성 문제를 완화할 수 있다.

수산자원평가모델 가운데 가장 간단한 잉여생산량 모델일지라도 모델의 모수들을 한꺼번에 추정하는 것은 쉽지 않다. 우리는 위의 문제를 극복하기 위해 모수의 사전분포를 결정하는 Froese et al. (2017)의 방법에 근거하여 r 과 K에 사전확률분포를 제공하였다. 이로서 흔히 사용되는 초기자원량을 환경수용력이 같다(B₁=K)는 가정을 사용하지 않고 초기자원량을 자유 모수로서 추정할 수 있었다. 또한 Froese et al. (2017)의 방법으로 r과 K의 사전분포의 형태를 결정하는 과정에서 사전확률분포 각각의 CV를 최대한 증가시킴으로서 사전확률분포의 객관성을 강화하고자 하였다. 본 연구에서 제시한 베이지언 상태 공간 잉여생산량 모델은 r과 K의 사전확률분포의 형태를 결정함에 있어서 Froese et al. (2017)의 방법에 근거하지만, 최종적으로 다루고자 하는 개체군의 자원복원력을 결정하는 재량권을 사용자에게 부여한다는 장점이 있다.

Punt and Hilborn (1997)은 베이지언 자원평가모델에서 모수의 사전확률분포를 가정할 때 상당한 주의를 기울일 것을 당부했다. 선행연구에 근거하여 잉여생산량 모델의 모수 r과 K에 대한 사전확률분포를 가정할 경우, 동일한 계군 또는 종일한 어종의 사전정보를 사용하는 것이 권장된다(Millar and Meyer, 2000). 그러나 동일 계군 또는 동일한 종으로부터 모델의 모수에 대한 사전정보를 획득하기가 어려운 경우가 존재할 수 있다. 그 예로 Jung (2019)은 베이지언 상태공간 잉여생산량 모델을 이용하여 한국 고등어 자원을 평가하는 연구에서 모수 r과 K 대한 사전정보를 얻을 수 없었다. 이에 Jung (2019)은 여러 선행연구를 참고하여 r 과 K대한 정보를 수집하고, 이 가운데 그럴듯한 정보를 취하여 r과 K의 사전확률분포를 가정하였다. 반면, 본 연구에서 사용한 방법은 평가하고자하는 개체군의 r과 K에 대한 사전정보가 부재하더라도 FishBase에 등록된 자원복원력 기준에 근거하여 어느 정도 객관적인 사전확률분포를 가정할 수 있다는 장점이 있다.

결론적으로 본 연구에서 사용한 방법은 모델의 모수에 가정하는 사전확률분포에 어느 정도의 객관성을 확보하면서 잉여생산량 모델의 모수들을 한꺼번에 자유모수로서 추정하게 한다는 점에서 유용하다. 따라서 어획량 자료와 CPUE자료가 확보된 수산생물 개체군을 베이지언 잉여생산량 모델을 이용하여 평가하는 방법으로서 권장할 만하다.