반응표면법을 이용한 아쿠아포닉스 전처리조 최적 운전 조건 평가 모델

Evaluation Model of Optimal Operating Conditions for Aquaponics Pretreatment Using Response Surface Methodology

Abstract

The aim of this research was to apply a method designed to derive the factors influencing total ammonia removal when operating an additional pretreatment system at Aquaponics. The Box-Behnken method, among response surface analysis methods was used to characterize and determine the optimal nitrification conditions when operating the pretreatment system. Among the mathematically and statistically calculated prediction equations, the total ammonia nitrogen concentration Y1 measured on day 8 was derived as Y₁=-195.8+2.23X₁+42.9X₂+47.5X₃+0.1856X₁²-1.380X₁X₂-1.770X₁X₃, and the time taken to reach the maximum total ammonia nitrogen concentration during the experiment period was derived as Y₂=271-5.04X₁+60.5X₂-64.8X₃+0.1654X₁²+6.54X₃²-0.600X₁X₃-9.00X₂X₃. The coefficients of determination of the regression models of Y₁ and Y₂ were 93.99% and 94.46%, respectively. The modified coefficients of determination were also high, at 89.48% and 88.91%, respectively. The prediction coefficients of determination of Y₁ and Y₂, were 70.68% and 62.11%, respectively, which was relatively lower than that of Y1, but still indicated a reliable prediction performance.

서론

세계인구는 1959년에 약 30억 명을 기점으로 2022년에는 약 80억 명으로 매우 빠른 속도로 증가하고 있다. UN (2022)에 따르면 세계인구는 2030년에는 약 85억명, 2050년에는 약 97억명을 넘어설 것으로 예상하고 있으며, 2100년에는 약 104억명으로 현재 인구수 대비 약 30% 증가할 것으로 예상하고 있다(UN, 2022). 또한 급격한 인구 증가로 인해 지속가능한 단백질 공급원에 대한 논의가 이루어지고 있으며, 현재 수산물은 인류에게 동물성 단백질의 약 15%를 공급하고 있어 식량보급에 있어 중요한 역할로 자리 매김하고 있다. 하지만 어족자원 고갈로 인하여 생산량의 감소 및 어획 정체 상황에 따라 양식산업이 발전하고 있으며, 이로 인해 어획량보다 양식생산량이 더 증가되고 있다(FAO, 2022). 국내 내수면 양식은 지속적인 투자로 과거 15년 동안 많은 양식장에서 순환여과양식시스템을 적용하고 있으며, 내수면 양식 총 생산량의 약 30% 이상을 차지하는 뱀장어 양식(KOSIS, 2022)이 가장 대표적으로 순환여과 양식시스템을 이용하고 있다. 그러나 국내의 경우, 순환여과양식시스템의 메커니즘, 운영기술 부족, 노후화 등으로 인하여 생산성이 낮아지고 있는데, 이는 시스템 내 수질관리 실패(Kim and Khang, 1982)를 주요한 요인으로 들 수 있다. Russo et al. (1981)는 수질 악화가 어체의 스트레스 증가와 이로 인한 다양한 질병 발생의 원인이 되며 각종 항생제와 소독약품의 남용으로 이어져 수산물을 이용하는 소비자들의 건강을 저해하고 식품안전성을 저하시킨다고 하였다. 이러한 수질 악화를 방지하기 위해 간헐적으로 대량의 환수를 실시하면서 환경으로 질소화합물(nitrogen compound)이 배출되어 환경오염을 일으킬 우려도 있다(Heinsbroek and Kamstra, 1990). 국내 내수면 양식의 고질적인 문제를 해결하기 위하여 현재 순환여과양식시스템과 아쿠아포닉스 시스템을 융합한 시스템 개발이 이루어지고 있다(FAO, 2014). 기존 내수면 양식시스템의 구조 변경 없이 양식장의 배출수를 아쿠아포닉스에 이용한 융복합 확장 시스템으로, 양식 현장의 구조개선을 꺼리는 어가들도 적용이 가능하다. 이를 통해 어류 생산시설과 아쿠아포닉스 생산시설이 별도로 운영되기 때문에 보다 안정적이고 추가로 식물을 재배하여 부가 수익을 창출할 수 있는 장점이 있다(Kloas et al., 2015). 그러나, 현재 순환여과시스템에서 배출수는 유수식 양식장의 배출수 보다 배출되는 양은 적지만 높은 농도의 질소 화합물을 배출하고 있다. 적당한 질소 화합물은 식물이 성장하는데 도움이 되지만, 높은 농도의 질소 화합물은 오히려 식물 성장에 방해요인이 된다(Wongkiew et al., 2017). 따라서 본 연구는 뱀장어 양식장의 생물 성장 데이터 및 생체량에 따른 순환여과양식 시스템에서 물리적 여과장치인 드럼 필터에서 발생하는 역세수(backwashing water) 및 침전조 배출수(wastewater)의 질소 부하량을 예측하고, 계산을 통하여 시스템에 전처리조를 추가로 도입하여 고형물을 제거 및 배출수 내 질소 화합물인 총 암모니아의 농도를 낮추고자 하였다. 아쿠아포닉스와 연계한 전처리조 설계 시 성능에 영향을 주는 요인에 대한 자료가 매우 부족하다. 따라서, 본 연구에서는 전처리조 시스템 운영에 미치는 영향 인자를 파악하기 위해 반응표면분석법(response surface methodology, RSM) 중 Box-Behnken법(Ferreira et al., 2007)을 이용하였다. 이 방법을 이용하여 전처리의 성능에 영향을 주는 인자 중 가장 중요한 총 암모니아에 대해 수온, 사료부하량, pH 가 미치는 영향과 상호관계를 평가하였다. 이를 통해 순환 여과양식시스템과 아쿠아포닉스가 결합된 융복합 시스템에 전처리조 시스템의 설계 및 운영의 기초 통계적 모델을 구축하고자 하였다.

재료 및 방법

실험 시료

경남 창원시 소재 국립수산과학원 첨단양식실증센터 뱀장어 아쿠아포닉스 시스템에서 채취한 배출수를 현장에서 직접 채취하여 냉장 상태로 이송한 다음, 500 L 원형 폭기조에서 약 4주간 batch test 및 conditioning 후 본 실험에 사용하였다.

실험 시스템 구성

실험을 위한 전처리조(pretreatment system) 시스템은 폭기용 수조, 온도 유지용 수조, 히터 및 온도조절기, 온도조절이 가능한 냉각기, 에어 브로어로 구성되었다. 폭기용 수조는 5 L를 수용할 수 있는 유리 원통수조이고, 온도 유지용 수조는 폭기용 수조 4개를 수용할 수 있는 약 100 L 유리 사각수조이다. 온도 유지용 수조를 따로 구성하여 온도조절 장치들을 직접적으로 폭기용 수조 안에 설치하지 않고, 온도 유지용 수조에 설치하고 수온을 유지하여 온도 유지 장치들이 폭기조에 영향을 주지 않게 하였다. 온도 유지용 수조에 물을 채운 뒤 1 kW 수중히터(DH-1000; Donghwa electronics, Busan, Korea)를 넣고 수온을 시험조건에 맞춰 조절 및 유지하였다. 수온 조건이 낮아 냉각이 필요한 실험구에는 냉각기(DBC-150; Daeil, Busan, Korea)를 추가로 설치하였다. 폭기는 분지 및 에어스톤을 이용하여 실험 구간별 동일하게 공급하였다.

실험 과정 및 수질 분석

실험은 시작 후 실험조건별 사료부하량을 맞춰준 뒤, 17일 동안에 전처리조의 수질 변화를 관찰하였다. 매일 14:00에 시료 채취 및 pH 조건 보정을 진행하였다. 폭기용 수조에 숙성된 물을 각 5 L씩 채워주고 실험 구간별 조건에 맞춰 수질을 조정해 주었다. 실험조건인 사료부하량에 사용된 사료는 뱀장어 시판 사료(골드일 2호; Purina, Cheongju, Korea)를 사용하였고, 사료 성분은 조단백 52% 이상, 조지방 12% 이상, 조섬유 6% 이하, 조회분 18% 이하, 수분 12% 이하로 조성되어 있다. 배출 고형물의 조건을 주기 위해 믹서기로 분쇄하여 사용하였다. pH는 NaHCO₃ (탄산수소나트륨) 및 HCl (염산)을 이용하여 유지해주었다. 매일 취수한 실험 구간별 시료의 총 암모니아성 질소(total ammonia nitrogen, TAN)는 Salicylate method (Hach method 8155)을 이용하여 분광광도계(DR 900; Hach, Loveland, CO, USA)로 측정하였다. 수온과 pH는 휴대용 수질측정기(HQ2200; Hach)를 이용하여 측정하였다.

분석 방법 및 통계 처리

전처리조의 총암모니아성 질소 제거 특성과 영향 인자들의 상호관계를 평가하고 최적 조건을 찾기 위해 다변량 기법 중 RSM을 이용하였다. RSM은 실험 데이터에 대한 다항식의 적합성에 기초한 수학적·통계적 기법의 집합으로, 반응 또는 관심 반응의 집합이 여러 변수에 의해 영향을 받을 때 이용될 수 있다(Bezerra et al., 2008). 본 실험에서는 RSM 중 대표적인 실험설계법인 Box-Behnken design (BBD)을 사용하였다. BBD에 선정된 각 3개의 독립변수는 전처리조 운영 중 TAN에 영향을 미치는 주요 인자인 수온 20–30°C, 사료부하량 0.5–1.5 mg/L, pH 6–7로 설정하고, X₁, 수온; X₂, 사료부하량; X₃, pH로 배열하여 실험을 설계하였다(Table 1).

Table 1. Box-Behnken design of pretreatment system

반응모델의 통계적 유의성 및 적합성을 검증하기 위해 Minitab (version 21; Minitab Inc., State College, PA, USA)을 이용하여 각 독립변수와 종속변수에 대한 분산분석(analysis of variance, ANOVA)과 회귀분석(regression analysis)을 실시하였다. 또한, 통계적으로 적합하지 않은 반응표면 모델의 적합성을 높이기 위해 P-value를 기준으로 유의하지 않는 항(P>0.1)을 삭제하여 모델 최적화를 진행하였다.

결과 및 고찰

Box-Behnken법을 통해 세 가지 독립 변수(수온, 사료부하량, pH)의 코드화된 조건(X₁, X₂, X₃)에서 실험을 수행한 결과, 전체 실험기간 동안 TAN (mg/L)의 농도는 점차 상승하다 하강하는 경향을 나타내는 것을 확인할 수 있었다. 8일차에 측정된 TAN의 농도(Y₁)와 실험 기간 중 최대 TAN 농도 도달까지 걸린 시간(Y₂)을 도출하였다(Fig. 1).

Fig. 1. (a) is a scatter plot for Y₁, and (b) is a scatter plot for Y₂.

수집된 데이터를 Minitab (version 21; Minitab Inc.)의 반응 표면설계분석을 이용하여 분석을 진행하였다. 회귀분석의 유의수준은 수집한 데이터의 양이 한정적이고 선행연구의 부재로 탐색적 데이터 분석이 필요하였기에 10%로 설정하였다. 회귀분석의 변수선택은 모델의 단순화 및 계산 효율성을 위해 후진제거법(Song and Kim, 2012)을 이용하여 효과적으로 불필요한 변수를 제거하였다. Table 1에서 제시한 변수들의 단위를 고려하지 않은 원래 변수를 적용한 경우와 단위가 다르기 때문에 이들 변수들을 코드화 시켜 통일된 변수로 적용하여 산출한 각 반응 model은 다음 식과 같다.

Parameter estimates from coded factor:

y₁=5.671-6.888x₁+4.188 x₂+1.650 x₃+4.64 x²₁-3.45x₁x₂-4.43x₁x₃

y₂=7.923-3.375 x₁+1.000 x₂-1.875 x₃+4.135 x²₁-1.635 x²₃-1.50 x₁x₃-2.25 x₂x₃

Parameter estimate from un-coded factor:

Y₁ = -195.8+2.23X₁+42.9X₂+47.5X₃+0.1856X²₁-1.380X₁X₂-1.770X₁X₃

Y₂ = 271-5.04X₁+60.5X₂-64.8X₃+0.1654X²₁+6.54X²₃-0.600X₁X₃-9.00X₂X₃

도출된 회귀모형을 통해 Y₁에는 X₁, X₂, X₃, X₁의 제곱, X₁과 X₂의 상호작용. X₁과 X₃의 상호작용이 유의한 영향을 미치며, Y₂에는 X₁, X₂, X₃, X₁의 제곱, X₃의 제곱, X₁과 X₃의 상호작용. X₂와 X₃의 상호작용이 유의한 영향을 미치는 것을 확인할 수 있다. Y₁과 Y₂의 결과는 선형적 관계뿐만 아니라 독립변수 간의 상호작용과 비선형 효과의 영향을 받는 것을 알 수 있다. 회귀모형의 적합성과 변수의 유의성을 검증하기 위해 분산분석을 진행하였다(Table 2). Y₁에 대한 회귀모형에서 모형은 데이터의 변동성 중 93.99%를 설명하고 있다. Y₁ 회귀모델에 대한 X₁의 기여도는 47.69%, F-통계량은 63.45로 Y₁ 결과에 수온 변수가 통계적으로 매우 유의미하게 영향을 미치는 것을 확인할 수 있으며, X₂의 기여도 역시 17.63%로 Y₁의 결과에 영향을 미치는 중요한 요인으로 분석된다. X₃의 기여도는 2.74%로 상대적으로 작지만, P-value가 10%인 신뢰수준보다 낮다는 점에서 X₃는 통계적으로 분명히 유의한 영향을 미친다고 해석된다. 또한, 적합도 부족의 기여도가 5.73%이고 P-value는 0.134로 모델이 데이터의 주요 패턴과 구조를 잘 파악하고 적합시키고 있음을 확인할 수 있다. Y₂에 대한 회귀모형에서 모형은 데이터의 변동성 중 94.46%를 설명하고 있다. Y₂ 회귀모형에 대한 X₁의 기여도는 37.97%, P-value는 0.000이며, X₁의 제곱 항의 기여도 25.11%, P-value 0.001로 수온 변수가 Y₂의 결과에 통계적으로 매우 유의한 영향을 미치는 것을 확인할 수 있다. X₃의 기여도는 11.72%로 X₁ 다음으로 높은 영향을 미치며, X₂와 X₃ 상호작용의 기여도는 8.44%로 X₃ 다음으로 높은 영향을 미친다. X₂와 X₃의 제곱 항, X₁과 X₃의 상호작용의 경우 기여도는 낮지만 각각의 P-value가 0.079, 0.056, 0.066으로 유의수준 10%에서 분명히 유의한 영향을 가지는 것을 알 수 있다. Y₂에 대한 회귀모형의 오차 중 적합도 부족은 5.27%, P-value는 0.121로 적합도 부족이 통계적으로 유의하지 않으며, 순수오차는 0.28%로 측정 오류나 데이터의 자연스러운 변동성이 낮다고 해석된다.

Table 2. Analysis of variance results for the regression models of response variables Y₁ and Y₂

회귀모형의 결과는 Table 3과 같다. Y₁과 Y₂의 회귀모형의 결정계수(R²)는 각각 93.99%, 94.46%로 변수들이 데이터의 변동성을 잘 설명하고 있음을 시사한다. 수정된 결정계수(R²_adj)역시 각각 89.48%, 88.91%로 높은 수준을 보여주고 있다. Y₁ 회귀모형의 표준오차(S)는 2.446, Y₂ 회귀모형은 1.379로 Y₁의 데이터가 Y₂의 데이터보다 모델 예측 값에서 더 큰 편차를 보일 것으로 해석된다. Y₁ 회귀모형의 예측오차제곱합(predicted residual error sum of squares, PRESS)은 233.29로 모형이 새로운 데이터에 대해 높은 오차를 가질 수 있음을 시사한다. 그러나, 예측 결정계수(R²_pred)는 70.68%로 상당한 수준의 예측 성능을 보일 것으로 기대된다. Y₂ 회귀모형의 PRESS은 90.941로 새로운 데이터에 대해 낮은 오차를 보일 것으로 기대되지만, 예측 결정계수(R²_pred)가 62.11%이라는 점에서 모델이 새로운 데이터의 변동성을 얼마나 잘 설명할 수 있을지에 대한 추가적인 검증이 필요하다.

Table 3. Regression model performance metrics for response variables Y₁ and Y₂

PRESS, Predicted residual error sum of squares.

모델이 데이터를 적절하게 예측하는지 확인하기 위해 uncoded 변수로 도출한 회귀모형의 식을 이용하여 예측 값과 잔차를 계산하였다(Table 4). Y₁의 실제 값과 예측 값 사이 잔차는 대부분의 실행에서는 ±2 이내의 오차를 보여주지만, 8번의 실행에서는 -4.3의 높은 잔차를 보였다. Y₂의 실제 값과 예측 값 사이 잔차는 모두 ±2 이내의 값을 가진다는 점에서 모형이 데이터를 적절하게 예측하고 있음을 시사한다. 잔차의 분포와 패턴을 시각화하기 위해 그래프로 나타내었다(Fig. 2).

Fig. 2. (a) is a residual analysis plot for Y₁, and (b) is a residual analysis plot for Y₂.

Table 4. Comparison of the observed and predicted Y₁ and Y₂

Fig. 2a와 Fig. 2b는 잔차를 시각적으로 평가하기 위해 Y₁과 Y₂에 대한 잔차 분석 결과를 정규확률도(normal probability plot), 잔차 대 적합값(residuals vs fitted), 히스토그램(histogram), 잔차 대 관측순서(residuals vs order)로 나타내었다. Fig. 2a의 정규확률도에서 데이터 대부분이 참조 선을 따라 분포하고 있으며 히스토그램이 정규분포의 형태를 보인다는 점에서 회귀모델의 적합성이 상당히 좋다고 평가할 수 있다. 잔차 대적합값과 잔차 대 관측순서 그래프에서 명확한 추세 혹은 패턴이 나타나지 않으므로 잔차가 독립적이라는 가정을 지지할 수 있다. Fig. 2b의 정규확률도에서는 모든 데이터가 참조선을 벗어나지 않고 주변에 밀집한 것을 확인할 수 있으며, 히스토그램 역시 중심점을 기준으로 좌우 대칭의 형태를 가지고 있다. 잔차대 적합값, 잔차 대 관측순서 그래프에서 명확한 패턴이나 추세는 나타나지 않고 있지만 적합값이 증가함에 따라 잔차의 분산이 약간 증가하였다.

2차원 등고선 그래프를 이용하여 Y₁과 Y₂에 미치는 X₁, X₂, X₃의 영향을 제시하였다(Fig. 3). Fig. 3a는 온도와 사료부하량이 Y₁에 미치는 영향을 나타내는데 수온이 낮아지고 사료부하량이 증가할수록 Y₁의 값은 높아지는 경향을 보인다. 특히 사료 부하량이 1.0 이상인 경우 더욱 명확한 경향을 보이며, 수온이 26°C 이상인 경우에는 Y₁이 감소하는 경향이 뚜렷하게 나타난다. Fig. 3b에서는 온도와 pH가 Y₁에 미치는 영향을 제시하며 수온이 높아지고 pH가 낮아질수록 Y₁은 낮아지는 추세가 나타났다. 등고선이 직선에 가까운 형태로 나타난다는 점에서 pH와 수온의 상호작용 효과가 크게 나타나지 않고 두 변수의 영향이 독립적임을 시사한다. Fig. 3c에서 사료부하량이 증가함에 따라 Y₁도 같이 증가하는 경향을 보이는데, 특히 사료부하량이 1.0 이상인 구간에서 더욱 확실한 경향이 나타난다. 등고선 간 거리가 먼 부분에서는 pH 혹은 사료부하량의 변화가 Y₁의 값을 크게 변화시키지 않을 것을 시사한다. Fig. 3d는 온도와 사료부하량이 Y₂에 미치는 영향으로 낮은 온도와 높은 사료부하량에서 Y₂의 값이 높게 나타나는데, 특히 사료부하량이 1.0 이상, 수온이 24–26°C 사이일 때 Y₂가 증가하는 경향이 뚜렷하게 보인다. 그래프의 중앙 부분은 등고선이 넓게 퍼진 형태를 가지는데, 이는 두 변수의 상호작용이 Y₂에 미치는 영향이 덜 민감하며 상대적으로 안정적인 범위임을 시사한다. Fig. 3e는 온도와 pH가 Y₂에 미치는 영향을 시각화한 그래프이다. 왼쪽 하단 영역은 등고선이 좁은 간격으로 밀집해 있다는 점에서, pH와 온도가 낮아질 경우 Y₂의 값이 빠르게 감소하는 경향을 가진다. Fig. 3f는 사료부하량과 pH가 Y₂에 미치는 영향을 나타내는데, 등고선이 가로 방향으로 넓게 퍼져 있어 사료부하량이 Y₂에 미치는 영향은 상대적으로 일정하게 유지되는 경향을 가지는 것을 알 수 있다. 다만, 사료부하량이 1.0 이상이고 pH가 낮은 상황에서는 Y₂값이 상대적으로 높게 나타나는 것을 확인할 수 있다.

Fig. 3. Series of contour plots (a)-(f) displaying the effects of X₁, X₂, X₃ on two responses, Y₁, and Y₂. Plots (a)-(c) detail these effects on Y₁, with (a) focusing on temperature and feeding load, (b) on Temperature and pH, and (c) on pH and feeding Load. Plots (d)-(f) correspond to the same interactions for Y₂.

본 연구는 아쿠아포닉스에서 전처리조(pretreatment)를 추가로 결합한 융복합 시스템에서 전처리조를 운영할 때 총 암모니아 제거에 영향을 미치는 인자를 도출하기 위해 실험을 수행하였다. 이 때 총 암모니아는 질산화 반응에 영향을 미치는 3개 요인인 수온, pH, 사료부하량에 따라 달라지므로 이들 특성을 해석하고 전처리조 운영 시 최적의 질산화 조건을 찾기 위해 RSM 중 박스-벤켄법(Box-Behnken method)을 이용하였다. 실험을 통한 최적화 도출 시 8일 차에 측정된 TAN의 농도(Y₁)에서는 주요 변수가 6개 정도 포함되어 있으며, 실험 기간 중 최대 TAN 농도 도달까지 걸린 시간(Y₂)에서는 변수가 7개 정도 포함되어 있다. 이들 변수들을 객관적으로 설계하여 변수간 관계 및 변수와 반응 간의 관계를 찾아 해석하고 최적의 조건들을 찾기란 어려운 작업이다. 본 연구는 변동성이 포함된 실제 아쿠아포닉스 시스템 배출수를 대상으로 수행하였으며, 초기 실험조건에서 일정한 조건을 유지하기가 매우 어려웠다. 전처리조에 의한 총암모니아 제거율에 대하여 반응표면분석을 적용하여 얻은 모델 예측식을 도출한 결과, model Y₁은 70.68%, Y₂는 62.11%로 Y₁은 상당한 수준의 예측 성능을 보일 것으로 보이지만, Y₂의 경우 그 보다 낮은 수준의 예측 성능을 보일 것으로 판단된다. 예측 값과 잔차를 계산하였을 때, Y₁의 실제 값과 예측 값 사이 잔차는 대부분의 실행에서는 ±2 이내의 오차를 보였고, Y₂의 실제 값과 예측 값 사이 잔차는 모두 ±2 이내의 값을 가진다는 점에서 모형이 데이터를 적절하게 예측하고 있음을 시사한다. 그러나, 이 모델이 실제 현장에서 다른 새로운 데이터의 변동성을 얼마나 잘 예측할 수 있을지에 대한 추가적인 실험 또는 검증 과정이 필요할 것으로 보인다.

이런 연구를 통해 향후 아쿠아포닉스에 새로운 수질 관련 설비 및 관련 공정을 설계하고자 할 때 과학적이고 체계적인 실험 계획을 통하여 커다란 설비투자 없이 공정의 핵심 인자를 규명할 수 있을 것이다. 또한, 다양한 인자들에게 영향을 미치는 관계를 모형과 함께 원인 분석하거나 최적화를 도출하여 실제 현장에 적용 시 발생하는 다양한 변수들을 더 정확하게 파악하고 제어한다면 위험 부담을 줄이고 공정의 효율도 극대화할 것으로 판단된다.