Title

A competing Markov model for cracking prediction on civil
structures

Author(s)Kobayashi, K.; Kaito, K.; Lethanh, N.

CitationTransportation Research Part B: Methodological (2014), 68:
345-362

Issue Date

2014-10

URL

http://hdl.handle.net/2433/191023

Right© 2014 Elsevier Ltd.; This is not the published version. Pleasecite only the published version.; この論文は出版社版でありません。引用の際には出版社版をご確認ご利用ください。

Type

Journal Article

Textversion

author

Kyoto University

A Competing Markov Model for Cracking Prediction on Civil Structures
K. Kobayashia , K. Kaitob , N. Lethanh∗,ca Graduate School of Management, Kyoto University, Kyoto 606-8501, Japanb Department of Civil Engineering, Osaka University, Suita, Osaka 565-0871, Japanc Institute of Construction and Infrastructure Management, Swiss Federal Institute of Technology (ETH), Zurich 8093, Switzerland

Abstract
Cracks on the surface of civil structures (e.g. pavement sections, concrete structures) progress in several formations andunder different deterioration mechanisms. In monitoring practice, it is often that cracking type with its worst damagelevel is selected as a representative condition state, while other cracking types and their damage levels are neglected inrecords, remaining as hidden information. Therefore, the practice in monitoring has a potential to conceal with a bias
selection process, which possibly result in not optimal intervention strategies. In overcoming these problems, our paper
presents a non-homogeneous Markov hazard model, with competing hazard rates. Cracking condition states are classified
in three types (longitudinal crack, horizontal crack, and alligator crack), with three respective damage levels. The dynamic
selection of cracking condition states are undergone a competing process of cracking types and damage levels. We apply anumerical solution using Bayesian estimation and Markov Chain Monte Carlo method to solve the problem of high-orderintegration of complete likelihood function. An empirical study on a data-set of Japanese pavement system is presented
to demonstrate the applicability and contribution of the model.
Key words: Bayesian estimation, Markov Chain Monte Carlo, Gibbs sampling, Statistical deterioration prediction,
Markov chain model, Infrastructure management, Pavement cracking processes

1. Introduction
Modeling deterioration of civil infrastructures with statistical approach has been significantly documented in recent
years. One of the great advantages of applying statistical deterioration forecasting models is of its capability to incorporate
uncertainties. The uncertainty in deterioration process of civil infrastructure can be reduced by making use of historicalmonitoring data. Among the statistical modeling models, models with Markov chain have been widely applied bothin research (Tsuda et al., 2006; Robelin and Madanat, 2007; Kobayashi et al., 2012b) and in practice (Thompson et al.,
1998; AASHTO, 2004). In Markov deterioration forecasting models, the deterioration process of civil infrastructures is
described by the transition probability among discrete condition states, which are deducted values or composite values
from performance indicators of civil infrastructures (e.g. Pavement service indicator (PSI) is a composite condition state,
which is evaluated using performance indicators such as cracking, roughness, etc. (Shahin, 2005).
In Pavement management systems (PMSs), quality of a road surface is generally valuated by its riding quality and theskid resistance. The riding quality and skid resistance are reflected and quantified by evaluating pavement performanceindicators such as cracking, rutting, and longitudinal profile (Fukuhara et al., 1990). Measured values of those indicators
are monitoring data, which are stored in the data bank of the PMS. Monitoring data are used for the evaluation of pavement
condition, deterioration forecasting, and supporting decision makers to select optimal intervention strategies (OISs). If
the monitoring data appears to possess measurement errors or bias selection of data, it is likely that accurate deterioration
forecasting is flaw, as a results, decision makers would end up with intervention strategies (ISs), which are not optimaland considerable amount of benefits incurred to stakeholders (e.g. owner, users, and the public) reduced (Kobayashi et al.,2012b).
The deterioration of pavement is a complex process and therefore a single performance indicator would not be possible
to perfectly capture the process (Dore and Zubeck, 2009; Kenneth, 2010) (e.g. riding quality of pavement is measured
not only by cracking, but also by rutting, and longitudinal profile). Even when considering a single performance indicatoralone, its deterioration mechanism and progress could be complicated. One of the typical example that shows the complex
deterioration concerning one performance indicator is the mechanism and development of cracks.
Cracking is one of the popular deterioration phenomena in the PMSs, especially in the cold regions (Dore and Zubeck,2009). Cracks appear on the road surface in different patterns and directions. Cracks are progressed and observed inhorizontal direction, longitudinal direction, or the combination of the twos. Reasons for crack initiation can be numerous∗ Corresponding author. Tel.: +41 446 333 083; Fax +41 762 072 105.

Email addresses: kobayashi.kiyoshi.6n@kyoto-u.ac.jp (K. Kobayashi), kaito@ga.eng.osaka-u.ac.jp (K. Kaito),
lethanh@ibi.baug.ethz.ch (N. Lethanh )Preprint submitted to Transportation Research Part B

May 14, 2013

(e.g. rain water infiltrated directly into pavement structures and causes the reduction in the cohesion of materials, orcontraction and expansion of materials due to high fluctuation of ambient temperature). If a small crack initiated on aparticular location of the pavement section, the overall deterioration process of a road section would significantly advancedwith a greater speed.
In practices, when percentage of cracks on a pavement section reaches to a certain alarming level (level at which riding
quality, skid resistance, and greater negative impact is anticipated, or the alarming level specified as technical standard andrequirement for interventions, etc.), interventions (e.g. crack sealing and crack filling) should be implemented to prevent
cracking from further advancement and also to ensure the riding quality in acceptable standards for users. However, as
earlier mentioned, cracking is a complex one. It progresses in various directions and possibily with different levels ofdamages (Kaito et al., 2007). Mechanism of horizontal cracking could be different from that of longitudinal cracking.
Moreover, treatments for different cracking types are not always the sames. Therefore, it is important to understand thedeterioration process of each crack types. This paper focuses on deterioration prediction of cracking, which is a good
representation of complex pavement performance indicators, with multiple types.
Thanks to the development of high-tech inspection and measurement devices for civil infrastructures in the past
decades, especially the application of high definition cameras and image processing techniques (Kawashima et al., 1984;
Fukuhara et al., 1990; Mohajeri and Manning, 1991; Wang and Smadi, 2011), a small crack, even with 1 mm width, canbe detected (Fukuhara et al., 1990). Basically, cracks are classified into three types: horizontal crack (or transverse crack),
longitudinal crack, and alligator crack (Rababaah et al., 2005; Kaito et al., 2007; Nakat and Madanat, 2008). Despite thefact that there are different types of cracks that can be recorded and stored in the data bank of the PMS, the use of crack isoften referred and represented to its general name as “crack” without much attention of the actual type of crack itself. Inmany practical cases, the percentage of crack on a road section is the representation of the percentage of a type of crack
at its worst deterioration condition. The representative crack and its corresponding deterioration condition is used directly for deterioration prediction using Markov model with single dimension of condition states (e.g. a single dimension
Markov model is formulated using only a vector of condition states).
There is a problem in practice that selection of representative crack to be used in Markov models is, in many cases,
driven by the bias or default assumption on which cracking types and their corresponding levels of deterioration. For
example, during the courses of inspections, engineers often select the cracking type with its worst deterioration condition
among other cracking types to be the representative of crack for the examined road section. Under this situation, there isa bias in the selection of crack for use in the Markov deterioration model with single dimension of condition states. The
practice of selecting the typical cracking type with its worst deterioration condition among other cracking types is not only
considered as bias selection process, but also a competitive selection process as well. For example, in the first inspection
time, horizontal crack is selected to be representative crack of the investigated road section. However, in the second
inspection time on the same road section (or other inspection on different road section), longitudinal crack is selected tobe the representative one.
To overcome the limitation of using Markov model with single dimension of condition states, Kaito et al. (2007)
developed a novel hierarchical Markov model (hereafter referred as HIMA model) and successfully tested it on a set ofcracking data in Japan. To the best of authors knowledge, the cited model is the first one in the field addressing the multiple
dimensions of deterioration and multiple condition states with Markov chain model. The HIMA model was developed asa novel extension of the multi-stage exponential Markov model (MUSTEM), which was previously developed also by thesame research group (Tsuda et al., 2006; Lethanh, 2009; Kobayashi et al., 2012a) at Kyoto University. The fundamental
difference between the HIMA model and the MUSTEM model is their dimensions of deterioration types and condition
states. In the MUSTEM model, only single dimension of condition stages is used. Whilst, in the HIMA model, the authors
proposed to use three dimensions of crack initiation with multiple discrete condition states.
In the HIMA model, cracks are classified into three types: horizontal crack, longitudinal crack, and alligator crack.
The deterioration of each cracking type is expressed as Markov transition probabilities (m.t.p) among its damage level.
Based on cracking types and their corresponding damage levels, a multi-dimension discrete condition states are defined(Table 1). The authors have successfully conducted an empirical study of the model on a set of cracking data of Japanesenational roads to prove the applicability and usefulness of the model. Moreover, the authors concluded that the model ismore effective under the condition of having sufficient amount of data concerning each cracking type and its damage level.
However, the HIMA model has a default assumption that no selection bias and measurement errors exist in monitoring
activities concerning the measurement of cracking types and their damage levels.
The present paper proposes a new deterioration forecasting model in the streamline of Markov deterioration forecasting models for infrastructure asset management. The model is a novel extension of the HIMA model, which expresses thecracking condition states on pavement structural surface in multiple directions. We discuss the bias in selecting cracking
condition states as a competitive process, which was not yet developed with the HIMA model, and propose a numericalsolution using Bayesian estimation and Markov Chain Monte Carlo (MCMC) simulation. To test our model, we apply anempirical study on a set of cracking data of national road system in Japan.
Following section gives a brief review of the HIMA model, which discussed cracking process and selection bias inmonitoring and management practices in the PMS. Section 3 summaries important formulation of the MUSTEM model,
which is considered as the core of the HIMA model and the proposed model in this paper. Section 4 includes the math-2

Table 1: Condition states concerning the cracking modes.
Damage levels
Cracking types (l)Longitudinal crack Horizontal crack
Alligator crack
i=0(0, 0)
i=1(1, 1)(1, 2)(1, 3)
i=2(2, 1)(2, 2)(2, 3)(3, 1)(3, 2)(3, 3)
i=3

ematical formulation of the new model. Numerical solution as estimation method for the proposed model is detailed in
Section 5. Section 6 includes an empirical example. The last section concludes the paper and gives recommendation forpractical use of the model and future study.
2. Cracking Process and Selection Bias
2.1. Cracking and deterioration condition states
Cracks appears due to the combination of various aspects and therefore being considered as a complex process. There
always exists a certain level of uncertainty in forecasting cracking condition states. In view of management, ambiguity incapturing the true cracking types and the damage levels could lead to non optimal intervention strategies. For example,
cracks are observed in longitudinal direction, horizontal direction, and in alligator crack. Depending on cracking types
and actual damage levels, interventions activities are then designated. However, interventions are arranged on regularbasis, and thus interventions cannot be spontaneously applied at right time and right place to heal all cracking types. In
addition, an interventions activity, designated for mending horizontal cracks, do not genuinely heal longitudinal cracks.
As a result, expected values of longitudinal cracks could be neglected. Thus, it is true to state that using a single cracking
type for making decision is not an optimal choice.
In the paper of Kaito et al. (2007), the authors developed the HIMA model to address the multiple progression of
cracking. In the cited paper, cracking process of road surface is described by means of condition states with two variables,
damage level i and cracking type l. Cracking type is expressed by variable l (l = 0, · · · , L), with l = 0 indicating no sign ofcrack. Cracking type l (l = 1, · · · , L) can be of longitudinal crack, horizontal crack, and alligator crack. Damage level isdefined as state variable i (i = 0, 1, · · · , I), representing deterioration level of respective cracking type. Damage level i = 0
when there is no sign of deterioration and damage level I is referred as “absorbing state”, which demands an immediate
intervention. Under this assumption, overall deterioration of crack can be relatively described by a pair of state variables,
which is defined as “pair condition state” (i, l). Condition state (0, 0) and (I, L) respectively refer to no sign of crack and“absorbing condition state” of a road section. Description of cracking types and damage levels are given in Table 1.
2.2. Selection bias and representation matter
In many PMSs, an entire road length is often divided into smaller sections (sk , k = 1, · · · , K) or divisions for conveniences of monitoring and management (Fig. 1). Within a road section, cracks might occur in several types (l), withdifferent damage levels (i). In practice, inspectors often decide to use the worst damage level, corresponding to a specific cracking type, to be the representative condition state. However, there is no guarantee which condition state is theworst. This practice exposes as a shortcoming that true information of other cracking types is omitted, remaining as
hidden information. For example, as shown in Table 1, the progress of crack on a road section is described by 3 condition
states: longitudinal crack (2, 1), horizontal crack (2, 2), and alligator crack (1, 3)). It is not certain at a point in time toverify which condition state is the best representative among these three. However, it is often the case that inspectors
select condition state (2, 2) (horizontal crack) to be the representative value for deterioration of cracks of that road section.
Given that reality, true information of other cracking types (longitudinal crack (2, 1) and alligator crack (1, 3)) is neglected
and not recorded in the database. This problem could result in non-optimal estimation results of deterioration and also
influence the true likelihood of the OISs.
The process to select representative condition states as mentioned is viewed as a bias selection process. The selection
bias is a type of systematic measurement errors, which deems unsolvable by using conventional errors elimination techniques (Kobayashi et al., 2012b). The bias selection process in the practice of inspections therefore targets only a typicalcracking type and its damage level. However, in our paper, we develop a model to capture the cracking progress of not
only typical crack but also other cracking types. The missing information of cracking types, which were not selected andrecorded, and their corresponding damage levels, is addressed. In the paper, condition states concerning missing cracking
types and their damage levels are regarded as “latent” or “hidden” variables.
The selection bias in defining a representative cracking type can be interpreted from Fig. 1. As shown in the figure, aroad section can be divided into smaller rectangular areas sk = (s1 , · · · , sK ). Measurement of cracks at any inspection time,in general, is carried out on each rectangular area sk . The results of different cracking types and their damage level are thenweighted (or aggregated) to become the means, which represent the overall cracking condition states of the road section.
3

Measured areas sk in a road section
s1

s2

s3

s4

Longitudinal
crack

Horizontal
crack

sk

sK

β

α

β'

Figure 1: Measured areas in a road section and the value distribution of cracking types

Conditionstatesinspection

inspection

Figure 2: Deterioration process

We denote the mean values of damage levels of longitudinal crack and horizontal crack as α and β, respectively. The value
distributions of longitudinal crack and horizontal crack around their means can be interpreted also from the figure. The
solid distribution curve represents the value distribution of longitudinal crack, and the two dotted lines represent value
distribution of horizontal crack at a inspection time. Due to uncertainty, at any inspection time, the means of distributionof horizontal crack β, can be either greater or smaller than the mean of longitudinal crack α.
The selection of damage level and typical cracking type is assumed as follows:• at any inspection time, if damage level of horizontal crack follows the pattern of solid distribution curve, which hasmean value of α, then α is selected as the representative damage level.
• On the other hand, if damage level of horizontal crack progresses as like in the dotted distribution curve, which hasmean value of β, then β is selected as the representative damage level.
As a matter of fact, the development of horizontal crack and longitudinal crack varies differently. Hence, the mean valueof horizontal crack (dotted distribution curves) is fluctuated, either with value of β or β′ . This type of selecting therepresentation of damage level is considered as a competing process among damage levels and types of cracking.
In order to address the competitive crack selection process, it is important to develop a new methodology that takes
into account the selection bias. The new methodology should be able to estimate damage levels of hidden cracking types,
especially under the situations of having only partial or incomplete monitoring data, which is quite popular in many PMSs.
3. The MUSTEM model
In this section, we summary the research of Tsuda et al. (2006); Kobayashi et al. (2012a), which proposes a modelingapproach to estimate the Markov transition probability (m.t.p) based on historical monitoring data. Careful readers arerecommended to refer to the original paper for greater details of the methodology.
It is assumed that the deterioration of a road section follows the path shown in Fig. 2. In the figure, time τ represents
the actual time on a calendar (hereafter referred as “time”). Deterioration of a road section can be represented by discretecondition state i(i = 1, .., J), with i = 1 as initial condition state (when structure is new) and i = J as absorbing condition
4

state. Time τA and τB are inspection times, while time τi is any arbitrary time in between. Duration between two inspection
times is Z. Given monitoring data of two inspection times τA and τB , the m.t.p is described as follows:
Prob[h(τB ) = j | h(τA ) = i] = pi j

(1)

Suppose that condition state changes from i to i + 1 at time τi (timing yC ). At that time, the duration of condition state
i can be expressed by the following equation: ζi = τi − τi−1 = yC . Assume that the duration ζi of condition state i is arandom variable, and is subjected to the probability density function fi (ζi ) and the distribution function Fi (ζi ). Here, thedomain of the duration ζi is [0, ∞). The following expression is defined from the definition of distribution function:Z yifi (ζi )dζi(2)
Fi (yi ) =0

The distribution function Fi (yi ) represents the cumulative probability of the change of the condition state from i to i + 1
in the period from the initial timing yi = 0 (time τi−1 ), at which the condition state has become i, to the timing yi (timeτi−1 + yi ). Accordingly, the probability F̃i (yi ) of remaining at condition state i from the initial timing yi = 0 to the sample
timing yi ∈ [0, ∞) can be expressed by the following equation, using the cumulative probability of the change of thecondition state from i to i + 1 until timing yi :Prob{ζi ≥ yi } = F̃i (yi ) = 1 − Fi (yi )(3)The conditional probability of the event that the road section remains in condition state i until timing yi and change tocondition state i + 1 in the period [yi , yi + ∆yi ) is defined as:
λi (yi )∆yi =fi (yi )∆yiF̃i (yi )(4)

The instantaneous rate λi (yi ) of the change in the condition state from i to i + 1 at timing yi is called a hazard function. By
using a hazard function suited for the assumed deterioration process, it is possible to describe the evolution of deterioration
over time.
Under the assumption that the Markov characteristics concerning the deterioration processes of the road section do
not depend on the history of deterioration and the hazard function is constant θi > 0, in another words, the hazard function
is independent of the timing yi , following equation is defined:
λi (yi ) = θi

(5)Using the hazard function λi (yi ) = θi , the probability of the even that condition state i remains over a duration yi is further
described as:" Z yi#F̃i (yi ) = exp −(6)
λi (u)du = exp(−θi yi )0

The survival probability function is identical to the transition probability pii when the duration yi equals to interval z. By
defining the subsequent conditional probability of condition state j to i, with respect to z, a general mathematical formulafor estimating the m.t.p pi j is formulated:
pi j (z) = Prob[h(τB ) = j|h(τA ) = i] =

j−1
θm
θm Y
exp(−θk z)θ − θk m=k θm+1 − θkk=i m=i m

j Y
k−1
X

where there are the following conditions: Q
k−1 θm
=1
at (k ≤ i + 1)
 m=i θm −θ
k
Q

j−1
θm


m=k θm+1 −θk = 1 at (k ≥ j)

(7)

(8)

In addition, with regard to the m.t.p from any condition state to absorbing condition state piJ , following equation is used:piJ (z) = 1 −J−1Xpi j (z) (i = 1, · · · , J − 1)

(9)

j=i

4. Formulation of the competing Markov hazard (COHA) model
4.1. Competing Markov transition probabilities
Cracks could appear on the surface of a road section after the road section is newly paved or intervened. In order touphold a certain service level and maintain riding quality, ISs need to be evaluated for future management. It is assumed

5

that immediately after an intervention, surface condition of the road section is renewed, with condition state (i, l) = (0, 0).
To visualize the deterioration of cracks, we define τ0 , with y = 0, as starting point immediately after an intervention.
Initial condition state is updated as (0, 0) at time τ0 . In next inspection time, also the time of next intervention (if required),passing duration is assumed as z, with y + z as time indication of when measurement being counted. The uncertainty ofdeterioration progress with respect to cracking is described by means of transition probability among condition states inthat period.
Condition states of a road surface in duration y is defined as state variable h(y) = (i, l) (i = 0, · · · , I; l = 0, · · · , L). Afterany arbitrary time duration y + z, state variable becomes h(y + z) = ( j, m) ( j = 0, · · · , I; m = 0, · · · , L). Thus, following
transition probability is defined:
Prob[h(y + z) = ( j, m)|h(y) = (i, l)] =Prob[h(y + z) = ( j, m), h(y) = (i, l)]= πil, jm (y, z).
Prob[h(y) = (i, l)](10)

The m.t.p πil, jm (y, z) is estimated based on the condition states of road sections, which are recorded as a result of thecompetitive selection among damage levels and cracking types. Therefore, we refer the m.t.p πil, jm (y, z) as competing
m.t.p. It is noted that the m.t.p of the MUSTEM model depends on transition duration z and time y. However, thecompeting m.t.p of the COHA is defined based only on the information of typical cracking type and damage level. It
should be also noted that the competing m.t.p might not satisfy the condition of homogeneous m.t.p with respect tohidden cracking type as earlier mentioned in section 2. Only if information of all condition states is available, had thefollowing competing m.t.p assured.
 π00 (y, z) · · · π0I (y, z) 


..
..
..
 .
(11)
Π(y, z) = 
.
.
.


0
· · · πII (y, z)
Where 0 is a block procession, and πi j (y, z) (i, j = 0, · · · , I) satisfies following conditions:
π00 (y, z) = π00,00 (y, z),
π0 j (y, z) = ( π00, j1 (y, z) · · · π00, jL (y, z) ), πi1, j1 (y, z) · · · πi1, jL (y, z) 


..
..
..
 .
πi j (y, z) = 
.
.
.


πiL, j1 (y, z) · · · πiL, jL (y, z)(12)Eq. (11) shows the competing m.t.p between two consecutive inspection times t = y and t + 1 = y + z. Probabilistically,
properties of the competing m.t.p vary upon different inspection period. If an intervention has not been applied, cracking
happens in natural process. Hence, condition πil, jm (y, z) = 0 (i > j) must be satisfied. Moreover, summation of all
P PL
πil, jm (y, z) = 1. Following conditions
properties in a row of transition probability matrix must be equal to 1 as Ij=i m=0
are obligated as rules of the competing m.t.p, with πIm,Im (y, z) = 1 for absorbing state of damage level I.

πil, jm (y, z) ≥ 0πil, jm (y, z) = 0 (i > j)
.
(13)


PI PL

π (y, z) = 1 j=iil, jmm=0

The competing m.t.p describes the transition among condition states of cracks based on the selection of typical cracking
type and damage level. The selection and estimation of competing m.t.p for each road section is mutually independent
from other sections.
In order to obtain the explicit mathematical formula for the competing m.t.p, we extend the formula in Eq. (7) for thecase of cracking. We assume l∗ from the range l (l = 1, · · · , L) as an example of cracking type. Transition concerning
damage level i is described as (0, 0) → (1, l∗ ) → · · · → (I, l∗ ). Hazard rate θil∗ is defined as transition from condition state(i, l∗ ) to (i + 1, l∗ ). For transition from initial condition state (0, 0) to (1, l∗ ), hazard rate is expressed as θ0l∗ . Remaining
duration in condition state (i, l∗ ) (i = 0, · · · , I −1), hereafter referred as “life expectancy” of condition state (i, l∗ ), is definedby means of a stochastic variable ζil∗ , with its probability density function fil∗ (ζil∗ ) and distribution function Fil∗ (ζil∗ ). In
the period (t,t + 1), the probability density that condition state (i, l∗ ) disappears at the end of the period is expressed byhazard function λil∗ (zil∗ ), with elapsed time zil∗ . As a result, the competing m.t.p is defined as:∗plii (z) = Prob[h(τB ) = (i, l∗ )|h(τA ) = (i, l∗ )] = exp(−θil∗ z).
∗

pli j (z) = Prob[h(τB ) = ( j, l∗ )|h(τA ) = (i, l∗ )] =

6

j Y
k−1
X

j−1
Y
θml∗
θml∗
exp(−θkl∗ z)
∗
∗
−
θ
θ
θ
kl m=k (m+1)l∗ − θkl∗
k=i m=i ml(14-a)(14-b)

Eq. (14-b) satisfies following conditions: Qk−1 θ ∗ml(k ≤ i + 1)
 m=i θml∗ −θkl∗ = 1
Q j−1
.


 m=k θ θml∗ ∗−θ ∗ = 1 (k ≥ j)(m+1)lkl

Using hazard rate θil∗ , we define RMD as an important management indicator (RMDi,l∗ stands for Expected Remaining
Duration of damage i of cracking type l∗ ). The indicator reflects the duration, in which damage level i survives, given a
condition that it has been observed in previous inspection time. Formula of RMDi,l∗ is actually analogous to the formulaof survival function F̃il∗ in infinite domain (Lancaster, 1990):
Z ∞
∗
(15)
F̃il (z|θil∗ )dz = θil−1∗ .
RMDil =
0

Expected life expectancy (denoted as ET jl∗ ) of damage level j (> 1), if considering all damage levels i(i < j) of cracking
type l∗ , is thus a summation of all transition duration from every damage level i.
ET jl∗ =jX

θil−1∗ .

(16)

i=1

4.2. The COHA model.
In previous section, l∗ is an example of cracking type. In this section, we consider all cracking types in the rangen (n = 1, · · · , L). The cracking process of an individual cracking type is independent from other cracking types. Selection
method for representative condition state follows two rules. The first rule applies to select the worst damage level among
cracking types. The second rule applies when road section exposes with identical damage level for all cracking types andcracking type with higher index value is selected. Rules to decide condition state (i, l) can be expressible as:(
i = max{in (n = 1, · · · , L)}.
(17)l = max{n|in = i (n = 1, · · · , L)}
Eq. (17) expresses recorded data, which is considered as a result of bias selection process. In addition, we denote in (n , l)as damage level of hidden cracking type n.
(in ≤ i (n = 1, · · · , l − 1).
(18)
in < i (n = l + 1, · · · , L)
Hazard rate θ00 , which infers the change of damage level i from 0 to 1, is defined in the whole range of cracking type
(1, n) (n = 1, · · · , L):
θ00 =

L
X

(19)

θ0n .

n=1

Condition state variables at time t = τA and t = τB are h(τA ) = (i, l) and h(τB ) = ( j, m) respectively. The transition of
condition states from (i, l) to ( j, m) in duration z = τB − τA is expressed as the m.t.p in following subsections:4.2.1. when i = j = 0
In the case that condition state change from (0, 0) to (1, n) (n = 1, · · · , L), hazard rate is defined in Eq. (19). The m.t.p
π00,00 (y, z) over duration z is then defined:
 L

 X

π00,00 (y, z) = exp −θ0,n z .
(20)
n=1

The competing m.t.p π00,00 (y, z) does not depend on time t, but on duration z.

4.2.2. when j , 0 and l , m
The competing m.t.p is defined similarly as in Eq. (10).
πil, jm (y, z) =Prob[h(y + z) = ( j, m), h(y) = (i, l)].
Prob[h(y) = (i, l)](21)

Formulation of the nominator in Eq. (21) can be described in three different scenarios as: 1) cracking type l, observed attime τA , is not selected as representative one at time τB ; 2) cracking type m, observed at time τB , is not appeared at timeτA ; 3) both scenarios 1 and 2 occurs simultaneously.

7

In scenario 1, at time τA , the probability of damage level i with cracking type l is pl0i (y). To satisfy the assumption thatat time τB , damage level jl is not chosen, following condition must be satisfied:(i ≤ jl ≤ j when l < m.
(22)i ≤ jl < j when l > m
Under this assumption, probability for scenario 1 to happen is defined as:
P1 (y, z) = pl0i (y)j̃lXplit (z).

(23)

t=i

Probability plit (z) is considered in the entire range of damage levels from i to t in duration z, with its detailed properties
presented in Eqs. (14-a) and (14-b). The affixing character j̃l has following meaning:
(
j
l<m
.
(24)
j̃l =
j−1 l>mIn scenario 2, at time τA , following condition must be satisfied to hold the assumption m , l, meaning that cracking typem is not chosen at time τA .
(im ≤ i when l > m.
(25)
im < i when l < mThus, the simultaneous occurrence probability for cracking type m in damage level j at time τB is formulated as:
P2 (y, z) =

ĩm
X

m
pm
0s (y)p s j (z).

(26)

s=0

In scenario 3, the damage level of cracking type n should satisfy((
jn ≤ j when m > nin ≤ i when l > n,, andjn < j when m < nin < i when l < n(27)

and following occurrence probability for scenario 3 is obtained:
P3 (y, z) =

L
Y

j̃n
ĩn X
X

pn0s (y)pnst (z).

(28)

n=1,,l,,m s=0 t=s

Q
Product sign n=1,,l,,m in Eq. (28) concerns cracking type n other than cracking types l and m. Therefore, the nominatorof Eq. (21) is further expressed as




j̃l
j̃n






ĩm
ĩn
L






X
X











 Y XX nmmlnl. (29)
p0s (y)p s j (z)p0i (y)p0s (y)p st (z)pit (z)
Prob[h(y + z) = ( j, m), h(y) = (i, l)] = 
















  s=0n=1,,l,,m s=0 t=s
t=i

The denominator of Eq. (21) concerns only the probability of condition state (i, l) at time τA . If transition from initial time
is considered, following probability is defined:
Prob[h(y) = (i, l)] = pl0i (y)

ĩn
L X
Y

pn0s (y).

(30)

n=1,,l s=0

Finally, the explicit form of competing m.t.p πil, jm (y, z) is obtained:

−1



j̃l
j̃n








ĩn
ĩm
ĩn
L
L











X l 
X m
Y X n  Y XX n
n
m
.
p
(y)p
(z)
p
(y)
p
(z)
p
(y)p
(z)
πil, jm (y, z) =








st
it
s
j
0s
0s
0s













 t=i
 s=0



n=1,,l s=0n=1,,l,,m s=0 t=s(31)

It is noted that the competing m.t.p depends not only on z but also y. As earlier mentioned in section 3, the m.t.p of the
MUSTEM model depends only on z (Eq. (7)).
4.2.3. when j , 0, l = m
This assumption occurs under two scenarios: 1) cracking type l is in damage levels i and j respectively at time τA andτB ; 2) at time τA , damage level of cracking type l is i. However, in time τB , damage level j is observed simultaneously
for cracking types l and m. To differentiate the description with the scenarios in previous section, the sign [′ ] is used infollowing paragraphs.
8In scenario 1′ , we apply Eqs. (14-a) and (14-b) for the probability of cracking type l, with its damage level i and j attime τA and τB respectively.
P1′ (y, z) = pl0 j (y)pli j (z).

(32)

In scenario 2′ , conditions in Eqs. (22) and (25) must be satisfied, and thus, the probability of scenario 2′ is defined:
P2′ (y, z) =

j̃n
ĩn X
L X
Y

pn0s (y)pnst (z).

(33)

n=1,,l s=0 t=s

Probability Prob[h(y) = (i, l)], which happens at time τA = y is also defined asProb[h(y) = (i, l)] = pl0 j (y)

ĩn
L X
Y

p0s (y).

(34)

n=1,,l s=0

To this point, an explicit mathematical form for competing m.t.p is specified:

−1


j̃n




ĩn
ĩn
L
L







Y X n  l Y XX n
n
.
p
(z)
p
(y)
p
(y)p
(z)
πil, jl (y, z) =




st
i
j
0s
0s








n=1,,l s=0n=1,,l s=0 t=s(35)

the competing m.t.p πil, jm (y, z) is expressible through the m.t.p pli j (z).
To this point, it is our target to estimate the model’s parameters, which are embedded in the likelihood function of them.t.p in Eq. (35), based on inspection data. One way to estimate the model’s parameters is using the maximum likelihood
estimation (MLE) method. However, It is likely impossible to use the MLE method as the likelihood function involves ahigh-order of integration, which limits the derivation for the Jacobian matrix (first order derivative) and Hessian matrix(second order derivative). To overcome this problem, in section 5, we propose a methodology using frontier stochastic
approach with Bayesian estimation and MCMC simulation.
4.3. Prediction of non-homogeneous Markov distribution.
The competing m.t.p has Markov property in its form. However, the competing m.t.p is not purely homogeneous
Markov process. In fact, it is regarded as non-homogeneous Markov process. To describe this fact, we express thecompeting m.t.p in three respective initial times yA , yB , and yC , with durations τA , τB , and τC respectively. Time interval
is denoted as z. If competing m.t.p satisfies homogenous Markov process, following equation can be derived:
Prob[h(yC ) = ( j, m)|h(yA ) = (i, l)]i
L XYProb[h(yC ) = ( j, m)|h(yB ) = (k, n)] · Prob[h(yB ) = (k, n)|h(yA ) = (i, l)].
=(36)n=1,,l s=0

On the other hand, if considering two time intervals 2z in Markov process, following equation cannot be derived for thecase of two time intervals 2z. In Eqs. (31) and (35), there is an evidence that competing m.t.p depends not only z but also
on y.
πil, jm (yA , 2z) =

L
I X
X

πil,kn (yA , z)πkn, jm (yB , z).

(37)

k=1 n=1

For any arbitrary time interval z, another explicit mathematical form should be developed in line with Eqs. (31) and (35).
We define probability Pil (y) for the event “condition state (i, l) (i = 0, · · · , I; l = 0, · · · , L) occurs at a certain time withelapsed time y”.




ĩl




 l
Y X np0i (y).
(38)p(y)
Pil (y) = 

0in





 n,l in =0

In addition, the probability that damage level i of any cracking type occurs in elapsed time y after any intervention isdefined as maximum damage level probability Pi (y).
P∗i (y) =LX

Pil (y).

(39)

l=1

9

Life expectancy ET ( j) of damage level j after any intervention is formulated as:
Z ∞ X
j−1
ET ( j) =yP∗i (y)dy.
0

(40)

i=0

5. Estimation Method
5.1. MCMC Method.
In statistic with Bayesian inference, prior and posterior probability are employed with aim to estimate the valuesof model’s parameters. However, in actual analysis, it is hard to define a prior probability distribution, even with asimple condition state hazard model (Ibrahim et al., 2001). Methods to overcome the problems in the assumption ofprior probability distribution often require numerical analysis with multi-dimensional integration, and thus remaining asa limitation in Bayesian estimation.
In recent years, an appealing solution to the problems in Bayesian estimation has been proposed, with the applicationof MCMC simulation. The MCMC simulation technique does not require a high level of derivative and multi-dimensional
integration of model’s objective functions (Robert, 1996). As a result, estimation results, in a great number of applied
statistic research, have been improved through the combination of Bayesian estimation and MCMC simulation.
In MCMC simulation, Gibbs sampling and Metropolis Hastings (Metropolis-Hastings or MH) techniques have been
extensively discussed (Robert, 1996). Reference to research on image restoration is a good example of MCMC simulation(Geman and Geman, 1984). Of that study, the algorithm of Gibbs sampling was used to estimate the posterior distribution
in Bayesian estimation. In MH law, the iterative parameter β is defined by repeatedly generating random numbers throughthe conditional probability density function.
Regarding application of Bayesian estimation and MCMC method in infrastructure management, the authors of thispaper has developed a hidden Markov model for elimination of selection bias Kobayashi et al. (2012b). The use ofBayesian estimation and MCMC method has showed a great advantage over the conventional MLE approach in the caseof having complete likelihood function with multiple integrations. Following sections detail our numerical solution to
overcome the challanges in estimation of the COHA model’s parameters.
5.2. Formulation
Two visual inspections are conducted for each section k of the entire road system (with K is the total number of roadsections). Condition states being observed at two inspection time τkA = yk and τkB = yk + zk are h(yk ) and h(yk + zk )respectively. Duration zk is time interval between two inspections times. Based on actual observed condition states, adummy variable δkil, jm is assumed, with following characteristics:(
1h(yk ) = (i, l) and h(yk + zk ) = ( j, m)kδil, jm =.
(41)
0 otherwise
To further describe the information of sample k, a characteristic vector xk = (x1k , · · · , xkM ), representing characteristickvariables influencing on cracking progress is also defined. For instance, characteristic variable xm(m = 1, · · · , M) is
referred as a characteristic variable with its numbering index m. Overall information on sample k can be summarized invector ξk = (δk , yk , zk , xk ), with δk = δkil, jm is a vector of dummy variable. We pay attention to damage level i and cracking
type l. The formulation of hazard rate θilk of sample k, with its characteristic variable xk = (x1k , · · · , xkM ), is then defined as:
θilk = exp(xk β′il ),(42)

with βil = (β1il , · · · , βilM ) is the line vector of unknown parameters βm
il (m = 1, · · · , M). The sign ′ indicates the vectortransposition. The m.t.p pli j of cracking type l is expressible as a function of zk , xk , and βil as pli j (zk , xk : βil ). As mentioned
in previous section, competing m.t.p πil, jm comprises of two components pli j (yk , xk : βil ) and pnij (zk , xk : βil ) (n , l). Thus,
it can also be described as πil, jm (yk , zk , xk : β), with β = (β01 , · · · , βil , · · · , βI−1L ).
5.3. Complete Likelihood Function
In the COHA model, competing m.t.p πil, jm is defined by means of condition states of both observed typical cracking
type and hidden cracking type. Hence, the simultaneous occurrence probability (likelihood) L(β : ξ̄) at the time of1
K
observing the sampling information ξ̄ = (ξ̄ , · · · , ξ̄ ) in entire K road sections can be defined as
L(β : ξ̄) =

K n
L Y
I Y
L Y
I−1 Y
Y
oδ̄k
πil, jm (ȳk , z̄k , x̄k : β) il, jm
i=0 l=0 j=i m=0 k=1

10




δ̄kil, jm −1   Y

ĩn
L 
L
L 
I Y
L X
L Y
K 
I−1 Y
Y



Y
Y





 Y

 

n
k 
p
(ȳ
)
=



 
0s









  
i=0 l=0 j=i k=1 m=0

n=1,,l s=0

m=0,,lδ̄kil, jm

pn0s (ȳk )pnst (z̄k )
n=1,,l,,m s=0 t=s
j̃n
ĩn X
X

δ̄kil, jlδ̄kil, jm  δ̄kil, jm   Y
j̃n
j̃l





ĩn X
ĩm
L X









k
X l k   X m k m k 
k n k 
n
(ȳ
)p
(z̄
)
p
{pli j (z̄k )}δ̄il, jl ,
(z̄
)
(ȳ
)p
(z̄
)
p
p







st

s
j
it
0s
0s














 t=i
  n=1,,l s=0 t=s

 s=0

(43)

where pli j (z̄k ) = pli j (z̄k , x̄k : βil ). Estimating a set of unknown parameters β to maximize the likelihood of function (43) is
our objective.
In freqenties statistics, the model’s parameters are often estimated by employing regression methods. The method of
Maximum Likelhihood Estimation (MLE) is among the most popular used in frequentist statistics. However, the MLE
method has its limitation when the objective function (or likelihood function) of a statistical model is in the form of highorder dimension and multiple integrations (Andrew et al., 2006; Dani and Hedibert, 2006; Jeff, 2006). This limitation
appears to be true with the COHA model since the likelihood function of the model is of high dimensions. In order toovercome the limitation of the MLE method, Bayesian statistic is recommended (Geman and Geman, 1984; Capper et al.,2005; Gamerman and Lopes, 2006; Kobayashi et al., 2012a,b).
In Bayesian estimation, Robert (1996) recommended to use complete likelihood function instead of conventional
likelihood function. To come up with the complete likelihood function in our model, we assume a latent variable skn torepresent the damage level of cracking type n, which occurs at time τkA . In another words, the damage level of road sectionk, with respect to cracking type l, is sk = (sk1 , · · · , ski , · · · , skL ) at time τkA . Damage level i has its potential range of values asfollows:(
0 ≤ skn ≤ i when n < l.
(44)
0 ≤ skn < i when n > lSupposing that latent variable vector s = (s1 , · · · , sK ) is measurable, the likelihood function in Eq. (43) can be further
defined: Lδ̄kil, jm −1
Y
L  Y
K
I Y
L Y
I−1 Y

Y

 


n
k 

p0sk (ȳ )
L̃(s, β : ξ̄) =




n

 
i=0 l=0 j=i k=1

m=0

n=1,,lδ̄kil, jm δ̄kil, jm
j̃n
j̃l




L
L 
 Y




X
X
Y
oδ̄kil, jm 
n










pnsk t (z̄k )
pnsk j (z̄k )plit (z̄k )








n
n








k
t=i
m=0,,l n=1,,l,,mt=sn

δ̄kil, jl

j̃n
j̃n


L

δ̄kil, jl Y XX n k lk. (45)
(z̄
)
p
p
(z̄
)


k
i
j




n=1,,l sk =0 t=sk sn t nnEq. (45) is referred as the complete likelihood function in Bayesian statistic. The complete likelihood function contains
a latent variable s. However, by assuming the posterior distribution of s through its prior distribution, it is possible togenerate the value of s. Following equation describes the posterior distribution of s.
Prob{skn = s|sk−n , ξ̄ : β} = P

L̃(skn = s, sk−n , β : ξ̄)ĩnkk
s=0 L̃(sn = s, s−n , β : ξ̄)=P

pn0s (ȳk , x̄k : β0s )
ĩn
n
k k
s=0 p0s (ȳ , x̄ : β0s )

.

(46)

5.4. Metropolis Hastings Algorithm.
In this study, we apply Metropolis Hastings (MH) algorithm (Metropolis et al., 1953; Hastings, 1970) to generate
the sample of unknown parameter β. In the algorithm, a conjugate distribution to be used as the prior distribution oftarget parameter needs to be defined. Sampled values of β are then generated using random walk process. In our study,
multi-dimensional normal distribution with mean 0 is used as the conjugate distribution.
2βilm(t) − βilm(t−1) ∼ N(0, (σmil ) ).

(47)

2
In Eq. (47), t is frequency of sample and standard deviation (σmil ) can be of any arbitrary set. The steps on generating
values of β using MH method with random walk process are briefly described in following steps:

Step 1 - Setting up initial value.
(0)
Value of standard deviation σm=il in the conjugate distribution (Eq. (47)) is an arbitrary set. The initial value of s
(0)
(1,0)
(K,0)
(s , · · · , s) is selected under the condition in Eq. (44). The initial value of β is also an arbitrary set. Through thegeneration of the MH algorithm, the influence of initial values on estimation results will gradually weaken. In our study,
sampling frequency t is equal to 1.

11

Step 2 - Sampling latent variable s(k,t) .
By applying Eq. (46), we can calculate new value of s(k,t) through using the obtained value s(k,t−1) and β(t−1) in step
5.4. The m.t.p pli j (ȳk ) is also estimated in this step.
Step 3 - Sampling unknown parameter β(t) .
The sample for unknown parameter is generated through the random walk process in the MH algorithm. To begin
M,(t−1)
with, an affixing character of each element of parameter vector β(t−1) = (β1,(t−1), · · · , βI−1L) is redefined and rewritten as00(t−1)
β(t−1) = (β(t−1)
,
·
·
·
,
β
).
Followingprocedures
explainthesampling
process:J
1
Step 3.1: s(k,t) and β(t−1) are assumed to be available from previous steps.
Step 3.2: To increase the loop iteration (frequency) t and to define sub-domain of damage level q:t,q t,q−1t,q−1 ′
β(t,q) = (βt,q).
1 , · · · , βq , βq+1 , · · · , β J(48)In addition, step-width vector ξtq = (0, · · · , 0, ξqt , 0, · · · , 0)′ , with element i in its value of ξqk , is defined. The approval of
expected (on average) step-width in random walk process follows normal distribution with mean 0.
ξqt ∼ N(0, (σq )2 ).

(49)

As a result, the standard deviation can be estimated:#"L̃(s(t) , β(t,q) , ξ̄),
1,α(t,q) = minL̃(s(t) , β(t,q−1) , ξ̄)(50)with (L̃) as the completed likelihood function.
Step 3.3: Random number u ∼ U(0, 1) is generated throught uniform distribution U(0, 1) in [0, 1]. Unknown parameterβ(t,q) is decided according to the following rules:

(t,q)

β

=

(

β(t,q) + ξtq
β(t,q)

i f u ≤ α(t,q).
otherwise(51)

The above-mentioned procedure is applied with q from q = 1 to q = L.
Step 4-Updating Parameter
Updated value β(t) is recorded at any iteration and compared for acceptance or rejection. If t ≤ t as t = t + 1, theprogram returns to Step 2. Otherwise, the algorithm stops.
It is likely that initial values of parameters remains only in the early stages of iteration process. When t increases in ahigh number, reaching a certain numbers t, we can define our preferable value of β(t) (t = t + 1, t + 2, · · · , t) as the outcome.
5.5. Statistical Inference on Posterior Distribution.
Statistical test for parameter β can be carried out based on the generated samples through the MCMC simulation.
However, in the simulation, probability density function π(β|ξ̄) cannot be considered as an analytical function. Therefore,
instead of using full parametric approach for statistical test, non-parametric approach is recommended Robert (1996).
According to MH method, among generated samples β(t) (t = 1, · · · , t), the first t samples are removed. A new set ofsamples is then defined as a replacement, with its subscriptions M = {t + 1, · · · , t). By applying this approach, the joint
probability distribution function G(β) is defined:
G(β) =#{β(t) ≤ β, t ∈ M}
,
t−t

(52)

where #{β(t) ≤ β, t ∈ M} is regarded as the total number of samples, from which logical expression β(t) ≤ β, t ∈ M issatisfied. Moreover, expected values of ζ̃(β) and standard covariance Σ̃(β) of the posterior distribution of β are defined
respectively as:ζ̃(β) = (ζ̃(β1 ), · · · , ζ̃(β J ))′ = 2··· σ̃ (β1 ) ..
..
Σ̃(β) =  .
.
σ̃(β J β1 ) · · ·

t
X
β(t)
1

t−t
t=t+1

σ̃(β1 β J ) 

..
 ,
.

σ̃2 (β J )
12

,···,

t
X
β(t) ′
J

t−t
t=t+1

,

(53-a)(53-b)

Table 2: Cracking Condition States.
Condition states
Description
0,0 (0)New condition states
1,1 (1)longitudinal crack Small crack
2,1 (2)Medium crack
3,1 (3)Large crack
1,2 (4)
Horizontal crack
Small crack
2,2 (5)
Medium crack
3,2 (6)Large crack
1,3 (7)
Alligator crack
Small crack
2,3 (8)
Medium crack
3,3 (9)Large crackNote) Figures in the parentheses shows measured condition states of cracks.

where
σ̃2 (βi ) =tX{β(t) − ζ̃(βi )}2
i

t=t+1

σ̃(βi β j ) =t−tt {β(t) − ζ̃(β )}{β(t) − ζ̃(β )}
X
i
j
i
j

t−t

t=t+1

(54-a),(54-b).

The credible interval of parameter β is examined and determined by using generated samples. For example, the 100(1 −εε2ε)% credible interval of parameter β is defined by using statistical sampling order βε , β j ( j = 1, · · · , J) with βε < β j < β j :
j

 (t)



#{β j ≤ β∗j , t ∈ M}βε = arg max
≤
ε
,




j

β∗j 
t−t
 (t)

∗∗




ε
 #{β j ≥ β j , t ∈ M}
β j = arg min
≤ ε
.



∗∗

βj 
t−t

j

(55-a)(55-b)

6. Empirical Study
6.1. Data and definition of condition states
We conducted an empirical application of the HIMA model using a representative set of cracking data of road sectionsin Japan. Data was recorded during the period from 1992 to 2004 and consisted of three consecutive inspection dataon same road sections. Total numbers of investigated samples are 2, 751, with each sample representing for an average
sectional length of 100 meters. Beside the percentage values of cracks measured in three inspection times, values ofvariables such as traffic volume, ambient temperature, thickness of road sections were also recorded.
In the study, cracks are classified into three types (l): longitudinal crack (l = 1), horizontal crack (l = 2), and alligatorcrack (l = 3). Damage levels i for each cracking type are also defined as discrete value in a range of [0, 3] (i(i = 0, ..., 3))according to its size (small, medium, and large). The condition states of road sections are described as the combination of
cracking types and their damage levels. As a result of combination, there are 10 condition states. The description of eachcondition state is shown in Table 2.
6.2. Transition of condition states
This section explains the transition among condition states defined in Table 2. It is assumed that at times τA and τB(refer to Fig. 1), condition states are (i, l) and ( j, m), respectively, with z as time interval. The absorbing condition state isdefined as a combination of (i = 3, l = 3), with its transition probability π33,33 = 1.
If at time τA , damage level is i = 2, the transition of condition state can be: 1) no further deterioration till time τB .
Condition state (2, l) (l = 1, 2, 3) remains in period z, with transition probability (π2l,2l ); 2) Damage level i = 2 remains induration z. However, there is a change in cracking type l. Transition probability of this event is π2l,2m (m = l + 1, · · · , 3);
3) Both damage level i and cracking type l change at time τB with transition probability π2l,3m (m = 1, 2, 3).
In case damage level at time τA is i = 1, transitions can be described as one of the events: 1) no change of either
damage level and cracking type from time τA to τB . The probability of the event is (π1l,1l ); 2) Cracking type changes
at time τB , with transition probability (π1l,1m (m = l + 1, · · · , 3)); 3) Damage level i changes one step at time τB , withtransition probability (π1l,2m (m = 1, 2, 3)); 4) Damage level i changes to i = 3 at time τB , with transition probability(π1l,3m (m = 1, 2, 3)).
13

Table 3: Model’s parameters βi,l , hazard rate θi,l , and duration RMDi,l of transition.
Condition states
Constant term
Traffic volume
Hazard rate Duration(i, l)
β1i,l
β2i,l
θi,l
RMDi,l(0,1)-2.7050.067
14.955(-2.795,-2.623)(1,1)-1.8780.153
6.538(-2.082,-1.697)(2,1)-0.5720.564
1.772(-0.895,-0.278)(0,2)-2.299
0.284
0.126
7.944
(-2.417,-2.186)(0.196,0.433)(1,2)-1.560
0.720
0.373
2.685
(-1.831,-1.360)(0.509,1.037)
(2,2)
0.227
1.255
0.797
(0.0886,0.364)(0,3)-5.8300.003(-6.233,-5.508)(1,3)-1.3670.255
3.925(-2.242,-0.565)(2,3)Note) Values in parentheses are lower bound and upper bound values of credible interval corresponding to 95% of its significant level.

Finally, in case condition state at time τA is (0, 0), the transitions among condition states can be: 1) no deteriorationfrom time τA to τB . The transition probability of the event is (π00,00 ); 2) Damage level changes one step, crack appears, withtransition probability (π00,1m (m = 1, 2, 3)); 3) Damage level changes two steps, crack appears, with transition probability(π00,2m (m = 1, 2, 3)); 4) Damage level changes three steps, crack appears, with transition probability (π00,3m (m = 1, 2, 3)).
6.3. Results
In the study, hazard rate θi,l of each condition state was assumed to be dependent only on traffic volume (Eq. (56)).
This assumption was acceptable as traffic volume is considered as one of the main factors causing deterioration of road
(Tsuda et al., 2006). As a matter of fact, crack initiation could be also related to other factors such as: ambient temperature,materials, and axes load. However, in our study, those factors can be considered as already incorporated in the constantterm of regression function. The denotation of constant term and traffic volume used in the study are x1k = 1 and x2k ,respectively.
The objective of estimation is to estimate the model’s parameters βi,l , which is often referred as unknown parametersin statistical models (Eq. (45)). Following equation describes the function form assumed for the hazard rate θi,l .
k
θ00
=

2
3 X
X
n=1 m=1kkβm0n xm and θin =

2
X

k
βm
in xm ,

(56)m=1

(i = 1, 2; n = 1, 2, 3; k = 1, · · · , 2751).
In the program (coded in FORTRAN 90), at each iteration, new values βi,l are recorded. To verify its values, likelihood
ratio R(β̂−m ) (m = 1, · · · , I) is examined. Values of β will be accepted as convergent values with a certain degree ofsignificance R(β̂−m ) = 2{ln L(β̂) − ln L(β̂−m )}, where β̂−m is a vector, in which, m element of β̂m is substituted for β̂ withits mean equals to 0. When |R(β̂−m )| ≥ 3.48 is observed, null hypothesis βm = 0 can be dismissed by a significant level
5%.
Results of estimation using the COHA model with MCMC simulation are shown in Table 3. In the table, values inparentheses are lower bound and upper bound values of the credible interval with 95% significant level, calculated byusing Eqs. (55-a) and (55-b). Values of β1il in Table 3 infers that other influencing factors to deterioration other than traffic
volume have considerably impacts on the hazard rate of condition state (i, l) = (1, 1),(2, 1),(2, 2), and (1, 3). Hazard rateof condition state (2, 3) is not obtained due to a reason that data on damage level 3 is not available. Estimation results
also highlight a fact that annual traffic volume has a significant impact on the deterioration of horizontal crack, with theirparameter values β202 and β212 equal to 0.284 and 0.720, respectively.
The two last right column of Table 3 show the values of hazard rate θi,l calculated based on Eq. (56) and the life
expectancy of condition state (or the duration of being in a condition state) calculated based on Eq. (15). It can beinterpreted from the values of hazard rates θ that longitudinal crack progresses significantly slower than horizontal crack
and alligator crack. Furthermore, the life expectancy of horizontal cracks becomes relatively short after longitudinalcracks occurs. The life expectancy of damage level 1 of alligator crack is not counted as a result of insufficient recorded
data. This problem might be due to the past interventions, with objectives to heal only for horizontal and longitudinal

14

Table 4: Markov transition probability according to cracking type
Damage levels
Cracking types
01
2
3
Longitudinal crack
00.935 0.060 0.004 0.001
1
0.0
0.858 0.108 0.034
0.0
0.0
0.569 0.431
2
3
0.0
0.0
0.0
1.0
Horizontal crack
0.882 0.098 0.013 0.007
01
0.0
0.689 0.170 0.141
2
0.0
0.0
0.285 0.715
3
0.0
0.0
0.0
1.0
Alligator crack
0.997 0.003 0.000 0.000
01
0.0
0.775 0.225 0.000
2
0.0
0.0
0.997 0.003
3
0.0
0.0
0.0
1.0
Note) Transition probabilities are in one year term.

Condition states
(0,0)
(1,1)
(1,2)
(1,3)
(2,1)
(2,2)
(2,3)
(3,1)
(3,2)
(3,3)

(0,0)
0.822
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0

Table 5: Competing m.t.p for pairs of condition states.
(1,1)
(1,2)
(1,3)
(2,1)
(2,2)
(2,3)
0.053 0.098 0.003 0.004 0.013 0.000
0.755 0.084 0.002 0.105 0.014 0.000
0.0
0.678 0.003 0.007 0.170 0.000
0.0
0.0
0.728 0.007 0.023 0.219
0.0
0.0
0.0
0.541 0.016 0.001
0.0
0.0
0.0
0.0
0.284 0.000
0.0
0.0
0.0
0.0
0.0
0.963
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0

(3,1)
0.001
0.034
0.002
0.003
0.422
0.001
0.004
0.971
0.0
0.0

(3,2)
0.006
0.006
0.140
0.020
0.020
0.715
0.029
0.029
1.000
0.0

(3,3)
0.000
0.000
0.000
0.000
0.000
0.000
0.004
0.000
0.000
1.0

Note) Transition probabilities are in one year term.

cracks. Overall, it takes about more than 21 years for a road section in perfect condition state to encounter heavily damagelevel of longitudinal crack and about 10 years to be in the last damage level (i = 3) of horizontal crack.
The transition probabilities among damage levels of respective cracking types are presented in Table 4. The propertiesof transition probabilities are estimated by means of expected hazard rates, which are calculated on average basis using
the MUSTEM model.
Using the Eqs. (20), (31), and (35), the properties of the competing m.t.p for each pair of condition state (i, l) are
estimated and shown in Table 5. The properties of the competing m.t.p in Table 5 are not purely homogenous Markov
properties. The properties reflect transitions among pairs of condition states in one year term after last interventions. Aninteresting finding is that transition probabilities of initial damage levels of either longitudinal crack and horizontal crack
are considerable high, reflecting the true cracking deterioration in the area of targeted road sections.
We present in Fig. 3 deterioration curves of respective cracking types. The horizontal axis indicates elapsed time (inyears). The blue lines demonstrate the change in value of transition probability, in which, damage levels of cracking types
remain in the same states. The two lines in the lower part of the figure are referred to same cracking types. It is observablethat there is a sharp decrease in the values of transition probability concerning damage levels of horizontal crack, while
it is relatively slow with longitudinal crack. This is corresponding to the results in Table 3. Thus, generally, suffice tosay that the overall cracking of road sections is dominated by deterioration of horizontal crack. In addition, correlationin term of deterioration curves between horizontal crack and longitudinal crack reveals a potential of horizontal cracking
influence over the longitudinal crack. Finally, the pink lines illustrate the change of transition probability among different
cracking types.
Fig. 4-a additionally gives information on the distribution of condition state (i, l) over the time span of 25 years.
Calculation for drawing this figure is based on obtained competing m.t.p Pil (y) through Eq. (38). Colored patterns in the
figure have their respective implications to cracking types and damage levels.
It can be seen from the figure that, within 3.5 years of usage, numbers of road sections in good condition account for
about 50% out of the total numbers. Within the same period, deterioration concerning horizontal crack is more intensethan longitudinal crack, with a percentage of 20 % in term of distribution. Approximately 8 years after the initial time,numbers of road sections with longitudinal crack are forecasted to reach 50%. The distribution of alligator crack in thisfigure is not significant observed. With this observation, we realize that horizontal crack is the dominant cracking type.

15

π 13,21
π 13,22

π 13,13
π 21,210.03

π 13,23
π 21,31

π 21,22
π 21,23

Transition probability

0.02
0.01
0
0

10

20

30

40

50

-0.01
-0.02
-0.03
-0.04
-0.05

Elapsed time (years)

Figure 3: Expected deterioration curves of cracks.

In practice, most of interventions have been carried out to recover the horizontal crack and longitudinal crack before theoccurrence of alligator crack.
100%

100%

(3,3)
(3,2)
(3,1)
(2,3)
(2,2)
(2,1)
(1,3)
(1,2)
(1,1)
(0,0)

80%

Percentage

70%
60%
50%
40%
30%

(3,3)
(3,2)
(3,1)
(2,3)
(2,2)
(2,1)
(1,3)
(1,2)
(1,1)
(0,0)

90%
80%
70%
Percentage

90%

60%
50%
40%
30%

20%

20%

10%

10%

0%

0%
1

5

9

13

17

21

25

1

Elapsed time (years)

5

9

13

17

21

25

Elapsed time (years)(a) - COHA model(b) - HIMA model
Figure 4: Condition states distribution

6.4. Comparison
Empirical applications on the same sort of cracking data were further carried out with the MUSTEM model and the
HIMA model for the comparison of estimation results. The MUSTEM model was used separately with each cracking
type. Whilst, the HIMA model was used with three cracking types at the same time like in the COHA model. Due to thedifferences of the three models, the MUSTEM model and the HIMA model consider only the observed sampling valuesof typical cracking type, hidden cracking types are not considered in these two models. Moreover, with the HIMA model,
it is not possible to construct deterioration overall curves since the HIMA model considers only each cracking type andits damage level.
Distribution of condition states over time using the HIMA model is shown in Fig. 4-b. This figure can be compareddirectly with the distribution of condition states distribution using the COHA model (Fig. 4-a). As can be seen fromcomparison between the two figures (Fig. 4-a and Fig. 4-b), the differences in the shapes and patterns corresponding to
the damage levels under two models are not noteworthy observed. However, there is a minor noticeable difference withregard to the probability distribution of horizontal crack (i.e. becomes smaller along with time axis). Reason for thedifference of condition state distribution of horizontal crack is believed due to the overestimation of horizontal cracking
type when using the HIMA model.

16

Elapsed time (Years)
0

5

10

15

20

25

Damage level

0

1

2

3
Longitudial crack - COHA modelLongitudial crack - MUSTEM model

Horizontal crack - COHA model

Horizontal crack - MUSTEM model

Figure 5: Comparison of deterioration curves between the MUSTEM model and the COHA model.

The comparison of estimation results between the MUSTEM model and the COHA model are shown in Fig. 5. In
either cases, there appear a substantial discrepancy. Deterioration speed obtained from the MUSTEM model is sharply
faster than that of the COHA model. These differences are due to the fact that the MUSTEM model uses only representing
sampling data, which are always in worse condition states. This problem is regarded as overestimation of true condition
states. Meanwhile, the COHA model has ability to reveal the hidden condition states of damage levels and cracking types.
Thus, its estimation results are greatly improved.
7. Conclusions
Deterioration of an infrastructure object can be represented by means of several performance indicators. In monitoring
activities, it is often the case that worse value of performance indicator is selected as representative condition state. Valuesof other performance indicators could be neglected, and thus remaining as hidden information in the data bank. A typicalexample of such a problem is with monitoring data of cracks. Cracks occur on the surface of a road section can be inmultiple directions (e.g. longitudinal crack, horizontal crack, and alligator crack). In many cases, only crack with worse
value (i.e. percentage of crack on road surface), is chosen to be representative crack. This selection bias could possibly
result in non-optimal decisions on chosing intervention strategy to be executed on road sections in order to provide
adequate level of service.
In order to tackle the problems, we developed a novel Markov hazard model that can be used for deterioration prediction of infrastructure system with more than a single performance indicator and with multiple condition states. The
model was formulated to mitigate the selection bias embedded in monitoring data. Precisely, in the model, it is assumed
that deterioration process can be of multiple types. Each type of deterioration is measured by use of damage levels. The
overall condition state of an infrastructure object is the combination of deterioration type and damage level. In addition, it
is also assumed that condition state is defined based on a competitive selection among damage levels and cracking types,
and therefore, the model’s name is a competitive Markov model.
In the paper, we also presented a novel numerical solution to obtain model’s parameters by using Bayesian estimation
method and Markov Chain Monte Carlo simulation. The model was then tested with a cracking dataset of road sectionsin Japan. It was found that horizontal crack is the dominant cracking type. Longitudinal crack progresses slower thanhorizontal crack and alligator crack.
Finally, it is concluded that the model is robust and the estimation results are significant. The model can be applied
not only for modelling crack initiation but also for other types of deterioration as well, e.g. the corrosion of rebars due toboth carbonition and chloride attack.
References
AASHTO, 2004. Pontis Release 4.0 Technical Manual. American Association of State Highway and Transportation Officials, Washington DC.
Andrew, G., John, B. C., Hal, S. S., Donald, B. R., 2006. Bayesian Data Analysis. Chapman and Hall/CRC.
Capper, O., Moulines, E., Ryden, T., 2005. Inference in Hidden Markov Models. Springer.

17

Dani, G., Hedibert, F. L., 2006. Markov Chain Monte Carlo-Stochastic Simulation for Bayesian Inference. Chapman andHall/CRC-Taylor and Francis Group.
Dore, G., Zubeck, H., 2009. Cold Regions Pavement Engineering. Mc Graw Hill.
Fukuhara, T., Terada, K., Nagao, M., Kasahara, A., Ichihashi, S., 1990. Automatic pavement-distress-survey system.
Journal of Transportation Engineering 116 (3), 280–286.
URL http://link.aip.org/link/?QTE/116/280/1
Gamerman, D., Lopes, H. F., 2006. Markov Chain Monte Carlo, Stochastic Simulation for Bayesian Inference. Chapman& Hall/CRC.
Geman, S., Geman, D., 1984. Stochastic relaxation, gibbs distributions and the bayesian restoration of images. IEEE
Transactions on Pattern Analysis and Machine Intelligence 6 (6), 721–741.
Hastings, W. K., 1970. Monte carlo sampling methods using markov chains and their applications. Biometrika 57 (1),
97–109.
URL http://biomet.oxfordjournals.org/content/57/1/97.abstract
Ibrahim, J., Ming-Hui, C., Sinha, D., 2001. Bayesian Survival Analysis. Springer Series in Statistics.
Jeff, G., 2006. Bayesian Methods: A Social and Behavioral Sciences Approach, 2nd Edition. Chapman and Hall/CRC
Taylor and Francis Group.
Kaito, K., Kumada, K., Hayashi, H., Kobayashi, K., 2007. Modeling process by cracking pavement deterioration hazard
index hierarchical model. Journal of Civil Engineering, JSCE (in Japanese) 63 (3), 386–402.
Kawashima, Y., Fukushima, D., Yasuo, M., 1984. A study on mechanism of cracking on asphalt pavement. Laboratory
Report in Fiscal 1983, 67–79.
Kenneth, D. K., June 2010. Network-level infrastructure managemenent using approximate dynamic program. Journal of
Infrastructure Systems 16 (2), 103–111.
Kobayashi, K., Kaito, K., Lethanh, N., April 2012a. A bayesian estimation method to improve deterioration prediction for
infrastructure system with markov chain model. International Journal of Architecture, Engineering and Construction
1 (1), 1–13.
Kobayashi, K., Kaito, K., Lethanh, N., May 2012b. A statistical deterioration forecasting method using hidden markov
model for infrastructure management. Transportation Research Part B: Methodological 46 (4), 544–561.
Lancaster, T., 1990. The econometric analysis of transition data. Cambridge University Press.
Lethanh, N., September 2009. Stochastic optimization methods for infrastructure management with incomplete monitoring data. Ph.D. thesis, Graduate School of Engineering, Kyoto University.
URL http://hdl.handle.net/2433/85384
Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., Teller, A. H., Teller, E., 1953. Equation of state calculations byfast computing machines. The Journal of Chemical Physics 21 (6), 1087–1092.
URL http://link.aip.org/link/?JCP/21/1087/1
Mohajeri, M. H., Manning, P. J., 1991. An operating system of pavement distress diagnosis by image processing. Transportation Research Record 1311, 120–130.
Nakat, Z., Madanat, S., September 2008. Stochastic duration modeling of pavement overlay crack initiation. Journal of
Infrastructure Systems (ASCE) 14 (3), 185–192.
Rababaah, H., Vrajitoru, D., Wolfer, J., 2005. Asphalt pavement crack classification: A comparision of ga, mlp, and som.
In: The Genetic and Evolutionary Computation Conference. Washington DC.
Robelin, C., Madanat, S., 2007. History-dependent bridge deck maintenance and replacement optimization with markov
decision process. Journal of Infrastructure System, ASCE 13 (3), 195–201.
Robert, C., 1996. Mixtures of Distributions: Inference and Estimation, in: Gillks, W.R., Richardson, S. and Spiegelhalter,
D.J. (eds.) Markov Chain Monte Carlo in Practice. Chapman and Hall.
Shahin, M., 2005. Pavement Management for Airports, Roads, and Parking Lots. Springer.

18

Thompson, P. D., Small, E. P., Johnson, M., Marshall, A. R., 1998. The pontis bridge management system. StructuralEngineering International 8 (4), 303–308.
URL http://www.ingentaconnect.com/content/iabse/sei/1998/00000008/00000004/art00018
Tsuda, Y., Kaito, K., Aoki, K., Kobayashi, K., 2006. Estimating markovian transition probabilities for bridge deterioration
forecasting. Journal of Structural Engineering and Earthquake Engineering 23 (2), 241–256.
Wang, K., Smadi, O., 2011. Automated imaging technologies for pavement distress surveys. in: Transportation research
circular no. e-c156. Tech. rep., Transportation Research Board of the National Academies.

19