Seismic vulnerability assessment of residential buildings using logistic regression and geographic information system (GIS) in Pleret Sub District (Yogyakarta, Indonesia)

The Southeast of Yogyakarta City has had the heaviest damages to buildings in the 2006 of Yogyakarta Earthquake disaster. A moderate to strong earthquake of 6.3 Mw shook the 20 km southeast part of the Yogyakarta City early in the morning at 5:54 local time. On top of extensive damage in Yogyakarta and Central Java, more than 5700 people perished; 37,927 people were injured in the collapse of more than 240,396 residential buildings. Furthermore, the earthquake also affected the infrastructure and local economic activities. The total damages and losses because of the earthquake was 29.1 trillion rupiahs or equal to approximately 3.1 million US dollar. Two main factors that caused the severe damages were a dense population and the lack of seismic design of residential buildings. After reconstruction and rehabilitation, the area where the study was conducted grew into a densely populated area. This urbanistic change is feared to be potentially the lead to a great disaster if an earthquake occurs again. Thus, a comprehensive study about building vulnerability is absolutely needed in study area. Therefore, the main objective of this study has been the provision of a probabilistic model of seismic building vulnerability based on the damage data of the last big earthquake. By considering the relationship between building characteristics, site conditions, and the damage level based on probabilistic analysis, this study can offer a better understanding of earthquake damage estimation for residential building in Java. The main findings of this study were as follows: The most vulnerable building type is the reinforced masonry structure with clay tile roof, it is located between 8.1-10 km of the epicentre and it is built on young Merapi volcanic deposits. On the contrary, the safest building type is the houses which has characteristics of reinforced masonry structure, asbestos or zinc roof type, and being located in Semilir Formation. The results showed that the building damage probability provided a high accuracy of prediction about 75.81%. The results explain the prediction of building vulnerability based on the building damaged of the Yogyakarta earthquake 2006. This study is suitable for preliminary study at the region scale. Thus, the site investigation still needs to be conducted for the future research to determine the safety and vulnerability of residential building.


Background
An earthquake is a sudden motion of the earth caused by the release of accumulated energy that mostly occurs within or along the edges of earth's tectonic plates. The released energy generates the seismic vibration radiating in the earth's body and perceived as an earthquake (Bath, 1979). This natural phenomenon strongly relates to the geological condition and the configuration of tectonic plates. The areas located at the edges of tectonic plates are more vulnerable to earthquake hazard. However, few earthquakes can also occur away from the tectonic plates boundary such as in the interior of a tectonic plate or intraplate regions (Stein & Wysession, 2003).
Indonesia is an archipelago country located between three active tectonic plates: The Pacific Plate, the Eurasian Plate, and the Indo-Australian Plate. Demets et al. (1994) and Schluter et al. (2002) stated that the Eurasian Plate is relatively more stable than the other two tectonic plates. The Pacific Plate is moving northward with an average velocity of 11-12.5 cm/year while the Indo-Australian Plate is moving 7.23 cm/year westward. According to U.S. Geological Survey (USGS), at least 14,000 earthquakes greater than 5.0 Richter occurred in this area between 1900 and 2009 with the biggest earthquake (9.1 MW) having occurred in the Andaman Sea, west coast of Aceh (northern Sumatra) on 26 December 2004. The Andaman earthquake also generated a destructive tsunami and caused massive casualties not only in Indonesia but also in other nations surrounding the Indian Ocean such as Sri Lanka, India, Thailand, Somalia, Myanmar, Maldives, etc. The other big earthquakes occurred in Indonesia after Aceh Earthquake 2004 were Nias Earthquake in 2005 (Mw = 8.7), Yogyakarta Earthquake in 2006 (Mw = 6.3), Tasikmalaya Earthquake in 2009 (Mw = 7.4), Padang earthquake in 2009 (Mw = 7.6), and Kebumen, Central Java Earthquake in 2014 (Mw = 6.1) (Irsyam, et al., 2010).
Yogyakarta earthquake on May 27th, 2006, was unexpected not only for Yogyakarta residents but also for all Indonesians. A moderate to strong, 6.3Mw earthquake hit Yogyakarta early in the morning at 5:54 local time. According to USGS, the epicentre was located 20 km Southeast of Yogyakarta City at geographic coordinates 7.9620 0 S, 110.4580 0 E. This earthquake caused extensive damage in Yogyakarta and Central Java. More than 5,700 people died, 37,927 people were injured, and 240,396 residential buildings were destroyed, and the local infrastructure and economic activities were largely disrupted. According to the damage assessment conducted in June 2006, the cost of damage and losses reached 29.1 trillion rupiahs (US$ 3.1 million). This made the Yogyakarta earthquake one of the worst disasters in Indonesia in the last ten years (BAPPENAS 2006).
Housing was the most severely impacted sector, with more than half of the total private building stocks of 2004 being affected and an estimated total loss of about 15.3 trillion rupiahs. The impacts on other sectors such as infrastructure, lifeline, and trade were widespread but relatively limited in severity. For example, the Pedan electricity substation was switched off directly after earthquake. Three sets of 500KV circuit breakers, five sets of 500KV and two sets of 500 KV/150KV electricity transformers, and a 500KV lighting arrester suffered significant damages. Social, education and productive sectors also experienced significant damage due to the earthquake. At least 2,155 school buildings, 17 private hospitals, 41 private clinics and 45 health clinics (PUS-KESMAS) were destroyed, and the same happened with as many as 2,201 religious facilities (20% of their total number in 2004). Additionally, a large number of enterprises, mostly small and mid-sized ones, such as shops and traders, were completely damaged. Furthermore, several main traditional markets such as Niten, Imogiri, Pleret and Piyungan were closed. The summary of estimated damages and losses in infrastructure and productive sectors are shown in Table 1 while the details about affected small and medium enterprises can be seen in Table 2.
Two primary factors that caused the severe damages were the high density of population in the affected areas (1,600 people per km 2 ) and the lack of seismic design of housing units (Elnashai et al., 2006). Before 2006, the house types in Yogyakarta and its suburbs, could be divided into three categories: 1) unreinforced masonry or URM (in older houses, i.e., pre-1990), 2) partially reinforced masonry (newer houses, post-1990) and 3) Source: (BAPPENAS 2006) traditional timber houses (Joglo). Most of the URM collapsed due to the lack of mechanical connection between roof, walls and floor. This structural failure was responsible for most of the deaths and injuries (Earthquake Engineering Research Institute (ERRI), 2006). In general, RM 1 and RM2 performed well during the 27 May 2006 earthquake, but several of them collapsed due to the poor connections between walls or columns and roof. Timber frame houses such as the Javanese traditional house or "Joglo" performed well because of their good connection between roof, columns, and floor. However, detailed investigation of 8 Joglo houses showed that 5 Joglo houses had collapsed while the rest of them were only slightly damaged. Based on thorough visual investigation, the damage of the Joglo houses can be classified into three main categories: columns-to leg connection; joint between main columns and beams; and roof construction and its attachment with the core structure (Prihatmaji et al. 2012). The 27 May 2006 earthquake strongly impacted the Pleret Sub District in Bantul District, 10 km southeast of Yogyakarta city. Approximately 8,309 building units collapsed and 579 people died (BPS-Statistics of Bantul Regency 2010). After the reconstruction and rehabilitation processes finished, Pleret Sub District has been growing into one of the most densely populated areas in the southeast part of Yogyakarta City. The population of Pleret Sub District in 2013 was 45,136 people with a population growth rate of 0.05 (BPS-Statistics of Bantul Regency 2014). The increasing numbers of population in Pleret Sub District will also be followed by a rising demand for housing. As the population increases, the number of structures at risk also increases. Moreover, most of the buildings in Pleret Sub District are built on the top of dense volcanic sediment (Volcanic deposit of Merapi Volcano) which can amplify the earthquake wave and will increase the surface tremor when an earthquake occurs (Daryono, 2011). This condition increases the probability of a great disaster. Moreover, Java islands, especially in the south part of Java, often experience moderate to strong earthquakes (M > 6.0) with 50-100 years of return period (Prihatmaji et al. 2012). Despite being experienced with the earthquakes, an integrated risk analysis of earthquake hazard for rural communities is still necessary in the research area. Therefore, the study of the comprehensive earthquake vulnerability of residential buildings based on the probabilistic model in geographic information systems (GIS) environment is needed in earthquake prone areas like Pleret Sub District. The buildings damage prediction can be obtained through the probabilistic analysis based on the building damage data of the 27 May 2006 earthquake. Additionally, GIS will provide a better explanation of building damage distribution in order to support the local government in disaster and risk management phase.
There was an immediate international response and reaction after the main shock at 05:54 local time at 27 May 2006. Approximately 100 response operations reacted quickly and brought to a rapid acquisition of satellite imagery to support the preliminary damage assessment process. A variety of damage maps were available in a week after the earthquake. Those damage maps mostly were produced by UN Institute for Training and Research (UNITAR)'s Operational Satellite Applications Programs (UNOSAT) (Kerle, 2010). The rapid ground survey was conducted by the Universitas Gadjah Mada (UGM) in Yogyakarta. As a result, an extensive database of building type, function and construction material of affected buildings were made available soon after the earthquake.
Probabilistic seismic hazard assessment (PSHA) and deterministic seismic hazard assessment (DSHA) have been growing rapidly both for governmental purposes and scientific research since the Yogyakarta earthquake in 2006. PSHA and DSHA help producing an earthquake hazard zonation in Yogyakarta and surrounding areas. However, the study about the damage patterns of Yogyakarta earthquake in research is limited and absolutely needed. A previous relevant study was conducted by Nurwihastuti, et al. (2014) which had the main objective to investigate the Yogyakarta damage pattern Source: (BAPPENAS 2006) through the geomorphological approach including investigation of surface and subsurface characteristics. The results showed that the most severe damage tended to occur in areas which had particular characteristics: a deep basement layer, low gravity anomaly, thick surface sediment, and unconsolidated surface material. Based on this result, a study about the relationship between building characteristics, site condition, and the damage level based on probabilistic analysis was regarded as important to give a complete understanding of earthquake damage pattern in Yogyakarta. The objectives of this study is to conduct a building vulnerability assessment based on probabilistic analysis and to estimate the vulnerability of buildings and population in the study area under different scenarios of population distribution.
The Pleret Sub District is predominantly located on two major geological units. The main part of Pleret Sub District (the capital sub district) is situated in a low-relief alluvial deposit while the East area of Pleret consists of an ancient volcanic deposit which is collectively referred to as the Semilir Formation and Nglanggran Formation. According to the geological map of Yogyakarta (scale 1:100,000) (Rahardjo 1995), those two geological units can be subdivided into several smaller formation units as follows (see Fig. 2).

Alluvium (Qa):
it generally consists of gravel, sand, silt and clay along the river.
Both the Nglanggran Formation (Tmn) and the Semilir Formation (Tmse) are the Tertiary volcanic deposits that were formed between late Oligocene and early Miocene, respectively (Mulyaningsih et al., 2011). Tmse consists of volcanic clastic materials with the pumice as the main material. This fragmental material has various grain size fragments starting from a very fine tuff until breccia pumice that has very coarse grains fragments (Yusliandi et al. 2013). The abundance of pumice fragments in Tmse indicates that Tmse is typically co-ignimbrite deposit. In term of volcanology sediment, the coignimbrite deposit can be classified as a volcanic material which was formed by close explosive eruption (Mulyaningsih et al., 2011); (Yusliandi et al. 2013); (Winarti 2015); (Bronto et al., 2009). On the contrary, Nglanggran Formation (Tmn) was formed as a result of an effusive eruption and deposited on the top of the Semilir Formation. Tmn consists of solid material of breccia, lava andesite, and basalt. This formation is wide-spread along the Baturagung Escarpment in the west part of Parangtritis and east part of Gunung Panggung. Similar to Semilir Formation, Nglanggran is lacking fossils. However, based on the foraminifera content that was found in the insertion of sandstone and claystone in the bottom-most layer, Nglanggran Formation is the middle Miocene deposits. These Miocene formations were buried by younger sedimentary rock including the Young Volcanic Deposits of Merapi Volcano (Qmi) and the Alluvium (Qa). The young volcanic deposits of Merapi Volcano consist of the Young Merapi Volcano sediment that was transported by several big rivers such as Opak River while the alluvium was formed through the denudation processes on the steep areas. The Qa and Qmi are characterized as dense soil located in the extensive flat land along the Opak River. In the eastern part of the research area, uplift and erosion have stripped away much of the cover rock, exposing the underlying rock, Nglanggran and Semilir Formation.
There are many geological faults formed in the research area, as a result of the plate movements along the subduction zone in the south part of the Java islands. One of them is known as Opak Fault which is often associated with the Opak River. Abidin, et al. (2009) concluded that this SW-NE normal fault is an active one. Another normal fault which has the same orientation with Opak Fault lies in the middle of the research area. This normal fault is known as Bawuran Fault. Both of them have the same movement, i.e., the west part of the fault line is moving downward while the east part of the fault line is moving upward. There are also two major strikes-slip faults that located in the research area, namely Bawuran-Cinomati (centre part) and Becucu-Tekek Fault (north part) (see Fig. 2). Both of them were formed later after the development of Opak Fault. In the post-stage of uplifting movement of Opak Fault, the strike-slip faults were created and trimmed horizontally the research area into north and south area. The north part is moving eastward while the south part is moving westward (Sanjoto 2004).
Located on the western flank of Baturagung Escarpment, The Pleret Sub District has various topographical conditions. The lowest point (39.9 m) is located on the alluvial plain along the Opak River while the highest point (344.7 m) is located at the summit of Baturagung Escarpment. The steepest slopes (>350) are found at the hilly area of Baturagung Escarpment, located on the east side of the study area.
Based on their genesis, the general landform of Pleret Sub District consists of three major landforms of structural, fluvial and denudation origin. The structural landforms can be recognized from the morphological uplift and depress in the east and the west part of the study area. The intensive denudation processes occur in the middle and upper slope of the escarpment and the hilly area of Semilir Formation, which has less vegetation due to the intensive mining activities. The fluvial landform containing alluvial plain is located along the Opak River to the west part of study area including the main part of the Pleret Sub District (see Fig. 3)

Historical and recent earthquake events
Yogyakarta Province that lies in the south part of Java island is one of the most seismically active regions in Indonesia. This seismic condition is greatly influenced by the subduction process of the Indo-Australian and Eurasian Plates. Additionally, an active normal fault is also located in the middle part of the study area. This normal fault, Opak Fault, is associated with Opak River with SW-NE direction. The uplift zone is located to the east of this fault and vice versa. (Rahardjo et al. 1995).
Twelve great earthquakes between 1840 and 2006 were historically or instrumentally recorded near the Yogyakarta and Central Java Province. Many of these earthquakes caused damage, or even worse, casualties (Table 3). Based on the historical source, the 7.2 Mw earthquake with the epicentre offshore, 120 km SE of Yogyakarta City Centre, severely damaged 2,200 houses in 1937. Also, in 1943, an earthquake with the epicentre near the previous earthquake caused a death toll of 213 and destroyed 15,275 houses. One of the most recent strong earthquakes that caused several damages in Yogyakarta and Central Java Province including study area was on May 27, 2006, at 5:54 a.m. local time, with a moment magnitude of 6.4 which occurred 20 km SE of Yogyakarta City. The earthquake epicentre's as estimated by the United State Geological Survey (USGS) was 7.962 0 S -110.458 0 E, and the focal depth was 10 km.
The tremors lasted for 52 s. The earthquake directly affected the Yogyakarta Province and Central Java area The damage distribution was well correlated with the epicentre and presumed fault region. All regions nearby the Opak fault suffered a high death toll and massive damage including the study area. Walter et al. (2007) and Daryono (2011) concluded that damage pattern was highly controlled by the amplification factors. The earthquake amplification phenomena occurred in Young Merapi Volcanic Deposits which is located along the Opak River. This area was a densely populated area. Several sub districts that were within the amplification zone and had severe damage were Imogiri, Bambanglipuro, Jetis, Pleret, Banguntapan, Prambanan, Gantiwarno, Wedi and Bayat Sub Districts.
The collapsed buildings due to the ground shaking were mostly non-engineering buildings consisting of one or two storey houses, shops, religious and school buildings. Most of the collapsed residential buildings were non-engineering masonry houses. They were one storey tall with unreinforced clay, brick or block masonry in cement or lime mortar and no particular connection frame between timber roof and the walls. According to Elnashai, et al. (2006) these houses commonly had 8-20 m 2 of plan dimensions and 2.5 to 3 m of average heights. They usually used half brick masonry infill walls as a reinforced concrete framing. Damage was also found in reinforced masonry buildings. However, it was only moderate. The earthquake also affected several commercial buildings in the research area. Most of them suffered from slight damage until collapse. Those buildings were engineered multi-storey reinforced concrete structures (Elnashai et al., 2006). Examples of the building damages can be seen in Fig. 4 below

An overview of existing building condition
The description of damaged buildings above may reveal that the buildings had been constructed with lack of anticipation on seismic events as well as quality of construction (Elnashai et al., 2006). Generally, buildings consist of non-structural and structural elements. Nonstructural elements refer to the components of the building that are not considered to support its selfweight as well as the external forces of the building. On the other hand, the structural elements refer to the components of the building that are carefully designed to hold and distribute the forces acting on the building and transfer every single force to the ground continuously without any significant deformation. Each part of the structural elements of the building collaborates with each other to form a structural system that guarantees the utility of a building. Rapid Visual Screening of Building for Potential Seismic Hazards (FEMA P-154) classifies building structure into 15 categories. Classification takes into account the combination between structural element material and structural load-bearing system. The building classification in this study is also in line with the FEMA P-154. Based on the previous research, (Aswandono, 2011 andSaputra, 2012), four categories of structure found in the study area, i.e.: URM -unreinforced masonry bearingwall, W 1 -light wood-frame residential and commercial buildings smaller than or equal to 5,000 ft 2~4 60 m 2 , RM 1 -reinforced masonry buildings with flexible floor and roof diaphragms, and RM 2 -reinforced masonry buildings with rigid floor and roof diaphragms (FEMA, 2015).
Within the URM building, the unreinforced fired clay brick masonry walls built in cement mortar act as the main load-bearing structure, while the wooden or bamboo roof as upper structure which also carries the load of clay tiles thereon. The upper structure is laid directly on the walls as the main structure with no special connection. This type of building works best in carrying gravity force, but it is very poor in resisting lateral forces. Consequently, URM building will experience severe damage or collapse under earthquakes' force (Wijanto & Sinha, 2003). URM is generally categorized as a nonengineered construction which is built under the traditional construction practice on an empirical basis.
Light wood-frame buildings (W 1 ) in the study area are represented by traditional buildings which are the result of a long empirical process of trial and error in establishing a unique type structure (Prihatmaji et al. 2012). Traditional wood-frame structures were designed to deal with environmental threats to provide safety and convenience conditions for the people inside. The building is mainly built with wood which relatively results in significant reduction on the weight of building compared to that of different types of material. Since the seismic forces are proportional to the total weight of the building, the heavier building will suffer from greater horizontal force which might cause damage and even the collapse of the building. Another essential characteristic of light wood-frame structure is the common use of special connections which provide both strength and flexibility in responding forces. The connections will assure that the structural elements will be kept connected under a certain magnitude of forces and structural displacement in respond of the forces. Both light-weight and wooden connections provide structural flexibility to deal with seismic forces. However, the special characteristics of W 1 -type still do not guarantee that the building would be safe under an earthquake shake. According to Prihatmaji et al. (2012), there are three common damages causing the total of the Joglo Structure. First, the failure of the connection between the wood columns and the foundation due to the lacked of columnfoundation especially in the side structure. Additionally, the column can easily slipped out from the bas stone when the column's leg decayed. Second, the failur of connection between lower and upper beam because the damage of both column and beam after receiving the lateral force of the earthquake. Third, the roof strucuture failure which is triggered by the detachment of the roof rafter from the main beam due to the deformation, unstable, or collpases of the outer structure.
Reinforced masonry buildings (RM) which use wood diaphragm and RM 2 which often use precast concrete  (Daryono, 2011) diaphragm, rely on the perimeter bearing wall to hold acting forces. Generally, reinforced masonry building types have relative stiffer structure which works well under moderate earthquake forces as far as the building code implementation and adequate construction practice and control are well executed (FEMA, 2015). However, the evidence of post-disaster building condition in the study area reveals that many buildings were not constructed to meet the requirement of earthquake-resistant building as regulated in several building code. The building cost is considered as the main problem for lowerincome people at the study area. A considerable number of buildings were constructed under poor construction practice by low-cost, incompetent traditional builders, without compliance with any safety rules and regulations and sufficient knowledge in construction science and earthquake engineering. In addition, the choice of lowquality material also significantly contributes to the poor building quality.

Demographic condition
According to the civil registration data of Bantul District, the population of Pleret Sub District in 2013 was 45,316 people with an average population growth of 1.99%. The sex ratio of Pleret Sub District was 100.34, which means that the numbers of females and males in the population are relatively balanced. Wonokromo and Pleret Villages are the densest areas in Pleret Sub District with the population density being 3,231 and 2,953 people per km 2 respectively. Whereas, the lowest population density (1,003 people per km 2 ) is Wonolelo Village, which is located in the eastern part of Pleret. The demographic structure in this sub district displays a near-stationary pyramid. The 2013 demographic data showed that this area has stable growth and was dominated by productive people between15 and 49 years of age. This age group is relatively less vulnerable than the other ones. On the other hand, approximately 43% of the total population in Pleret Sub District is more susceptible to earthquake casualties due to their ages (>14 and < 49 years old) (Fig. 5).
Based on the latest uploaded data (Wonokromo Pleret Village, 2014;Segoroyoso Village, 2014;Bawuran Village, 2014;and Wonolelo Village 2014), there are five major types of occupation in Pleret Sub District, i.e., the casual worker, student, unemployment, entrepreneur and farm worker. Almost 20% of the total population in Pleret Sub District are casual workers, 13.55% are students, 13.39% work as entrepreneurs, 11.37% are farm workers, and 13.10% are unemployed. The detail of occupation types in Pleret Sub District can be seen in Table 4. This type of data, i.e. type of population occupation of each village, are imperative in vulnerability studies especially for modelling human vulnerability based on spatiotemporal distribution (Freire et al. 2013).

Method
Three main analyses were applied in this study; first, the visual interpretation of geological features and land use; second, the probabilistic analysis of building damage; and third, the population distribution analysis. The Landsat 7 ETM+ and Quickbird 2012 satellite imagery were utilized to generate the geology and land use map in the areas of interest. The probabilistic analysis was done by using the building damage dataset which is obtained from the rapid survey of earthquake damage after the Yogyakarta earthquake, 2006. In addition, to disaggregate the numbers of population into particular land use units, this study used the data of local livelihood obtained from the local government at village level. The final results of this study is the multi vulnerability model which was combination between the damage probability of building block and the population distribution model. The multi-vulnerability model was visualised in a map form with the scale of 1: 30,000.

Visual interpretation
A series of visual interpretation elements such as colour/tone, size, shape, texture, pattern, shadow, and association were used to distinguish the types of land use and the geological features. The first stage is the land use data extraction through the visual interpretation of Quickbird imagery. The USGS land use and land cover classification system was used to divide the land use into several major groups. Those major groups are urban or built-up land; agricultural land; rangeland; forestland; water; wetlands; barren land; and managed wetland. The main difficulties of land use interpretation in the research area are to differentiate the residential and non-residential buildings such as schools, hospitals, governmental offices, mosques, light industrial buildings, and traditional markets. The residential buildings usually have a regular shape (rectangle) and the building size lies between 21 m 2 and 200 m 2 . The school buildings have irregular shapes (e.g. shape of letter "O", "L", "U", "T", and "H") and are mostly located near playgrounds or open areas. The governmental offices and hospitals are usually comprised of several smaller building units with regular shapes, close to each other, and have the same colour of roof material. Light industrial buildings and traditional markets usually have a large building size (>200 m 2 ) and they use asbestos or zinc roofs which is indicated by the white or grey roof colour in Quickbird Imagery. The mosques in the research area have a specific roof shape ("limasan"), and are located near the public cemetery (Fig. 10). The commercial strip development refers to the commercial activity developed along the main road such as shops, retail stores, fast food services, gas stations, etc. In general, the commercial strip has the same characteristics as the residential buildings.
The key criterion to differentiate the commercial strip from the settlement is the proximity to the main road. Most shops, retail stores, fast food services and other similar goods and services are located along the main road and vehicular transportation road. The next step is to identify the geological characteristics through visual interpretation of Landsat 7 ETM+. The visual interpretation elements of colour/tone, shape, pattern, texture, and association were used to extract the geological information. Moreover, a 1: 100,000 geology map of Yogyakarta and a topographic map with a 12.5 m contour interval were also used to support the interpretation process. The geological features were obtained from the visual interpretation process of Landsat 7 ETM + imagery with RGB colour composite of 4,5,7. The Band 4 or near infrared (0.7-0.9μm) is good for determining water or land surfaces because almost all radiation in this wavelength range is absorbed by water. The band 5 or middle infrared (1.55-1.75 μm) is very sensitive to moisture and very suitable to monitor vegetation. The band 7 or middle infrared (2.08-2.35 μm) is good for soil and geological mapping. By using the elements of interpretation and topographic information such as relief and slope, a detailed unit of geological unit can be produced to complement the geology map. Finally, a fieldwork has been conducted to justify qualitatively the accuracy of geological interpretation. At least 53 observation locations have been identified and verified in term of the geological unit characteristic. A physical characteristic identification of the outcrops or surface sediment have been conducted in the field. The rock and formation classification referred to the Geological map of Yogyakarta scale 1: 100,000.

Probabilistic analysis of building damage
A statistical analysis was applied to calculate the probability of building damage based on the damage building inventory data. An ordinal logistic regression was used to generate a probabilistic model of building damages. This modelling method was used because the dependent variable (damage level) is a 3-level ordinal data. The number "1" means not damaged, 2 is moderately damaged, and 3 is collapsed. These damage levels were obtained from rapid damage assessment after the Yogyakarta earthquake 2006 (Kerle, 2010). The independent variables were the structure types (reinforced masonry, unreinforced masonry, and wood); roof material (asbestos or zinc, cement tile, clay tile, concrete slap, and thatch); distance (within 8 km, 8.1-10 km, 10.1-12 km, and 12.1-15 km); geology (Qa, Qmi, Tmn, and Tmse). However, because of the low frequency, thatch roof variable was taken out from the analysis. The summary of data used in the statistical analysis is provided in Table 5 below.
Before running the above ordinal logistic regression model in SPSS software, a benchmark or base value was defined for each variable (Leech et al. 2005). The building structure in this case is represented by three variables namely Structure = 0; Structure = 1; and Structure = 2. It is necessary to define one variable to be a base value among those three variables in order to help with the interpretation of the results. Reinforced masonry (Structure = 2) is defined as the base value for variable 1 and 2. Roof = 3 is defined as the base value for variables 4, 5, and 6. Distance = 4 is defined as the base value for variables 8, 9, and 10 and the last benchmark is Geology = 3 for variables 12, 13, and 14. The last step was to implement the probabilistic model into the 2012 building foot print data and to model the spatio-temporal patterns of residential building damage. A dasymetric mapping was applied in this model to disaggregate the numbers of population in space and time. The summary of the data used in this study can be seen in the Table 6 below.

Population distribution analysis
The population distribution in this study refers to the population distribution based on the community habits during the weekdays and holidays. The basic technique used was dasymetric mapping. The land use and the type of occupation were used to increase the accuracy of population distribution model. The adopted concept of dasymetric mapping can be seen in the Fig. 6 below.
The dasymetric model was conducted in several steps; first, the land use data was divided into binary values of population distribution (non-populated area and populated area). Second, the population density of each land use unit was estimated. Third, a dynamic population distribution based on the community habits was made. The latter is very important, because human distribution in space varies largely during daily life (Freire et al. 2013).
Several equations are used to disaggregate the numbers of population   The general framework of the research can be seen in Fig. 7 below.

Minimum, maximum, and average scenarios
In order to generalize the building damage probability value at building blocks level, three scenarios (minimum, maximum, and average scenarios) of building damage probability were applied. The illustration of converting the damage probability of each building unit into the damage probability of building block can be seen in the Fig. 7.
The predicted damage category of each building unit was calculated using the equation resulted from the Afterwards, the probabilities of not damaged, moderately damaged, and collapse were calculated using the equations 5, 6, and 7). Then, the separation of building land use from the others was done with the "select by attribute" tools in ArcGis 9.0. Based on the modified Anderson system 2002 (U.S. Geological Survey, 2007), the building land use consists of commercial strip development, educational institution, government centre, hospital, light industrial, mosque, residential area (high, medium, and low density), rural, single unit residential, and traditional market.

Estimating multi-vulnerability (Building damage category and population distribution)
The multi-vulnerability of buildings was determined by adding the building damage category (average scenario) and population distribution class. There are three classes of the numbers of population staying in the particular land use units. Class 1 refers to numbers of population below 100 people, class 2 to numbers of population between 100 and 200 people, and class 3 to numbers of population above 200 people. This classification follows the rule of natural break technique in GIS reclassify tools. The multi-vulnerability value will range between 2 and 6. The values of 2 and 3 are defined as low vulnerability, the value of 3 and 5-6 as medium and high vulnerabilities, respectively. The overall research framework can be seen in the Fig. 8 below.

Visual interpretation of Landsat 7 ETM+ imagery to improve geological map of research area
Based on the visual interpretation of Landsat 7 ETM+, the research area was classified into four geological units. These units have unique characteristic of colour and tone, relief, and texture. Those geological units are Semilir Formation (Tmse), Nglanggran Formation (Tmn), Alluvium (Qa), and Young Merapi Volcano deposit (Qmi). The total area of Tmse in east part of Pleret Sub District is about 945.03 Ha (39.67%). Tmse is characterized with the domination of light blue colour and some random of reddish brown spot. Tmse has very rough texture, especially in the east part area, and is located between undulating slope and step slope areas (4 0 < slope < 35 0 ). Different from Tmse, Tmn has a striking red colour with a few random of light blue spots. This unit is located in the eastern part of the research area and most of Tmn is in the flat area with the highest elevation being 334.7 m above sea level. This area extends to approximate 82.16 Ha. Qa and Qmi are very difficult to distinguish in the interpretation process, because both of them have similar characteristics of morphology, colour or tone, texture, and relief. Therefore, the field investigation was conducted to distinguish both geological units and to check the visual interpretation results. There are some fundamental differences between the geological units derived from visual interpretation and the 1:100,000 geological map of Yogyakarta. Based on the latter, the border between Tmse and Tmn is located along the foot slope of Baturagung Escarpment near the Guyangan Villages, but we found that the border is located in the upper slope of Baturagung Escarpment near the Dlingo Villages. Then, we also found that there are some isolated hills in the north part of research area which is belong to the Tmn units. The main difference was also found in the extents of Qmi. Based on the visual interpretation, Qmi is spread along the left and the right side of Opak River, while the alluvium only occupies the small part of flat area near the hilly area of Semilir Formation. The summary of interpretation results and the differences between both geological maps are shown in Table 7 and Fig. 9.
The rapid field identification was conducted to determine the border between Qmi and Qa. Two approaches were used to determine both geological units. The first approach was to locate the boulder location. Qmi is dominated by fluvial process and produces very well sorted sediment. The upper layer consists of very fine sediment, while the lower layer consists of rough sediment. The boulder indicate that the particular area is influenced by the colluvium sediment which is characterized as bad sorted or mixture sediment from the hilly area nearby. The second approach is the clay content of both Qmi and Qa. Qmi tends to have less clay content being the Merapi Volcano sediment, while the Qa tends to have more clay content as it is the denudation material of the weathering process in the hilly surrounding areas. Therefore, one of the objectives of the rapid field investigation was to locate the traditional brick factory along the Opak River. The brick makers tend to use the best soil which has less clay rather than soil which has more clay. The reason is that the brick will crack during the heating process if there are a lot of clay contents inside the soil. Therefore, the location of traditional brick makers can be used to indicate the Qmi area.
Based on the field observation and investigation, the geological unit's characteristic in the field is 100% the same as the information given in the improved geological map of the research area. The field documentation of rapid field investigation, is shown in Fig. 10.

Land use interpretation
Based on the Land use-Land cover classification system, NJDEP modified Anderson System 2002, there are 29 land use units in the research areas. Those land use units cover 10 units of level 1 i.e., residential, commercial, industrial, transportation, other built up area, agricultural area, confined feeding operation, forestland, stream & canal, and wetland. The land use classification show that the research area is dominated by the shrub land which cover approximately of 42.36% of the total area. The wetland agricultural area or paddy field only cover 25.34% of the total area, while the residential area (high, medium, and low density) cover approximately 3.79%, 5.70% and 3.01%, respectively. Most of them are distributed in the west and middle part of research area which are the Qmi and Qa. Detail information is shown in Fig. 11, and Table 8, while the land use map is shown in Fig. 12.

Statistical results
The main objective of using statistical analysis in this study was to derive a model showing the relationship between the level of damage in residential houses and their characteristics. This study applied the ordinal logistic  (Nurwihastuti et al., 2014). Based on the Table 11, there are three categories of predicted response value namely, Y1, Y2 and Y3. Y1 refers to the not damaged category, Y2 refers to the moderately damaged category, and Y3 refers to the collapsed category. These categories were determined based on the threshold value in Table 11. Therefore, the predicted response values are defined as follow.
The Y*i value can be calculated by applying the regression equation (equation 4) to the dataset of building footprint. There are 17,512 data of building units with 33 combination of building structure, roof material, distance from the epicentre, and the geological type. The example of the calculation can be seen in Table 12.
Based on the Table 12, specific combinations of building attributes (structure, roof, distance, and geology) give a particular damage category (Y*i value). For instance, combination 1 has 22 buildings with the same combination, i.e., wood structure, reinforced masonry structure, located between 8.1 and 10 km from the epicentre, and located in Semilir Formation. This combination gives a value of Y*i equal to 0.915 which means that this building tends to have a low level of damage or not damaged. The probability of collapse can be also calculated by put the Y*i into general equation of ordinal logistic regression. The probability of damage of building attribute number 1 is as follows: -Not damaged probabilistic (Y1) Therefore,  Table 13 shows the complete results of probability calculation of 33 buildings attribute combination in research area. Based on the calculation of building damage probability, the safest building type is the combination of reinforce masonry, asbestos or zinc roof, located between 10.1 and 12 km from the epicentre, and being located on Semilir Formation, with a probability of not damaged equal to 0.85 or probability of collapse equal to 0.07. Only five buildings or 0.028% of the total buildings in the research area have this building attribute combination. The most vulnerable buildings in the research area are the ones with a combination of reinforced masonry structure, clay tile roof, being located between 8.1 and 10 km from the epicentre, and being located on young Merapi volcanic deposits. This building attribute combination gives a value of collapsed probability equal to 0.84 or not damaged probability equal to 0.07. The number of the buildings with a vulnerable combination is 45 buildings or 0.25% of the total buildings. Most of the buildings in research area have a combination of reinforced masonry, clay tile roof, are located more than 12 km from the epicentre, and are located on young Merapi volcanic deposits (37.71%). This building attribute combination gives the collapse probability value of 0.67 and it is predicted as highly damaged building or collapsed.
The building foot print model of building damage probability is shown in the Fig. 13. The scenario of minimum, maximum, and average damage was also applied to generalise the building foot print model into the building block probability. The conversion method follows the illustration shown in Fig. 7 above. We found that the pattern of probability value between minimum, maximum, and average scenarios were slightly different. Based on the minimum scenario, the building block located on the western part of Opak river has a probability of collapse between 0.22 and 0.81, while based on the maximum and average scenarios, probability of collapse is the same, about 0.29-0.81. The rest of the building located in eastern part of Opak river has collapsed probability value about 0.07-0.84; 0.08-0.84; and 0.08-0.84 of minimum, maximum, and average scenarios, respectively. The complete map of the collapsed probability and the predicted building damage categories under minimum, maximum, and average scenarios are shown in the Fig. 14 and Fig. 15.

Population distribution
The population distribution analysis was applied to disaggregate the number of people at Sub District level into a particular land use types within the area. The temporal analysis of population distribution based on the types of occupation was also conducted in this stage. The main concept of this analysis is to breakdown the population data using dasymetric technique. The main equation used follows the equations 1, 2 and 3 above. The characteristic of livelihood of each village was also used to support the analysis. The result shows that the top 5 occupation types in research area are casual worker, student, unemployment, entrepreneur, and farmer. The percentage of people   Based on the population percentage calculation of each land use type, during the working hour the population is scattered into three major land use types such as     The final result, multi-vulnerability of building blocks, is obtained by combining the population distribution model and the predicted building damage categories. The results ( Fig. 20 and 21) show that the western and middle part of the research area are dominated by moderate and high levels of vulnerability, while the eastern parts are dominated by low levels of vulnerability. The vulnerability patterns of working hours and holiday times are not considerably different. However, the vulnerability level of some settlement blocks are increasing during night and holiday times.

Model validation
Several validation tests have been conducted in this study. First, a validation of geological unit interpretation, second, a validation of building structure, and last, a validation of the predicted building damage categories. The validation of geological units shows that 72 fieldwork locations have the same geological characteristic with the information of geological map used in this study ( Fig. 9  and 10). This result indicate that the interpretation is accurate and it can be used for further analysis.
The building structure validation using confusion matrix analysis shows that the visual interpretation of building structure is highly accurate. The accuracy level of interpretation from the total of 332 buildings sample is 95.18% (Table 14). This result also indicate that the interpretation result of building structure is accurate and it can be used for further analysis.
The last part of the validation process was to calculate the accuracy level of the ordinal logistic regression model. The main calculation was to compare and calculate the matching value between the response variable (building damage) and the predicted damage categories. The result shows that, the ordinal logistic regression model was able to predict precisely the response variable for about 5,795 out of 7,645 building damage data points or about 75.80% of the total. This means that the model can predict accurately the building damage. Additionally, this model can be applied to a dataset of another area of interest. Table 15 below shows an example of the model validation process.

Discussion
The visual interpretation results of Landsat 7 ETM+ produced a 1:40.000 geological map which is slightly different from the existing geological map produced by (Rahardjo et al. 1995). The main distinctions are 1) the border between Semilir Formation and Nglanggran Formation in the Baturagung escarpment, 2) the border  between young volcano deposits of Merapi Volcano and alluvium, 3) and the isolated hill member of Nglanggran formation (G.Gelap) which does not exist in the geological map of Yogyakarta. Based on the interpretation results and field work validation, we found that the border between Semilir and Nglanggran Formation is located on the upper slope of Baturagung, while according to the Geological map of Yogyakarta, the border is located in the middle slope of Baturagung. We found also two isolated hills in the west part of Opak river namely Gunung Gelap. This result is very similar to the Sanjoto (2004) results. They found also two isolated hills in the same location which consist of volcanic breccia, flow breccia, agglomerate, lava and tuff members of Nglanggran Formation (Fig. 9). By combining the Yogyakarta map with the geological map resulted from interpretation process, a better undertsanding of the local geological characteristics can be obtained, which is important for further analysis. The Yogyakarta earthquake in 2006 caused a major impact to the local infrastructure, especially in the housing sector. Scientists believe that the damage pattern is controlled by certain factors. Qualitative, semi-qualitative, and quantitative methods have been developed to model the controlling factors. Some studies have been conducted to get a better understanding about the damage pattern caused by the Yogyakarta earthquake in 2006, such as the study conducted by Nurwihastuti, et al. (2014). The study concludes that there is connection between earthquake damage pattern, geomorphological characteristics and sub-surface characteristics. They found that the severe damage tends to occur in areas that have deep basement, low gravity anomaly, thick sediment, and unconsolidated material, while the slight damages tend to occur in areas that have shallow basement, high gravity anomaly, thin sediment, and consolidated material. These results are also in accordance with the previous study conducted by (Daryono, 2011). They found that the high vulnerability of seismic index (kg), tends to give a major impact on buildings and infrastructures. This area covers the west part of research area which consists of young volcanic deposits of Merapi volcano and alluvium.
The aforementioned studies allowed for a better understanding on how to characterize the pattern of damage caused by The Yogyakarta earthquake in 2006. However, the scope of those studies, explain only the physical characteristics. Therefore, the probabilistic modelling for estimating earthquake vulnerability was conducted to provide further understanding of the problem. For this, we used an integrated analysis methodology using remote sensing, GIS analysis and statistical modelling. A combination of physical and building characteristics was used as the independent variables. These were the type of building structure, roof type, distance to the epicentre, and geological characteristics. As mentioned in the results chapter, the model was able to explain the probability of building damage with a 75.80% success rate. This model shows that the combination of reinforced masonry structure, clay tile roof, being located between 8.1 and 10 km from the epicentre, and being located on young Merapi volcanic deposits defines the most vulnerable type of building in the study area.
Most of the existing reinforced masonry buildings in the research area consist of three main structure components such as the building foundation, reinforced concrete columns, and roof structure. Sometimes the buildings have no ground beam and roof beam, which implies lack of connection between roof, concrete columns, and roof structure. Moreover, they have no diagonal concrete columns which are useful to reduce the tensile stresses during a seismic vibration. The use of clay tile in this case, will increase also the probability of a collapse. As shown in the Table 11, clay tile roof has a value of 1.634, which means that clay gives the higher probability of a collapse (equal to 1.634) compared to the roof benchmark (concrete slap). In fact, it has a heavier weight (1.49 kg per square meter) than asbestos and concrete tile or cement tile. The distance from the epicentre also contributes to the probability of collapse. According to the model, the further from the epicentre, the lower the probability of damage.
The young volcanic deposits of Merapi Volcano give also the highest contribution to the probability of building collapse among the other geological units. This unit is characterized as deep basement, thick sediment, unconsolidated material and having abundant shallow groundwater. The area has high possibility of ground amplification occurrence (Daryono, 2011) (Nurwihastuti et al., 2014) (Koseki, et al., 2007) (Pandita et al. 2016). Parallel with the result of Nurwihastuti, et al. (2014), most collapsed buildings due to The Yogyakarta earthquake were located in fluvial landforms, in this case the fluvial plain that consist of Young Merapi volcano deposits. The probabilistic model also shows that the usage of ordinary reinforced masonry and clay tile is not suitable in this area. it estimates that the combination of structure and roof material types in this zone gives a probability of collpase between 0.51 and 0.84. Thus, it is recommended to build a house by using the building structure of wood and light steel roof or reinforeced masonry with a good connection between foundation, concrete column, and roof structre.
However, the combination of reinforced masonry, asbestos zinc, and being located in Semilir Formation gives the lowest probability of collapse in the model. This is because the Semilir Formation is a massive volcanic rock which consists of interbedded tuff breccia, pumice breccia, dacite tuff, andesite tuff, and tuffaceous clay stone. This rock formation can reduce the ground amplification caused by an earthquake. In line with Daryono (2011) results, the lowest value of seismic index vulnerability is located in the structural and denudational landform member of Semilir Formation.
This probabilistic method is very suitable for developing countries like Indonesia where a national standard of building code and map of building damage prediction are less available. It provides a comprehensive and simple probabilistic modelling of building damage which can be used easily by scientists and decision makers at both local and regional levels to improve the current disaster and risk management reduction programs. By using this method, an accurate results of preliminary damage prediction up to 75% level of correctness can be provided with relatively limited time and cost.
There are several limitations that need to be taken into consideration for future research and improvement of the model. Within this study, only few variables were used in the ordinal logistic regression. Only about 33.20% of the model can explain the variance. There still is about 66.80% that might be can be explained by other factors. However, the model fit information (Table 9) showed the significant result. It means, this model is appropriate and can be used for further analysis, although the model only can explain of 33.20% of variance. Thus, the advanced study with more additional independent variable is needed to improve the results and to gain higher level of variance explanation. In term of engineering purposes, a detailed investigation of residential building is absolutely need to be conducted to determine the safety factor and vulnerability of the buildings. Based on the main requirement on how to build the safer houses, which was recommended by the Indonesian Ministry of Public Works, there are some aspects that determine the vulnerability of the buildings. First, a good quality of construction material. Second, appropriate structural existence and dimension. Third, a good connection within the main structural elements, and fourth, a good processing quality. Additionally, for the Joglo houses, the quality of the timber (age, type of the material, quality of the carpentry, and maintenance) also determines the level of earthquake vulnerability. Thus, for further research the combination of this probabilistic model and the detailed site investigation are very important to support the existing building code.
The others research that can be carried out in the future is the study about the outcrop or geological mapping activities. By improving the geological map of Yogyakarta 1: 100,000; the better understanding of geological characteristic and surface material behaviour during the earthquake can be obtained. Moreover, the outcrop study is able to reconstruct the minor fault which was unidentified from the geology map of Yogyakarta. By revealing the configuration of this minor faults, the collateral damages due to amplification of the fault can be avoided.
Further research also needs to be done to model the population distribution in space and time. This study only used a percentage analysis to disaggregate the population data at a Village level into particular land use types. This model gives only a general description about how the population spreads in time based on their livelihood type. A good local knowledge and local community habit understanding are absolutely necessary to generate this model. Therefore, primary data of population distribution are very important to explain the patterns of local community behaviour. However, this model provides a valuable information for the local government to support decision making concerning earthquake hazards and risk management.

Conclusion
Earthquake hazard becomes a serious problem in Java and particularly in Southern Java as a result of the subduction zone along the south part of the area, about 320 km south of Yogyakarta. In addition, there is a growing demand for residential buildings due to the rapid population growth. This condition might lead to a higher level of vulnerability to earthquake disaster. This phenomenon has triggered scientists to further study and model the corresponding earthquake hazards, vulnerability, and make risk assessments to help reduce disaster risks. Several  studies have been conducted to analyse the damage patterns caused by The Yogyakarta earthquake in 2006. By combining the geophysical and building characteristics, this study provides an accurate prediction about 75.81% of building damage probability. Based on the results of an ordinal logistic regression model, the reinforced masonry house with clay tile roof, it is located between 8.1 and 10 km from the epicentre of the 2006 earthquake and was built on young Merapi volcanic deposits has higher probability of building damaged. On the contrary, the combination of reinforced masonry, asbestos or zinc roof, being located in Semilir Formation and far (>12 Km) from the 2006 epicentre has lower probability of buildings damaged. Additionally, the results explain the estimation or prediction of building vulnerability based on the building damaged of the Yogyakarta earthquake 2006. This study is also suitable for preliminary study at the regional scale. Thus, the site investigation still needs to be conducted for the future research to determine the safety and vulnerability of residential building.