単位 日语量あたりの大きさ 什么意思

一般化状態空間モデルで漁業動態を記述する ―マサバ努力量管理効果の定量評価 (PDF Download Available)
See all >23 References
19.96Fisheries Research Agency28.45Fisheries Research AgencyAbstractシミュレーションモデルを用いて様々な資源管理の方策を評価するmanagement strategy evaluation(MSE)は,水産資源や野生生物の管理を主な目的の一つとする水産資源学や応用生態学において重要なツールとなりつつある.ここでは,マサバの努力量管理のMSEにおいて統計モデルが用いられた一つの研究例を紹介する.この研究は,漁獲努力量(操業日数)と漁獲量の関係を確率モデルで表し,管理によって獲り控えられた漁獲分が個体群動態を通してその後の資源回復にどれだけ寄与したかをシミュレーションにより評価したものである.結果として,日々の漁獲量と出漁隻数を表す時系列データに潜む自己相関構造と,資源の増減に反応した漁業者の努力量の変化が努力量管理の効果に大きな影響を与えていることが明らかになった.通常のMSEでは努力量と漁獲量の関係(漁業動態)が単純な線形関係で表現されることが多いが,本研究の成果は漁業動態を実データに即した統計モデルで表現することの重要性を示している.今後,資源管理の分野の中で統計モデルが積極的に使われ,より実際的な管理方策の評価に繋がることを期待する.Discover the world's research14+ million members100+ million publications700k+ research projects
統計数理(2016)第64 巻第1号59–75c?2016 統計数理研究所特集「生態学における統計モデリング」[研究詳解]
一般化状態空間モデルで漁業動態を記述する—マサバ努力量管理効果の定量評価市野川 桃子+?岡村 寛+(受付 2015 年7月16 日;改訂 12 月7日;採択 2016 年1月4日)要旨シミュレーションモデルを用いて様々な資源管理の方策を評価する management strategy eval-uation(MSE)は,水産資源や野生生物の管理を主な目的の一つとする水産資源学や応用生態学において重要なツールとなりつつある.ここでは,マサバの努力量管理の MSE において統計モデルが用いられた一つの研究例を紹介する.この研究は,漁獲努力量(操業日数)と漁獲量の関係を確率モデルで表し,管理によって獲り控えられた漁獲分が個体群動態を通してその後の資源回復にどれだけ寄与したかをシミュレーションにより評価したものである.結果として,日々の漁獲量と出漁隻数を表す時系列データに潜む自己相関構造と,資源の増減に反応した漁業者の努力量の変化が努力量管理の効果に大きな影響を与えていることが明らかになった.通常の MSE では努力量と漁獲量の関係(漁業動態)が単純な線形関係で表現されることが多いが,本研究の成果は漁業動態を実データに即した統計モデルで表現することの重要性を示している.今後,資源管理の分野の中で統計モデルが積極的に使われ,より実際的な管理方策の評価に繋がることを期待する.キーワード:資源管理,水産資源,管理方策評価,マサバ,漁業動態.1. はじめに:水産資源学の中での統計モデルの利用水産資源の管理は水産資源学という一つの学問分野の中で取り扱われるが,野生生物資源を管理するという点で応用生態学の中の一つの分野とも言える.但し,水産資源の管理を目的として行われる解析や研究には三つの特徴的な点がある.ひとつめは,陸上の生物と違い,管理の対象となる水産生物を目で見て個体数の計数をすることが困難という点である.この問題に対して,水産資源学の分野では個体群動態モデルを漁業データや調査データにあてはめて個体数を推定する「資源評価」が昔から行われてきた(Hilbornand Walters, 1991; Quinn and Deriso, 1999).そして,資源評価手法は最尤法の枠組みを取り入れることで 1990 年代から飛躍的に発展した(Quinn, 2003).2000 年代になると,漁獲量?漁獲物の体長組成?相対資源量指数(対象資源の相対的な資源量を示す指数)といった複数の質の異なるデータから得られる尤度を統合して資源評価を行う統合型資源評価(integrated analysis)が盛んに利用されるようになった(Maunder and Punt, 2013; 市野川 他, 2015b).統合型資源評価には未だ多くの課題があるものの,個体群動態や漁業データの収集のプロセスをより柔軟に表現でき,それらにまつわる不確実性を資源評価結果の不確実性として定量評価することができる有+水産総合研究センター 中央水産研究所:〒 236–8648 神奈川県横浜市金沢区福浦 2–12–4
60 統計数理 第 64 巻第1号2016図1.コンピュータシミュレーションを用いて効果的な管理方策を探索する MSE(ManagementStrategy Evaluation)の概念図.漁業資源管理に関わる 1?~5?のプロセス全てをコンピュータ内にシミュレーションモデルとして再現し,様々な管理方策を試す.個体群動態や資源評価,管理の実施について複数のシナリオを用意し,どのシナリオでも平均的にうまく働く管理方策を選ぶことで,不確実性に対して頑健な管理方策を選ぶことができる.用な手法である.同様の考え方で個体数推定を行うモデルは個体群生態学の分野でも integratedpopulation model(Schaub and Abadi, 2011)と呼ばれ,利用されるようになっている.水産資源学におけるふたつめの特徴は,対象とする生物の豊度に関する情報が漁業を通して得られる(漁業データ)ことが多いという点である.漁業データの例としては,漁獲量や漁獲物の体長組成,CPUE(catch per unit of effort,単位努力量あたりの漁獲量)などがある.しかし,科学的に設計された調査や実験と異なり,いつ?どこで漁獲が行われるか,つまり,いつ?どこで漁業データが得られるかは海況?漁業者の行動?社会経済的な要因等に影響される.そのため,漁業データには欠損やサンプル数の偏りの問題が常に付きまとう.この問題に対処するため,水産資源学の分野では比較的早い時期から様々な統計的手法が応用されてきた.特に,漁業から得られる CPUE データにおけるサンプルサイズの時空間的な偏りを統計的手法によって補正し,相対資源量指数を抽出するプロセスは CPUE 標準化と呼ばれている.一般化線形モデル?一般化線形混合モデル?一般化加法モデル?機械学習法などの様々な統計学的手法が CPUE標準化の際に用いられている(Maunder and Punt, 2004).漁業データ解析や資源評価で用いられる統計解析手法については,本誌でより詳しく総説されている(岡村?市野川, 2016).本稿では,水産資源管理に向けた研究の中にある三つめの特徴?問題点とそれに対する近年のアプローチ,さらに,その中で統計モデルが効果的に使われた一つの研究例を紹介したい.三つめの特徴?問題点とは,漁業という産業が関わるため,様々な資源管理方策が提案されたとしても,複数の管理方策を実験的に導入して比較したりすることが難しいという点である.そして,この問題を克服するツールとして,対象とする資源や生態系の動態とそれに対する評価?管理プロセスをシミュレーションモデルとしてコンピュータ上に再現して効果的な管理方策を探索する手法(management strategy evaluation,MSE,図 1)が近年注目を浴びている.現実の世界で複数の管理方策を一つ一つ試すのは不可能でも,コンピュータ上に再現した資源に対してなら何度でも,どのような方策でも,失敗を恐れずに試すことができるのである.
一般化状態空間モデルで漁業動態を記述する 61MSE の端緒は持続的な商業捕鯨のための管理方策を探索する目的で国際捕鯨委員会の科学委員会が開発したものである(Punt and Donovan, 2007).商業捕鯨が再開されていない現在,そこで開発された管理方策は実用に至っていないが,その後,MSE の手法は様々な漁業資源管理で利用されることとなった.MSE が資源管理に実際に貢献した代表的な例としてはミナミマグロの管理が挙げられる.MSE を用いて開発されたミナミマグロの管理方策は 10 年以上の年月をかけて完成し,管理に適用され,ミナミマグロ資源の回復に貢献している(Hillary et al., 2015;黒田 他, 2015).また,我が国の水産資源管理においても,データ不足のせいで資源量が推定できない資源には,簡単な MSE を用いて決定された管理方策が適用されている(Ohshimo andNaya, 2014; 市野川 他, 2015a).MSE は水産資源管理の中で今後より一層重要なツールになっていくと期待されており,陸上の野生生物管理への応用も提案されている(Bunnefeld et al., 2011).MSE では,対象とする資源で既に得られている漁獲データを再現するようにパラメータを調整される(conditioning)ことが多いが,この際に尤度が用いられ,統計的な手法が使われる.また,仮定が異なる複数の個体群動態シミュレーション結果を尤度で重みづけして統合するといったことも行われている.しかし,MSE の中で研究がまだ十分に進んでいない分野がある.それは,漁獲に対する努力量と漁獲量の関係(漁業動態)のモデル化である。資源管理は漁獲量の制限だけでなく,漁獲努力量を制限すること(努力量管理)によっても行われることがある.具体的には,禁漁区の設置や漁獲日数の制限などである.このような管理手法の効果を MSE で定量化する場合,例えば,漁獲日数の削減がどの程度の漁獲量の削減に繋がるかをモデル化する必要がある.また,管理の導入や資源の増減によって漁業者の行動が変わり,努力量が変化する可能性がある.例えば,資源の増加によって漁業者のモチベーションが上がり,より長い時間操業するようになる,といった行動の変化があるかもしれない.このように,漁業者の漁獲努力量は資源管理の導入やその他の様々な要因によって変化し,それによって漁獲量も変化する.これらのプロセスをここでは漁業動態と呼ぶことにする.そして,漁業動態のモデル化は努力量管理の効果を MSE で定量化する場合に重要になってくる(Bunnefeld et al., 2011).しかし,今まで,実際を反映した漁業動態が MSE に導入され,その影響が評価された事例はほとんどなかった.本稿で紹介する研究例は,統計モデルを用いて実際の漁業動態(特に,努力量の時間的な分布と漁獲量の関係性)を確率モデルとして明示的に記述し,漁業動態のモデル化によって漁業管理の有効性がどの程度変わるかを定量化したものである.この例を通して,今後の MSE 研究の中で,漁業動態を明らかにするために統計モデルを利用すること,そして,様々な管理方策を定量的に評価することの重要性を示したい.2. サバの管理と統計モデル:太平洋マサバ資源に対する休漁管理効果の評価2.1 背景水産資源の管理手法には大きく分けて三つのタイプがある(Pope, 2009).一つは技術的コントロールと呼ばれ,漁船の大きさや漁具の形状といった魚を漁獲する手法そのものに制限をかける方法である.二つ目は努力量管理(インプット?コントロール)と呼ばれ,漁船の操業時間や日数?操業場所等を制限する方法である.三つ目は漁獲量管理(アウトプット?コントロール)と呼ばれ,設定された漁獲量以上に漁獲を許さないという方法である.全ての水揚げ量を逐一把握する必要があることから,漁獲量管理は資源管理の手法として最もコストがかかる.一方,技術的コントロールや努力量管理は比較的実施しやすいという利点があり,世界や日本で古くから利用されてきた(e.g., Worm et al., 2009; Makino, 2011; Melnychuk et al., 2013).日本では,広域に分布し,漁獲量が多い 7魚種(マイワシやマサバ等)を対象として漁獲量管理が
62 統計数理 第 64 巻第1号2016実施されているが,その他の多くの魚種の管理は努力量管理によって実施されている(Makino,2011).これは,努力量管理が漁獲量管理に比べて実施しやすく,古くから使われてきた管理手法であるために,利害関係者間で管理方策に関する合意を得やすいことによると思われる.しかし,努力量管理や技術的コントロールは間接的にしか漁獲量を制限できないため,これらの管理が実際にどのくらいの管理効果を持つかを定量的に評価することが困難という欠点がある.MSE の枠組みの中で努力量管理の効果を定量化しようとする際には,管理の実施によって削減された努力量が実際にどのくらいの漁獲量の削減に結びつくかという漁業動態を明示的にモデル化してやる必要がある.しかし,実際の漁業動態のモデル化は一般に難しい.通常の努力量管理は漁獲日数や操業時間等を制限することによってなされる.しかし,同じやり方で1日操業をしたとしても,資源量や海況?その日の魚の分布といった様々な要因によって 1日あたりの漁獲量は変動するだろう.そのため,漁獲日数の 1日の削減がどのくらいの漁獲量の削減に繋がるかは,これらの要因を考慮して予測する必要がある.また,管理の導入によって漁業者の行動が変化し,制限された努力量が他の部分の努力量の増加によって補填されてしまうことがある(Salas and Gaertner, 2004; Fulton et al., 2011).もし,1日分の操業を削減する,としか決められていない場合,削減された 1日分を取り戻すため,他の日により長い時間操業するようなことが起きるかもしれない.このような漁業者の行動の変化で有名な例は,禁漁区の設置によって禁漁区周辺の漁獲努力量が増加するような現象である(Kellner et al., 2007; vander Lee et al., 2013).適切な努力量管理の評価のためには,漁業者の行動の変化を考慮して漁業動態をモデル化する必要がある.このような評価の困難さが原因となり,日本で古くから実施されてきた努力量管理の効果は今まで定量的に評価されることがなかった.また,世界で実施されている漁業管理に目を向けても,漁獲量管理や技術的コントロールといった他の管理方策に比べて努力量管理の有効性はあまり評価されていない(Worm et al., 2009).しかし,本研究で紹介する Ichinokawa et al.(2015)は,漁業動態を記述した統計モデルと個体群動態モデルを組み合わせることにより,太平洋で漁獲されるマサバ資源(マサバ太平洋系群)に対して実施された努力量管理の効果を定量化することに成功した.マサバ太平洋系群に対する努力量管理は 2003 年から実施され,投網時刻や投網回数の制限に加え,1日のサバの総水揚げ量が閾値(2000~3000 トン)以上あった翌日を臨時の休漁日とすることで操業日数が削減されてきた.本研究が解析の対象とした 2003 年10 月から 2009 年6月までには,のべ 131 日が臨時休漁となった.本資源に対する公の管理としては,漁獲枠以上の漁獲を許さない漁獲量管理が1997 年から実施されていた.しかし,上限となる漁獲量の総枠は実漁獲量よりもはるかに高く(Swartz and Ishimura, 2014),漁獲枠が実質的な意味をなしていなかったため,2004 年から臨時休漁による努力量管理も併せて導入されたという経緯がある.では,これらの管理の実施に対してマサバの資源量はどのように変動したのだろうか? 図 2に,マサバ太平洋系群の推定資源量(川端 他, 2013)を示した.この推定によると,マサバの資源量は 1980 年以前に 300 万トン以上あったが,1980 年代に減少し始め,1980 年代後半の低い加入尾数(0歳魚の生き残り数)と1990 年代の高い漁獲率によって 1990 年代から 2000 年代前半に10–20 万トンにまで減少した.しかしその後,資源量は増加傾向となった.資源管理の歴史と併せて図 2を見ると,漁獲量管理が開始された 1997 年以降も資源量は低迷したままだったのに対し,努力量管理が開始された 2004 年以降は資源が回復しており,あたかも努力量管理がマサバ資源を回復させたような印象を受ける.しかし,加入尾数を見ると(図 2右),2004 年は特に加入尾数が多い年でもあった.加入尾数は親の量だけでなく,環境条件の自然な変動に大きく左右されるため,2004 年以降の資源回復は,2004 年にたまたま環境条件が良く,多い加入尾数が得られたことが原因だと見る人もいるかもしれない.マサバの資源回復が管理のおかげなのか,それとも単なる偶然なのか,科学的に答えること
一般化状態空間モデルで漁業動態を記述する 63図2.マサバ太平洋系群の資源量(左)と加入尾数?親魚資源量あたりの加入量(RPS)?漁獲率(右).漁獲率は灰色の実線で示している.Ichinokawa et al.(2015)の Fig. 2 を改変.ここでの年は漁期年で表し,例えば 2004 年7月から翌年の 6月までを 2004 年の漁期年としている.が求められていた.そこで,臨時休漁を用いたマサバの努力量管理の効果を定量的に評価するため,次の 2つの問いに答えることにした.一つ目は漁業動態に関するもので,臨時休漁による 1日分の努力量の削減がどのくらいの漁獲量の削減に繋がるか?という問いである.この問いに答えるため,努力量管理が実施されていた 2004 年から 2008 年までの日々の漁獲データ(いつサバ操業が行われるか,サバ操業が行われる場合何隻の船が水揚げし,漁獲量がどのくらいか)を様々な共変量(資源量?月?自己相関構造)で説明する確率モデル(一般化状態空間モデル)を作成した.特に,臨時休漁分を補填しようとして休漁日の翌日の努力量や漁獲量が多くなる可能性が指摘されていたため,共変量の一つに「前日が臨時休漁日であるかどうか」を加え,そのようなことが実際に起こっていたかを統計的に検討できるようにした.二つ目の問いは,獲り控えられた漁獲分が実際に資源回復にどの程度貢献したか?という問いである.これは,漁業動態に加えて個体群動態を同時に考えた場合にどうなるか?という問いでもある.これに答えるため,前述の解析で作成した確率モデルを用いて臨時休漁によって獲り控えられた漁獲分を予測し,さらに,個体群動態モデルを用いて獲り控えられた漁獲分がサバの資源量の増加にどのくらい寄与したかを予測した(図 3).この簡易的な MSE を用いて,臨時休漁を実施した場合としなかった場合で資源量や漁獲率を比較し,臨時休漁の効果を定量化した.研究の詳細は Ichinokawa et al.(2015)に原著論文として発表されているので,ここではこの研究における 2つの重要な成果,(1)確率モデルによって記述された漁業動態を MSE の中で利用することで努力量管理の効果を定量化できること,(2)現実のデータに即した漁業動態が一般化状態空間モデルによって記述され,その結果,特に漁獲量や努力量の自己相関構造と漁業者の行動の変化に関するモデルの違いが管理効果の評価を大きく変えること,の 2点に絞って内容を紹介したい.
64 統計数理 第 64 巻第1号2016図3.漁業動態モデルと個体群動態モデルを組み合わせたシミュレーションの概要.漁業動態モデルは漁獲量を予測するためのモデルで,実データをもとにパラメータ推定した一般化線形モデルまたは自己相関を考慮した一般化状態空間モデルで記述される.また,ここで示している個体群動態モデルは概要の説明のために単純化したもので,実際には年齢別モデルを使用している.図4.巻き網漁業によって漁獲されたマサバ太平洋系群の日別の漁獲量(灰色の棒グラフ)とその日の水揚げ隻数(棒グラフ上の数字)の例.点線は臨時休漁の閾値.Ichinokawa et al.(2015)の Fig. 3 を改変.2.2 漁業動態を記述する統計モデルと個体群動態モデルのカップリング2.2.1 漁業動態モデル図4に,太平洋の三陸沖から銚子沖にかけて漁獲された日別のサバ類漁獲量と水揚げ隻数のデータの例を示した.この図から,漁獲量の閾値を超えた日の翌日は臨時休漁になるが,必ずしも臨時休漁になっていない日もあること,臨時休漁以外でも悪天候や市場休場等の理由により漁獲が全く行われない日もあること,漁獲がある場合でも日々の漁獲量,水揚げ隻数は大きく変動していることがわかる.つまり,日々の漁獲は様々な確率プロセスの組み合わせで成り立っていると考えることができる.ここではそれを以下のように定式化した.まず,ある t日において,その日が臨時休漁日となるかどうか?が最初に決められる.そこで,t日が臨時休漁である場合に 1,そうでない場合に 0となる確率変数 rtを考え,それは
一般化状態空間モデルで漁業動態を記述する 65表1.休漁?漁獲の有無?漁獲量?隻数モデルに導入した共変量.rt~Bernoulli(qt)に従うとした(休漁モデル).次に,t日が臨時休漁でない日でも,他の要因によってサバ漁獲が行われない場合がある.そこで,その日のサバ漁獲の有無を,サバ漁獲がある場合に 1,ない場合に 0となる確率変数 st~Bernoulli(ut)によって表現した(漁獲の有無モデル).さらに,臨時休漁でなく,かつ,サバ漁獲が行われる日の漁獲量 ctをct~Gamma(μt,δ)とした(漁獲量モデル).これより,ある t日の漁獲量 Ctは,これらの確率変数の積(2.1) Ct=(1-rt)stctとして表現できる.それぞれの確率変数 rt,st,ctで用いられるパラメータ qt,ut,μtは(2.2) f(θt)=α+βTxtのように,共変量 xtとそれらに対する係数パラメータ βと,切片 αの線形結合として表現した.ここで θtはパラメータ qt,ut,μtのどれかで,fはリンク関数を示す.リンク関数は,qt,utでlogit,μtでlog リンク関数を用いた.パラメータ α,βは漁獲データ(図 4)を用いて最尤法により推定した.それぞれのモデルに導入する共変量は,今年または去年のマサバの資源量(cbiom または pbiom),月(MON),前日が臨時休漁である(pC),前日にサバ漁獲がない(pF),その日の水揚げ隻数(L),臨時休漁となる漁獲量と前日の漁獲量の差の絶対値(D),前々日に漁獲量が閾値を超え,かつ,前日が臨時休漁でない(ppC)を考えた(表 1).pC 効果は「臨時休漁の実施を補填するために臨時休漁の翌日の努力量や漁獲量が多くなる可能性」を検討するために導入した共変量である.また,水揚げ隻数(L)はそれ自体も様々な要因に影響され,日によって大きく変動する.そこで,t日の水揚げ隻数も,Lt-1~Bin(NLt-1,lt)に従う確率変数として扱った(隻数モデル).ここで,NLtは,当該海域でまき網操業の許可を持つ漁船団の総数で,Ltで使われるパラメータ ltはlogit リンク関数を仮定した式(2.2)により定式化し,日別の水揚げ隻数データよりパラメータ推定を行った.それぞれの共変量をどのモデルに導入するかについては,それぞれの応答変数の特徴から可能性のある共変量を事前に選択した(表 1).その後,
66 統計数理 第 64 巻第1号2016図5.時系列データに対する GASSM とGLM の予測の違い.GASSM では自己相関構造を状態変数として推定することができるため,漁獲データが存在しない臨時休漁日の漁獲量や努力量を高い精度で予測することができる.また,自己相関を考慮するかどうかによって,臨時休漁日の翌日の漁獲量?努力量に対する解釈が GASSM とGLM で大きく変わってくる.bias-corrected Akaike Information Criterion(AICc)を用いて,累積 Akaike weight が90%以上となる候補モデルを用いてモデル平均を行った(Burnham and Anderson, 2002).但し,図 4のデータをよく見ると,日々の漁獲量や水揚げ隻数の変動は完全にランダムでなく,漁獲量や水揚げ隻数が多い日が連続して続くような傾向が見受けられた.サバ漁場となる東北沖においては数日~数十日スケールで好漁場が形成されることが知られており,その間には特に漁獲量や出漁隻数が多くなっているものと思われる.そこで,漁獲量モデルと隻数モデルについては式(2.2)を式(2.3)のように拡張することで,自己相関構造を考慮できるようにした.(2.3) f(θt)=α+βTxt+yt,yt=ρ1yt-1+ρ2yt-2+?t,?~N(0,σ2)ここで ytは1日前(t-1)と 2日前(t-2)の状態を反映する状態変数となっており,自己相関係数ρ1,ρ2が大きいほど前日?前々日の状態の影響を強く受けることを示す.ρ1,ρ2は,ytが定常状態となるように |ρ1|&1, ρ1+ρ2&1, ρ2-ρ1&1の制約条件をつけた.式(2.3)のモデルは,観測できない状態変数 ytと共変量に対するパラメータ βを同時に推定するものとなっており,ここでは自己相関を考慮した一般化状態空間モデル(Generalized Autoregressive StateSpace Model, GASSM)と呼ぶことにする.一方,自己相関を考慮しない式(2.2)は通常の一般化線形モデル(Generalized Linear Model, GLM)に対応する.時系列データに対するそれぞれのモデルの予測の概念は図 5に示した.2.2.2 漁業動態と個体群動態モデルのカップリング以上より,個々の漁獲イベントを示す確率モデル(休漁?漁獲の有無?漁獲量?隻数モデル)の積としてまき網船の漁業動態を確率プロセス(式(2.1))で記述することができた.次に,この漁獲プロセスがマサバ個体群にどのような影響をあたえるかを調べるため,漁業動態と個体群動態モデルをカップリングさせたシミュレーションを行った(図 3).マサバの個体群動態モデルはマサバの資源評価(川端 他,2013)で用いられているものと同じ年齢別モデルを用いた(Ichinokawaet al., 2015)が,ここでは年齢構造を含まない単純化した式で概要を説明する.まず,y年のマサバの個体数(Ny)は前年の漁獲による死亡係数(Fy-1)と自然の原因による死亡係数(M,0.4 を仮定)により減耗する一方で,新規の加入 Ryを得る(式(2.4)).(2.4) Ny=Ny-1exp(-Fy-1-M)+Ry
一般化状態空間モデルで漁業動態を記述する 67ここで,漁獲死亡係数 Fyは総漁獲量(TC y)を導く下式(Pope, 1972)を解くことで求められる.(2.5) TCy=Ny(1 -exp(-Fy)) exp(-0.5M)WyWyはy年におけるマサバの平均体重である.TC yは年間の総漁獲量であるから,式(2.1)を用いてランダムに発生させた日々の漁獲量を 1日目から 365 日目まで積算することで TC yを決めることができる.つまり,TCyが与えられれば Fyが決まり(式(2.5)),Fyが決まれば翌年の資源尾数 Ny+1 を決定できる(式(2.4)).このように,確率モデルで表された日々の漁獲量の予測値を個体群動態モデルに導入することで,今年の漁獲量の多い?少ないが,翌年の資源尾数にどのくらい影響を与えるかを見積もることができる.さらに,翌年の個体数に年齢別の平均体重を乗じて翌年の資源重量を算出し,それを翌年の漁獲量(Ct)を予測する際の共変量として用いることで,前年の漁獲量の多い?少ないが,資源量の増減を介して翌年の漁獲量にフィードバックされることとなる.シミュレーションでは,式(2.1)から発生させた漁獲量を用いた休漁管理がある場合のシナリオと,式(2.1)で rt≡0とした休漁管理がない場合のシナリオの両方を実施した.そして,休漁ありに対する休漁なしの資源量?漁獲量?漁獲率(資源量に対する漁獲量の比)の比を算出し,それを休漁の実施による管理効果とした.具体的には,漁獲量と漁獲率については,シミュレーションを実施した 2004 年から 2008 年までの平均値の比に,資源量については 2009 年当初の資源量の比に注目した.シミュレーションは休漁あり?なしのシナリオを 1000 回ずつ繰り返した.2.3 管理効果に影響を与える要因2.3.1 結果の概要漁獲の有無モデル,隻数モデル,漁獲量モデルで推定されたパラメータをそれぞれ表 2–4 に示した(休漁モデルの結果は省略した).ここではモデル平均によって推定されたパラメータに加え,感度分析で用いるため,AICc 最小モデル?AICc 最小モデルに「前日が臨時休漁(pC)」の効果を加えたモデル?AICc 最小モデルから資源量(cbiom,pbiom)の効果を除いたモデルの結果も示した.また,隻数モデルと漁獲量モデルについては,自己相関を考慮しない単純な GLMを用いた場合のモデル平均と AICc 最小モデルの結果も示した.モデル平均から得られたパラメータを用いて臨時休漁のあり?なしでシミュレーションした結果(ベースケース)を図 6に,ベースケースの結果に加えて,自己相関構造や資源量の効果を除いたモデルを組み合わせてシミュレーションした感度分析の結果を表 5にまとめた.まず,ベースケースにおけるシミュレーションの結果(図 6)を概説する.漁獲量については,本来なら臨時休漁となる日に漁獲を行うことができるため,休漁なしの場合の
年の漁獲量は休漁ありの場合の漁獲量を上回った(図 6(a)).しかし,最初の 3年により多く漁獲した分,休漁なしの場合の総資源量は休漁ありの場合に比べて少なくなった(図 6(c)).サバ漁獲日,水揚げ隻数,1日あたりの漁獲量は総資源量に対して正の関係があると推定された(表 2–4)ため,総資源量の減少は漁獲量?漁獲日数の減少に繋がり, 年の休漁なしの漁獲量は休漁ありの漁獲量よりも少なくなった.結果として,2009 年当初の資源量は休漁ありの場合がなしの場合に比べて 1.4 倍と見積もられた(表 5).つまり,休漁による努力量管理は,5年間でマサバの資源を 1.4 倍に増やすのに確かに貢献したことが示されたのである.但し,年の総資源量は休漁がなかった場合でも増加していた.つまり, 年の総資源量の増加は 2004 年の大きな加入尾数にも依存しており,臨時休漁による管理だけが貢献していたわけではないことも同時に示された.2.3.2 自己相関構造の影響以上が結果の概要であるが,様々な統計モデルを試した感度分析により,自己相関を考慮し
68 統計数理 第 64 巻第1号2016表2.漁獲の有無モデルで推定された係数と AICc.ΔAICc はAICc 最小モデルとの AICcの差.表中の “**” は,モデルに導入されているが推定された係数を省略していることを示す.“-” は,そのモデルに導入されていないか,または,計算不可であることを示す.表3.隻数モデルで推定されたパラメータと AICc.ΔAICc はAICc 最小モデルとの AICcの差.表中の “**” は,モデルに導入されているが推定された係数を省略していることを示す.“-” は,そのモデルに導入されていないか,または,計算不可であることを示す.表4.漁獲量モデルで推定されたパラメータと AICc.ΔAICc はAICc 最小モデルとの AICcの差.表中の “**” は,モデルに導入されているが推定された係数を省略していることを示す.“-” は,そのモデルに導入されていないか,または,計算不可であることを示す.
一般化状態空間モデルで漁業動態を記述する 69図6.休漁あり(白抜きの箱型図)と休漁なし(灰色の箱型図)の場合の年間漁獲量(a),漁獲率(b),総資源量(c),親魚資源量(d)のシミュレーション結果.白抜きの菱形は 2012 年の資源評価で推定された実際の値.モデル平均によって推定されたパラメータを用いた(表 5, S1)場合の結果.実際には休漁がなされているので,白抜きの箱型図と菱形が近いほどシミュレーションが現実を良く表現していることになる.Ichinokawa et al.(2015)のFig. 7 を改変.ない単純な GLM と自己相関を考慮する GASSM とでは管理効果の評価結果が大きく変わることが明らかとなった.隻数?漁獲量モデルでは GLM での AICc 最小モデル(表 3のb6,表 4のc6)よりも GASSM での AICc 最小モデル(表 3のb2,表 4のc2)で AICc が大幅に小さくなっているため,GASSM はGLM よりもデータをより良く説明できるモデルであるのが明らかである.実際,モデル平均で推定されたパラメータでランダムに出漁隻数(図 7)と漁獲量(図 8)を発生させると,GASSM によるモデルのほうが観測値をより良く説明できている.そして GLM とGASSM ではパラメータの推定値に大きな違いがあった.隻数?漁獲量モデルで推定された「前日が臨時休漁(pC)」の係数は,GLM でそれぞれ 0.70,0.50(モデル平均)だったのに対し,GASSM では 0.09 と0.02(モデル平均)と非常に小さくなったのである.また,GASSMによる AICc 最小モデルにはこれらの効果が含まれていない.pC は「臨時休漁の実施を補填するために臨時休漁の翌日の努力量や漁獲量が多くなる可能性」を検討するための共変量であった.そのため,GLM とGASSM でpC に対するパラメータの推定値が大きく異なることは,休漁の実施に対する補填行動があるかどうかという仮説に対して GLM とGASSM が全く異なる結論を示していることを意味している.単純な GLM を用いた場合の結論は「臨時休漁の翌日に漁獲量も出漁隻数も大幅に増加する」である一方,GASSM を用いた場合の結論は「そのようなことはないか,増加したとしてもその程度は小さい」である.但し,AICc の差から考えると,後者の仮説の方がより尤もらしい.
70 統計数理 第 64 巻第1号2016表5.様々な統計モデルを組み合わせた場合のシミュレーション結果.「漁獲量の比」,「漁獲率の比」は, 年の休漁ありの場合の平均漁獲量または平均漁獲率を休漁なしの場合の値で割ったもの.「資源量の比」は2009 年当初の休漁ありの場合の資源量を休漁なしの場合の値で割ったもの.5%,50%,95%は,休漁あり?なしシナリオの各 1000回ずつの確率的シミュレーションの結果からランダムに 1組結果を抽出し,その比をとったときのパーセンタイル.モデル平均を組み合わせた結果(S1)と ΔAICc の合計値が小さい S2,S3 がより尤もらしい結果で,その他のシナリオ(S4–S6)は感度分析や他の管理方策を試すために実施したもの.単純な GLM を用いた場合に GASSM と全く異なる結論が得られたのは何故なのだろうか? それは GLM が時系列データの自己相関構造を考慮できないからである.観測値に対する GASSMの予測力の高さや GASSM とGLM の間の AICc の差からわかるように,日々の漁獲量と出漁隻数の時系列データには明らかな自己相関構造があった.そして,臨時休漁は漁獲量がある閾値を超えた場合(出漁隻数が多く,漁獲量も多い日)の翌日に実施されることになっている.そのような条件の下で日々の漁獲量や努力量に強い自己相関構造があれば,臨時休漁の翌日もやはり出漁隻数や漁獲量は多くて当然というのが GASSM による解釈である(図 5).しかしこの自己相関構造を考慮できない GLM では,「臨時休漁の翌日」の漁獲量?出漁隻数は GLM による期待値よりも常に高いところにあるとなり,pC 効果を過大に見積もる結果となる.これらの解釈の違いにより,GLM をもとにしたシミュレーションから見積もられた管理効果の大きさは GASSM の結果と大きく異なった(表 5).GASSM をもとにしたシミュレーション結果からは 2009 年当初におけるマサバ資源量は休漁がなかった場合の約 1.4 倍と見積もられたのに対し(表 5, S1),GLM で推定されたパラメータをもとにしたシミュレーションの結果(表 5,S5)では,休漁ありの場合となしの場合で資源量にほとんど差が見られなかった.これは,臨時休漁の翌日に大幅に漁獲量と努力量が増えると GLM が予測し,臨時休漁によって削減された漁獲量がこれによって補填されてしまったためである.GASSM によるデータの解釈が真実に近いとするなら,漁獲が多かった翌日を臨時休漁とする管理方策は,獲り控えの効果が特に大きい日に特定して休漁を実施することになり,ランダムに(または一定間隔で)休漁日を設定するよりも少ない休漁日数で大きな効果を上げる方策と考えることができる.実際,本シミュレーションにおいて,臨時休漁の発生の割合を変えずに休漁日の配置をランダムにした場合(表 5, S6)には,臨時休漁をベースケースと同程度の日数
一般化状態空間モデルで漁業動態を記述する 71図7.2004 年から 2008 年までの日別のサバ水揚げ隻数の観測値(白丸)と隻数モデルによる1000 回のシミュレーション結果の 90 パーセンタイル(灰色).(a)自己相関を考慮したGASSM をもとにしたシミュレーション結果.(b)自己相関を考慮しない GLM をもとにしたシミュレーション結果.Ichinokawa et al.(2015)の Fig. 6, Fig. C1 を改変.図8.2004 年から 2008 年までの日別のマサバ漁獲量の観測値(白丸)と漁獲量モデルによる1000 回のシミュレーション結果の 90 パーセンタイル(灰色).(a)自己相関を考慮したGASSM をもとにしたシミュレーション結果.(b)自己相関を考慮しない GLM をもとにしたシミュレーション結果.Ichinokawa et al.(2015)の Fig. 6, Fig. C1 を改変.だけ実施しているにも関わらず,休漁がないときの資源量に対する休漁ありの資源量比が 1.19と,ベースケースの 1.4 よりも小さかった.2.3.3 資源量の増減に対する努力量?漁獲量の変化マサバの努力量管理では,自己相関構造の他にもう一つ,その効果に大きな影響を与えていた要因があった.それはマサバの資源量(cbiom,pbiom)の効果である.これらの効果は漁獲の有無モデル(表 2)?隻数モデル(表 3)?漁獲量モデル(表 4)の全てで比較的大きな係数が推定された.漁獲の有無モデルと隻数モデルは,いつ?どのくらいの数の船が出漁するか,といった漁業者の漁獲努力量を示すモデルである.特にこれらのモデルで資源量の係数が正の大きい値
72 統計数理 第 64 巻第1号2016として推定されたことは,マサバの資源量が増えると漁業者の漁獲努力量も増える傾向があることを意味している.では,このような資源量と漁業者の努力量との関係を見落としてしまった場合,管理効果の評価はどのくらい変わるだろうか? そこで,努力量に関するモデル(漁獲の有無モデルと隻数モデル)に資源量の効果を入れないでパラメータ推定したモデル(表 2, a4; 表3, b4)を用いて,休漁あり?なしのシミュレーションを行った(表 5, S4).その結果,2009 年当初の資源量は休漁なしの場合の資源量の約 1.6 倍となり,ベースケースで推定された 1.4 倍よりも管理効果が高く見積もられることとなった.資源量の効果を除いたモデルと AICc 最小モデルの AICc の差は32.9 と大きいため,このシナリオは現実的でない.しかし,管理効果をより高めるための方策を検討する際,このような検討結果が役に立つ.つまり,臨時休漁による管理と同時に,資源量が増加しても出漁隻数や出漁日数が増加しないように上限を設ける管理を取り入れることで,管理効果をより高められることが期待されるのである.3. まとめ以上,実際に行われた努力量管理の効果評価の中で,現実に即した漁業動態をモデル化しているかどうかによって管理効果の評価が大きく変わった例を紹介した.特にこの例では,漁獲量や水揚げ隻数の時間的な自己相関構造と,管理の実施による漁獲量の増減が個体群動態を介して翌年の漁業者の努力量に影響を与えるフィードバック構造が管理効果に大きな影響を与えていることが明らかとなった.漁業動態が努力量管理に影響を与える過程は対象とする資源や漁業によって大きく異なるだろうから,本研究の結果を単純に一般化することはできない.管理効果を適切に評価するためには,努力量管理の効果に影響を与え得る漁業動態をそれぞれの資源ごとに漁業データと統計モデルを用いて明らかにしていく必要がある.このような試みを様々な資源で行うことにより,努力量管理が漁業動態を通じて漁獲量?資源量に与える過程をどのような形で一般的にモデル化すればよいかが明らかになってくることが期待される.そうすれば,MSE の枠組みの中で努力量管理の効果を定量的に評価することがより一般的になり,今まであまりその効果が評価されてこなかった努力量管理が再度注目されることになるかもしれない.水産資源の管理手法には,今,様々なものが提案されている(Hilborn and Hilborn, 2012).特に近年では,漁獲枠を個人や船団に割り当てる個別割り当て方式(individual quota, IQ)が漁獲にかける無駄な投資を抑制し,利益を重視した経済的な漁獲を行うインセンティブを与える管理方式として注目されている.また,移動の少ない底生性の資源の管理においては,地元の漁業者に漁場の排他的な利用権を与える地域漁業権(territorial user rights to fish, TURF)が有効に働くことが多くの実例とともに示された.海洋保護区(marine protected area, MPA)は,生態系や生物多様性の保全という目的に対して有効な手段であることが多くの研究から科学的に証明されており,多くの国で海洋保護区の設置が義務付けられるようになっている.これらの管理手法は特に近年注目を浴びているものであるが,有効な管理方策というのは管理される資源の目的や特性によって異なる.世界的に有名な水産資源学者アナ?パルマ教授は 2014 年の米国水産学会の招待講演で「資源管理に特効薬はない.常に全ての手札を使って戦え.」というメッセージを発した.このメッセージには,IQ やTURF など,特に有効に働くことが証明された管理方策が近年もてはやされる傾向にあるが,これらの管理手法が全ての資源に通用するわけではないこと,そして,資源の特性に合わせて利用可能な管理方策を全て検討し,最適なものを選択すべきだ,という意味が込められている.本稿で取り上げた休漁を用いた努力量管理はIQ?TURF?MPA などと比べると時代遅れの感があるかもしれない.しかし,比較的広い範囲
一般化状態空間モデルで漁業動態を記述する 73を回遊し,散発的に漁場が形成されるマサバの管理においては,漁獲量?努力量の自己相関構造を利用した休漁による管理が効果的に働いたことが本研究によって示された.他の多くの管理手法においても,流行廃りにこだわらず,用いられている管理方策の有効性を統計モデルやMSE 等の手法を用いて科学的に評価していくことが重要である.それにより管理方策の手札が増え,将来,より多くの水産?野生生物資源が適切に管理されることに繋がるであろう.参 考 文 献Bunnefeld, N., Hoshino, E. and Milner-Gulland, E. J. (2011). Management strategy evaluation: A pow-erful tool for conservation?, Trends in Ecology and Evolution,26, 441–447.Burnham, K. P. and Anderson, D. R. (2002). Model Selection and Multimodel Inference:APracticalInformation-theoretic Approach, 2nd ed., Springer, New York.Fulton, E. A., Smith, A. D. M., Smith, D. C. and van Putten, I. E. (2011). Human behaviour: The keysource of uncertainty in fisheries management, Fish and Fisheries,12, 2–17.Hilborn, R. and Hilborn, U. (2012). Overfishing :What Everyone Needs to Know, Oxford UniversityPress, Oxford, New York.(市野川桃子,岡村寛 訳 (2015).『乱獲—漁業資源の今とこれから』,東海大学出版部,神奈川.)Hilborn, R. and Walters, C. J. (1991). Quantitative Fisheries Stocks Assessment :Choice, Dynamics,and Uncertainty, Chapman and Hall, New York, London.Hillary, R. M., Preece, A. L., Davis, C. R., Kurota, H., Sakai, O. and Itoh, T. (2015). A scientificalternative to moratoria for rebuilding depleted international tuna stocks, Fish and Fisheries(in press).市野川桃子,岡村 寛,黒田啓行,由上龍嗣,田中寛繁,柴田泰宙 (2015a). 管理目標の数値化による最適なABC 算定規則の探索,水産学会誌,81, 206–218.市野川桃子,北門利英,竹内幸夫 (2015b). 統合型資源評価モデル stock synthesis の検討会を開催,水産学会誌,81, 756–761.Ichinokawa, M., Okamura, H., Watanabe, C., Kawabata, A. and Oozeki, Y. (2015). Effective timeclosures: Quantifying the conservation benefits of input control for the Pacific chub mackerelfishery, Ecological Applications,25, .川端 淳,渡邊千夏子,本田 聡,久保田洋 (2013). 平成 24 年度マサバ太平洋系群の資源評価,平成 24 年度我が国周辺水域の漁業資源評価 第一分冊(水産庁増殖推進部?水産総合研究センター 編), 133–166,水産庁増殖推進部?水産総合研究センター,東京.Kellner, J. B., Tetreault, I., Gaines, S. D. and Nisbet, R. M. (2007). Fishing the line near marinereserves in single and multispecies fisheries, Ecological Applications,17, .黒田啓行,境磨,高橋紀夫,伊藤智幸 (2015). TAC を算定する新しいアプローチ:ミナミマグロの管理方式の開発と運用,水産海洋研究,79, 297–307.Makino, M. (2011). Fisheries Management in Japan:Its Institutional Features and Case Studies,Springer, Dordrecht, New York.Maunder, M. N. and Punt, A. E. (2004). Standardizing catch and effort data: A review of recentapproaches, Fisheries Research,70, 141–159.Maunder, M. N. and Punt, A. E. (2013). A review of integrated analysis in fisheries stock assessment,Fisheries Research,142, 61–74.Melnychuk, M. C., Banobi, J. A. and Hilborn, R. (2013). Effects of management tactics on meetingconservation objectives for western north american groundfish fisheries, Plos One,8.Ohshimo, S. and Naya, M. (2014). Management strategy evaluation of fisheries resources in data-poorsituations using an operating model based on a production model, Ja pan Agr ic ul tura l Res earch
74 統計数理 第 64 巻第1号2016Quarterly,48, 237–244.岡村 寛,市野川桃子 (2016). 水産資源学における統計モデリング,統計数理,64(1), 39–57.Pope, J. G. (1972). An investigation into the accuracy of virtual population analysis using cohortanalysis, Research Bulletin of the International Commission for North Atlantic Fisheries,9,65–74.Pope, J. (2009). Input and output controls: the practice of fishing effort and catch management in re-sponsible fisheries, A Fishery Manager’s Guidebook (eds.K.L.Cochrane,S.GarciaandFoodand Agriculture Organization of the United Nations), 220–252, Wiley-Blackwell, Chichester,West Sussex, Ames, Iowa.Punt, A. E. and Donovan, G. P. (2007). Developing management procedures that are robust to uncer-tainty: Lessons from the International Whaling Commission, ICES Journal of Marine Science,64, 603–612.Quinn, T. J. (2003). Ruminations of the development and future of population dynamics models infisheries, Natu ral Re sou rce Mod el ing ,16, 341–392.Quinn, T. J. and Deriso, R. B. (1999). Quantitative Fish Dynamics, Oxford University Press, NewYor k .Salas, S. and Gaertner, D. (2004). The behavioural dynamics of fishers: Management implications, Fishand Fisheries,5, 153–167.Schaub, M. and Abadi, F. (2011). Integrated population models: A novel analysis framework for deeperinsights into population dynamics, Journal of Ornithology,152, 227–237.Swartz, W. and Ishimura, G. (2014). Baseline assessment of total fisheries-related biomass removal fromJapan’s Exclusive Economic Zones: , Fisheries Science,80, 643–651.van der Lee, A., Gillis, D. M., Comeau, P. and Hurley, P. (2013). Fishing the line: Catch and effortdistribution around the seasonal haddock (Melanogrammus aeglefinus) spawning closure on theScotian Shelf, Canadian Journal of Fisheries and Aquatic Sciences,70, 973–981.Worm, B., Hilborn, R., Baum, J. K., Branch, T. A., Collie, J. S. and Costello, C. (2009). Rebuildingglobal fisheries, Science,325, 578–585.
Proceedings of the Institute of Statistical Mathematics Vol. 64, No. 1, 59–75 (2016) 75Modeling of Fishery Dynamics with Autoregressive State-space Modelsfor Quantifying Management Effectivenessin the Pacific Chub Mackerel FisheryMomoko Ichinokawa and Hiroshi OkamuraNational Research Institute of Fisheries Science, Fisheries Research AgencyManagement strategy evaluation (MSE) has increasingly become one of the most im-portant tools for natural resource managements in the applied ecology and fishery sciences.This paper briefly introduces the concept of MSE and then reviews the study in which sta-tistical models describing fishery dynamics work efficiently in MSE. The study quantifiedmanagement strategy of effort control actually applied to the purse seine fishery catch-ing Pacific chub mackerel, by coupling population dynamics simulation and statisticalmodels based on fishery data. The statistical models, generalized autoregressive state-space models, explicitly describe the relationships among fishing effort, daily catches, andbiomass of chub mackerel as stochastic processes. The study revealed two important fac-tors that affect management effectiveness: the autoregressive processes hidden in the dailycatch and effort data and the fisher’s behavioral change in response to increase of stockbiomass. While most MSE applications tend to simplify relationships between fishing ef-fort and catches in a linear manner without consideration of the fisher’s behavior, theactual fishery dynamics considering those factors estimated from appropriate statisticalmodels would advance MSE and improve future wildlife and fishery managements.Key words: Fishery dynamics, management strategy evaluation, chub mackerel, effort control.
CitationsCitations0ReferencesReferences23BookJan 1992Rev Fish Biol FishArticleFull-text availableJun 2015Fish FishArticleJan 2015ECOL APPLShow moreProjectArticleJanuary 2015 · Ecological Applications · Impact Factor: 4.09Restricting human access to a specific wildlife species, community, or ecosystem, i.e., input control, is one of the most popular tools to control human impacts for natural resource management and wildlife conservation. However, quantitative evaluations of input control are generally difficult, because it is unclear how much human impacts can actually be reduced by the control. We present a... [Show full abstract]ArticleJanuary 2016水産資源学は,生態学とは異なる側面を持っている.実験が難しく,主要なデータが漁業からのものであるという点で,不確実性が大きく,バイアスの混入もしばしば見られる.そのような問題に対処するため,古くから統計モデルの活用が積極的に進められてきた.水産資源学は,大きく分けて,資源評価と資源管理からなる.資源評価,資源管理で使われる統計モデルの概要を紹介する.生態学と水産資源学で使われる統計モデルには多くの共通点があり,統計モデルを媒介として,水産資源学と生態学の協調?融合が進むことを期待する. ArticleMay 2017 · NIPPON SUISAN GAKKAISHI · Impact Factor: 0.15ArticleMarch 2015 · NIPPON SUISAN GAKKAISHI · Impact Factor: 0.15This study explored optimum sets of coefficients used in the management procedure (MP) currently applied to the fishery stocks in Japan when estimates of the total stock biomass are not available. The MP determines annual allowable biological catch (ABC) every year using the trend of relative stock abundance and the estimated current stock status categorized into &high&, &middle&, and &low&.... [Show full abstract]ArticleDecember 2016 · ICES Journal of Marine Science · Impact Factor: 2.38We present the first quantitative review of the stock status relative to the stock biomass (B) and the exploitation rate (U) that achieved the maximum sustainable yield (MSY) (B MSY and U MSY , respectively) for 37 Japanese stocks contributing 61% of the total marine capture production in Japan. B MSY and U MSY were estimated by assuming three types of stock-recruitment (S-R) relationships and... [Show full abstract]Data provided are for informational purposes only. Although carefully collected, accuracy cannot be guaranteed. Publisher conditions are provided by RoMEO. Differing provisions from the publisher's actual policy or licence agreement may be applicable.This publication is from a journal that may support self archiving.Last Updated: 26 Sep 17}

我要回帖

更多关于 単位取得退学 的文章

更多推荐

版权声明:文章内容来源于网络,版权归原作者所有,如有侵权请点击这里与我们联系,我们将及时删除。

点击添加站长微信