ISO/TR 16250:2013

TECHNICAL ISO/TR
REPORT 16250
First edition
2013-07-15
Road vehicles — Objective rating
metrics for dynamic systems
Véhicules routiers — Mesures pour l’évaluation objective des
systèmes dynamiques
Reference number
©
ISO 2013
© ISO 2013
All rights reserved. Unless otherwise specified, no part of this publication may be reproduced or utilized otherwise in any form
or by any means, electronic or mechanical, including photocopying, or posting on the internet or an intranet, without prior
written permission. Permission can be requested from either ISO at the address below or ISO’s member body in the country of
the requester.
ISO copyright office
Case postale 56 • CH-1211 Geneva 20
Tel. + 41 22 749 01 11
Fax + 41 22 749 09 47
E-mail copyright@iso.org
Web www.iso.org
Published in Switzerland
ii © ISO 2013 – All rights reserved

Contents Page
Foreword .iv
Introduction .v
1 Scope . 1
2 Terms and definitions . 1
3 Symbols and abbreviated terms . 1
3.1 General abbreviated terms . 1
3.2 General symbols and subscripts . 2
3.3 CORA . 2
3.4 EARTH and EEARTH . 3
3.5 Model reliability metric . 4
3.6 Bayesian confidence metric . 5
3.7 Overall ISO rating . 5
4 General requirements to the data . 5
5 CORA metric . 6
5.1 Corridor rating . 6
5.2 Cross-correlation rating . 8
5.3 Step-by-step procedure .10
6 EARTH metric .11
6.1 EARTH phase score .12
6.2 EARTH magnitude score .13
6.3 EARTH slope score .14
6.4 Overall EARTH score .15
6.5 Step-by-step procedure .15
7 Model reliability metric .16
8 Bayesian confidence metric .16
9 ISO metric .18
9.1 CORA corridor method .18
9.2 EEARTH method .18
9.3 Calculation of the overall ISO rating .23
9.4 Meaning of the objective rating score .24
10 Pre-processing of the data .24
10.1 Sampling rate.25
10.2 Filtering .25
10.3 Interval of evaluation .25
11 Limitations .26
11.1 Type of signals .26
11.2 Metrics validation .26
11.3 Meaning of the results .26
11.4 Multiple responses .27
Annex A (informative) Child restraint example .28
Annex B (informative) Sled test example .46
Annex C (informative) Case studies .51
Bibliography .65
Foreword
ISO (the International Organization for Standardization) is a worldwide federation of national standards
bodies (ISO member bodies). The work of preparing International Standards is normally carried out
through ISO technical committees. Each member body interested in a subject for which a technical
committee has been established has the right to be represented on that committee. International
organizations, governmental and non-governmental, in liaison with ISO, also take part in the work.
ISO collaborates closely with the International Electrotechnical Commission (IEC) on all matters of
electrotechnical standardization.
The procedures used to develop this document and those intended for its further maintenance are
described in the ISO/IEC Directives, Part 1. In particular the different approval criteria needed for the
different types of ISO documents should be noted. This document was drafted in accordance with the
editorial rules of the ISO/IEC Directives, Part 2. www.iso.org/directives
Attention is drawn to the possibility that some of the elements of this document may be the subject of
patent rights. ISO shall not be held responsible for identifying any or all such patent rights. Details of
any patent rights identified during the development of the document will be in the Introduction and/or
on the ISO list of patent declarations received. www.iso.org/patents
Any trade name used in this document is information given for the convenience of users and does not
constitute an endorsement.
The committee responsible for this document is ISO/TC 22, Road vehicles, Subcommittee SC 10, Impact
test procedures, and SC 12, Passive safety crash protection systems.
iv © ISO 2013 – All rights reserved

Introduction
Computer-Aided Engineering (CAE) has become a vital tool for product development in the automobile
industry. Various computer programs and models are developed to simulate dynamic systems. To
maximize the use of these models, their validity and predictive capabilities need to be assessed
quantitatively. Model validation is the process of comparing CAE model outputs with test measurements
in order to assess the validity or predictive capabilities of the CAE model for its intended usage. The
fundamental concepts and terminology of model validation have been established mainly by standard
[6]
committees including the United States Department of Energy (DOE), the American Institute of
[1]
Aeronautics and Astronautics (AIAA), the Defense Modeling and Simulation Office (DMSO) of the US
[5]
Department of Defense (DOD), the American Society of Mechanical Engineers Standards Committee
[2]
(ASME) on verification and validation of Computational Solid Mechanics, Computational Fluid
[3] [4][22][23]
Dynamics and Heat Transfer, and various other professional societies.
One of the critical tasks to achieve quantitative assessment of models is to develop a validation metric
that has the desirable metric properties to quantify the discrepancy between functional or time history
[7][19][20]
responses from both physical test and simulation result of a dynamic system. Developing
quantitative model validation methods has attracted considerable researchers’ interest in recent years.
[12][13][14][18][20][21][26][28][29][32]
However, the primary consideration in the selection of an effective
metric should be based on the application requirements. In general, the validation metric is a quantitative
measurement of the degree of agreement between the physical test and simulation result.
In this Technical Report, four state-of-the-art objective rating metrics are investigated and they are:
[10][30][31]
CORrelation and Analysis (CORA) metric, Error Assessment of Response Time Histories
[28][34] [18][27][35] [14][16][36]
(EARTH) metric, model reliability metric, and Bayesian confidence metric.
Multiple dynamic system examples for both tests and CAE models are used to show their advantages
and limitations. Further enhancements of the CORA corridor rating and the development of an Enhanced
Error Assessment of Response Time Histories (EEARTH) metric are proposed to improve the robustness
of these metrics. A new combined objective rating metric is developed to standardize the calculation
of the correlation between two time history signals of dynamic systems. Multiple vehicle safety case
studies are used to demonstrate the effectiveness and usefulness of the proposed metric for an ISO
Technical Report.
TECHNICAL REPORT ISO/TR 16250:2013(E)
Road vehicles — Objective rating metrics for dynamic
systems
1 Scope
This Technical Report specifies a method to calculate the level of correlation between two non-ambiguous
signals. The focus of the methods described in this Technical Report is on the comparison of time-history
signals or functional responses obtained in all kinds of tests of the passive safety of vehicles and the
corresponding numerical simulations. It is validated with signals of various kinds of physical loads such
as forces, moments, accelerations, velocities, and displacements. However, other applications might be
possible too, but are not in the scope of this Technical Report.
2 Terms and definitions
For the purposes of this document, the following terms and definitions apply.
2.1
filtering
smoothing of signals by using standardized algorithms
2.2
goodness or level of correlation
similarity of two signals
2.3
interval of evaluation
time domain that is used to calculate the correlation between two signals
2.4
rating
rating score
calculated value that represents a certain level of correlation (objective rating)
2.5
sampling rate
recording frequency of a signal
2.6
time sample
pair values (e.g. time and amplitude) of a recorded signal
2.7
time-history signal
physical value recorded in a time domain; those signals are non-ambiguous
3 Symbols and abbreviated terms
3.1 General abbreviated terms
CAE Computer-Aided Engineering
CORA CORrelation and Analysis
DTW Dynamic Time Warping
EARTH Error Assessment of Response Time Histories
EEARTH Enhanced Error Assessment of Response Time Histories
SME Subject Matter Expert
3.2 General symbols and subscripts
C, C(t) analysed signal (CAE signal)
T, T(t) reference signal (test signal)
t time signal (axis of abscissa)
∆t interval between two time samples
t time zero of an event (e.g. test, crash, impact, etc.)
t starting time of the interval of evaluation
start
t ending time of the interval of evaluation
end
N total number of sample points (e.g. time steps) between the starting time,
t , and ending time, t
start end
Ν all natural numbers without zero
>0
3.3 CORA
Z CORA rating
CORA
Z corridor rating
Zt() corridor rating at time t (curve)
Z cross-correlation rating
Z phase-shift rating
2a
Z size rating
2b
Z shape (progression) rating
2c
w weighting factor of the corridor rating, Z
Z1 1
w weighting factor of the cross-correlation rating, Z
Z2 2
w weighting factor of the phase-shift rating, Z
Za2 2a
w weighting factor of the size rating, Z
Zb2 2b
w weighting factor of the shape rating, Z
Zc2 2c
k exponent factor for calculating the corridor rating between the inner and outer corridors
Z1
k exponent factor for calculating phase-shift rating, Z
Za2 2a
k exponent factor for calculating size rating, Z
Zb2 2b
k exponent factor for calculating shape rating, Z
Zc2 2c
T absolute maximum amplitude of the reference signal, T
norm
a relative half width of the inner corridor
2 © ISO 2013 – All rights reserved

b relative half width of the outer corridor
δ half width of the inner corridor
i
δ half width of the outer corridor
o
δ ()t lower/upper inner corridor at time t (curve)
i
δ ()t lower/upper outer corridor at time t (curve)
o
D coefficient of the allowable lower limit of the phase shift
min
D coefficient of the allowable upper limit of the phase shift
max
F sum of the square of the area for the time-shifted evaluated curve, C
C
F sum of the square of the area for the reference curve, T
T
INT percentage of the minimum remaining overlapping time of the reference and
min
evaluation curves after time shift
m shift of a signal along the axis of abscissa
m minimum m shift of a signal
min
m maximum m shift of a signal
max
n number of samples
n time step shifted to get the maximum cross correlation
min
ρ cross correlation
ρ()m cross correlation at shift, m
δ phase-shift time at the maximum cross correlation, ρ
δ lower limit of CORA phase shift
min
δ upper limit of CORA phase shift
max
3.4 EARTH and EEARTH
E overall EARTH score
E
E EARTH magnitude score
M
E EARTH phase score
P
E EARTH slope (topology) score
S
w weighting factor of the magnitude score, E
M M
w weighting factor of the phase score, E
P P
w weighting factor of the slope score, E
S S
k exponent factor for calculating the magnitude score, E
M M
k exponent factor for calculating the phase score, E
P P
k exponent factor for calculating the slope score, E
S S
ε EARTH magnitude error
mag
ε EARTH slope error
slope
*
ε maximum allowable magnitude error
M
*
ε maximum allowable percentage of time shift
P
*
ε maximum allowable slope error
S
Ct() mean value of CAE curve
ts ts
C , Ci() truncated and shifted CAE curve
ts+d ts
C derivative CAE curve, C
ts+w ts
C warped CAE curve, C
ts++dw ts
C derivative warped CAE curve, C
Tt() mean value of test curve
ts ts
T , Tj() truncated and shifted test curve
ts+d ts
T derivative test curve, T
ts+w ts
T warped test curve, T
ts++dw ts
T derivative warped test curve, T
ρ maximum cross correlation of all ρ ()m and ρ ()m
E L R
ρ ()m cross correlation — signal is moved to the left
L
ρ ()m cross correlation — signal is moved to the right
R
di(, j) local cost function to perform the dynamic time warping
m time steps moved to evaluate the EARTH phase error
n number of time shifts to get, ρ
ε E
3.5 Model reliability metric
T absolute maximum amplitude of the reference signal, T
norm
a reliability target
b threshold factor of the reliability assessment
ε lower bound of the threshold interval
L
ε upper bound of the threshold interval
U
L
ε lower bound of the Bayesian interval hypothesis in probabilistic principal
Φ
component analysis space (PPCA)
U
ε lower bound of the Bayesian interval hypothesis in PPCA space
Φ
r model reliability
P cumulative probability
4 © ISO 2013 – All rights reserved

Φ , Φ()t pn× reduced data matrix
3.6 Bayesian confidence metric
A constant vector
B Bayes factor for multivariate case
iM
δ likelihood function
ε predefined threshold vector
Φ
H null hypothesis
H alternative hypothesis
K confidence of accepting the model
κ measure of confidence within an interval
*
Λ variance of error variable, ε
i
*
μ p mean values obtained from Φ
Φ
N normal distribution
N(,ρΛ) normal distribution of δ with mean vector, p , and variance matrix, Λ
π prior probability of hypothesis
ρ prior mean, δ
*
Σ variance matrix of Φ
Φ
Φ , Φ()t difference curve between test curve, T , and CAE curve, C
f()δ prior density function of δ
3.7 Overall ISO rating
R combined rating of EEARTH and the CORA corridor method
E EEARTH rating score
Z CORA corridor rating ( ZZ= )
w weighting factor of the EEARTH rating, E
E
w weighting factor of the CORA corridor rating, Z
Z
r rank of the sliding scale of the ISO metric
SC ()r lower threshold of rank, r
lower
SC ()r upper threshold of rank, r
upper
4 General requirements to the data
The metrics described in this Technical Report require non-ambiguous curves (e.g. time-history curves).
Furthermore, it is required that the reference curve, T(t), and the evaluated curve, C(t), are both defined
between starting time, t , and ending time, t . Both curves shall have the same number of sample
start end
points, N, with a constant time interval, ∆t, within the evaluation interval.
5 CORA metric
[10][30][31]
The objective evaluation metric called CORA — correlation and analysis — uses two
independent sub-ratings, a corridor rating, and a cross-correlation rating to assess the correlation of
two signals. The rating structure of CORA is shown in Figure 1.
Figure 1 — CORA rating structure
The corridor and cross-correlation ratings are used to compensate each other’s disadvantages, and the
CORA rating tool is trying to separate an engineer’s knowledge from the objective rating metric by
using external parameters. However, it is possible to fine-tune the evaluation to the specific needs of the
applications by adjusting those metric parameters to reflect the SME’s knowledge of the applications.
The corridor rating, Z , calculates the deviation between both curves with the help of user-defined or
automatically generated corridors. The cross-correlation rating, Z , analyses specific curve
characteristics, such as phase shift, Z , size, Z , and shape of the signals, Z . The rating results
2a 2b 2c
range from “0” (no correlation) to “1” (perfect match). The influence of the sub-ratings on the global
rating is adjusted by user-defined weighting factors. Formulae (1) and (3) show how to calculate the
CORA rating by using weighting factors [see Formulae (2) and (4)]. Details of each sub-rating are
introduced in the following subsections.
Zw=⋅Zw+⋅Z (1)
CORA ZZ11 22
ww+=1 (2)
ZZ12
Zw=⋅Zw+⋅Zw+⋅Z (3)
22Za 22aZ bb22Zc 2c
www++ =1 (4)
Za222Zb Zc
5.1 Corridor rating
The corridor rating calculates the deviation between two signals by means of corridor fitting. The two
sets of corridors, the inner and the outer corridors, are defined along the mean curve. If the evaluated
curve (e.g. CAE curve) is within the inner corridor bounds, a score of “1” is given, and if it is outside the
outer corridors, the rating is set to “0”. The assessment declines from “1” to “0” between the bounds of
6 © ISO 2013 – All rights reserved

inner and outer corridors resulting in three different rating zones as shown in Figure 2. This transition
is user-defined. The compliance with the corridors is calculated at each specific time, t, and the final
corridor rating, Z , of a signal is the average of all ratings, Zt(), at specific times, t.
1 1
[10]
Figure 2 — Rating zones of the corridor method (corridors of constant width)
[17]
The philosophy is to use a narrow inner corridor and a wide outer corridor. It limits the number of
“1” ratings to only good correlations and gives the opportunity to distinguish between poor and fair
correlations. If the outer corridor is too narrow, too many curves of a fair or moderate correlation would
get the same poor rating of “0”, like signals of almost no correlation with the reference. The width of the
corridors can be adjusted in order to reflect the specific signal characteristic, and it can be constant for
the whole duration of the dynamic responses or vary at the different time steps.
This Technical Report applies the most common approach of using the constant corridor widths for the
whole duration of the dynamic response. The parameters a and b define the relative half width of the
0 0
inner and the outer corridors. Both shall be between “0” and “1”, and a must be less than b . The
0 0
absolute half widths of both corridors are defined as the product of relative half width and the absolute
maximum amplitude, T , of the reference signal, T. Formula (5) shows the calculation of T .
norm norm
TT=maxmin() ,max()T (5)
{}
norm
The absolute half width of the inner corridors (absolute distance from reference signal to outer bounds
of the inner corridors) is defined by Formula (6). The calculation of the absolute half width of the outer
corridors [see Formula (7)] is similar to that of the inner corridors.
δ =⋅aT 01≤≤a (6)
in00orm
δ =⋅bT 01≤≤baand on00orm 00
Based on these definitions, the upper and lower bounds of the inner corridors are defined by Formula (8)
and the upper and lower bounds of the outer corridors are defined by Formula (9).
δδ()tT=±()t (8)
ii
δδ()tT=±()t (9)
oo
Formula (10) shows the calculation of the corridor rating for the correlation between the reference signal,
T, and the analysed signal, C, at each evaluation time, t. If the absolute difference between the signals T and
C is less than the half width of the inner corridors, δ , then the rating is set to “1”. The rating is calculated
i
by Formula (10) when the absolute difference between both signals is in between δδ≤−Tt() Ct() ≤ . If
io
the absolute difference between both signals is greater than the half width of the outer corridors, δ , then
o
the rating is set to “0”. The parameter k assesses the location of the analysed signal within the outer
Z1
corridor and it applies the appropriate penalty on the rating score. A linear (k = 1), quadratic (k =2
Z1 Z1
), cubical (k =3), or any other regression relationship can be defined accordingly.
Z1
 1 if Tt()− i

k
Z1
 
 δ −−Tt() Ct()
o
Zt()= k ∈Ν (10)
 

1 Z10>
 
δδ−
oi
 

0 if Tt()−C(()t >δ
o

The final corridor rating, Z , is calculated by averaging all single time step ratings, Zt(), as shown in
1 1
Formula (11). The parameter N represents the total number of sample points (e.g. time steps) between
the starting and ending times of the interval of evaluation.
t
end
Zt()
∑ 1
tt=
start
Z = (11)
N
One of the advantages of the corridor rating is the simplicity and the clearness of the algorithm. It reflects
criteria which are used intuitively in engineering judgment. Sometimes, this simplicity may be the
disadvantage of the method. For example, a small distortion of the phase can lead to a undesirable rating.
5.2 Cross-correlation rating
The cross-correlation rating may compensate for the disadvantages resulting from the corridor rating
by analysing the characteristics of signals. Three sub-ratings (phase, size, and shape) with individual
weighting factors are implemented.
The calculated maximum cross correlation is the base of the analysis of phase shift, size, and shape
of the signals.
5.2.1 Maximum cross correlation
In general, the cross-correlation metric moves the test signal, T, by multiples of ∆t in relation to the CAE
signal, C. The cross-correlation value, ρ()m , is calculated at each shifted state. The time shift with the
maximum ρ is the base of the calculation of the cross-correlation rating.
Formula (12) shows the calculation of the cross correlation at each time shift, m. Curve C is moved by
multiples mm∈(, m ) of Δt between the minimum and the maximum shift [see Formulae (13) and
minmax
(14)]. The range of the time shift is limited by the parameter INT . Therefore, the signals T and C are
min
8 © ISO 2013 – All rights reserved

at least overlapping within the interval INTt⋅−()t . The parameter n in Formula (12) is not
min endstart
constant but is reduced to n [see Formula (15)].
min
n−1
Ct((++mi))⋅⋅ΔΔtT()ti+⋅ t
∑ startstart
i=0
ρ()m = −≤11ρ≤ (12)
n−1 n−1
2 2
Ct((++m iit))⋅⋅ΔΔTt()+⋅it
∑∑start start
i=0 i=0
()INTt−⋅1 ()−t
min endstart
m = 01< min min
Δt
()1−⋅INTt()−t
min endstart
m = 01< max min
Δt
nI=⋅NT n (15)
minmin
As described above, the cross correlation ρ is the maximum of all ρ()m [see Formula (16)].
ρρ=max{ ()m } (16)
5.2.2 Phase-shift rating
The phase-shift rating requires the parameters D and D to limit the phase shift. Both are defined
min max
within (,01). The thresholds of the phase shift are calculated by Formulae (17) and (18).
δ =⋅Dt()− minmin endstart min
δ =⋅Dt()− maxmax endstart max
The phase-shift rating is calculated by Formula (19) at the maximum cross correlation, ρ , and the
corresponding time shift, δ . The parameter k describes the decline of the rating between “1” and
Za2
“0”. A linear (k =1 ), quadratic (k = 2 ), cubical (k =3 ), or any other regression relationship
Za2 Za2 Za2
can be defined accordingly.
 1 if δδ
<
min

k
Za2
 

δδ−
 max
 
Zk= ∈Ν (19)

22b Za >0
 
δδ−
maxmin

 

0 if δδ>
max


5.2.3 Size rating
The size of the signals is analysed by comparing the area below the two curves after the phase shift.
It is a necessary evaluation but it may not be sufficient to evaluate the overall level of correlation. For
instance, the area below a signal with high and narrow peak could be identical to the area of a curve
with low but wide peak. The size method would evaluate this example with “1” although the shape of the
signals is completely different.
n
Ct()++δΔit⋅
∑ start
F
C i=1
= (20)
n
F
T
Tt()+⋅itΔ
∑ start
i=1
k
 Zb2
 
F
C

if FF>
 
TC
 F
 T 
Z = k ∈Ν (21)

22b Zb >0
k
Zb2
 
F
T
 

F
 C 

The analysis of size is done by comparing the square of the areas between the curves and the time axis.
Due to equidistant supporting points of both signals, the area can be represented by Formula (20).
Finally, the size rating is calculated by Formula (21). The meaning of k is similar to k (see 5.2.2).
Zb2 Za2
5.2.4 Shape rating
As shown in Formula (22), the shape rating, Z , is derived from the maximum cross correlation, ρ ,
2c
described in 5.2.1. The meaning of k is similar to k (see 5.2.2).
Zc2 Za2
k
Zc2
Zk=+()ρΝ1 ∈ (22)
22c Zc >0
5.3 Step-by-step procedure
The signals shall be pre-processed as described in Clause 10. After preparing the signals for the analysis
and defining the interval of evaluation, the maximum absolute amplitude, T , of the reference signal,
norm
T, shall be determined within this interval. It is used to calculate the inner and outer corridors. The
actual corridor and cross-correlation assessment can be executed within the defined interval. The
overall rating ranges between “0” and “1”. A score of “1” does not mean that both signals are identical.
Solely, their correlation is mathematically perfect within the defined tolerances.
To summarize, the following step-by-step procedures shall be followed to calculate CORA rating:
1) Pre-process both signals according to Clause 10.
2) Calculate T within this interval by using the reference signal.
norm
3) Calculate the inner and the outer corridors.
4) Calculate the corridor rating, Zt(), at every specific time, t, within the interval of evaluation.
5) Calculate the corridor rating, Z , based on Zt() and the number, N, of time samples.
1 1
6) Calculate the maximum cross correlation, ρ , between T and C.
7) Calculate phase rating, Z .
2a
8) Calculate size rating, Z .
2b
9) Calculate shape rating, Z .
2c
10) Calculate the cross-correlation rating, Z .
10 © ISO 2013 – All rights reserved

11) Calculate the overall CORA rating, Z .
CORA
6 EARTH metric
Another objective rating metric called Error Assessment of Response Time Histories (EARTH) was developed
[28][34]
for dynamic system applications. Figure 3 shows the structure of the EARTH error measures.
[28]
Figure 3 — EARTH error measures structure
The EARTH metric is divided into two categories: global response error and target point response
error. The global response error is defined as the error associated with the complete time history with
equal weight on each point. The three main components of the global response error are phase error,
magnitude error, and topology (or so-called slope) error. The target point error is defined as the error
associated with a certain localized phenomenon of interest, such as peak error and time-to-peak error.
The target point error represents the characteristic of a part of the time history, but does not indicate an
overall performance of the entire time history. In addition, the target point error is generally application-
dependent; hence, it is not the focus of this Technical Report.
Quantifying the errors associated with these features of phase, magnitude, and topology (slope)
separately is challenging because there are strong interactions among them. For example, to quantify
the error associated with magnitude, the presence of a phase difference between the time histories may
[25]
result in a misleading measurement. A unique feature dynamic time warping (DTW) of the EARTH
metric is used to separate the interaction of phase, magnitude, and topology (slope) errors. DTW is an
algorithm for measuring discrepancy between time histories. It aligns peaks and valleys as much as
[9]
possible by expanding and compressing the time axis according to a given cost (distance) function.
Figure 4 — Workflow of the revised EARTH metric
Since the ranges of the three errors are quite different and there is no single rating that can provide a
quantitative assessment alone, the initial EARTH metric study employs a linear regression method to
combine the three errors into one score. Numerical optimization method is employed to identify the
linear coefficients so that the resulting EARTH rating can match the SME’s ratings closely for a specific
[8][24]
application. However, the resulting linear combination of the EARTH metric is mainly numerical-
based and application-dependent; hence, it may not be scalable to other applications. In order to provide
more intuitive ratings while maintaining the advantages of the original EARTH metric, a revised EARTH
[28]
metric is proposed. Figure 4 shows the flowchart of the revised EARTH metric, and the details of the
algorithms are described in the following subclauses.
6.1 EARTH phase score
The EARTH phase score, E , is used to measure the phase lag between the two time histories. The
P
*
maximum allowable percentage of time shift is ε , and it is predefined. In this step, the initial curve, C,
P
is shifted left then right one step at a time to the original test data, and the cross correlation between the
12 © ISO 2013 – All rights reserved

truncated test curve and shifted and truncated C are calculated until reaching the maximum allowable
*
time-shift limits ε ⋅−()tt .
Pend start
When the initial curve, C, is moved to the left by m time steps, the number of overlap points of the two
time histories after time shift mt⋅Δ is reduced to Nm− and the corresponding cross-correlation value,
ρ ()m , is calculated by Formula (23).
L
nn−1
[(Ct((++mi))⋅−ΔΔtC()tT)(⋅+()ti⋅−tT(t))]
∑ startstart
i=0
ρ ()m = (23)
L
n−1 n−1
2 2
[(Ct ++()mi ⋅−ΔΔtC)(tT)] ⋅+[(ti⋅−tT)(t))]
∑ start ∑ start
i=0 i=0
When the initial curve, C, is moved to the right by m time steps, the number of overlap points after time
shift mt⋅Δ is reduced to Nm− and the corresponding cross-correlation value, ρ ()m , is calculated by
R
Formula (24).
nn−1
[(Ct()+⋅itΔΔ−⋅Ct())(Tt((++mi))⋅−tT(t))]
∑ startstart
i=0
ρ ()m = (24)
R
n−1 n−1
2 2
[(Ct +⋅itΔΔ)(−⋅Ct)] [(Tt ++()mi ⋅−tT)(t))]
start start
∑ ∑
i=0 i=0
The maximum cross correlation, ρ , is the maximum of all ρ ()m and ρ ()m . The number of the time-
E L R
shifting steps that yields the maximum cross correlation, ρ , is defined as the EARTH phase error, n .
E ε
ts
The corresponding shifted and truncated CAE curve C is recorded as C and the corresponding
ts
truncated test curve is recorded as T .
The EARTH phase score, E , is calculated by Formula (25). The best EARTH phase score is “1”, which
P
means there is no need to shift CAE curve to reach the maximum cross correlation between the initial
test and CAE curves. If the time shift, n , is equal to or greater than the maximum allowable time-shift
ε
*
threshold ε ⋅N , then the EARTH phase score is “0”. In between, the EARTH phase score is calculated by
P
a regression method; it is either linear (k = 1 ), quadratic (k = 2 ), or cubical (k = 2 ).
P P P
10if n =

ε

k
P
*
 
 ε ⋅−Nn
P ε
E =   k ∈ 12,,3 (25)
{}

P P
*
 
ε ⋅N

 P 

*
0 if nN≥⋅ε
ε P

6.2 EARTH magnitude score
The magnitude error is a measure of discrepancy in the amplitude of the two time histories. The
magnitude error is defined as the difference in amplitude of the two time histories when there is no time
lag between them. Before calculating the magnitude error, the difference between the time histories
caused by error in phase and topology (slope) are minimized by using dynamic time warping. The local
cost function, di(, j) , for DTW used considers both the distance and the slope as shown in Formula (26).
Once the local cost matrix is built, the algorithm finds the alignment path which runs through the low-
ts
cost areas on the cost matrix. This alignment path defines the corresponding elements of both Ci()
ts
and Tj() that will lead to the minimum accumulated cost function.
ts tss
   
dC ()i dT ()j
ts ts 22
di(, jC)(=−()iT ()jt)(+−t ) ⋅   −  (26)
ij
   
dt dt
   
tt= tt=
i j
ts
The derivatives at each data points of the truncated test curve, T , and the shifted and truncated CAE
ts ts+d
curve, C , are first calculated, which yields the truncated test slope curve, T , and the shifted and
ts+d ts
truncated CAE slope curve, C . Next, DTW is performed to the truncated test curve, T , and the
ts
shifted and truncated CAE curve, C , which results in the truncated and warped test slope curve,
ts+w ts+w
T , and the shifted, truncated, and warped CAE slope curve, C . The EARTH magnitude error,
ε , is calculated by Formula (27).
mag
ts++wtsw
CT−
ε = (27)
mag
ts+w
T
*
Formula (28) is used to calculate the EARTH magnitude score, E , where ε is the maximum allowable
M M
magnitude error, and k defines the order of the regression. The best EARTH magnitude score is “1”,
M
which means there is no difference in the amplitudes after phase shift and dynamic time warping. If the
EARTH magnitude error, ε , is equal to or greater than the maximum allowable magnitude error
mag
*
threshold, ε , then the EARTH magnitude score is “0”. In between, the EEARTH magnitude score is
M
calculated by regression method.
10if ε =

mag

k
M
*
 

εε−
 Mmag
 
Ek= ∈ 1,223, (28)
{}

M M
*
 
ε

M
 

*
0 if εε≥
 magM

6.3 EARTH slope score
The topological error is a measure of discrepancy in topology (slope) of the two time histories. The
topology of a time history is defined by the slope at each point. In order to ensure that the effect of global
ts ts
time shift is minimized, the slope is calculated from the time-shifted histories T and C . Thus, by
ts+d
taking the derivative at each point, the derivative time-shifted histories, represented by T and
ts+d ts+d ts+d
C , are obtained. Dynamic time warping is performed on T and C , which results in the
ts++dw
truncated and warped test slope curve, T , and the shifted, truncated, and warped CAE slope
ts++dw
curve, C . These two curves are then used to calculate the EARTH slope error, ε , by
slope
Formula (29).
ts++dw ts++dw
CT−
ε = (29)
slope
ts++dw
T
*
Formula (30) is used to calculate the EARTH slope score, E , where ε is the maximum allowable slope
S S
error, and k defines the order of the regression. The best EARTH slope score is “1”, which means there
S
is no difference between the two slope curves. If the slope error, ε , is equal to or greater than the
slope
14 © ISO 2013 – All rights reserved

*
maximum allowable slope error, ε , then the EARTH slope score is “0”. In between, the EARTH slope
S
score is calculated by regression method.
 10if ε =
slope

k
S
*
 

εε−
 Sslope
 
Ek= ∈ 12,,3 (30)
{{}

S S
*
 
ε

S
 

*
0 if εε≥
 slopeS

6.4 Overall EARTH score
The above three EARTH scores are combined into one EARTH score by using weighting factors shown
in Formulae (31) and (32).
Ew=⋅Ew+⋅Ew+⋅E (31)
EP PM MS S
ww++w = 1 (32)
PM S
The goal is to translate the initial three EARTH errors into a standardized score between “0” and “1”. Due
to the similarity with the cross-correlation rating in the CORA metric, the consistent weighting factors
with the CORA cross-correlation rating are used. In addition, the parameters of the EARTH metric, such
as thresholds for phase, magnitude, and topological errors, are defined by matching the EARTH rating
[28]
with the SME’s knowledge using a validation database and optimization technique.
6.5 Step-by-step procedure
The following step-by-step process shall be followed to calculate the EARTH rating:
1) Pre-process both signals according to Clause 10 (T and C).
2) Calculate the phase error in terms of time steps, n , by maximizing cross correlation.
ε
3) Calculate the phase score, E .
P
ts ts
4) Calculate the shifted and truncated time history curves, T
and C .
ts+d ts+d
5) Calculate the slope curves of the shifted and truncated time history curves, T and C .
6) Perform dynamic time warping to the shifted and truncated time history curves to generate the
ts+w ts+w
shifted, truncated, warped time history curves, T and C .
ts+w ts+w
7) Calculate the magnitude error, ε , between T and C .
mag
8) Calculate the magnitude score, E .
M
9) Perform dynamic time warping to the slope curves to generate the shifted, truncated, warped slope
ts++dw ts++dw
time history curves, T and C .
ts++dw ts++dw
10) Calculate the slope error, ε , between T and C .
slope
11) Calculate the slope score, E .
S
12) Calculate the overall EARTH rating, E .
E
7 Model reliability metric
[14]
A model reliability-based validation metric was developed for dynamic system applications. Figure 5
shows the illustration of this metric. The difference between CAE and test curves is taken as the
validation feature [see Formula (33)]. The threshold factor, b, is defined by the SME’s experience; the
lower and upper bounds of the threshold interval, εε, , are defined as the product of the threshold
 
 LU 
factor, b, and the absolute maximum amplitude, T , of the reference sig
...