볼테라 급수는 일반적으로 비선형 동적 시스템의 응답을 식별하기 위해 널리 사용된다. 볼테라 급수는 응답의 차수에 기반한 볼테라 커널을 가지고 있으며, n차 볼테라 커널은 n차 충격응답함수(impulse response function)와 같으며 n차 볼테라 커널과 n차 입력의 컨볼루션으로 n차 응답이 구해진다. 시스템의 응답은 모든 차수의 응답의 합으로 이루어지며 단일 입력 u(t)하에 단일 출력 y(t)는 볼테라 급수로 식 (1)과 같이 나타낼 수 있다.

h_n은 n차 볼테라 커널이다. 볼테라 급수의 차수는 N에 의해 결정되며 시스템이 선형일 때, 즉 N=1일 때 식 (1)은 식 (2)와 같이 1차 볼테라 급수로 표현된다.

만약 시스템이 유한한 기억 효과(Memory effect)를 갖는 인과(causal) 시스템이면 연속적인 볼테라 커널을 볼테라 계수로 이산화하여 식 (3)과 같이 나타낼 수 있다.

M₁는 시스템의 기억 길이(Memory length)를 의미하며 Δτ₁는 샘플링 시간 간격을 의미한다. 본 연구에서는 시스템이 선형이라는 가정하에, 계측을 통해 얻어진 데이터를 이용하여 식 (3)에 주어진 1차 볼테라 커널 h₁(τ₁)을 확률적으로 식별하고자 하였다. 볼테라 커널의 효과적인 식별을 위해 먼저 커널을 적절한 기저 함수로 급수 전개하여 계산에 소요되는 과다한 시간 및 과적합 등의 문제를 해결하고자 하였다.

볼테라 커널은 다양한 기저 함수를 이용하여 급수 전개할 수 있는데, 본 연구에서는 자유감쇠신호의 추정에 적합한 라게르 다항식을 적용하였다 (Schetzen, 1980). 라게르 다항식은 직교하는 성질을 가지고 있으며 충격응답함수와 같이 시간에 따라 감쇠하는 함수를 추정하기에 적합하다. 라게르 다항식은 식 (4)와 같으며 라게르 다항식의 직교관계는 식 (5)와 같다 (Israelsen & Smith, 2014).

라게르 다항식의 모수 a는 시간축에 대한 스케일 팩터(Scale factor)이며 1/a 값에 비례하여 t축 방향으로 확대 변환된다. 추정하는 시스템의 특성에 따른 적당한 a값을 선정하는 것은 추정도에 큰 영향을 미치는 요소이다. Fig. 1은 a = 0.1일 때 라게르 다항식을 시간축에 따라 도시한 것으로 시간에 따라 감쇠하는 임의의 자유감쇠신호를 추정하기에 적합한 형상을 가지고 있다. 볼테라 급수를 라게르 다항식으로 급수 전개했기 때문에 볼테라 계수는 라게르 다항식의 계수를 구함으로서 얻을 수 있으며, 본 연구에서는 시스템이 가지는 불확실성을 고려하기 위해 라게르 다항식의 계수를 랜덤 변수로 취급하고, 그들의 확률분포를 계산하고자 하였다.

Fig. 1 
The first five Laguerre functions

베이지안 추론은 베이즈 정리에 기반한 추론기법으로, 시스템을 특정 짓는 변수들의 확률 분포를 추정하는 방법 중 하나이다. 베이지안 추론에 따르면, 식 (6)에 표현된 바와 같이, 시스템의 확률변수 β에 대한 사후 확률분포(posterior probability distribution)는 확률변수 β의 사전 확률분포(prior probability distribution)와 계측을 통해 얻어진 데이터 D의 우도함수(likelihood function)의 곱으로 표현이 가능하다.

여기서, p(β)는 확률변수 β에 대한 사전 확률분포이며, p(D|β)는 우도 함수로서 확률 변수 β에 대해 데이터가 발현될 가능도를 의미한다. p(D)는 사후 확률 분포 p(β|D)의 전체 확률이 1이 되도록 하는 정규화 상수이다.

베이지안 선형 회귀모델은 베이지안 추론에 기반하여 선형회귀 방정식의 계수들을 확률적으로 구하는 모델이다. 베이지안 선형 회귀모델에서는 입력 값 행렬이 X, 출력 값 벡터를 Y, 라게르 다항식의 계수를 β, 출력의 에러를 ε이라 하고, 에러 ε은 평균이 0이고 분산이 σ² 인 정규 분포를 따른다고 가정한다. 식 (7)은 선형회귀모델의 방정식을 행렬로 표현한 것이다 (Stevens et al., 2020).

일반적으로 베이지안 선형 회귀모델에서는 β와 σ²의 사전 확률분포를 정규 역감마(Multivariate normal inverse gamma)분포로 가정하여 사후 확률분포를 추정한다. 정규 역감마분포는 확률 변수 β의 다변수 정규 분포와 확률 변수 σ²의 역 감마 분포의 결합 확률 분포이다. 또한 정규 역감마분포는 켤레 사전 확률 분포로서 사전 확률분포와 사후 확률분포가 같은 종류의 확률분포를 갖는 특성이 있다. 일반적으로, 사후 확률분포는 해석적으로 구하거나 계산하기에 어려움이 있지만 켤레 사전 확률분포의 경우 해석적으로 구할 수 있어 매우 유용하다.

사전 확률분포의 평균을 μ_β, 분산을 σ² V_β, 형태 모수(shape parameter)를 A, 척도 모수(scale parameter)를 B라 하면 사전 확률분포는 식 (8)과 같다.

켤레 사전 확률분포를 가정했으므로 사후 확률분포 또한 정규 역감마분포로 식 (9)와 같이 계산된다. 사후 확률분포의 초모수(hyper-parameter)인 μ_β*, V_β*, A*, B*는 식 (10)과 같이 구할 수 있다.

사후 확률분포는 사전 확률분포가 측정된 데이터로 인해 업데이트되어 생성되는 확률분포이기 때문에 사전 확률분포가 균일 확률분포(uniform probability distribution)이거나 측정 데이터가 많을수록 사후 확률분포의 초모수는 측정 데이터에 지배적으로 결정된다.

사후 확률분포가 주어져 있을 때 새로 측정 및 관측된 상황을 고려하면 사후 확률분포를 이용하여 시스템의 응답 시계열을 확률론적으로 추정할 수 있으며 이를 사후 예측분포(Posterior predictive distribution)라 한다. 새로운 입력 행렬 X~가 들어왔을 때 사후 예측분포는 식 (11)과 같이 계산되며 다변수 t 분포(Multivariate t distribution)를 따른다.

실험 데이터는 400,000 DWT VLOC의 1:60 척도의 모델 기반 실험 수행을 통해 측정되었다(Kim & Park, 2015). JONSWAP 스펙트럼에 기반한 불규칙파랑과 모형선의 수직 굽힘 모멘트가 측정되었으며 실험에 사용된 모형선은 Fig. 4에 주어진 것과 같다.

Fig. 4 
Model ship

모형선은 경하 상태(Ballast condition)로 해상 상태에 따라 주어진 파랑 조건하에서 선수 방향으로 파랑을 받으며 운행되었다. 유의 파고, 첨두 주기에 따른 여러 조합의 해상 상태 조건하에 Fig. 4(b)의 Cut 3 위치에서 수직 굽힘 모멘트가 측정되었다. 실험에 사용된 불규칙파는 식 (14)와 같은 JONSWAP 스펙트럼에 기반하여 생성되었다.

측정된 수직 굽힘 모멘트 시계열은 파랑 주파수 응답 부분(wave frequency part)와 진동 응답 부분(vibration part)을 포함하고 있으며 파랑 주파수 응답 부분만 본 연구에 사용되었다. 파랑 주파수 부분은 저 주파수 통과 필터(low frequency pass filter)를 통해 2H_z 이상의 신호는 필터링 되었으며 모든 데이터는 실선 스케일의 데이터 스케일을 갖는다.

모형선의 선속은 해상 상태에 따라 다르게 운행되었으며 이로 인해 시스템이 변화할 수 있기 때문에 경하 상태 하에 비슷한 선속을 가진 데이터를 같은 그룹으로 묶어 시스템을 식별하였다. 6가지 해상 상태는 JONSWAP 스펙트럼의 모수 γ, 유의 파고 H_s, 첨두 주기 T_p에 따라 Table 1과 같이 분류하였다.

파랑의 시계열 데이터를 입력, 모형선의 수직 굽힘 모멘트 시계열 데이터를 출력 시계열로 이용하여 사후 확률분포를 추정하였다. 라게르 다항식의 계수와 응답 시계열의 오차에 대한 사전 정보, 즉 사전 확률분포가 없기 때문에 사전 확률분포가 균등확률분포(uniform distribution)에 가깝게 초모수(hyper-parameter) 값을 설정하였다. 이를 위해 사전 확률 분포의 라게르 계수의 평균을 μβ=0, 분산을 V_β=10¹¹으로 설정하였으며, 식 (10)을 통해 알 수 있듯이 IG분포의 초모수가 사전 확률분포의 가중치에 비해 측정된 데이터에 의해 지배적으로 결정되도록 A=0, B=0 으로 설정하였다. 하나의 해상 상태의 시계열 데이터는 약 40 min의 시간동안 측정되었으며 Table 1의 각 그룹마다 사후 확률분포 추정 및 검증을 수행하였다. 사전 확률분포를 업데이트하기 위해 각 그룹에서 T_p=14.8 sec, 즉 No.2, No.5에 해당하는 해상 상태의 파랑 및 수직 굽힘 모멘트 시계열이 사용되었다. Group 1의 경우 No.1, No.3, Group 2의 경우 No.4, No.6에 해당하는 해상 상태의 파랑 및 수직 굽힘 모멘트 시계열은 추정된 사후 확률분포의 검증에 사용하였다.

Fig. 5는 No.2, No.5의 해상 상태의 데이터로 추정한 각각의 사후 확률 분포의 평균으로부터 계산한 1차 볼테라 커널이다. No.2와 No.5 해상 상태에서의 실험의 선속이 3 kts 만큼의 차이가 있기 때문에 두 그룹 간의 1차 볼테라 커널이 시간축에서 이동되어 나타남을 확인할 수 있다.

Fig. 5 
First order Volterra kernels calculated from the mean of the posterior probability distribution

Fig. 6은 Group1, Group2 각각의 사후 확률 분포의 평균으로 계산한 1차 볼테라 커널을 푸리에 변환(Fourier transfrom)하여 구한 RAO(Response amplitude operator)를 나타내며 일반적인 수직 굽힘 모멘트 RAO의 형상을 보여주고 있다. Group 1, Group 2에서의 수직 굽힘 모멘트의 1차 고유진동수는 각각 실험 선속에 따라 조우주파수 값으로 5.4, 5.2 rad/sec이며 Fig. 6에 보여지듯이 ω = 0.5 rad/sec 부근에서 가장 큰 응답 보여준다.

Fig. 6 
Response amplitude operator

Fig. 7은 모든 해상 상태에 대해 사후 확률분포로부터 추정한 수직 굽힘 모멘트 시계열을 나타내고 있다. Fig. 7(a), (b),(c)는 Group 1에 대한 시계열 데이터로서, (b)는 학습 데이터로 사용된 응답 시계열과 사후 확률분포로 추정한 시계열을, (a),(c)는 검증 데이터로 사용된 No.1, No.3의 해상 상태 시계열과 사후 확률분포로 추정한 시계열을 비교하여 나타내고 있다. 이를 통해 학습에 사용된 데이터에 대한 학습 정도 및 학습에 사용되지 않은 해상 상태의 시계열 데이터에 대한 추정도를 알 수 있다. Fig. 7(d),(e),(f)는 Group 2에 대한 시계열 데이터로서, Group 1과 마찬가지로 (e)는 학습 데이터로 사용된 시계열과 사후 확률분포로 추정한 시계열을, (d),(f)는 검증 데이터로 사용된 시계열과 사후 확률분포로 추정한 시계열을 비교하여 보여주며 Fig. 7의 추정 시계열은 사후 확률분포의 평균으로부터 추정한 시계열, 사후 예측분포의 ±1 표준편차, ±2 표준편차에 해당하는 영역을 나타낸다. Table 2는 모든 해상 상태에 대해 측정 데이터와 사후 확률분포의 평균으로부터 추정한 시계열 간의 결정계수 값을 보여준다.

Fig. 7 
Time series comparison between the measured data and the prediction of posterior