Research Paper

Journal of the Computational Structural Engineering Institute of Korea. 31 August 2023. 233-242
https://doi.org/10.7734/COSEIK.2023.36.4.233

ABSTRACT


MAIN

  • 1. 서 론

  • 2. Bouc-Wen Hysteresis 모델

  • 3. 분산점 칼만필터(Unscented Kalman Filter)

  •   3.1 분산점 칼만필터 알고리즘

  •   3.2 다중 분산점을 이용한 기법 제안

  • 4. 수치 예제 적용

  •   4.1 분산점 칼만필터 적용 및 결과

  •   4.2 다중 분산점 칼만필터 적용 및 결과

  • 5. 결 론

1. 서 론

칼만 필터(Kalman Filter, KF)는 측정 가능한 데이터를 이용하여 대상 시스템의 상태(state)를 편향되지 않고(unbiased), 최소의 분산을 가진(Minimum variance), 최적의 추정값을 도출하기 위한 방법론으로 동적 시스템의 상태 추정 및 예측, 적응 제어 등 다양한 공학적 분야에서 사용되고 있다(Kalman, 1960; Lee and Song, 2019; Nguyen et al., 2022). 칼만 필터 적용을 위해서는 대상 동적 시스템을 State-Space 모델로 표현해야 하며, 이때 대상 시스템이 비선형 거동을 보일 경우, 확장 칼만 필터(Extended KF, EKF) 혹은 분산점 칼만필터(Unscented KF, UKF)가 널리 사용된다. EKF는 테일러 근사(Taylor expansion)를 이용한 1차 선형 근사 기법이고, UKF는 시그마포인트(Sigma- points)를 이용한 비선형 근사 기법인 Unscented transformation을 통해 2차 선형 근사가 가능하다. 추정하고자 하는 상태 벡터 크기에 따라 UKF가 다소 많은 계산량을 요구하지만, 미분 값을 구하기 쉽지 않거나 벡터 사이즈가 매우 크지 않은 일반적인 경우 정확한 비선형 거동 모사를 위해 EKF보다 UKF를 사용하는 것이 효과적이다.

최근 UKF를 이용하여 시간이력에 따라 구조 특성이 변하는(Time-variant) 동적 시스템의 상태를 실시간으로 추적하는 연구가 수행되어 왔다(Lee and Song, 2019; Nguyen et al., 2022; Papakonstantinou et al., 2020). 이러한 연구들은 가용한 데이터, 예를 들어 측정된 가속도 이력을 기반으로 대상 시스템의 구조적 정보들(강성, 감쇠, 모델 특성 등)을 역추정하는 System Identification(SI)을 수행하기 위한 목적으로 UKF를 사용하였다. 칼만 필터의 추정 성능을 최적화하기 위해서는 기본적으로 자체 파라미터인 두 가지 인공오차 Artificial-process- noise와 Artificial-measurement-noise를 적절히 설정하는 튜닝(Tuning) 과정을 거친다. 일반적으로 두 가지 중 한 값은 주어진 문제에 따라 가정되고, 다른 나머지 값을 조정하여 최적의 추정값을 도출한다. 앞서 언급된 연구들에서는 최적의 상태 추정값을 얻기 위해 시행착오(Trial & Error) 방법 혹은 전역 최적화를 이용하여 인공오차를 제시하였다. 그런데 UKF는 추가적으로 시그마포인트의 분포를 결정하는 스케일링 파라미터(Scaling parameter)인 α를 선정해야 하는 어려움이 있다.

두 가지 인공오차와 스케일링 파라미터가 적절히 설정될 경우 UKF는 대상 시스템의 상태를 정확하게 추정할 수 있지만, 실제 상황에서 세 가지 파라미터를 적절하게 결정하는 것은 매우 어렵다. Lee와 Song(2019)은 제안된 칼만 필터 튜닝을 위해 전역 최적화 기법인 PSO(Particle Swarm Optimization)을 사용하였고, 국소해(Local minimum)를 피하기 위해 목적함수의 입력변수 경계조건(Boundary Condition)을 변경해가며 PSO를 여러 번 수행하는 방안을 제안하였다. 하지만, 해당 방안은 전역최적화 특성상 매우 많은 계산량을 요구하며, 추가적으로 스케일링 파라미터를 선정하기 위해서는 기하급수적으로 계산량이 증가할 수밖에 없다. 이러한 칼만 필터 튜닝 문제를 완화하기 위해서는 설정해야 하는 파라미터의 수를 줄여야 하는데, 칼만 필터의 종류와 상관없이 State-Space 모델에서 필요로 하는 두 가지 인공오차는 필수적으로 결정되어야 하는 파라미터들이다. 따라서, 시그마포인트 스케일링 파라미터 선정을 튜닝 과정에서 제외하는 것이 가장 현실적이며, 이를 위해 본 연구에서는 다양한 스케일링 파라미터에 해당하는 시그마포인트 세트를 동시에 사용하는 UKF를 제안하였고 이를 통해 튜닝 문제를 완화하였다.

2장에서는 비선형 동적 시스템을 모사하기 위해 선정된 Bouc- Wen 모델에 대해 간략히 설명하고, 칼만 필터 적용을 위한 State-Space 모델을 정의한다. 3장에서는 앞서 정의된 비선형성을 갖는 State-Space 모델을 이용하여 상태 추정을 하기 위한 UKF 알고리즘을 설명하고, 본 연구에서 제안한 다중 시그마포인트 세트를 사용하는 방법을 제안한다. 4장에서는 수치 예제인 급격한 구조손상 시나리오를 가정하고 다양한 노이즈 레벨에 대해 기존 UKF 및 제안된 UKF를 수행한다. 최종적으로 수치 예제 결과를 통해, 제안된 방안의 정확성(Accuracy)과 강인성(Robust)을 검증한다.

2. Bouc-Wen Hysteresis 모델

구조공학분야에서 Bouc-Wen 모델은 시간이력해석 및 Random Vibration 해석분야에서 가장 많이 사용되는 이력현상(Hysteresis) 모델 중 하나로, 1967년 Bouc에 의해 제안되었고 1976년 Wen에 의해 일반화되었다. 이 모델은 효율적이고 부드러운 이력현상 곡선을 표현할 수 있는 것으로 알려져 있다(Song and Der Kiureghian, 2006; Vaiana et al., 2018). Bouc-Wen 모델은 단 하나의 보조 비선형 방정식을 이용하여 이력현상 거동을 효율적으로 계산가능한 장점이 있으며, Strength Degradation, Stiffness Degradation, Strain-Hardening, Pinching-Effect 등 다양한 거동 현상과 비대칭 이력현상도 표현 가능하다(Baber and Noori, 1984; Baber and Wen, 1981; Foliente et al., 1996; Noori et al., 1986; Wang and Wen, 1998).

본 연구에서는 단자유도 시스템(SDOF Oscillator, Fig. 1)의 거동 모사를 위해 Bouc-Wen모델을 적용하였으며, 이때 지진 하중을 받는 시스템의 평형 방정식은 다음과 같다(Fujimura and Der Kiureghian, 2007; Yi et al., 2017).

(1)
mu¨(t)+cu˙(t)+k0[αBWu(t)+(1-α)z(t)]=F(t)

여기서, u¨(t),u˙(t),u(t)는 각각 가속도, 속도, 변위를 의미하며, m,c,k0는 질량, 감쇠, 강성을 의미하며 각각 227kg, 25.9652 N・s/m, 29700N/m을 사용하였다. Rayleigh 감쇠를 사용하였으며, 이때 고유주기(Natural Period) T = 0.48초이다. αBW는 항복 후 강성비를 나타내며 본 연구에서는 0.1을 사용하였다. F(t)=-mu¨g(t)는 지진 하중을 의미한다. 본 연구에서는 PEER (Pacific Earthquake Engineering Research Center)에서 제공하는 지진파 중, 1971년 San Fernando에서 발생한 규모 6.61 지진을 선정하여 사용하였다. 마지막으로, z(t)는 다음 Bouc-Wen 모델을 따르는 보조 변수이다(Fujimura and Der Kiureghian, 2007; Yi et al., 2017).

(2)
z˙(t)=u˙(t)A-[γsign(z(t)u˙(t))+η]z(t)n

위 식에서 z(t)u(t)를 제외한 나머지 파라미터들은 이력현상 곡선의 형상을 조정하는 값들이다. 본 연구에서는 A = 10, n = 4로 가정하였고, γ=η=1/(2uyn)이며, 항복 변위 uyn는 0.0202 값을 사용하였다. Fig. 2는 위와 같이 설정된 단자유도 시스템에 선정한 San Fernando 지진이력을 재하했을 때, 이력현상 곡선을 도출한 그래프이다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2023-036-04/N0040360403/images/Figure_jcoseik_36_04_03_F1.jpg
Fig. 1.

SDOF oscillator considering bouc-wen model

https://cdn.apub.kr/journalsite/sites/jcoseik/2023-036-04/N0040360403/images/Figure_jcoseik_36_04_03_F2.jpg
Fig. 2.

Hysteretic loop: SDOF oscillator considering bouc-wen model

다음장에서는 Bouc-Wen 모델로 모사된 비선형 거동을 보이는 단자유도 시스템의 상태 추정을 하기 위해 UKF를 적용한다. 이를 위해 State-Space 모델을 정의하고 UKF 알고리즘을 소개한다.

3. 분산점 칼만필터(Unscented Kalman Filter)

앞서 2장에서 소개된 비선형 동적 시스템의 상태를 적절히 추정하기 위해서는 비선형성을 고려할 수 있는 필터가 사용되어야 하며, 본 연구에서는 UKF를 적용하였다. UKF를 적용하기 위해서는 주어진 문제에 적절한 State-Space 모델을 구축해야 하며, 해당 모델은 다시 System-process 모델과 Measurement- output 모델로 구성된다(Chatzi and Smyth, 2009; Lee and Song, 2019; van der Merwe and Wan, 2001). 이는 일반적으로 다음과 같이 시간 축으로 이산화된 식으로 표현된다.

(3)
xt+1=f(xt)+wt,wt~N(0,Qt)
(4)
zt+1=h(xt+1,ut+1)+vt+1,vt+1~N(0,Rt+1)

위 식에서 xtut+1는 각각 강성 및 감쇠 계수와 같이 추정하고자 하는 변수로 이루어진 상태 벡터와 Input 벡터이며, QtRt는 각각 시스템 인공오차(Artificial-process-noise)와 측정 인공오차(Artificial-measurement-noise)의 공분산을 의미한다. wtvt+1은 평균이 0인 벡터이고 각 정의된 인공오차 공분산으로부터 생성된 가우시안(Gaussian) 랜덤 변수이다. 최적의 성능을 보이는 UKF를 위해서는 적절한 값의 QtRt를 설정해야 하며, 본 연구에서는 시행착오 방법을 이용하여 도출하였다. 마지막으로 함수 f(·)h(·)는 각각 비선형 예측 모델과 보정 모델을 대표하는 함수를 의미한다. 함수 f(·)의 경우 강성 및 감쇠 계수와 같은 시스템 파라미터들의 시간에 따른 관계는 정해진 수식이 없으므로 xt+1=xt+wt와 같이 정의가 가능하다. 반면 함수 h(·)는 업데이트된 시스템 파라미터 xt+1를 기반으로 동적 시스템의 거동을 계산하는 함수이며, 본 연구에서는 수식 (1)(2)를 기반으로 뉴마크(Newmark) 동해석 방법론을 사용하여 가속도를 계산하는 함수로 정의한다. 다음으로는 구축된 State-space 모델을 이용하기 위한 본래의 UKF 알고리즘을 소개한다.

3.1 분산점 칼만필터 알고리즘

칼만필터는 주어진 동적 시스템의 거동을 모사하는 State- space 모델을 이용하여 시간-업데이트(Time-update)와 측정-업데이트(Measurement-update)가 순차적, 반복적으로 수행되면서 시스템의 상태(state) 이력을 추정하게 된다(Kalman, 1960). 이때 시간-업데이트를 수행하기 위해 State-space 모델의 System- process 모델이 사용되며, 확률론적 관점에서 예측(Prediction)이라고 불리기도 한다. 측정-업데이트의 경우 Measurement- output 모델을 이용하여 수행되며 예측 과정에서 갱신된 시스템의 상태를 측정 정보를 기반으로 수정하는 과정이기 때문에 보정(Correction) 과정으로 불린다. 즉, 칼만 필터는 예측과 보정 과정이 순차적으로 반복되는 연속적 베이지안 업데이팅(Sequential Bayesian Updating)으로 이해할 수 있다.

UKF는 비선형성을 고려하기 위해 Unscented Transformation 기법을 도입하였고(Fig. 3), 이는 시그마포인트(Sigma-points)를 이용한 State-Space 모델의 비선형성을 2차근사화하는 것으로 알려져 있다(van der Merwe and Wan, 2001). Fig. 4는 UKF의 일련의 과정을 도식화하였으며, x^tt,Pttxx는 각각 추정하고자 하는 상태 벡터의 시간 t에서의 평균 벡터와 공분산을 의미한다. UKF는 시그마포인트를 계산하여 예측과 보정 과정을 수행하며, 다음은 예측 과정의 시그마포인트 계산 과정이다(Lee and Song, 2019).

(5)
Χtt[0]=X^ttΧtt[i]=X^tt+((nX+λ)PttXX)i,i=1,...,nXΧtt[i]=X^tt+((nX+λ)PttXX)i-nX,i=nX+1,...,2nX

위 식에서 Xtt[i]는 시그마포인트를 의미하며 총 2nX+1개가 생성된다. 𝜆는 α2(nX+κ)-nX로 정의되며, 𝛼와 𝜅는 스케일링 파라미터이다. 𝜅는 뒤틀림(Skewness), 첨도(Kurtosis)와 같은 통계적 고차 모멘트(High-moment)를 고려하기 위한 파라미터로 일반적으로 0으로 고정되며, 보통 𝛼값에 따라 시그마포인트 분포가 결정된다. 행렬 제곱근은 Cholesky Factorization을 통해 계산되며(van der Merwe and Wan, 2001), (·)i는 괄호 안 행렬의 i번째 열 벡터를 의미한다. 계산된 시그마포인트와 System-process 모델인 식 (3)을 이용하여 다음과 같이 시간-업데이트가 수행된다.

(6)
Xt+1t[i]=f(Xtt[i]),i=0,...,2nXXt+1t=i=02nXWmiXt+1t[i]Pt+1tXX=i=02nXWci[Xt+1t[i]-X^t+1t][Xt+1t[i]-X^t+1t]T+Qt

WmiWci는 가중치를 의미하며, Wm0=λ/(nX+λ),Wm0+(1-α2+β),Wmi=Wci=1/2(nX+λ)로 정의되며, 𝛽는 상태 벡터 사전분포(Prior)의 가우시안(Gaussian) 형태와 관련된 계수로써 2가 사용되었다(van der Merwe and Wan, 2001). 최종적으로 계산된 X^t+1tPt+1tXX는 시간-업데이트 과정에서 예측된 상태 백터의 평균과 공분산 값이며, 이 값들은 보정 과정에서 측정된 정보를 기반으로 갱신된다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2023-036-04/N0040360403/images/Figure_jcoseik_36_04_03_F3.jpg
Fig. 3.

Schematic diagram of unscented transformation

https://cdn.apub.kr/journalsite/sites/jcoseik/2023-036-04/N0040360403/images/Figure_jcoseik_36_04_03_F4.jpg
Fig. 4.

Algorithm of unscented kalman filter

보정을 위해 예측 과정과 마찬가지로 다음과 같이 시그마포인트를 계산한다. 이때, 예측을 수행하는 과정이 아니므로 스케일링 파라미터들은 생략된다.

(7)
Xt+1t[0]=X^t+1tXt+1t[i]=X^t+1t+Pt+1tXXi,i=1,...,nXXt+1t[i]=X^t+1t+Pt+1tXXi-nX,i=nX+1,...,2nX

위와 같이 시그마포인트를 계산한 후, Measurement-output 모델인 식 (4)를 이용하여 다음과 같이 대상 동적 시스템의 거동 값을 예측하고, 이를 토대로 두 가지 공분산을 계산한다.

(8)
ςt+1t[i]=h(Xt+1t[i]),i=0,...,2nXz^t+1t=i=02nX=Wmiςt+1t[i]νt+1=zt+1-z^t+1tPt+1tzz=i=02nX=Wciςt+1t[i]-z^t+1tςt+1t[i]-z^t+1tT+Rt+1Pt+1txz=i=02nX=WciXt+1t[i]-x^t+1tςt+1t[i]-z^t+1tT

ςt+1t[i]는 예측된 상태벡터의 시그마포인트인 Xt+1t[i]로부터 계산된 동적 시스템의 거동과 관련된 벡터이며, 본 연구에서는 예측된 가속도에 해당하는 시그마포인트이다. z^t+1t는 예측가속도 시그마포인트의 평균 값이며, zt+1는 현 시간에서의 계측된 가속도 u¨t+1를 의미한다. 그리고 두 값의 차이인 υt+1은 이론적 잔차(Innovation)로, 최종적으로 상태벡터를 업데이트할 때 사용된다. Pt+1tzzPt+1txz는 각각 예측된 가속도의 공분산과 예측 상태벡터와 가속도 벡터의 교차 공분산(Cross- Covariance)이다. 이 두 값으로부터 다음 수식과 같이 칼만 게인(Kalman Gain)이 계산된다.

(9)
K=Pt+1txzPt+1tzz-1x^t+1t+1=x^t+1t+Kνt+1Pt+1t+1xx=Pt+1txx-KPt+1tzzKT

K는 칼만필터에서 가장 중요한 칼만 게인 값이며, 이를 통해 식 (9)와 같이 현 시간에서의 예측과 보정 과정이 수행되어 산정된 최종 상태 벡터(x^t+1t+1)와 공분산(Pt+1t+1xx)이 계산된다.

3.2 다중 분산점을 이용한 기법 제안

대상 동적 시스템의 상태를 정확히 추정하는 최적의 UKF로 튜닝(Tuning)하기 위해서는 식 (6)(8)에 표현되어 있는 State-Space 모델에서 두 가지 인공오차 공분산을 의미하는 QtRt+1에 적절한 값을 결정해야 한다. 추가적으로, UKF의 핵심인 Unscented Transformation을 수행하기 위해 계산되는 시그마포인트의 분포를 결정하는 스케일링 파라미터 𝛼를 선정해야하는 문제가 남아있다. 결론적으로 총 세 가지의 UKF의 자체 파라미터를 설정해야 하는데, 앞의 두 공분산은 State- space 모델의 정의이기 때문에 튜닝 과정에서 제외 불가하다. 반면에 시그마포인트의 분포를 결정하는 𝛼는 Unscented Transformation 도입으로 인해 추가된 파라미터이다. 본 연구는 적절한 인공오차 공분산 값을 선정함과 동시에 튜닝 문제의 어려움을 완화하기 위해서 스케일링 파라미터 𝛼를 튜닝 과정 대상에서 제외할 수 있는 방안을 제시하고자 한다.

시그마포인트의 분포는 nX,κ,α 세 가지 파라미터에 의해 결정된다. nX는 주어진 동적 시스템에서 추정하고자 하는 변수의 수이며, 문제 정의에 의해 결정되는 값이다. 두 번째로 𝜅는 3.1절에서 설명한바와 같이 통계적 고차모멘트를 고려할 경우에 사용되는 값이며 일반적으로 0으로 사용된다. 따라서, 시그마포인트의 분포는 𝛼에 의해 주로 결정되며, State-Space 모델의 비선형성을 정확히 표현하기 위해서는 적절한 값의 선정이 필수적이다. 최근 UKF의 계산 성능 고도화 연구를 수행한 Papakonstantinou 등(2022)에서 인공오차 값을 고정하고, 𝛼값을 변화하였을 때 상태 추정 정확도가 달라지는 결과를 제시하였다.

적절한 𝛼란 Unscented Transformation이 대상 State-Space 모델의 비선형성을 잘 고려할 수 있는 시그마포인트 분포를 갖게 하는 값이며, 일반적으로 𝛼값의 범위는 0.01부터 1 사이로 제안된다(Candy, 2006). 하지만, 문제에 따라서 더 작은 값, 예를 들어 0.001이 사용되었을 때 주어진 동적시스템의 상태를 더욱 적절히 추정할 수 있음이 확인되었다(Papakonstantinou et al., 2022). Fig. 5의 그래프에서 Mean(빨간점)과 Covariance(검은 실선)는 각각 상태 벡터의 평균과 공분산을 의미하며, 파란 삼각형은 𝛼 = 1, 0.5, 0.1, 0.01값에 따른 시그마포인트 세트를 의미한다. 𝛼값이 작아짐에 따라 시그마포인트가 평균 근처에서 분포함을 확인할 수 있다. 이때 모든 적용예제에서 항상 잘 작동하는 시그마포인트 분포는 존재하지 않으며, 각 경우마다 적절한 스케일링 파라미터를 설정함에 따라 정확한 상태 이력 추정이 가능할 수도 있다. 하지만 무수히 많은 𝛼에 대해서 시행착오 방법 혹은 전역최적화와 같이 시간을 필요로 하는 작업을 수행하는 것은 실제 적용에 있어서 매우 비효율적이다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2023-036-04/N0040360403/images/Figure_jcoseik_36_04_03_F5.jpg
Fig. 5.

Sigma-points’ distribution based on each scaling parameter

따라서 본 연구에서는 여러 𝛼값으로 계산된 시그마포인트 분포 세트를 동시에 사용하는 것을 제안하며, 이를 다중 분산점 칼만필터(UKF using Multiple Sigma-Points, UKF-MSP)로 명명한다. Fig. 6은 𝛼 = 1, 0.5, 0.1, 0.05, 0.01일 때 시그마포인트 분포들을 표현한 것이며, 단일 시그마포인트 분포를 사용하는 경우와 달리 상태 벡터의 공분산의 통계적 특성을 자세히 파악하게 된다. 한편 각 시그마포인트는 한번의 구조 계산을 의미하며 기존에 비해 계산량이 증가되지만, 수 백개의 포인트를 사용하는 Ensemble-KF(Thurin et al., 2019), 수 만개의 포인트를 사용하는 Particle Filte(Yi and Song, 2018)들에 비해서는 여전히 적은 계산량을 요구하며, 다양한 𝛼값을 조사하는 과정이 생략된다는 이점이 있다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2023-036-04/N0040360403/images/Figure_jcoseik_36_04_03_F6.jpg
Fig. 6.

Multiple sigma-points’ distribution about 5 scaling parameters

UKF-MSP를 수행하기 위해 기존 UKF 수식 일부가 수정되어야 한다. 다음 식은 단일 시그마포인트 벡터를 계산하는 식 (5)를 간략히 표현한 것이다.

(10)
Xtt[·][αj]=Xtt[0]Xtt[0]±(nX+λ(αj))Pttxx

위 식에서 j=1,...n이며, n은 사용된 시그마포인트 벡터의 개수로 본 연구에서는 3이 사용되었다. 다음은 수정된 예측 과정이다.

(11)
Xt+1t[i][αj]=fXtt[i][αj],i=0,...,2nX,j=1,...,nx^t+1t=1nj=1ni=1ni=02nXWm[i][αj]Xt+1t[i][αj]Pt+1txx,[αj]=i=02nXWc[i][αj]Xt+1t[i][αj]-x^t+1tXt+1t[i][αj]-x^t+1tTPt+1txx=1nj=0nPt+1txx,[αj]+Qt

여기서 Pt+1txx,[αj]를 계산하는 이유는 다음 보정 과정에서 시그마포인트를 계산할 때 사용하기 위함이며 식 (10)과 같이 간략히 표현할 수 있다.

(12)
Xt+1t[·][αj]=Xt+1t[0]Xt+1t[0]±Pt+1txx,[αj],j=1,...,n

마지막으로 수정된 보정 과정은 다음과 같다.

(13)
ςt+1t[i][αj]=hXt+1t[i][αj],i=0,...,2nX,j=1,...,nz^t+1t=j=1ni=02nXWm[i][αj]ςt+1t[i][αj]νt+1=zt+1-z^t+1tPt+1tzz=j=0ni=02nXWc[i][αj]ςt+1t[i][αj]-z^t+1tςt+1t[i][αj]-z^t+1tT+Rt+1Pt+1txz=j=1ni=02nXWc[i][αj]Xt+1t[i][αj]-X^t+1tςt+1t[i][αj]-z^t+1tT

최종적으로 기존 UKF 알고리즘을 나타낸 Fig. 4에서 식 (5)(6) 대신 식 (10)(11)을, 식 (7)(8) 대신 식 (12)(13)을 사용하여 UKF-MSP를 수행한다.

다음으로 2장에서 설명한 단자유도 Bouc-Wen 모델에 대해 지진으로 인한 구조손상 시나리오를 가정하고, 각 시나리오에 대해 본래의 UKF와 제안된 UKF-MSP를 수행하여 그 결과를 비교하도록 한다.

4. 수치 예제 적용

앞서 2장에서 소개된 Bouc-Wen 모델과 San Fernando 지진에 대하여 Table 1과 같이 두 가지 시나리오를 가정하고, 각 시나리오에서 측정된 가속도 이력 정보를 이용하여 대상 시스템의 강성 변화 이력을 추정하였다. 이때 측정된 가속도 이력에 여러 신호 대 잡음비(Signal-to-Noise, SNR)를 추가하여 UKF와 UKF-MSP의 결과들을 도출하고 비교하였다. 두 필터들을 수행하기 위한 초기 설정으로 x^0=k0,P0xx=diag(k0·10-16)을 사용하였고, 측정 인공오차의 공분산은 diag(Rt)=RMS(u¨0:dt:T)·10-2로 가정하였다. 여기서 diag(·)RMS(·)는 각각 행렬의 대각성분과 Root-Mean-Square를 의미한다.

Table 1

Sudden damage scenarios for san fernando earthquake

Damage Description
Scenario 1 Stiffness (k) decreases suddenly by 40% at 4.96 sec
Scenario 2 Stiffness (k) decreases suddenly by 20% at 12.63 sec

4.1 분산점 칼만필터 적용 및 결과

UKF-MSP의 타당성을 검증하기 전에 기존 UKF를 수행하였다. 각 시나리오의 측정된 가속도에 SNR = 50, 30, 13dB을 가지는 백색 잡음(White Noise)을 생성하여 추가하고 Fig. 7에 나타내었다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2023-036-04/N0040360403/images/Figure_jcoseik_36_04_03_F7.jpg
Fig. 7.

Noisy-acceleration when SNR = 50, 30, 13dB

Fig. 8은 시나리오1에 대해 𝛼 = 1을 사용하여 추정한 결과이다. 위에서부터 SNR = 50, 30, 13dB일 때의 추정 결과이며, 각 결과는 서로 다른 시스템 인공오차 공분산 Qt로부터 도출된 최적의 UKF 결과이다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2023-036-04/N0040360403/images/Figure_jcoseik_36_04_03_F8.jpg
Fig. 8.

Estimation of UKF: scenario 1, 𝛼 = 1

본 연구에서는 시행착오 방법을 통해 Qt를 조정하여 식 (8)(13)에서 서술된 칼만필터의 부산물인 이론적잔차(Innovation)의 시간이력, ν0:dt:T의 RMS를 최소화(Minimum)하는 방식으로 최적의 UKF를 도출하였다. Table 2는 최적의 UKF일 때의 Qt와 이론적 잔차의 RMS의 값을 나타낸다.

Table 2

Optimal Qt with min-RMS of innovation: scenario 1

Scenario1 diag(Qt) RMS of ν0:dt:T
𝛼 = 1 50dB k0·10-2.8 0.0002
30dB k0·10-3.6 0.0071
13dB k0·10-6.6 0.2110
𝛼 = 0.1 50dB k0·10-3.2 0.0031
30dB k0·10-3.6 0.0120
13dB k0·10-6.9 0.2248

Table 2의 결과를 보면, 오차의 크기가 커질수록 ν0:dt:T의 RMS 값이 증가함을 알 수 있으며, 동시에 최적의 UKF를 위한 Qt값이 다름을 확인할 수 있다. 또한 시그마포인트 분포를 결정하는 𝛼가 달라질 경우에도 최적의 결과를 도출하기 위해서는 새로 튜닝과정을 거쳐야만 한다. 특히 오차가 30dB에서 13dB로 커질 때 Qt의 지수가 상당히 다른 값을 가지게 되며 동시에 ν0:dt:T의 RMS 값도 크게 변화한다. Fig. 9는 𝛼 = 0.1일 때, SNR = 50, 30, 13dB의 결과이다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2023-036-04/N0040360403/images/Figure_jcoseik_36_04_03_F9.jpg
Fig. 9.

Estimation of UKF: scenario 1, 𝛼 = 1

Table 2와 같이 Table 3은 시나리오2에 대해 최적의 UKF일 때의 Qt와 이론적 잔차의 RMS의 값을 서술하였다. Fig. 1011은 각각 𝛼 = 1, 0.1일 때의 최적 UKF의 결과이다.

Table 3

Optimal Qt with min-RMS of innovation: scenario 2

Scenario 2 diag(Qt) RMS of ν0:dt:T
𝛼 = 1 50dB k0·10-3.4 0.0002
30dB k0·10-4.1 0.0039
13dB k0·10-5.0 0.1716
𝛼 = 0.1 50dB k0·10-3.0 0.0013
30dB k0·10-4.1 0.0068
13dB k0·10-5.0 0.1720

https://cdn.apub.kr/journalsite/sites/jcoseik/2023-036-04/N0040360403/images/Figure_jcoseik_36_04_03_F10.jpg
Fig. 10.

Estimation of UKF: scenario 2, 𝛼 = 1

https://cdn.apub.kr/journalsite/sites/jcoseik/2023-036-04/N0040360403/images/Figure_jcoseik_36_04_03_F11.jpg
Fig. 11.

Estimation of UKF: scenario 2, 𝛼 = 0.1

시나리오1, 2 모두 𝛼 = 1일 때 급격한 변화를 빠르게 추정해내는 것을 확인할 수 있으며, 𝛼 = 0.1일 때 상대적으로 추정값의 변동성(Fluctuation)이 큰 것을 확인할 수 있다. 시나리오2의 경우 시나리오1에 비해 상대적으로 작은 강성 감소가 발생했기 때문에 오차의 크기가 커질수록 추정값의 불안정성이 증가하는 모습을 보인다. 특히, 측정 오차 13dB 결과에서 𝛼 = 1일 때 20~25초 사이 추정값은 정해에서 상당히 벗어나며, 𝛼 = 0.1일 때 25초 이후 추정값이 발산함을 확인할 수 있다.

4.2 다중 분산점 칼만필터 적용 및 결과

다음으로 𝛼 = 1, 0.5, 0.1을 사용하는 UKF-MSP를 수행하였다. Fig. 12는 시나리오1의 결과이며 이때 사용된 시스템 인공오차는 diag(Qt)=k0·10-6.8로 설정되었다. 이 공분산 값은 측정오차 13dB일 때의 최적화 값이며, 이 값을 다른 두 오차의 경우에도 사용하였다. 기존 UKF의 결과인 Fig. 8, 9과 비교해 보면 불규칙한 변동성이 거의 사라짐을 확인할 수 있다. 또한 실제 가속도 이력을 측정할 때 예상되는 가장 큰 오차에 대한 최적의 시스템 인공오차를 도출한 후, 해당 값을 그 외 작은 측정 오차때에 사용하여도 상당히 안정적으로 시스템 변수 이력을 잘 추정할 수 있음을 결과로부터 확인할 수 있다. Fig. 13은 다른 두 가속도 측정값에 대한 최적 인공오차를 설정했을 때의 결과이다. 측정오차 50, 30dB일 때 각각 diag(Qt)=k0·10-4,k0·10-5로 설정하였고, Fig. 12의 결과와 비교하여 급격한 변화를 상대적으로 빠르게 추정함을 확인할 수 있다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2023-036-04/N0040360403/images/Figure_jcoseik_36_04_03_F12.jpg
Fig. 12.

Estimation of UKF-MSP: scenario 1

https://cdn.apub.kr/journalsite/sites/jcoseik/2023-036-04/N0040360403/images/Figure_jcoseik_36_04_03_F13.jpg
Fig. 13.

Estimation of UKF-MSP: scenario 1 when optimal Qt for SNR = 50, 30dB

Fig. 14는 기존 UKF 𝛼 = 1일 때, 13dB에 최적화된 Qt를 이용하여 측정오차 50dB일 때 적용한 결과이다. UKF-MSP와 달리 지진 시간 이력내에 정해에 도달하지 못함을 확인할 수 있다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2023-036-04/N0040360403/images/Figure_jcoseik_36_04_03_F14.jpg
Fig. 14.

Estimation of UKF: scenario 1, 𝛼 = 1, optimal Qt when SNR = 13dB

Fig. 15는 시나리오2에 대한 UKF-MSP의 결과로 이때 사용된 시스템 인공오차는 diag(Qt)=k0·10-7.0로 설정되었으며, 이는 측정 오차 13dB일 때 최적화된 값이다. 앞선 시나리오1과 마찬가지로 작은 측정오차 경우들에 대해서도 정해를 안정적으로 추정함을 확인할 수 있었고, Fig. 16과 같이 각각 측정오차에 최적화된 인공오차를 사용할 경우 결과가 개선됨을 확인할 수 있다. 반면, Fig. 17은 13dB에 최적화된 Qt를 이용한 UKF 𝛼 = 1일 때의 결과이며 발산함을 확인하였다.

https://cdn.apub.kr/journalsite/sites/jcoseik/2023-036-04/N0040360403/images/Figure_jcoseik_36_04_03_F15.jpg
Fig. 15.

Estimation of UKF-MSP: scenario 2

https://cdn.apub.kr/journalsite/sites/jcoseik/2023-036-04/N0040360403/images/Figure_jcoseik_36_04_03_F16.jpg
Fig. 16.

Estimation of UKF-MSP: scenario 2 when diag(Qt)=k0·10-4,k0·10-4.5 for SNR = 50, 30dB

https://cdn.apub.kr/journalsite/sites/jcoseik/2023-036-04/N0040360403/images/Figure_jcoseik_36_04_03_F17.jpg
Fig. 17.

Estimation of UKF: scenario 2, 𝛼 = 1, optimal Qt when SNR = 13dB

가정한 두 시나리오에서 기존 UKF는 각 스케일링 파라미터, 𝛼에 따라, 그리고 측정 오차의 크기에 따라 적절한 인공오차 설정이 필요하였으며, 특히 작은 강성 감소 시나리오인 2에서는 추정 값에서 심한 변동성을 확인하였다. 반면에 UKF- MSP는 두 시나리오 모두 안정적인 추정을 성공함과 동시에 스케일링 파라미터 선정이 필요 없으며, 인공오차 값에 대한 추정 성능의 민감도가 크게 감소하여 최적의 인공오차 설정을 하지 않음에도 불구하고 정해에 가까운 추정이 가능하였다. 이때 작은 측정 오차를 갖는 가속도 측정값들에 대한 해석결과에서 급격한 변화를 빠르게 파악하지 못하는 결과처럼 보이지만, 이는 최적화된 인공오차를 사용하지 않았기 때문이며, 이들에 대해서도 최적의 인공오차를 설정하여 급격한 변화를 잘 추정함을 확인하였다. 추가적으로, 13dB의 결과에서 마지막 수렴 값이 다소 정해에서 벗어나 있지만, 지진이 끝난 직후 잠시 동안 구조물이 거동할 때의 추가 측정값이 활용될 수 있다면 개선될 수 있을 것으로 예상된다.

5. 결 론

본 연구에서는 비선형 동적 시스템을 표현하기 위해 Bouc- Wen 모델을 사용하여 단자유도 시스템을 가정하였다. Bouc- Wen 모델은 효율적이고 부드러운 지진 이력현상 곡선을 표현할 수 있으며, Bouc-Wen 모델이 적용된 비선형 동적 시스템의 상태를 추정하기위해 UKF가 사용될 수 있다. 하지만, UKF를 적용하기 위해서는 시그마포인트 분포를 결정하는 스케일링 파라미터와 두 가지 인공오차 공분산을 결정해야 하는 튜닝(Tuning) 문제가 남아있다. 칼만필터 적용시 필요로 하는 State- Space 모델의 정의에 따라 두 인공오차는 필수적으로 결정해야하는 값인 반면에, 스케일링 파라미터는 UKF의 Unscented Transformation 도입으로 인해 추가된 파라미터이다. 본 연구에서는 다중 시그마포인트 분포를 사용하는 UKF-MSP를 제안하여 스케일링 파라미터 선정 문제를 피하고자 하였고, 예제를 통해 우수성을 증명하였다.

제안된 기법은 칼만필터의 고질적 문제인 튜닝 문제를 완화하기 위해 고안되었다. 그런데 수행된 예제를 통해 제안된 기법이 본 목표였던 스케일링 파라미터 설정 완화뿐만 아니라, 다른 인공오차 설정에 대한 민감도를 감소시키는 효과가 있음을 확인하였다. 이는 UKF-MSP를 실제 상황에 적용시 추가적 튜닝과정 없이 동적 시스템의 상태를 빠르게 파악할 수 있는 큰 장점이 될 수 있다. 향후 연구에서는 Bouc-Wen 모델을 기반으로 복잡한 비선형성을 보이는 대형 구조물 혹은 설비를 모사하는 다자유도 동적 시스템에 대해 확장 적용이 이루어져야 할 것이다. 이때 일반적으로 알려진 스케일링 파라미터 범위에서 각각 작은 값과 큰 값을 요구하는 시스템들을 구축하여 UKF-MSP 거동을 파악한다면 더욱 적절한 MSP 설정 및 사용법을 제시할 수 있을 것으로 예상된다. 이와 같이 향상된 UKF- MSP를 통해 Bouc-Wen 모델 기반의 복잡 동적 시스템의 특성들을 안정적으로 파악 가능하다면, Reduced-Order 모델 구축이 가능하며 취약도 산정과 같은 연구에 활용될 수 있을 것이다. 또한 이중 칼만필터(Dual-Kalman Filter)와 같이 매시간 스텝마다 다중 시그마포인트 분포의 가중치를 설정할 수 있는 연구가 진행된다면 대상 동적 시스템의 급격한 변화에 더욱 잘 적응할 수 있는 기법이 개발될 수 있을 것으로 기대된다.

Acknowledgements

이 논문은 한국건설기술연구원 2023년 구조연구본부 목적형 R&R 과제 “국민 안전과 건강한 인프라 환경을 위한 지속 가능한 인프라 구조기술 연구”에 의하여 연구되었습니다.

References

1
Baber, T.T., Noori, M.N. (1984) Random Vibration of Pinching Hysteretic Systems, J. Eng. Mech., 110(7), pp.1036~1049. 10.1061/(ASCE)0733-9399(1984)110:7(1036)
2
Baber, T.T., Wen, Y.K. (1981) Random Vibration Hysteretic, Degrading Systems, J. Eng. Mech. Div., Am. Soc. Civil Eng., 107(6), pp.1069~1087. 10.1061/JMCEA3.0002768
3
Candy, J.V. (2006) Bayesian Signal Process: Classical, Modern, and Particle Filtering Methods, John Wiley & Sons, Inc.: Hoboken, NJ, USA.
4
Chatzi, E.N., Smyth, A.W. (2009) The Unscented Kalman Filter and Particle Filter Methods for Nonlinear Structural System Identification with Non-Collocated Heterogeneous Sensing, Struct. Control & Health Monit., 16(1), pp.99~123. 10.1002/stc.290
5
Foliente, G.C., Singh, M.P., Noori, M.N. (1996) Equivalent Linearization of Generally Pinching Hysteretic, Degrading Systems, Earthq. Eng. & Struct. Dyn., 25(6), pp.611~629. 10.1002/(SICI)1096-9845(199606)25:6<611::AID-EQE572>3.0.CO;2-S
6
Fujimura, K., Der Kiureghian, A. (2007) Tail-Equivalent Linearization Method for Nonlinear Random Vibration, Probab. Eng. Mech., 22(1), pp.63~76. 10.1016/j.probengmech.2006.08.001
7
Kalman, R.E. (1960) A New Approach to Linear Filtering and Prediction Problems, J. Basic Eng., 82(1), pp.35~45. 10.1115/1.3662552
8
Lee, S.-H., Song, J. (2019) Regularization-Based Dual Adaptive Kalman Filter for Identification of Sudden Structural Damage Using Sparse Measurements, Appl. Sci., 10(3), p.850. 10.3390/app10030850
9
Nguyen, V.H., Lee, S.-H., Lee, J.H. (2022) Application of the Unscented Kalman Filter to Estimate the Material Properties of a Layered Half-Space, Int. J. Appl. Mech., 15(2), p.2350005. 10.1142/S1758825123500059
10
Noori, M.N., Choi, J., Davoodi, H. (1986) Zero and Nonzero Mean Random Vibration Analysis of a New General Hysteresis Model, Probab. Eng. Mech., 1(4), pp.192~201. 10.1016/0266-8920(86)90012-3
11
Papakonstantinou, K.G., Amir, M., Warn, G.P. (2020) A Scaled Spherical Simplex Filter (S3F) with a Decreased n + 2 Sigma Points Set Size and Equivalent 2n + 1 Unscented Kalman (UKF) Accuracy, Mech. Syst. & Signal Process., 163(2022), p.107433. 10.1016/j.ymssp.2020.107433
12
Song, J., Der Kiureghin, A. (2006) Generalized Bouc-Wen Model for Highly Asymmetric Hysteresis, J. Eng. Mech., 123(6), pp.610~618. 10.1061/(ASCE)0733-9399(2006)132:6(610)
13
Thurin, J., Brossier, R., Metivier, L. (2019) Ensemble-based Uncertainty Estimation in Full Waveform Inversion, Geophysical J. Int., 219(3), pp.1613~1635. 10.1093/gji/ggz384
14
Vaiana, N., Sessa, S., Marmo, F., Rosati, L. (2018) A Class of Uniaxial Phenomenological Models for Simulating Hysteretic Phenomena in Rate-Independent Mechanical Systems and Materials, Nonlinear Dyn., 93(3), pp.1647~1669. 10.1007/s11071-018-4282-2
15
Van der Merwe, R., Wan, E. (2001) The Square-Root Unscented Kalman Filter for State and Parameter-Estimation, In Proceedings of the 2001 IEEE International Conference on Acoustics, Speech, and Signal Processing, Proceedings (Cat. No.01CH37221). Salt Lake City, UT, USA, 7-11 May 2001, 6, pp.3461~3464.
16
Wang, C.-H., Wen, Y.K. (1998) Reliability and Redundancy of pre-Northridge Low-Rise Steel Building under Seismic Excitation, Rep. No. UILU-ENG-99-2002, Univ. Illinois at Urbana-Champaign, Champaign, Ill.
17
Yi, S.-R., Song, J. (2018) Particle Filter Based Monitoring and Prediction of Spatiotemporal Corrosion Using Successive Measurements of Structural Responses, Sens., 18(11), p.3909. 10.3390/s1811390930428593PMC6263897
18
Yi, S.-R., Wang, Z., Song, J. (2017) Bivariate Gaussian Mixture- based Equivalent Linearization Method for Stochastic Seismic Analysis of Nonlinear Structures, Earthq. Eng. Struct. Dyn., 47(3), pp.678~696. 10.1002/eqe.2985
페이지 상단으로 이동하기