A Bayesian State-space Production Assessment Model for Common Squid Todarodes pacificus Stock Caught by Multiple Fisheries in Korean Waters

한국 해역의 살오징어(Todarodes pacificus) 개체군 자원평가를 위한 베이지안 상태공간 잉여생산량 모델의 적용

Abstract

Given data about the annual fishery yield of the common squid Todarodes pacificus, and the catch-per-unit-effort (CPUE) data from multiple fisheries from 2000-2018, we applied a Bayesian state - space assessment model for the squid population. One of our objectives was to do a stock assessment, simultaneously incorporating CPUE data from the following three fisheries, (i) large trawl, (ii) jigger, and (iii) large purse seine, which comprised on average a year about 65% of all fisheries, allowing possible correlations to be reflected. Other objectives were to consider both observation and process errors and to apply objective priors of parameters. The estimated annual exploitable biomass was in the range of 3.50×10⁵ to 1.22×10⁶ MT, the estimated intrinsic growth rate was 1.02, and the estimated carrying capacity was 1,151,259 MT. Comparison with available results from stock assessment of independently analyzed single fisheries revealed a large difference from the estimated values, suggesting that stock assessment based on multiple fisheries should be performed.

서론

한국해역의 살오징어(Todarodes pacificus)는 다양한 어구를 통해 어획되는 대표적인 복수어업 대상종 중 하나이다. 현재 총 41개의 연근해 어업 중 21개 어업에서 살오징어를 어획하고 있으며, 6개 어업(근해채낚기 어업, 대형 트롤 어업, 대형 선망 어업, 동해구 중형 트롤 어업, 근해자망 어업, 그리고 쌍끌이 대형저인망 어업)이 살오징어 TAC (total allowable catch) 대상어업에 참여하고 있다. 그런데 2000년대 초반 이후로 한국 해역의 살 오징어 연간 어획량은 크게 감소하고 있는 실정이다. 2003년 약 23만톤에 이르던 살오징어 어획량이 2015년에는 약 15만 톤으로 급격히 감소하였으며 이후, 지속적인 감소를 거쳐 2018년의 어획량이 약 4만톤으로 기록된 사실은 살오징어 자원의 심각성을 시사하는 바이다(KOSTAT, 2020). 하지만 이러한 상황에도 불구하고 살오징어 자원평가에 대한 연구는 현재까지 부재한 상황이다. 현존하는 수집 자료(e.g., 연령자료, 직접 자원조사자료, 등)도 단기간에 국한되어 있거나 부재하며, 단일 어업이 아닌 복수 어업을 고려한 자원평가 연구가 상당히 제한적인 상황이다. 이런 상황 속에서 살오징어를 어획하는 여러 어업들의 CPUE (catch per unit effort)자료가 수집되었다는 점은 잉여생산량 모델의 적용 가능성과 복수 어업을 고려한 자원평가의 활용 가능성을 제시한다. 현재까지 잉여생산량 모델을 적용한 대부분의 자원평가 연구를 살펴보면 단일어업의 CPUE자료만을 이용하여 자원평가를 수행했다(Hilborn and Walters, 1992; Millar and Meyer, 2000; Punt and Szuwaiski, 2021; Hyun and Kim, 2022). 하지만 실제로는 여러 어업에서 어느 한 개체군을 어획하기 때문에 복수 어업을 고려한 자원평가 연구가 필요하다(Pascoe et al., 2020). 이런 상황 속에서 최근에 잉여생산량 모델을 바탕으로 복수 어업을 고려한 자원평가 연구가 소개된바 있다(Hyun, 2018; Winker et al., 2020). 하지만 Hyun (2018) 에서는 자원량의 과정 오차를 별도로 고려하지 않았고, Winker et al. (2020)의 경우, 여러 어업의 CPUE를 서로 독립으로 가정하였으며, 어업 선택성을 반영하였으나 입력 값으로 고려했다는 단점을 지니고 있다.

따라서 본 연구는 크게 두 가지 목적에 기반을 두고 있다. 첫째로는 베이지안 상태공간 잉여생산량 모델에 복수어업을 고려한 살오징어 개체군 자원평가를 실시하는 것이다. 둘째로는 복수어업을 고려했을 때와 단일어업만을 고려했을 때의 자원평가 결과를 비교 및 분석하여 복수어업을 고려한 자원평가의 필요성을 증명하는 것이다. 본 연구에서는 기존 잉여생산량 모델에 상태 공간 구조를 적용함으로써 추정하고자 하는 자원량의 과정 오차와 자료의 관측 오차를 동시에 고려하였다. 그리고 여러 어업의 CPUE자료가 다변량 정규분포를 따른다고 가정하며 복수어업 CPUE간의 상관관계를 추가적으로 고려하였다. 그런데 TMB R package를 바탕으로 상태공간 잉여생산량 모델을 적용했을 때 수치 최적화가 제대로 이루어 지지 않았기 때문에 본연구에서는 정보적 사전분포를 바탕으로 베이지안 방법론을 적용하였다. 그리고 한국해역의 살오징어의 경우, 어업 선택성에 관한 정보가 부재한 상황이다. 국립수산과학원으로부터 제공받은 여러 어업의 CPUE자료는 상업 어업으로부터 수집된 자료이며 본 연구에서는 어업 선택성을 별도로 고려하지 않았다.

베이지안 상태공간 잉여생산량모델을 적용한 한국해역의 살오징어 자원평가는 두 가지 방식으로 실시하였다. 첫 번째는 복수어업을 고려한 방식으로써 세 가지 어업의 CPUE자료를 베이지안 상태공간 잉여생산량 모델에 적용하여 모수를 추정하였다(이하 Case I). 두 번째는 단일어업만을 고려한 방식으로써 각 어업별 하나의 CPUE자료를 베이지안 상태공간 잉여생산량 모델에 적용하여 모수를 추정하였다(이하 Case II). Case II의 경우 적용하는 단일어업에 따라 3개의 케이스(case)로 다시 구분하였다: i) 근해채낚기 어업만을 고려한 방식(이하 Case II-1), ii) 대형선망 어업만을 고려한 방식(이하 Case II-2), 그리고 iii) 대형트롤 어업만을 고려한 방식(이하 Case II-3). 베이지안 상태 공간 잉여생산량모델을 바탕으로 복수어업을 고려한 살오징어 자원평가(Case I)를 주요 결과로서 도출하였으며 복수어업을 고려한 경우(Case I)와 단일어업만을 고려한 경우(Case II) 의 자원평가 결과를 비교해보았다.

재료 및 방법

살오징어 자료

통계청에서 2000년부터 2018년까지 수집한 한국해역의 연도별 살오징어 총 어획량 자료(KOSTAT, 2020)와 동기간의 국립수산과학원에서 수집한 CPUE자료를 본 연구에 이용하였다. 한국 해역에서 살오징어를 어획하는 세 가지 어업의 CPUE자료를 사용하였으며, 구체적으로 2000년부터 2018년까지 살오징어 총 어획량의 약 33.5%를 차지하는 대형트롤(large trawl, tw) 어업과 약 27.3%를 차지하는 근해채낚기(jigger, jig) 어업, 그리고 약 4.2%를 차지하는 대형선망(large purse seine, ps) 어업의 CPUE자료를 사용하였다. 해당 세 종류의 어업만을 고려한 이유는 해당 어업이 2018년을 기준으로 살오징어 TAC 대상 어업에 참여하고 있음과 동시에 CPUE자료의 수집이 가능하였기 때문이다.

본 연구의 대상종인 살오징어는 우리나라 전체 해역, 동중국해, 그리고 일본 북해도 이남까지 북서 태평양 전체를 대상 시스템으로 분포하고 있다. 하지만 본 연구에서는 한국 해역에서 수집한 어획량, CPUE자료만을 연구에 이용하였기 때문에 대상 해역을 한국 해역으로 한정하였다. 그리고 국립수산과학원에서 수집한 여러 어업의 CPUE자료에는 조업 해역에 대한 정보, 그리고 어업 선택성에 대한 정보가 부재하기 때문에 어업별 살오징어에 대한 공간 영역을 우리나라의 전체 해역으로 가정하였다.

국립수산과학원에서 수집한 CPUE (근해채낚기 CPUE, 대형선망 CPUE, 대형트롤 CPUE)자료는 직접자원조사자료가 아닌 상업 어업으로부터 수집된 자료이기 때문에 상대적 개체군의 크기를 나타내는 CPUE의 단위와 추세가 어업마다 서로 상이함을 볼 수 있다(Fig. 1b-1c).

Fig. 1. Annual yield of the common squid from 2000 to 2018 (panel a), the jigger fishery CPUE (dashed line with dots in panel b), the purse seine fishery CPUE (dashed line with open circles in panel b), and the trawl fishery CPUE (dashed line with open circles in panel c).

상태공간 잉여생산량모델

본 연구에서는 다양한 형태의 잉여생산량 함수 중, Schaefer 함수(Schaefer, 1954)를 적용하였다. Pella and Tomlinson함수 (Pella and Tomlinson, 1969)도 있지만, 모수를 하나 더 추정해야 한다는 단점이 있기에 Schaefer함수를 적용하였다. 그리고 기존 잉여생산량모델에 상태공간 구조를 적용함으로써, 자료의 관측 오차와 자원량의 과정 오차를 동시에 고려하였다.

\(B_{1}=b K \exp \left(\varepsilon_{1}^{p}\right), \text { where } \varepsilon_{1}^{p} \sim N\left(0, \sigma_{p}^{2}\right)\)       (1)

\(B_{t+1}=\left[B_{t}+B_{t}\left(1-\frac{B_{t}}{K}\right) r-Y_{t}\right] \exp \left(\varepsilon_{t+1}^{p}\right),\)

\(\text { where } \varepsilon_{i+1}^{p} \sim N\left(0, \sigma_{p}^{2}\right), \text { and } t \geq 1 \)       (2)

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

여기서 t는 연도, Y_t는 연도별 어획량, I_t는 연도별 상대적 개체군 크기, r은 본원적 성장률, K는 최대 환경수용력, b는 초기 어획 가능 자원량에 최대 환경수용력을 나눈 계수, \(\varepsilon_{t}^{p}\)는 과정 오차, \(\varepsilon_{t}^{o}\)는 관측 오차, q는 어획능률계수, B_t는 연도별 어획가능 자원량을 의미한다. 본 연구에서 사용하는 연도별 CPUE자료는 상업 어업에서 수집된 자료이며 어업 선택성이 반영된 자료이기 때문에 B_t를 연도별 총 자원량이 아닌, 연도별 어획가능 자원량으로 정의하였다. 연도별 CPUE자료 I_t와 연도별 어획가능 자원량 B_t는 q 상수에 대하여 비례관계를 가진다.

본 연구에서 이용된 모든 표기기호에 대한 자세한 설명은 Table 1에 별도로 정리하였으며, 몇 가지 가정을 기반으로 모델을 표기하였다. 우선 과정 오차와 관측 오차가 정규분포를 따르며 평균이 0이고, 분산이 각각 \(\sigma_{p}^{2}\), \(\sigma_{o}^{2}\) 를 가진다고 가정하였다. 그리고, 수치 최적화 모수 추정을 위해 과정 오차의 분산과 관측 오차의 분산이 서로 동일한 값(\(\sigma_{p}^{2}\)=\(\sigma_{o}^{2}\) )을 가지도록 가정하였다.(Auger-Méthé et al., 2016). 식 (1)-(2)에서는 과정 오차를 고려함으로써 B_t를 고정 모수가 아닌, 상태 변수로 가정하였다. 초기 어획가능 자원량 B₁의 경우, 식 (1)과 같이 상수 b에 최대 환경수용력 K와 함께 과정 오차를 고려한 상태변수로 가정하였다. 식 (3)과 같이 CPUE자료는 어획가능 자원량 B_t에 비례하며 관측오차를 추가적으로 고려하였다(Polacheck et al., 1993; Winker et al., 2018).

Table 1. Notations used in this paper. Bold font is used for a vector to distinguish it from a scalar.

복수어업을 고려한 베이지안 상태공간 잉여생산량 모델 (Case I)

Case I은 복수어업을 고려한 방식으로써, 세 가지 어업의 CPUE자료를 이용하여 베이지안 상태공간 잉여생산량 모델에 적용한 방식이다. 식 (1)-(3)에서의 가정을 동일하게 유지하며 과정 방정식과 관측 방정식을 다음과 같이 표기하였다.

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

\(\text { where } \varepsilon_{t+1}^{p} \sim N\left(0, \sigma_{p}^{2}\right), \text { and } t \geq 1\)       (5)

\(i_{f, t}=q_{f} b_{t} \exp \left(\varepsilon_{t}^{f}\right), \text { where } \varepsilon_{t}^{f} \sim n\left(0, \sigma_{f}^{2}\right)\)       (6)

여기서 f는 어업의 종류, \(y_{t}^{\text {total }}\)는 모든 어업에서의 살오징어 연도별 총 어획량, \(\varepsilon_{t}^{f}\)는 각 어업별 CPUE 자료의 관측 오차, q_f는 각 어업별 어획능률계수를 의미한다.

살오징어 총 어획량자료와 세 가지 어업의 CPUE자료를 이용하여 어업 선택성이 반영된 어획가능 자원량 B_t를 추정하였다. 식 (6)과 같이 각 어업별 CPUE는 B_t에 직접 비례하며 각각의 어획능률계수 q_f와 관측오차 \(\varepsilon_{t}^{f}\) 를 가진다고 가정하였다. 세 가지 어업의 상관성을 추가적으로 고려하고자 각 어업별 CPUE자료를 확률벡터로서 고려하며 다변량 정규분포를 따른다고 가정하였다. 수치 최적화 및 모수 추정을 위해 과정 오차의 분산과 각 어업별 관측 오차의 분산이 동일한 값\(\left(\sigma_{p}^{2}=\sigma_{f}^{2}\right)\)을 가진다고 가정하였다(Auger-Méthé et al., 2016).

상태변수 B_t의 단위를 낮추면서 수치최적화의 안정성을 향상시키기 위해 B_t를 최대 환경수용력 K로 나누었다(P_t=B_t/K). 그리고 양변에 log를 취함으로써 logP₁, logP_t가 다음의 식 (7), (8) 과 같이 정규분포를 따르도록 식을 재구성하였다. 그리고 각 어업별 CPUE 자료의 경우 식 (9)와 같이 확률벡터로서 고려하며 다변량 정규분포를 따른다고 가정하였다.

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

\(\begin{array}{r} \log P_{t+1} \sim N\left(\log \left[P_{t}+r P_{t}\left(1-P_{t}\right)-Y / K\right], \sigma_{p}^{2}\right) \\ \text { where } t \geq 1 \end{array}\)       (8)

\(\mathrm{I} \sim N(\mathbf{u}, \Sigma), \text { where } \mathrm{I}=\left(\log \left(I_{i g, t}\right), \log \left(I_{p s, t}\right), \log \left(I_{t w, t}\right)\right)\)       (9)

식 (9)에서 아래 첨자 jig, ps, tw은 각각 근해채낚기 어업, 대형선망 어업, 대형트롤 어업을 의미한다. 다변량 정규분포의 평균은 식 (10)과 같이 표기하였으며, 각 어업별 상관계수를 고려한 분산 공분산 행렬은 다음의 식 (11)와 같이 표기하였다.

\(\mathbf{u}=\left(\log \left(q_{j i g} K P\right), \log \left(q_{p s} K P_{t}\right), \log \left(q_{t w} K P_{t}\right)\right)\)       (10)

\(\sum=\left[\begin{array}{ccc} \sigma_{j i g}{ }^{2} & \sigma_{j i g} \sigma_{p s} \rho_{j i g, p s} & \sigma_{j i g} \sigma_{t w} \rho_{j i g, t w} \\ \sigma_{j i g g} \sigma_{p s} \rho_{j i g, p s} & \sigma_{p s}{ }^{2} & \sigma_{p s} \sigma_{t W} \rho_{p s, t w} \\ \sigma_{j i g} \sigma_{t w} \rho_{j i g, t w} & \sigma_{p s} \sigma_{t w} \rho_{p s, t w} & \sigma_{t w}{ }^{2} \end{array}\right]\)       (11)

분산 공분산 행렬 중, 어업별 상관계수 ρ_jig,ps, ρ_jig,tw, ρ_ps,tw는 표준화한 어업별 CPUE자료를 기반으로 입력 값으로써 고려하였다. 그 값은 각각 ρ_jig,ps=0.473 (P<.000), ρ_jig,tw=0.208 (P<.000), ρ_ps,tw=0.767 (P<.000)이다.

사전 분포

상태공간 잉여생산량 모델을 적용했을 때 수치 최적화가 제대로 이루어 지지 않았기 때문에 정보적 사전분포를 바탕으로 베이지안 방법론을 도입하였다. 본 연구에서는 추정하고자 하는 모수 중, 본원적 성장률 r과 최대 환경수용력 K에 정보적 사전분포를 적용하였다. Froese et al. (2017)에서는 FishBase의 자원회복력(resilience) 정보를 활용하여 모수 r과 K에 정보적 사전분포를 주는 방법을 고안하였는데, Froese et al. (2017)의 방법을 적용하여 모수 r과 K에 정보적 사전분포를 적용했다. 한국해역의 살오징어는 자원 회복력이 높은 종으로써 Fish Base의 자원회복력은 high로 표시되어 있다. 이 정보를 바탕으로 r과 K 모수에 정보적 사전 분포를 설정할 수 있었다(Table 2). 모수 r 과 K가 로그 정규분포를 따른다고 가정하였으며, Froese et al. (2017)에서의 방법을 활용하여 각각의 최빈값을 1.06, 855, 264 로 적용하였다. 그리고 사전분포 불확실성을 낮추기 위해 변동계수(coefficient of variation)는 1.0 (=100%)으로 고려하였다. 그 외, 나머지 모수에는 무정보적 사전분포를 적용했다.

Table 2. Prior probability density functions (PDFs) on parameters

Note: CV, coefficient of variation.

본 논문에서 적용한 사전분포를 독립으로 가정하여 결합 사전 분포(joint prior distribution)를 다음의 식 (12)와 같이 표현하였다.

\(\operatorname{Pr}(r, K)=\operatorname{Pr}(r) \operatorname{Pr}(K)\)       (12)

우도 함수와 베이지안 방법론

과정오차 방정식과 관측오차 방정식에 대한 우도함수가 상호 독립으로 가정하며 다음의 식 (13), (14)과 같이 표시하였다.

\(L\left(b, K, \sigma_{p}^{2}, r, P \mid Y\right)=\frac{1}{\sqrt{2 \pi \sigma_{p}}} \exp \left[-\frac{\left(\log P_{2000}-\log (b)\right)^{2}}{2 \sigma_{p}^{2}}\right]\)

\(\times \prod_{2000}^{2018} \frac{1}{\sqrt{2 \pi \sigma_{p}}} \exp \left[-\frac{\left(\log P_{t+1}-\log \left[P_{t}+r P_{t}\left(1-P_{t}\right)-\frac{Y_{t}}{K}\right]\right)^{2}}{2 \sigma_{p}}\right]\)       (13)

\(L\left(q_{j i g}, q_{p s^{\prime}} q_{t w^{\prime}} K, \sigma_{j i g}{ }^{2}, \sigma_{p s}{ }^{2}, \sigma_{t w}{ }^{2}, P \mid \rho_{j i g, p s^{s}} \rho_{j i g, t w^{\prime}} \rho_{p s, t w}, \mathrm{I}\right)=\)

\(\times \prod_{2000}^{2018}(2 \pi)^{-\frac{3}{2}}\left|\sum\right|^{-\frac{1}{2}} \exp \left[-\frac{1}{2}(\mathrm{I}-\mathrm{u})^{\prime} \sum^{-1}(\mathrm{I}-\mathrm{u})\right]\)       (14)

여기서, Y=(Y₂₀₀₀, Y₂₀₀₁, …, Y₂₀₁₈), 그리고 I=(log(I_jig,t),log(I_ps,t),log(I_tw,t))이다.

식 (13)는 과정 방정식의 모수(b,r,K,\(\sigma_{p}^{2}\),P)를 추정하기 위한 우도 함수, 식 (14)는 관측 방정식의 모수(q_jig,q_ps,q_tw,K,\(\sigma_{j i g g}^{2}, \sigma_{p s^{2}}^{2}, \sigma_{t w}^{2}\),P)를 추정하기 위한 우도 함수이다. 우도 함수에 대한 식 (13), (14)의 결합 우도 함수(joint likelihood function)는 다음의 식 (15)와 같다.

\(\begin{aligned} &L\left(b, K, r, q_{j i g} q_{p s}, q_{t w^{2}} \sigma_{j i g g}{ }^{2}, \sigma_{p s}{ }^{2}, \sigma_{t W}{ }^{2}, \sigma_{p}{ }^{2}, \mathrm{P} \mid \rho_{j i g, p s s} \rho_{j i g, t w}\right. \\ &\left., \rho_{p s, t w}, \mathrm{I}, \mathrm{Y}\right)= \\ &L\left(b, K, \sigma_{p}{ }^{2}, r, \mathrm{P} \mid \mathrm{Y}\right) \cdot L\left(q_{j i g g} q_{p s} q_{t w}, K, \sigma_{j j g}{ }^{2}, \sigma_{p s}{ }^{2}, \sigma_{t w}{ }^{2}, \mathrm{P}\right. \\ &\left.\mid \rho_{j i g, p s}, \rho_{j i g, t w}, \rho_{p s, t w} \mathrm{I}\right) \end{aligned}\)       (15)

식 (16)과 같이 결합 사후 분포는 베이즈 정리에 따라 결합 사전 분포와 결합 우도 함수에 비례한다.

\(L\left(\theta \mid \rho_{j i g, p s}, \rho_{j i g, t w}, \rho_{p s, t w}, \mathrm{I}, \mathrm{Y}\right) \propto \int f(\mathbf{P}, \mathrm{I}, \mathrm{Y} \mid \theta) \cdot \operatorname{Pr}(r, K) d P\)       (16)

TMB와 상호 연동되는 tmbstan R package를 이용하여 MCMC (markov chain monte carlo) 샘플링을 실시하였으며, 총 3,000,000개의 표본을 추출하였다. 이 중 초기 15,000개의 표본을 번인(burn-in) 처리하였으며 추출된 표본 간의 자기 상관(auto correlation)을 제거하기 위해 매 100번째 표본을 추출하였다(thinning). 그 결과, 총 15,000개의 표본을 이용하여 사후분포를 형성하였다.

추정하고자하는 모수(b,r,K,q_jig,q_ps,q_tw,\(\sigma_{j i g}^{2}, \sigma_{p s}^{2}, \sigma_{t w}^{2} \sigma_{p}^{2}\))의 사후 분포가 형성된 후, 사후분포의 수렴 여부 및 자기 상관의 정도를 진단하는 과정을 실시하였다. 본 연구에서는 3가지 방법을 바탕으로 사후분포의 수렴여부 및 자기 상관의 정도를 진단하였다(Hyun et al., 2005): (i) 라프테리와 루이스 통계(Raftery and Lewis statistic), (ii) 레그 1 자기 상관(lag 1 auto correlation), (iii) 시계열 표준 오차(time series standard error)와 나이브 표준 오차(naïve standard error)간의 비율. R 소프트웨어 CODA package를 이용하여 진단을 수행하였다. 일반적으로 라프테리와 루이스 통계의 의존성 인자(dependence factor)값이 5보다 작거나, 레그 1 자기 상관의 결과값이 0에 가까울 때, 그리고 시계열 표준 오차와 나이브 표준 오차간의 비율 결과값이 1에 근사하면 추정된 사후분포들이 수렴한다고 진단할 수 있다.

단일어업만을 고려한 베이지안 상태공간 잉여생산량 모델 (Case II)

복수어업을 고려한 자원평가(Case I) 결과와 비교하기 위해, 단일어업만을 고려한 자원평가(Case II)를 수행했다. Case II는 단일어업만을 고려한 방식으로써 각 어업별 하나의 CPUE자료를 이용하여 베이지안 상태공간 잉여생산량 모델에 적용한 방식이다. Case II-1은 근해채낚기 어업의 CPUE자료만을 적용한 방식이고, Case II-2는 대형선망 어업의 CPUE자료만을 적용하였다. 마지막으로 Case II-3는 대형트롤 어업의 CPUE자료만을적용하였다. 식 (1), (2), (3)을 목적함수로 설정하며, MCMC 샘플링을 실시하였다. 총 3,000,000개의 표본을 추출하였으며 사후분포에 수렴하지 못하는 초기 15,000개의 표본을 번인(burn- in) 처리하였다. 그리고 추출된 표본 간의 자기 상관을 제거하기 위해 매 100번째 표본을 추출하였다. 그 결과 총 15,000개의 샘플을 이용하여 사후분포를 형성하였다.

본 모델에 각 어업별 하나의 CPUE자료를 적용하여 추정된 살오징어 자원평가 결과들을 비교 및 분석해보며 어떤 문제점이 존재하는지, 또한 차이점은 무엇이 있는지 살펴보았다.

결과

복수어업을 고려한 자원평가(Case I)

MCMC 샘플링을 통해 추정하고자 하는 모수(b,r,K,q_jig,q_ps,q_tw,\(\sigma_{j i g}^{2}, \sigma_{p s}^{2}, \sigma_{t w}^{2} \sigma_{p}^{2}\) )의 사후분포를 도출하였고(Fig. 2), 각 모수의 점 추정치(사후분포의 최빈값)와 95% 신용구간을 Table 4에 정리하였다. 각 모수의 점 추정치 결과를 간략히 살펴보면 본원적 성장률 r과 최대 환경수용력 K는 각각 1.02, 1,151,259톤으로 추정 되었으며 추정된 K는 2002년 최대 어획량을 기록했던 233,254톤의 약 5배에 달한다. 각 어업별 어획능률계수의 점 추정치를 살펴보면 q_jig는 3.29×10^-8, q_ps는 1.23×10^-6, q_tw는 3.19×10^-6로 추정되었으며, 단위는 각각 [[q_jig]]=hook^-1, [[q_ps]]=haul^-1, [[q_tw]]=haul^-1이다. 그리고 추정된 과정 오차의 분산 \(\sigma_{p}^{2}\), 어업별 관측 분산 \(\sigma_{f}^{2}\) 는 0.41로 동일하다. MSY, B_MSY, H_MSY는 각각 293,421톤, 575,629톤, 그리고 0.51으로 추정되었다. 추정하고자 하는 모수에 대한 사후분포의 수렴 여부 및 자기 상관의 정도를 진단한 결과는 다음과 같다(Table 3). 진단 결과, 모든 모수에 대해 사후 분포의 수렴성을 확인하였다.

Fig. 2. Posterior distributions of the respective parameters under Case I. The curve in panel a indicates the informative prior of r, and that inpanel b represents the informative prior of K.

Table 3. Diagnostic tests for Markov Chain Monte Carlo samples for parameters (b,r,K,\(\sigma_{p}^{2}, \sigma_{j i g}^{2}, \sigma_{p s}^{2}, \sigma_{t w}^{2}\),q_jig,q_ps,q_tw) from four different cases: Case I (when multiple fisheries CPUE data were used), Case II-1 (when only jigger fishery CPUE data were used), Case II-2 (when only large purse seine fishery CPUE data were used), Case II-3 (when only large trawl fishery CPUE data were used

The dependence factor of Raftery-Lewies statistics (DF), lag-1 autocorrelation, and the ratio of the naïve standard error to the time series standard error were checked.

Table 4. Point estimates of parameters and their 95% credible intervals from four different cases

Case I (when multiple fisheries CPUE data were used), Case II-1 (when only jigger fishery CPUE data were used), Case II-2 (when only large purse seine fishery CPUE data were used), Case II-3 (when only large trawl fishery CPUE data were used). The 95% credible intervals of estimates are provided in parentheses, respectively.

추정된 어업별 연도별 CPUE (predicted CPUE)는 근해채낚기 어업, 대형선망 어업, 대형트롤 어업의 CPUE 자료를 모두 반영한 추정치로서, 어업별 CPUE 자료의 추세를 잘 반영하였다(Fig. 3). 하지만 추정된 어업별 연도별 CPUE를 보면 어업별로 단위 및 크기는 서로 다르지만 같은 추세를 가지는 점을 볼 수 있다.

Fig. 3. Predicted annual CPUE under Case I. In panels a-c, dots denote CPUE data from jigger, purse seine, and trawl, respectively. The respective solid line in each panel indicates the corresponding predicted CPUE, and the dashed lines in each panel indicate the 95% credible intervals.

Case I에서 살오징어 어획가능 자원량을 추정한 결과, 한국해역의 살오징어 어획가능 자원량은 시간이 지남에 따라 크게 감소한 것으로 평가되었다(Fig. 4). 특히 2017년의 어획가능 자원량은 443,911톤으로 추정되며 15년 전인 2002년 자원량의 약 36% 수준으로 나타났으며, B_MSY보다 낮게 추정되었다.

Fig. 4. Estimated annual biomass and management references under Case I. In panel a, dots are annual yield data of the common squid, and the dashed line is the maximum sustainable yield (MSY). In panel b, the solid line is the estimated annual biomass with carrying capacity (dotted line), and the dashed line is biomass that produces the MSY (B_MSY). In panel c, the solid line is the estimated annual harvest rate, and the dashed line is harvest rate at the MSY (H_MSY).

단일어업만을 고려한 자원평가(Case II)

Case I과 동일한 방법으로 MCMC 샘플링을 통해 추정하고자 하는 모수 (b,r,K,q_jig,q_ps,q_tw,\(\sigma_{j i g}^{2}, \sigma_{p s}^{2}, \sigma_{t w}^{2} \sigma_{p}^{2}\))의 사후분포를 도출하였으며, Case II에서의 각 모수별 점 추정치와 95% 신용 구간은 Table 4에 별도로 정리하였다. Case II-1, Case II-2, 그리고 Case II-3에서의 각 모수별 점추정치 결과를 간략히 살펴보면 본원적 성장률 r은 각각 0.42, 0.74, 0.58으로 추정되었으며, 최대 환경수용력 K는 1,576,995톤, 817,025톤, 1,127,473톤으로추정되었다. 어획능률 계수의 점 추정치를 살펴보면 Case II-1 에서는 2.59×10^-8, Case II-2에서는 2.98×10^-6, 그리고 Case II-3에서는 4.12×10^-6으로 추정되었다. 그리고 과정 오차의 분산, 관측 오차의 분산 는 각각 0.03, 0.10, 0.10로 추정되었다.

Case II에서 추정하고자 하는 모수에 대한 사후분포의 수렴 여부 및 자기 상관의 정도를 진단한 결과는 다음과 같다(Table 3). 진단 결과, 모든 모수에 대해 사후 분포의 수렴성을 확인하였다.

Case II-1, Case II-2, 그리고 Case II-3에서 추정된 어업별 연도별 CPUE (predicted CPUE)와 CPUE 자료간 비교 결과는 다음과 같다(Fig. 5). 전체 연도에서 어업별 CPUE자료가 추정된 어업별 CPUE의 95% 신용구간 안에 모두 포함되는 것으로 나타났다.

Fig. 5. Predicted annual CPUE from three different cases: Case II-1 (panel a), Case II-2 (panel b), Case II-3 (panel c). In panels a-c, dots denote CPUE data from jigger, purse seine, and trawl, respectively. The respective solid line in each panel indicates the corresponding predicted CPUE, and the dashed lines in each panel indicate the 95% credible intervals.

Case II에서의 각 어업별 추정된 살오징어 어획가능 자원량과 관리 기준점은 상당한 차이를 보였다(Fig. 6). 일단 Case II-1의 경우, 추정된 살오징어 어획가능 자원량은 최소 6.71×10⁵톤, 최대 3.18×10⁶톤으로 변동이 심한 모습을 보여주며 시간이 지남에 따라 적정 수준의 어획가능 자원량을 유지하는 것으로 평가되었다. Case II-2의 경우, 살오징어 어획가능 자원량은 최소 69,458톤, 최대 8.78×10⁵으로 추정되며 Case II-1보다 더욱 변동이 심한 모습을 보여줬다. 그리고 2014년 이후로 살오징어 어획가능 자원량이 붕괴되는 형태를 볼 수 있었다. Case II-3의 경우, 최소 1.39×10⁵톤, 최대 9.52×10⁵톤으로 추정되며 2015년 이후로 어획가능 자원량이 B_MSY보다 낮게 추정되며 살오징어 어획가능 자원량이 급격히 감소하는 형태를 볼 수 있었다.

Fig. 6. Estimated annual biomass and management references from three different cases: Case II-1 (panels in the first row), Case II-2 (panels in the second row), Case II-3 (panels in the third row). In panels a-c, dots are annual yield data of the common squid, and the dashed line is the maximum sustainable yield (MSY). In panels d-f, the solid line is the estimated annual biomass with carrying capacity (dotted line), and the dashed line is biomass that produces the MSY (B_MSY). In panels g-i, the solid line is the estimated annual harvest rate, and the dashed line is harvest rate at the MSY (H_MSY).

고찰

Case I을 통해 추정된 연도별 어획가능 자원량을 살펴볼 때, 한국 해역의 살오징어 어획가능 자원량은 2002년부터 급격하게 감소하였으며 2017년 이후로 B_MSY보다 낮게 추정되었다 (Fig. 4). 그리고 2016년, 589, 024톤에 이르던 살오징어 어획 가능 자원량이 2018년에는 350,360톤으로 약 41%감소하며, 살오징어 어획량 또한 121,691톤에서 46,274톤으로 급감한 사실은 살오징어 어업의 심각성을 시사하는 바이다. 그런데 매년 살오징어 총 어획량이 최대 지속가능한 어획량보다 낮게 기록되었고, 어획률이 H_MSY보다 낮게 추정되었음에도 불구하고, 어획 가능 자원량이 매년 감소하는 이유로는 어획강도의 변화 외에 다른 원인이 적용한 것으로 보인다. 일각에서는 이러한 어획 변동의 원인으로 중국어선의 불법조업 및 근해채낚기 어업과 트롤 어업간의 공동조업이 주요 원인이라고 발표했다(Kim and Shin, 2019). 그리고 기후체제 전환과 같은 기후변화와 해류의 변화가 살오징어 산란장 환경을 변화시키며 살오징어 자원량 변동에 큰 영향을 주었다고 보고한 바 있다(Sakurai et al., 2000; Sakurai et al., 2002; Choi et al., 2008; Kim et al., 2018).

본 연구에서는 복수어업으로 근해채낚기 어업, 대형선망 어업, 대형트롤 어업을 고려하였다. 세 가지 어업의 CPUE자료는 상업 어업에 의해 수집된 자료로서 어업 의존적인 자료이다. 어획능률계수 q의 경우 상수로 가정하였기 때문에 각 어업별 CPUE자료를 어획능률계수 q와 함께 어획가능 자원량에 비례하는 값으로 가정할 수 밖에 없었다. 어획능률계수 q를 상수가 아닌 시 변동성 모수로 고려하며 어업 선택성에 따라 변한다고 가정한다면 각 어업별 CPUE자료는 시 변동성 어획능률계수 q_t 와 함께 총 자원량에 비례하는 값으로 가정할 수 있다(Winker et al., 2020).

살오징어를 주 어획종으로써 조업하는 근해채낚기 어업과 대형 트롤 어업과는 달리, 대형선망 어업은 살오징어를 부수 어획물로 조업하는 어업이다. 그럼에도 불구하고 대형선망 어업 CPUE자료를 이용한 이유는 수집된 대형선망 CPUE자료 역시 나머지 두 어업과 동일하게 한국 해역의 총 살오징어 어획 가능 자원량의 상대적 지표로써 이용 가능하기 때문이다. 그뿐만 아니라, 2018년 기준 대형선망 어업의 살오징어 어획 비중은 약 4.2%로서 큰 비중을 차지하기 때문에, 한국해역의 살오징어 평가에 중요한 어업으로 적용 될 것으로 판단한다.

Case II에서 Case II-1, Case II-2, 그리고 Case II-3의 자원량 추정 결과를 살펴보면 2015년 이후로 살오징어 어획가능 자원량은 지속적으로 감소하고 있음을 추정할 수 있었다. 하지만 각 어업별 모수(b,r,K,q_jig,q_ps,q_tw,\(\sigma_{j i g}^{2}, \sigma_{p s}^{2}, \sigma_{t w}^{2} \sigma_{p}^{2}\) )의 점 추정치 및 어획가능 자원량 추정치를 구체적으로 살펴보면, 각 어업별 상당한 차이가 존재함을 알 수 있다. 이에 단일어업만을 고려한 자원평가(case II)에서는 관리 기준점을 명확히 확립하기 어렵다는 문제점을 도출할 수 있었다. 그러므로 살오징어 자원의 정확한 자원평가를 위해서는 단일어업이 아닌 복수어업을 고려한 자원평가가 필요함을 시사하는 바이다.

우리나라 주변 해역에 서식하는 살오징어는 동중국해, 동해, 서해 등 한국을 넘어 일본 EEZ 뿐만 아니라 북서 태평양 전체를 왕래하는 자원이다(Sakurai et al., 2013). 그런데 수천 킬로미터를 왕래하는 살오징어의 공간 시스템 아래에 아주 작은 부분에 지나지 않는 한국 해역만을 고려한 사실은 본 연구의 문제점 중 하나이다. 향후 살오징어 자원의 정확한 자원평가를 위해서는 닫힌 시스템을 확립할 필요성이 있고, 한국해역 밖에서 또는 밖으로 유입, 유출되는 개체군들을 고려할 필요성이 있다. 그리고 중국, 일본, 그리고 러시아와 같은 인접국가의 살오징어 조업정보 및 생태자료 등을 함께 고려한 공동 자원평가 연구가 필요하다.

본 연구는 복수 어업(근해채낚기 어업, 대형선망 어업, 대형트롤 어업)에서의 살오징어 자료(어획량과 CPUE자료)를 이용하여 베이지안 상태공간 잉여생산량 모델에 적용하며, 살오징어 어획 가능 자원량을 추정하였다. 비록 복합적인 가정을 적용하여 모수를 추정했기 때문에 실제 살오징어 어획가능 자원량을 정확히 나타낸다고 할 수는 없다. 하지만 복수어업을 고려했을 때와 단일어업만을 독립적으로 고려했을 때의 자원평가 결과를 비교 및 분석하며 복수어업을 고려한 자원평가의 필요성을 증명한 연구라는 점에서 의의를 가진다.