ISSN 1880-8468
Technical Report ofInternational Development Engineering
国際開発工学報告

TRIDE-2013-02

March 1, 2013

Abstracts of Bachelor Theses
Presented in March 2013

Department of International Development Engineering,
Graduate School of Science and Engineering,Tokyo Institute of Technology
http://www.ide.titech.ac.jp/TR

Preface

Bachelor theses of Department of International Development Engineering, Tokyo Institute ofTechnology were presented successfully on August 3, 2012, March 1, 2013. This technicalreport consists of the abstracts of those theses.

ⅰTechnical Report of International Development Engineering
TRIDE-2013-03b
Table of Contents(Completing in September 2012)
E-Learning System Based on Wisdom of Crowds on Android and its
Vietnamese Version........................................................................................................ Huynh Dao Hoang

1

Titan Aerocapture for Saturn Orbit Insertion.................................................................................................................. Seiya FURUKI

3

Freight Mode Choice from Vietnam to Thailand................................................................................................................. Vu Mai Phuong

5(Completing in March 2013)
Fundamental Study on Proposed Efficient Repair with Combination ofElectro-chemical Methods.......................................................................................................... Daisuke KUNIOKA

7

Adhesion Phenomenon between Rigid Body and Tapering Elastic Beam.................................................................................................................... Boqing Wang

9

The Research of the Reinforcement Learning for Humanoid
Robot’s walk control......................................................................................................... Soichi UCHIYAMA

11

Application of TIV for Multi-scale Atmospheric Flows........................................................................................................................ ｱｰﾘﾌ ﾑﾀｽﾞ

13

発電デバイスの最大出力点を追従するエネルギーハーベスティング用昇圧コンバータの試作...................................................................................................................... 浅見立志

15............................................................................................................. Ryuuhei KUDOU

17

Financial Analysis of Vertical Integration and Privatization for Airports
in Hokkaido
Geometrical Shape of Arch Formed by Collapse of Undercut Slope................................................................................................................ Hidehira HIRAI

19

Numerical Analysis of Natural Convection in an Arch Vesselby a Spectral Finite Difference Scheme.................................................................................................................. Keita MITANI

ⅱ

21

ドップラーライダーで観測された流れ場の分類と特性の把握...................................................................................................................... 有場次郎

23................................................................................................................. Huiting WANG

25

CoOx-CeO2 Mixed Oxide Catalysts for N2O Decomposition
Synthesis of Zeolites from Lake Sludge as Ammonium Adsorbent Materials............................................................................................................... Wenjing WANG

27

衛星雲画像を用いた海風前線の推定および都市が海風に与える影響.................................................................................................................. 大久保洸平

29Measurement of Audible Sound due to Crushing Progression of Snail Fossils....................................................................................................................... Wenjie WU

31

Image Coding with 1D Discrete Wavelet Transform using..................................................................................................................... Shiyu ZHOU33

A Pareto Frontier Analysis of Household Budget and Environment for Energy
Saving Goods.............................................................................................................. Ryoya SUEHARA

35

Study on the Propagation of the Tidal Waves in the Mekong River.......................................................................................................... Yuto TAKAHASHI

37

Production of Monomers for Biodegradable Plastics with Lignin as Raw
Material............................................................................................................... Ngoc Linh TRAN

39

Tsunami Vulnerability Assessment in Tokyo Bay Using a Numerical
Simulation............................................................................................................ Shouhei NAGASE

41

Carbonation Phenomena and Corrosion Behavior of Stainless Steel Bar in the
Mortar Mixed with Seawater......................................................................................................... Tomohiro NAGATA

43

Image Compression by using Vector Quantization and Vector-Embedded
Karhunen-Loeve Transform.................................................................................................................... Kiung PARK

45

Development of Millimeter-wave Radio Channel Sounder............................................................................................................ Pham Kinh HUNG

47

Identification of Riboflavin Producing Yeast in the Fermentation Process for
Puer Tea............................................................................................................... Ryutaro BESSHO

ⅲ

49

大域的射影変換を用いたパターンマッチングに関する研究.............................................................................................................................. 苗春亭

51

Effect of Beam’s Shape on Adhesion between a Multi-beam Structure andRough Surface........................................................................................................... Takehiko YAZAKI

53

Influences of Atmospheric Conditions to EHD Thruster for Planetary
Explorations.........................................................................................................................Tianyue LI

55

Comparative Analysis of Aviation Passenger Attributes in Vietnam..................................................................................................................... Le Yen Lan

57

Prediction of Indoor MIMO Channel for Specific Systems using Ray Tracing
Simulation...............................................................................................................チョウ

ⅳ

カツシ

59

E-Learning System Based on Wisdom of Crowds on Android and
its Vietnamese Version
Students Number: 08 20565 Name: Huynh Dao Hoang Nam Supervisor: Yukihiko Yamasita
July 31, 20121

Introduction

E-learning (consisting of electronic learning) comprises all forms of electronically supported learningand teaching. The information and communication
systems, whether networked or not, serve as specificmedia to implement the learning process. In recent
years, with the increase of the importance with elearning, not only in Japan or Vietnam, but also
in many countries, more and more people study bye-learning. However we have been facing with large
problems in the creation or management of teaching materials.
The tobe (Text Oriented Bi-stream Explanation)system may solve those problems with wisdom of
crowds. However, PC is necessary to learn by thetobe system. As a result, we can not use it when
we are on the way to school or company.
I focus on the increase of use of the tablet andsmartphone and the highest share of Android OS.
N. Takemura made a tobeML player on Android
devices [1]. In this paper, I improve the system bymultilingualization and pointing function to highlight an important area.

2

Figure 1: Expositor

3

Tobe Player on Android

There are two main reasons that I use Android devices for my research.
1. It is the most popular mobile OS in the world.
Its share 79.6 % in Japan[2] and 28.7 % in
Vietnam[3].
2. Its applications can be developed by Java.

Table 1:Worldwide Share of OS Smartphone in 2012Year[4]

No
1
2
3
4
5
6

Tobe System

Company
AndroidiOS
Symbian
RIM
Bada
Microsoft
OthersTotal

4Q12 Units
81,067.4
33,120.5
12,466.9
9,939.3
3,842.2
2,712.5
1,242.9
144,391.7

Share(%)
56.1
22.9
8.6
6.9
2.7
1.9
0.9
100.0

4Q11 Units
36,350.1
16,883.2
27,598.5
13,004.0
1,862.2
2,582.1
1,495.0
99,775.0

Share(%)
36.4
16.9
27.7
13.0
1.9
2.6
1.5
100.0

Tobe system is an e-learning system for which anyone can make the contents. The contents aredescribed using tobeML (tobe markup language).
The process of creating and studying contents is as
In order to realize e-learning in mobile environfollows:
ment, I develop a tobeML player on Android devices.
1. Creators all over the world describe thetobeML documents using text editor on the
3.1 Text To Speechweb.
We use Text To Speech (TTS) to generate speech
2. The tobeML documents are stored in thefrom text. This system sends text data to Google
server.
Translate, and gets mp3 data from there. We can
3. Users receive the tobeML documents from the play and save that mp3 date on the device. Furtherserver, and load them into Expositor (Figure more, we can replay the data even when we cannot1), that is a lecture style player of tobeML doc- connect device to Internet. TTS in Google Translate can deal with a lot of languages. Android hasuments.
an own TTS system. However, its quality is not so
The format of tobeML documents is based on high and the number of languages is fewer. ThereXML. A virtual teacher in Expositor can point an fore, I used Google Translate for my system. In
important area in text or figure and read the ex- my research, I conducted experiments in English,
planation.
Japanese, and Vietnamese.
11

3.2

Structure of the player

1. Install this application.

In Figure 2, we can review the structure of the main
window of player.

2. Set the resources (texts, texts for speech, images, coordinations of lines)
3. Start and choose a language, and a content.
4. Push the button to play, back, or skip.
5. Push ’Language’ button if you want to change
the language.
6. Delete mp3 file if you want to renew the content of that file.
In the language selection, if you click “OK” button without selection, “Please select a language” isdisplayed in the bottom.

Figure 2: Display of this System

1. Image Area: Show the image to explain the 4 Evaluation Experimenttext.
I had taken the survey to evaluate the proposed
2. Teacher: Virtual teacher can point an imporsystem by 8 subjects. Result is as follows:tant area in image.
• Good point:
3. Text Area: Show the text.
– With using sounds, pointing and multi4. Play Button: Start speech.
lingualization, it is easy understanding.
– We can change language at the same
slide.

5. Back Button: Return to the previous slide.
6. Next Button: Skip to the next slide.

• Points to improve:

7. Language Button: Change the language.

– Pronunciation of Vietnamese is not sogood.

In this program, we use some important methods
and resources as follows:– The size of text or image should be increased.

• onCreate(): This method is called at first when
the application starts.

– The system should proceed to the next
slide automatically.

• onClick(): This method is called when a button is clicked.
5

Conclusions• getHttpRequestFile(): This method is called I have discussed on tobeML player on Android thatwhen this application sends the text to Google is an e-Learning system using wisdom of crowds,
Translate. Its return is an mp3 file.
and I have developed it on Android. With using
images, texts and sounds, and pointing an impor• setLocale(): This method is called to specify
tant area, it enabled easy understanding comparing
the locale of the response.
to text only or English only system. We totally can• invalidate(): If this method is called, onDraw use this system not only for e-learning but also formethod will be called again. In this applica- the manual of products, sharing of the information,
tion, it is used to reload the language in the or distance training.
same page when we want to change language.
• Toast.makeText(): A toast notification is aReferences
message that pops up on the surface of the[1] N. Takemura, “E-Learning System based on wisdom of crowds on Android”,
window.
Master thesis, Tokyo Tech, 2012 .

• Canvas.Drawline(): Draw a line segment withthe specified start and stop points in (x,y) coordinates, using the specified paint. In this application, I use to draw the pointer of teacher.
The process to use the proposed system is as follows:22[2]

Gartner News, http://www.gartner.com/,[3]

VietNam marketing, http://marketing.24h.com.vn,[4]

ITMedia, http://www.itmedia.co.jp/,[5]Google Translate, http://translate.google.co.jp/,[6]Android Learning, http://developer.android.com/,[7]

Jakarta Project, http://www.jajakarta.org/.

Titan Aerocapture for Saturn Orbit Insertion
Student Number: 08B21381 Name: Seiya Furuki

Supervisor: Yoshihiro Mochimaru

1. Introduction

2. Flight trajectory calculation method

Saturn and its moons have many mysterious points.
Particularly, Titan, one of Saturn’s moons, is the only one ,except for Earth, that has a dense atmosphere and a dynamic
climate system in our solar system.
In order to observe a planet from an orbit around it, aspacecraft is required to reduce its speed from hyperbolic.
Conventionally, thrusters are used to reduce the speed. However,
the method requires much extra propellant. Since specificimpulse of the chemical rocket is not so large, the ratio of
propellant mass to total mass becomes large, which leads toreduction of the payload mass.
For planets having atmosphere, it is possible to use a techniquecalled aerocapture, shown in Fig. 1. Aerocapture is a techniqueto reduce a spacecraft speed by aerodynamic drag. This
technique could drastically reduce to required propellant mass.
However, the spacecraft during atmospheric entry receives
severe aerodynamic heating. Moreover, the range of entry angle
for aerocapture is very narrow, and accurate control of flighttrajectory is required.
As for Saturn, it is possible to reduce a spacecraft speed by
Titan aerocapture as shown in Fig. 2. Saturn aerocapture fororbit insertion around Saturn requires an acceleration to preventthe spacecraft from atmospheric entry. Titan aerocapture fororbit insertion around Saturn does not requires an the maneuver.
This study compares the Titan aerocapture with Saturnaerocapture for an orbit insertion around Saturn.

2.1 Trajectory during atmospheric entry
Forces on the spacecraft are modeled as shown in Fig. 3.
Equations of motion are given bydr v sin 
(1)
dt

d v
 cos 
dt r

(2)

dv
1
1
 g sin   v 2
dt
2


(3)

d  g v     cos 
dt  v r 

(4)



m
CD S

(5)

where 𝑣 is the velocity of spacecraft, 𝜃 is the flight path angle,
𝑟 is the distance from the center of planet, 𝜌 is the atmosphericdensity, 𝑚 is total mass of the spacecraft, 𝑆 is the frontal area,
𝐶𝐷 is drag coefficient, and 𝛽 is the ballistic coefficient. The
equations of motion are calculated by using Runnge-Kutta
method.

probe
𝜃Interplanetary orbit

Target orbit

Entry

φ
Atmosphere

𝑃𝑙𝑎𝑛𝑒𝑡Acceleration

SaturnFig. 3 Trajectory calculation model during atmospheric entry
Exit

2.2 Aerodynamic heating estimation
The heating rate 𝑞̇ W/m2 are given from the folloing
equations
Saturn aerocapture

Fig. 1 Aerocapture
Target orbitq  0.25 
Exit
Titan’s orbit

Saturn v3(6)

Atmosphere

Titan aerocapture

Titan

Titan


RN

q  1.35  10 4 

Entry


RN

 v 3.04(7)

where 𝑣 is the velocity of spacecraft, 𝜌 is the atmosphericdensity, and 𝑅𝑁 is the nose radius.

Magnified view
Interplanetary orbit

Fig. 2 Titan aerocapture for orbit insertion around Saturn33. Flight trajectory calculation results
3.1 Saturn aerocapture for Saturn orbit insertion
Figure 4 shows the peat heat flux for various entry velocities
and entry angles. The ballistic coefficient is10.0 kg/m2 . The
aerocapture is not available at entry conditions in black area in
Fig. 4.

Over deceleration

Over deceleration
Under deceleration

Under deceleration
Fig.7 Peak heat flux at Titan aerocapture (𝛾 =180.0 deg)

Table 1 Comparison between Saturn and Titan aerocapture
SaturnTitan
0.36~11.0
0.03~8.28
Aerodynamic
heating
MW/m2
MW/m2
Aerocapture
0.26~0.73 deg
0.25~3.2 deg
corridorFig. 4 Peak heat flux at Saturn aerocapture

3.2 Titan aerocapture for Saturn orbit insertion
An angle between the velocity vector of a probe and therevolution velocity vector of Titan is 𝛾. Figures 5-7 show thepeak heat flux for various entry velocities and entry angles. The
ballistic coefficient is 10.0 kg/m2 . The aerocapture is notavailable at entry conditions in black area in Figs. 5-7. The
aerodynamic heating increases with 𝛾 . The range of entry
corridor increases with 𝛾.
Table 1 shows comparison between Saturn and Titan
aerocapture for orbit insertion around Saturn.

4. Conclusion
Titan aerocapture for Saturn orbit insertion has advantages,
compared with Saturn aerocapture, in the aerodynamic heating
during atmospheric entry and in the range of atmospheric entry
corridor. Additionally, Titan aerocapture for Saturn orbit
insertion does not require an acceleration at apoapsis, so that the
spacecraft does not need a propulsion system.

Over deceleration

References[1] Hiroshi Hirayama, “The slowdown of Saturn probe using the
Titan atmosphere”, Bachelor Thesis, 1991
[2] Eiko Nakazawa, “Feasibility Study on Low Bllistic
Coefficient Aerocapture Vehicle with Membrane Decelerator for
Planetary Exploration”, Master Thesis, Tokyo University, 2003
[3] Kumiko Nakamura, “Feasibility of the low ballisticcoefficient aerocapture satellite in Saturn”, Master Thesis,
Tokyo University, 2006
[4] M. E. Tauber, and J. V. Bowles, “Use of Atmospheric
Braking During Mars Missions”, Journal of Spacecraft andRockets, vol. 27, Sept-Oct, 1990, p. 514-521[5] Toshio Kawasaki, “Thieory and technology or astronautics”,
chijinnsyokann, 1986

Under deceleration

Fig.5 Peak heat flux at Titan aerocapture (𝛾 =0.0 deg)Over deceleration

Under deceleration

Fig.6 Peak heat flux at Titan aeroocapture (𝛾 =90.0 deg)4

Freight mode choice from Vietnam to Thailand
School register: 09B20278, Name: Vu Mai Phuong, Supervisor: Dr. Shinya Hanaoka
December, 2011was collected from Vietnam generaladministration of customs. It has 14,639 samples andeach sample includes: custom, border, mode and item.
In Vietnam, products are declared in the custom,transported through the border and to Thailand.

1. INTRODUCTION
ASEAN is one of the largest export markets ofVietnam. According to Vietnam general statistics(2011), about 14% of export turnover is generated
from ASEAN export.
Sea transport is the most popular transport totransfer cargo from Vietnam to Thailand. Road
transport does not gain sufficient attraction. One
reason of low popularity of road transport is crossing
third country border. This situation is changing withthe opening of East-West Economic Corridor(EWEC). EWEC connects Vietnam, Laos, Thailand
and Myanmar (Figure 1).

(2) Estimation
A SP survey with assumption that there are tworoutes to transport a container cargo from Hanoi (orNorthern Vietnam) to Bangkok. Two routes arediffered in terms of mode (road and sea), shipment
expense per FEU including customs, average
shipment time and customs procedure. Actual, it takes
3 days and 3800$/FEU to transport cargo by road, 10days and 2000$/FEU to transport cargo by sea. Costand time of these routes of SP survey are designed
based on these values.
Contact with 16 companies and request them tochoose one route which they prefer to transport.
(3) Qualitative survey
Interview 14 companies to answer these below
questions:(a) In mode choice from Vietnam to Thailand, is there
any specific factors more than time and cost?
(b) Why do you choose road rather than sea?
(c) Why do not you choose road rather than sea?
(d) How to make shippers change their choice fromsea transport to road transport?
(e) What kind of goods do you transport by mainly
road or mainly sea? Why?
(f) If you choose road transport, which route do you
choose?
(g) Is there the route of road transport that cross theborder of Lang Son? Why?

EWEC

BK

Figure 1 Vietnam-Thailand route
The cargo transport cost by road from Hanoi,
Vietnam to Bangkok, Thailand is more than twice astransported by sea mode. However, road transport is 3
times faster than sea transport. Therefore, seatransport will retain its advantage in bulk transport,
but because of rapidity, expectation for road is high.
This research aims to find out the reasons why
shippers and logistics companies use route between
road and sea, through data analysis, estimation
designing mode choice model with stated preference(SP) data and qualitative survey.

3. RESULTS(1) Data analysis results
Table 2 shows the ratio of samples of export fromeach region to Thailand by different transports in
December, 2011.
Table 2 Export from Vietnam to Thailand
North Middle
South
Total
Road 26.8% 10.4%
0.2%
8.8%
Sea
64.8% 71.9%
87.6%
80.2%
Air
8.2%
13.1%
10.9%
10.0%

2. DATA COLLECTION(1) Data analysis
The export statistics from Vietnam to Thailand in5

lead time, high frequency, flexibility the time
shipment. Besides, they can delivery in many places.
(c) However they have also many reasons to notchoose road rather than sea. For example: highexpense, complex of customs procedure, lack ofsafety, one-way cargo, not good quality road(especially Laos), custom (accustomed to sea
transport), unstable political situation in Thailand.
Shipper and logistics companies regard roadtransport as an alternative to other transports. The
reasons of using one route between road and sea is: atthat moment, shippers focus on what kind of factor.
For example, in emergency case, to keep up progress,
they choose road transport or even air transport
regardless of very high cost. Other reasons thatmentioned such as: type of goods, spot of starting
point and destination, amount of cargo.
(d) To make shippers change their choice from sea
transport to road transport, many companies answer
that the most important factor is reduce the cost ofroad transport. It is also necessary to simplify customs
procedure, improve infrastructure, persuade shipper(through advertisement, increase performance).
(e) Types of goods that they transport by road arefarm product (need fast transport), car, motorbike(high price product). Almost goods can be transported
by sea.
(f) Relate to borders, according to companies,
currently, mainly used borders are Cau Treo borderand Lao Bao border. Logistics companies are familiar
with which customs and border police, they will
choose that border. Beyond that, they also consider
based on spot of starting point and destination.
Because newly completed, Cha Lo border is not yet
known, but it is expected to this border will develop
from now on.
(g) Border of Lang Son: Many cargos from China aretransported to Thailand via by road transport.

The ratio of export by sea transport is biggest(about 80% and 65% for northern) and the ratio ofexport by road transport is small (about 9% and 27%
for northern).
When companies use road transport, the mainborders of middle region are: Cau Treo (HaTinh), Cha
Lo (QuangBinh) and Lao Bao (Quang Tri). Notice tocargo of northern and middle region, from data table,
we can see that the ratio of transport that go throughCha Lo border is 88%. The ratio that goes throughLao Bao and Cau Treo is very few (only 6% for eachborder).
Data table also shows types of goods thattransported by road. They’re farm product, clothes,
paper product, plastic product (e.g. shoe, bag), metal
product (e.g. accessory), watch, motorbike.
(2) Estimation of mode choice model
Make utility functions as below:
V_R = α +βC * R_Cost + βT * R_Time (1)V_S = βC * S_Cost + βT * S_Time(α, βC, βTare are coefficients to be estimated.

(2)

R_Cost, R_Time, S_Cost, S_Time which are the twovariables of the model associated with each alternative.)

Model was executed by open source softwarebiogeme. We get the following output:
Table 3 Results of mode choice model
Name
Value
t value Standard error-2.2700
-0.59
3.820
α
-0.0053
-5.58
0.001
βC
-0.6800
-1.15
0.593
βT
Null log-likelihood-88.728
Final log-likelihood-37.709
Rho-square
0.575
Adjusted rho-square
0.541
Number of samples
128
From the value of Rho-square, we can judge thatmodel can explain the real mode choice well.
The sign of value of βC and βT is minus. It means
there will be a negative impact to the selection, whenvalue of cost or time is increased.
Absolute value of t of βC is greater than 1.96. Itmeans that cost is a significant variable. Absolutevalue of t of βC is smaller than 1.96, so it means thattime is not significant.

4. CONCLUSION
The most important factor of choosing route
between road and sea is cost. In current situation,
shipper and logistics companies choose road transport
only in emergency case and regard road transport asan alternative to other transports. Improvement of costand procedure will bring high expectation for road.

(3) Results of qualitative survey
REFERENCES(a) In general, cost is the most important factor andfollowed by customs procedure, safety andinfrastructure.
(b) Reasons of choosing road rather than sea are quick

1) Import – export department [Accessed 1. 2012].
2) General administration of customs [Accessed 2. 2012].
3) General statistics office [Accessed 6. 2012].
4) Ministry of the industrial and commercial [Accessed 6].

6

Fundamental Study on Proposed Efficient Repair
with Combination of Electro-chemical Methods
ID: 04-09083

Name: Daisuke KUNIOKA Supervisor: Prof. Nobuaki OTSUKI

1. Introduction

Table.1: Types of specimen

Steel bars in reinforced concrete (RC) areeasily corroded by the influence of carbonationand chloride ion. It can cause cracks, and if cracks
are not repaired, it leads to further structural
deterioration. In early stage of deterioration,
there are no cracks on the surface of RC, butcauses of corrosion are contained in it. They haveto be removed from RC. Therefore, different
repairing method should be applied according to
the presence of cracks.
Realkalization,
Desalinatin
andElectrodeposition are the kind of electro-chemical
method. Realkalization is the method to recover
pH of RC. Desalination is the method to remove
inner chloride ion. Electrodeposition is themethod to close cracks and coat surface withelectrodeposits. Those three methods are based on
the same electro-chemical system, which is theelectrification between inner steel bars andexternal anode with solutions. Hence, if they are
combined, it is expected that repair of concrete isachieved efficiently by using the common
equipments. Also, by using the combination ofelectro-chemical methods, it is expected that theeffect of each method will be highly enhanced.
However, the removing effect of corrosion factors
is not clarified enough. Hence, the objectives ofthis study are as follows.
1) To clarify influence of solution type on theremoval effect of corrosion factors such as
chloride ion or pH degradation.
2) To propose the most efficient electro-chemical
combination method against chloride attackand carbonation.

Causes ofcorrosion
Carbonation
Chloride attack

Chloride(4kg/m3)

Chloride(8kg/m3)

2.2 Evaluation
Effectiveness of repair was evaluated fromfour points of view.
① Characteristics of electrodeposits② Removing effect of causes of corrosion③ Improvement of resistance to permeation④ Protecting effect against corrosion of steel bars
2.3 Pattern of combination
Table.2 shows the solution used for eachmethod; and Table.3 shows the pattern of
combination. Pattern 1, 2 and 3 were applied forboth specimens with crack and specimens without
crack. Pattern 4 and 5 were applied only forspecimen with crack. Electrification was
conducted for one week in each step. After the
repair, carbonated specimens were exposed tocarbonation, and chloride attacked specimens
were submerged into 3.5% NaCl aq, so that toevaluate the resistance to permeation.
Table.2: Solutions used for each method
MethodRealkalization
Desalination
Electrodeposition

2. Experimental procedure
2.1 Specimen
Mortar specimens in prism shape(40mm×40mm×160mm) were used in thisresearch. Steel bar embedded in mortar was
located at 15mm from mortar surface. Table
1 shows types of specimen.

Presence of cracks
withoutwithcrackscracksCarbonated
CarbonatedSolutions
K2CO3 aq
Li3BO3 aq
Mg(CH3COO)2
aq

Concentration
0.5(mol/L)
2N
0.1(mol/L)

Table.3: Patterns of combination
Step 1
Step 2
Step 3(for a week)(for a week)(for a week)
Pattern.1
ED
Pattern.2
Re or De
Pattern.3
Re or DeED
Pattern.4
ED
Re or De
Pattern.5
ED
Re or DeED
ED: Electrodeposition
Re: Realkalization
De: Desalination

7

15

Re
Re→ED
ED→Re
ED→Re
→ED

5

20

Position of
Steel bar in mortar
15

ED

10

20

0

Carbonation depth（mm）

Carbonation depth（mm）

Position of
Steel bar in mortar

Carbonation depth（mm）

20

ED
Re

10

Re→ED
5

0
0

7
14
21
Charging period （day）

15

After completion of repair

5

0

0

7
14
21
Charging period（day）

Fig.2 Increase of carbonation depth
by carbonation after repairing(W/C=0.55, with crack, carbonation)
12

20

After submersion

Concentration of chloride ion（kg/m3）

Carbonation depth（mm）

After carbonation for one week

10(a) with crack(b) without crack
Fig.1 Changes of carbonation depth during repair(W/C=0.55, carbonation)

15

Carbonation depth
before repairing

After carbonation for one week
After the completion of repair

Carbonation depth
before repairing

10

5

10

Initial Chloride

Before submersion

8
6
4
2
0

0

ED

Re

Re→ED

Untreated

Fig.4 Increase of chloride ion
by submersion test after repairing(W/C=0.55, with crack, chloride attack)Fig.3 Increase of carbonation depth
by carbonation after repairing(W/C=0.55, without crack, carbonation)

against carbonation; that means solution type andpattern of combination has very little influence on
the removing effect.
The resistance to permeation of chloride ion
was improved by repair compared with untreatedspecimen as shown in Fig 4. Especially, by applying
combination pattern, the resistance became higher.
These tendencies were also observed in thecase of specimen without crack.

3. Experimental result
3.1 Carbonation
Carbonation depth decreased regardless thepattern of combination or solution type as shown in
Fig 1. This tendency was also observed in the caseof specimen without crack. It means that the
realkalizing effect is almost the same in eachpattern regardless the presence of crack.
Fig 2 shows the increase of carbonation
depth by carbonation after repairing of specimens
with crack, and Fig 3 shows that of specimens
without crack. The resistance to permeation of CO2
is highly improved by repair. Especially, in the caseof specimens with crack, when the pattern “ED→Re→ED” is applied, the resistance became highest.
In the case of specimens without crack, pattern “Re→ED” obtained highest resistance. These results
are shown in Table 4.
3.2 Chloride attack
The removing effect of chloride ion had
almost same tendencies as realkalizing effect
Table.4: Most effective combination for each corrosion type
Presence of cracks
Causes ofcorrosion
without crackwith crack
Carbonation
Re → ED
ED → Re → ED
ED → De
Chloride attack
De → EDor
ED → De → ED

4. Conclusions
1)

2)

The influence of solution type on the removal
effect of corrosion factors such as pH
degradation or chloride ion is small and thiseffect is followed by charging periods.
The most effective combinations of electrochemical repair method for carbonation canbe proposed as listed in Table 4. For chloride
attack, combination pattern is more effective
than each single electro-chemical method.

5. Reference1)

8

Hitoshi AKIYAMA: An Experimental Studyon Closure of Cracks and Removal of
Chloride Ion by Combining Electrodeposition
and Desalination; bachelor thesis, Tokyo
institute of technology, 2007.2

Adhesion Phenomemon between Rigid Body and Tapering Elastic Beam
Student number: 04B04387 Name: Boqing Wang Supervisor: Kunio TAKAHASHI
Neutral axis of beam is considered as straight line,
thus moment of intertia of area can be defined byLinear beam theory[3] as

1 Introduction
Geckos have micro,nano-hierarchical hairstructure, i.e. setas and spatulas, on their feet.
They use them to climb wall and move on
ceiling[1]. Therefore the adhesion mechanism of
gecko’s foot hair is important for the design of
grip-and- release applications.
In order to analyse the adhesion mechanism of
the hair structure, Sekiguchi assumed the hair tobe an elastic beam of cuboid shape [2], andconsidered elastic energy, surface energy to
discuss the adhesion force between the elastic
beam and rigid body. But, in actuality, gecko’s
foot hairs are not cuboid, therefore, analysis of
the effects of the beam shape on adhesion issignificant forunderstanding adhesion
mechanism of the hair structure.
In this study, a model of tapering elastic beam
was constructed. Elastic energy, surface energy
and shape of the beam are considered to analyse
the phenomemon of adhesion between rigid bodyand the elastic beam.

(2)
By using eq. (2), the relation between length andmoment of intertia of area can be drawn as Fig2.
(b). By using eq. (1) we can know the shape of thebeam as Fig2. (c) when n=0,n=1,n=3.

Fig.1 :Model of the contact between anelastic beam and a rigid body

2 Model of tapering elastic beam
contact with rigid body
Figure 1 shows that an elastic beam contacts
with rigid body at static condition. Slope at root of
the beam is defined as
andis infinitesimal.
There are two kinds of contact at this model: line
contact and area contact Fig.1 (a), (b). Force acted
on the beam is defined as . Displacement of thebeam is defined as
and the origin ofisdefined at the surface of the rigid body. Length of
non-adhesion area of the beam is defined as incase of area contact.
Length, width and height of elastic beam are L,
W and H, the moment of intertia of area is I Fig2.
(a). Height and moment of intertia of area at theroot of beam isand . The surface of beam
which contacts with rigid body and the surface ofrigid body are flat. In order to get a general result,
H is defined as(1)

Fig.2: Image of model of tapering elastic beam

3 Relation between adhesion force
and displacement of beam
At line contact Fig2. (a), by using flexure of beam
by Linear beam theory[3] and geometric relation
between beam and rigid body, relation between
displacement and force can be calculated as

9

(3)

(a)

Normalized force(4)
At area contact Fig2. (b), elastic energy of beamand surface energy can be defined as
andas,:Young’s modulus,:bending moment,:
work of adhesion per unit area.Total energy ofthis model can be defined as
settled at(5)Thus, by numerical analyzing eq. (5), relation
between Displacementand Length of nonadhesion area can be calculated. Whenthe area contact will become line contact. Then,by using eq. (3), (4) and flexure of non-adhesionarea and geometric relation, relation between
displacement and force can be calculated and therelation curve can be drawn as Fig3.
From Fig3, we can know if
becomes greater,
the adhesion force becomes smaller, adhesion
happens with smaller displacement, and when
, line contact does not happen.

n)/1

0

2

(d
/

3
pe:(
Slo

0.1

(L
sin








0.2

-0.1
-0.2

n=1.5
n=0.5

-0.3

n=0.99

n=0

-0.4
-1.2

Normalized force fz/12EI0sin/(Lcos

Length of non-adhesion area
minimum of total energy0(WL2/(2EI))1/2/tan=1.5

(1/Where,,.
When the slope at the tip of the beam equals ,the line contact will become area contact. At thismoment displacement and force isf z/12EI0sin/(Lcos

0.3

-1
-0.8
-0.6
-0.4
-0.2
Normalized displacement d/(Lsin00.05(WL2/(2EI))1/2/tan=1.5
0

-0.05

(b)

-0.1

-0.15

3.9
-0.2

3.01
n=2.01-0.25-2.4-2.2 -2 -1.8-1.6-1.4-1.2 -1 -0.8-0.6-0.4-0.2 0

Normalized displacement d/(Lsin
Fig.3.1:The relation between normalized force and normalized
displacement when (WL2/(2EI))1/2/tan=1.5

4 Conclusion
Fig.4 Image of beam with s-shaped bending

In this study, we can know if
becomes greater,
the adhesion between rigid body and elastic beam
is likely to occur, but the adhesion force becomes
smaller. The condition for which the linear beam
theory can be applied in adhesion between elastic
beam and rigid body theory is when the relation
between displacement and Length of nonadhesion area of beam satisfies.

References[1]Adhesive force of a single gecko foot-hair, K.
Autumn, et al., Nature Vol. 405, 681-685, 2000
[2]Journal of Adhesion Science and Technology,
Y.sekiguchi,vol.26,No.23,December.2012,2615-26
26[3] 材料力学最新機械工学シリーズ１４ 渥美光
鈴木幸三 三ヶ田賢次 森北出版

10

The Research of the Reinforcement Learning for Humanoid
robot’s walk control
Student Number：07-04050 Name：Soichi Uchiyama

1

Introduction

Humanoid robots can easily fit into people’s daily
lives comparing to other types of robots. In recent years,
the robots are used in many fields. To adapt to complex and diverse environments, a number of researchers
study how to get a decision policy for the robot by its
self-learning. The reinforcement learning is one of approaches, and many studies have been done. In thelearning stage, the reinforcement learning is not necessary to teach the action in advance and the environment. Only the compensation for its action is necessary.
Therefor learning can be done autonomously so that it
is very effective.
In this paper, in order to make a humanoid robot adapt
to its many degrees of freedom, I develop a simulator
for a humanoid robot to be able to acquire the walking motion by using Q-learning with the dynamic structuring of exploration space based on genetic algorithm(QDSEGA). QDSEGA is an approach which combines
evolutionary computing and reinforcement learning.
Fig.1 shows the outline of QDSEGA. First, the features of the robot are incorporated into the learning device so that the learning device can limit the search spaceand extract a subset of the search space. Next, the extracted search space is adapted to the reinforcement learning. Then, the subset is rebuilt with the knowledge
gained from the reinforcement learning. By repeating
this process, an effective subset and a policy for the reinforcement learning are obtained. Usefulness of simulation using QDSEGA is indicated by Ito [4].
I apply this algorithm to the walking behavior acquisition problem for a humanoid robot, and verify the usefulness of results of the simulation.

2

Super visor：Yukihiko YAMASHITA

humanoid robot, Q-learning is difficult to adapt because
the state space and freedom of action are too large.
The classifying system is trained by genetic algorithm(GA). If we try to learn a policy in a large space, suchas that of humanoid robots, it is difficult to acquire of aneffective policy becase of the difficulties of maintaining
consistency between the set of states indicated by thecondition and the set of possible state transitions by theaction.
QDSEGA approach can solve the problems. QDSEGA is designed for complicated systems with large
action-state space like a robot with the redundant degrees of freedom.
.
.
.
.
.

Fig. 1: Outline of QDSEGA
.
.
.
.
.
.
.
.

Learning Algorithm

Q-learning is a typical approach in the reinforcement
learning. Q-learning has a theoretical framework regarding to optimazation. If the degree of freedom ofaction is small, Q-learning is very effective. However,for the problem of multi-degree-of-freedom, such as a

Fig. 2: Learning algorithm of QDSEGA

11

Fig.2 shows the flowchart of learning algorithm of
Action acquisition process using reinforcement learnQDSEGA. The learning algorithm is composed of a com- ing algorithm usually needs very long running time. In
bination of action generating dynamics with GA and be- addition, as this paper, for learning of waking behavhavior acquiring dynamics with Q-learning. Actions are ior acquisition, wide range position measurement is rerecorded as a gene and generated by GA. First, the ini- quired. Therefore, even if we produce a real robot, wetial population is generated randomly. Q-table is gen- need to use its kinetic simulator. In this paper, I built
erated from the elements of the actions and the possi- a simulator for a humanoid robot by using open-source
ble states. Q-learning is performed using this Q-table 3D physics engine called Open Dynamics Engine(ODE).
(Fig.3). If Q-learning is done sufficiently, Q-learning Usefulness of simulation using the ODE is indicated byexits the learning stage and the fitness of each action Endo [5]. In this paper, I don’t use an actual robot, andevaluates based on the Q-value. According to the eval- verify by using the simulator.
uation, parents are selected. Then, a new set of individFig.4 shows the screen of the simulator.
uals is created from the parents by crossover and mutation. By repeating this process, a set of effective actions
4 Conclusion
and a policy are obtained.
.
I made a walking simulator, and its possibility. For.
the future works it is necessary to set an appropriate ac.
tion space against the explosion of the number of di.
mension, generate a Generate a value function with amulti-modal, and apply the results to real robots.

References[1] Kazuyuki Ito, "The Robot Intelligence -the evolutionary computing and the reinforcement learning," Ohmsha, Tokyo, 2007 (in Japanese)

Fig. 3: Example of Q table

3[2] Kosei Demura, "Simple! Practice! The robotic simulation," Morikita publishing firm, Tokyo, 2007 (inJapanese)

Simulator
.
.
.
.
.

[3] "Open Dynamics Engine"http://www.ode.org/[4] Kazuyuki Ito, Fumitoshi Matsuno, "Control ofReal Robot Using QDSEGA -Extension of Layered
Structure and Application to Snake-like Robot"
JRSJ Vol.21 No.5, 2003 (in Japanese)[5] Gen Endo, Keisuke Arikawa, Shigeo Hirose,"Quantitative Evaluations of Open Dynamics Engine as a Research Tool -Comparison of JointTorque and Power Using Quadruped Robot TITANVIII-," 28th RSJ’s academic meeting, 2010 (inJapanese)

Fig. 4: Screen of Simulator[6] Masatoshi Yamaguchi, Tadahiko Murata, "The Effectiveness of Neighboring Crossover for Controlling Multi-Legged Robot Using QDSEGA" 16th
RSJ’s intelligent system symposium, 2006 (inJapanese)

I apply this algorithm to the walking behavior acquisition problem for a humanoid robot, and verify the usefulness of results of the simulation. In addition, by introducing a multi-modal reward, the acquired behaver is
optimized.

12

Application of TIV for multi-scale atmospheric flows
学籍番号：08B00019
1

氏名：アーリフ ムタズ

はじめに慣性を持ち，大きな表面温度変動が期待できるためである．表 1 に各種観測概要について記す．

TIV (Thermal image velocimetry) はサーモカメラ
で撮影された熱画像の時間的な変位を追跡すること
で，物体表面近傍の速度場を同定する手法である．
本手法の適用例はこれまでに屋外大気環境下における建物鉛直壁面において精度検証が行われているのみであり，今後様々な場への応用が期待される[1]．
本研究ではその手始めとして，（１）様々な空間ス
ケールの流れの観測と，（２）複雑形状に沿った流れの計測についての応用を検討した．（１）に関して，本手法の測定範囲はカメラの撮影範囲とほぼ等価であるため，数 km に渡る広範な領域から，数 m 以下の領域への適用も原理的には可能であるが，その技術的側面について実観測に基づいた検討を行った。
（２）の複雑形状への応用については主にその方法論について記述する．

2

表 1：サーモカメラの撮影領域と時空間解像度
L [m]
ds [mm/pixel]
F [Hz]
ケース
50
山
～1000
～1000
100
人工芝
～100
～100
30
30
壁
20ｘ3
1.8
100
円柱
0.3ｘ0.5
1.8
166
平板
0.58ｘ0.46
L: domain size, ds: resolution, f: sampling frequency
(a)

測定手法

TIV はサーモカメラで撮影した一連の連続熱画像
データから，表面近傍の風速場を定量的に推定する。
測定原理は PIV(Particle image velocimetry)[2]とほぼ同じであるが，粒子画像の代わりに熱画像を利用する
ため，場の不均一性に起因した温度むらを，ハイパ
スフィルターを用いて除去する操作が必要である[3]．
上記の解析を行う上で，3 つのパラメータを設定する必要がある．それは，画像相関を行う検査領域の大きさ S，ハイパスフィルターの平均化時間 Tave，
そして比較する 2 枚の画像の時間間隔 dt である．本
研究では，これらのパラメータの撮影領域依存性（画像解像度依存性）について検討し，パラメータを決定する指標を作ることを試みる．

3

指導教官：神田 学

(b)

観測概要
本観測では，空間スケールの大きく違う領域を選んだ。その中から紹介するのは（１）人工芝の運動場（約 100 m 四方），（２）屋外に曝したスタイロフ
ォーム平板（約 50 ㎝四方），（３）鋸山山麓（樹冠部，

図 1： 移流速度ベクトル分布

約 1 km 四方）、
（４）スタイロフォーム円柱（約 30(a) 人工芝グラウンド, (b)スタイロフォーム平板㎝四方），の 4 つであり，これに既往の研究で行われ4

た（５）建物壁面観測結果を加えて比較を行う．
（１），

スケール依存性（２），（５）は平坦な場であるが，その空間スケー
ルが大きく異なる．
（３），（４）は平面ではなく立体

平面の 2 つのケースと壁面のケースについて，各面の流れの計測を実施する．

種解析パラメータを比較する。図 1 はグラウンドとこれらの領域の共通項として，撮影面が小さな熱

スタイロフォーム平板の観測結果であり、流れの様

13

6

子が風速ベクトルにより表現されている．適切なパ結論

ラメータの選択は目視に基づいており，誤ベクトル
スケールが違う領域や、複雑な幾何形状を有する面に対して TIV を適用し、その適用可能性について検討した。今回テストした 1m 弱の平板及び、100m程度のグラウンドでは本手法を適用することが可能であることを定性的に確認した。また空間スケールの異なる 3 つの観測の比較から、時間フィルターの大きさや、解析画像間隔が撮影画像の空間解像度に大きく依存することを示した。その一方で、画像相関を行う検査領域の大きさについてはほぼ一定値であることが分かった。この結果は他の実験を行う際の、パラメータの設定基準に応用することができると考える。複雑形状面への応用については、流れの可視化のみ示した。今後、テクスチャーマッピング
による物体表面への投影により、物体面に沿ったベクトルの計算などができるものと考えられる。

の尐なさや実際の風向との一致などを判断基準に行った．
次に各観測例における最適な解析パラメータについて見てみる（図 2）．主なパラメータとしては，画像相関を検定する領域の大きさ𝑆，時間フィルターの大きさ𝑇𝑎𝑣𝑒 ，画像間の時間幅𝑑𝑡などがある．まず𝑆に関してはスケールに依らずピクセル数で決まって
おり，スケールにはあまり依存しない傾向が見られた．つまり画像相関という統計手法の制約で決まる．
𝑇𝑎𝑣𝑒 ，𝑑𝑡については明確なスケール依存性があり，それは流れ場の速度と画像の空間解像度によって主に決められるためである．

(a)

図 2: 画像解像度に対する(a)平均化時間 Tave，(b)画像相関の時間間隔 dt，(c)検査領域サイズ S 依存性(b)5

複雑な地形への応用

次に複雑な物体表面周りの流れの計測について，適用可能性を検討した．図 3 は鋸山とスタイロフォー
ム円柱の観測結果である．解析パラメータは図 2 の結果を参考に，表 2 のように決定した．図 3 は時間フィルターを施した温度変動の画像であり、それぞれ鋸山（図 3a）とスタイロフォーム円柱（図 3b）の観測結果である。鋸山の場合は画像に移る前後の山の空間的な連続性が無いため、温度画像をパッチ状に切り分けて解析を行った。その結果流れに沿って温度分布が引き伸ばされている様子が見て取れる。
速度ベクトルを算出する際、平坦と異なり 2 次元で表示される速度ベクトルを 3 次元形状の面にそれぞれ投影することで、実際のベクトルが算出可能であると考えられるが、これについては今後の課題で
ある。

図 3

変動温度

(a)鋸山, (b)スタイロフォーム円柱
コンターは表面温度分布を表し、赤い方が高温である

参考文献[1] 神田学,稲垣厚史,久米村秀明:“サーモカメラを用
い た 壁 面 近 傍 風 速 の 計 測 手 法 開 発 ”, な が れ ,Vol.31, pp.491-493, 2012
[2] Adrian RJ (1991) Particle-imaging techniques forexperimental fluid mechanics. Annu Rev Fluid Mech
23: 261-304[3] Christen A, Meier F, Scherer D (2012)
High-frequency fluctuations of surface temperaturesin an urban environment. Theor Appl Climatol 108:301-324

表 2 ：TIV の各種解析パラメータS [pixel]
Tave [sec]dt [sec]
ケース
21
6
0.8
山
17
0.5
0.01
円柱
S: interrogation window size, Tave: averaging time,
dt: timespan between two images

14

Prototype of boost converter for energy harvestingwith automatic Maximum Power Point searchingStudent number: 08B00597 Name: Ryushi ASAMI Supervisor : Kunio TAKAHASHI
defined as follows[2].

1. Introduction
Energy harvesting, which means gathering unused tinyenergy from environments such as lights, vibration, or
heat, is mainly used as power source of low powerdevices. For example, it is used for driving wirelessswitch [1], and expected for further application such as
driving wireless sensor networks.
An optimized storage method of storing tiny electricenergy using boost converter, was proposed and studied
by Tanaka [2], from which he could achieved storing
efficiency of 90%. However, tiny electric energy issensitive to fluctuation of environmental factors.
Therefore, searching Maximum Power Point (MPP) ofpower source is an essential function needed for theforementioned boost converter.
The purpose of this study is to make a prototype of theboost converter to storage electric energy from solar cellto capacitor, taking into account the searching of MPP.

1𝑑gt = 1 + log
𝑇̃𝑔̃CD =𝑉C + 𝑉D𝑉MP

𝑅
𝑇𝐿
𝑅 = 𝑅MP +𝑅L +𝑅C ≈ 𝑅MP +𝑅L +𝑅DS
We have to know 𝑉MP , 𝑅MP , and capacitor voltage 𝑉C inreal time to apply optimized duty ratio.
3. Prototype and experiment
3.1 MPP searching method
We introduce the MPP searching circuit including
power source as shown in Fig. 3.

Fig.3: MPP searcher
At first, the capacitor has low voltage, the solar cellstarts charging it from low voltage. The Voltage ofcapacitor increases, and the voltage of the solar cellincreases consequently. We can observe a wide range of
I-V characteristic of solar cell and decide where MPP is.
3.2 Circuit
We setup the circuit for experiment as shown in Fig.4.

𝐼MP

Maximum Power Point

Power[W]

Solar cell 's current I [A]

2
1
1 2
+√( −
) (1 − exp(−𝑇̃)) + exp(−𝑇̃)
2 𝑔̃CD()

𝑇̃ =2. Theory
2.1 Power source’s model
In this study, we use a solar cell as a power sourcewith the I-V characteristic and P-V characteristic
shown as Fig.1 [3]. The largest output point is called
Maximum Power Point, its voltage and current isdefined as 𝑉MP and 𝐼MP respectively. 𝑅MP is defined𝑉
as MP .

IMP[A]11( −) (1 − exp(−𝑇̃))
2 𝑔̃CD:I[A]:Power[W]
VMP[V]
Solar cell's voltage V[V]Fig.1: Theoretical characteristics of solar cell
2.2 Circuit model
Figure2 shows the boost converter circuit and its
equivalent model. It is assumed that power source
works nearby MPP, and 𝑉C is constant during one
switching period.

Fig.4: Prototype circuit
The specification of the circuit components are as
follows;
Boost converter
L: 0.47[mH]𝑅L : 0.21[Ω]Inductor
𝑅DS : 78[mΩ]MOS-FET
Diode
𝑉D : 0.32[V] 𝑅D :0.20[Ω]Capacitor
C : 0.22[F] 𝑅C : 6.0[mΩ]

Fig.2: Boost converter circuitand equivalent circuit model
Boost converter works by switching FET on and off.
The ratio of FET on-time in one switching period isdefined as duty ratio.
2.3 Formula
Optimized duty ratio 𝑑gt , maximizing efficiency, is

MPP searcher
MOS-FET’s D-S on resistance(x 3)
Shunt resistance
Capacitance15

0.1[Ω]
1.00[Ω]1000[μF]A computer is used to search P-V characteristics
shown as Fig.5, measure 𝑉C , calculate 𝑑gt , and output
voltage to PWM circuit to control FET Switching.

3.5 Discussion
From Fig6, charging efficiency of proposed method,
searching MPP and in real time is almost equal to thatof MPP parameters fixed. It can be said MPP searchingis almost exactly.
From Fig7, efficiency of proposed method is 135% ofthat of MPP fixed. It can be said proposed method has
advantage when MPP changed.
On the other hand, efficiency of proposed method is
80~90% of theoretical value. Figure.8 shows behavior ofsolar cell when charging. It shows that output voltageof solar cell doesn’t follow the theory after FET’sswitching. It is not clear that switching noises effectefficiency, but if solar cell didn’t work enough whenswitching noises occur, differences between theory andexperiment are reasonable. To discuss this defect, wehave to make the circuit not to be influenced from any
external factor like electromagnetic waves.

0.14
0.12
MPP

Power[W]

0.1
0.08
0.06
0.04
0.02
0
0

0.5
1
1.5
Solar cell's voltage V[V]2

Fig.5: Measurement of P-V characteristics
Measuring 𝑉C is done every 100ms. Measuring 𝑉MP
and 𝑅MP is done every 5000ms. One searching MPP
needs to separate the solar cell from boost converter
less than 30ms. Measuring 𝑉C , calculating 𝑑gt , andoutputting don’t need to separate
3.3 Experiment
We adjusted the solar cell as 𝑉MP ＝1.9[V], 𝑅MP = 24[Ω]and charged capacitor from 2V to 6V by this circuit.
After that, we gave to the computer fixed parameter
𝑉MP ＝1.9[V], 𝑅MP = 24[Ω], and charged capacitor again.
Next, we adjusted the solar cell as 𝑉MP ＝1.7[V], 𝑅MP =32[Ω] and experiment for both parameters settings
again.
3.4 Experimental results
Figure.6 and 7 show the experimental results

Solar cell's output voltage
Vout[V]FET off

Capacitor's voltage VC[V]4
C = 0.22[F]Switching Freq = 100kHz
VMP=1.9[V],RMP=24[]:realtime MPP searching:fixed VMP=1.9[V],RMP=24[]:theoretical optimized value
0

10

20
Time[s]Capacitor's voltage VC[V]654
C = 0.22[F]Switching Freq = 100kHz
VMP=1.75[V],RMP=32[]:realtime MPP searching:fixed VMP=1.9[V],RMP=24[]:theoretical optimized value

2
0

20

40

2
1.8
Solar cell's
output voltage
Vout[V]

1.6
1.4
0

2

4
6
Time[s]8

10References[1] EnOcean. http://www.enocean.com/home/[2] Masanori Tanaka ”Optimization of control
parameters of a boost converter for energy harvesting”
Journal of Physics: Conference Series 379 (2012)[3] Kiyoshi Takahashi “太陽光発電”(1980)
ISBN： 9784627740303(in Japanese)

30

Fig.6: Experimental result

3

2.2

4. Conclusion
We proposed method of the boost converter to search
MPP and measure 𝑉C in real time. The proposed method
can increase harvest when MPP varied. This MPP
searching method is applicable other power sources likea thermoelectric element and also have room for
improvement in deciding process of MPP.
To make more harvest, we need to identify what factor
make this method less efficient. The switching noises areworth to investigate.

5

2

FET driving
signalFig.8: Switching noises

6

3

FET on

FET on60

Time[s]Fig.7: Experimental result
16Financial Analysis of Vertical Integration and Privatization for Airports
in Hokkaido
Students Number ：08-08759
1

Name：Ryuuhei Kudou

Introduction

Supervisor：Shinya Hanaoka

figure 1, Okushiri, Rishiri and Rebun airport are

Air transport is very important for Hokkaido

located in isolated island. Chitose airport has no

because of distance between cities, and distance toprivate schedule air flight service, hence thisurban area outside Hokkaido. But the estimation byresearch excludes this airport from analysis.

Hokkaido prefecture government says many of

Methodology and Results

airports in Hokkaido are operated at a loss and3

deficit in these airports are paid out of government's

3.1 Calculate items

general account. To reduce the deficit covering by

This analysis is based on form of operating profit

government's general account and to operate more

and loss; P/L, and deal with subsidies to containing

flexibly, it is needs to introduce private capital andnothing, only fuel tax of plane, or containing all

operating know-how by privatization.

subsidies.

I

calculated

about

after

verticalOn the other hands, there is no estimation of

integration (VI) and 2 types of privatization, Saleeach airport’s balance until recently, and Airside;

type, and Concession. Table 1 is operating P/L of VI,

facilities of airports is operated by public like The

in Sale type, and in Concession.

Government, prefectures, and cities, in contrast,

Table 1, P/L in BI, Sale type, Concession

Landside like airport buildings is operated byVI

private capital with public, or only private. In

Sale type

Concession

Airside income
Landside income
Labor cost(A)Airside income
Landside income
Labor cost(A)are made much account for airport operation. So it isAirside income
Landside income
Labor cost(A)need for vertical integration of Airside and Landside

Rent for lot and building（A） Rent for lot and building（A）

abroad cases, the commercial revenues by Landside

before privatization in Japanese airports.
This

researchaims

tocalculatefinancialbalances of airports in Hokkaido after vertical
integration, and to estimate profit or loss of airports
after privatization.
2

Subsidy for location
Depreciation cost（A）

Fixed assets taxes
Depreciation cost（A）

Other costs(A)Labor cost（L）
Depreciation cost（L）
Other costs(L)Other costs(A)Labor cost（L）
Depreciation cost（L）
Other costs(L)Other costs(A)Labor cost（L）
Other costs(L)Operating profit and loss Operating profit and loss Operating profit and loss

Airports in Hokkaido
There are 14 airports in Hokkaido as shown in

P/L with fuel tax

P/L with fuel tax

P/L with fuel tax

P/L with all subsidies

P/L with all subsidies

P/L with all subsidies

Figure 1. 4 airports are operated by the central

There is no difference in P/L about airport

government (●), 2 airports by SDF (Self Defense

operation of after VI and Sale type because theForce) (☆), 6 airports by Hokkaido prefecture (▲),

operator and owner of airport translate private toand 2 airports by their cities government (○). In

public at same time. But profit on sales will be paidby new owner to current owner, public. ‘Subsidy forlocation’ is like fixed assets tax on government
properties. Subsidy for location will change to fixed
assets tax, because owner changes public to private,
but sum won’t change.
In case of Concession, operator is private, butowner remain public as VI. Therefore, the expenses
that the owner must pay will be paid by owner of
airport but not operator. The expenses are ‘rent for
lot

and

building’,‘subsidy

forlocation’,‘depreciation cost’ of Airside and Land.

Figure 1, Airports in Hokkaido

17

and

3.2 Items of P/LLabor cost (L) ＝((number of employee * employee’s

Airside income is assumption of landing fee andsalary + (number of full-time officer + number ofrent from lot and building. I use ordinary expenses

part-time officer*0.163) * officer’s salary) * 1.139

of airport buildings for ‘Landside income’. Because I

Other cost (L) is cost of goods sold and selling andofadministration costs otherwise labor costs andAsahikawa and Obihiro, I use the data of airports

depreciation costs. But it is difficult to get them, so I

whose histories of building runway are similar touse ordinary cost of buildings minus interest

these airports. Other cost (A) is assumption of

expenses, labor cost, and depreciation cost. The ratecannot

get

data

of

depreciation

cost

(A)

building cost, cost of supplies cost for environmentalof fuel tax is different as the operator, so I revise themeasure, and so on. Because I cannot get labor cost

rate as operator gets all of fuel tax in the airport. I(L), I use number of employee and officer. Thesecan get other costs.

number multiple by their salaries and add legal
3.3 Results

welfare expenses. I deal with part-time officer as
0.163 of full time officer according to the rate ofTable2 shows the results of calculation based onsalary in incorporated administrative agency.

Table1 about P/L in VI, Sale type privatization, andCalculation formula is follow;

Concession with conditions of subsidies about 13
airports in Hokkaido.
Table 2, P/L in BI, SALE type, Concession

OpertorAirportThe government
Hokkaido prefectureCity
SDF
New Wakk Kushi Hako Okad Mema Nakas Mema Okus Rishi Reb Asahi Obihi
Chitose anairodate ama nbets hibets nbetsu hiriri
un kawa ro million

Operating P/L

2,087-897 -1,048

-817

-343

-449

-340

-347

-277

-281

-6

-343

-440

P/L with fuel tax

4,595

-845

-813

-279

-343

-286

-308

-347

-277

-281

-6

584

543

P/L with all subsidies

6,915

-796

-577

233

-341

-268

-299

-346

-267

-281

-6

856

928

Operating P/L

2,087-897 -1,048

-817

-343

-449

-340

-347

-277

-281

-6

-343

-440

P/L with fuel tax

4,595

-845

-813

-279

-343

-286

-308

-347

-277

-281

-6

584

543

P/L with all subsidies

6,915

-796

-577

233

-341

-268

-299

-346

-267

-281

-6

856

928

Operating P/L

4,679

-388

-237

506

-169

-116

-201

-167

-83

-218

-5

289

-22

P/L with fuel tax

7,187

-336

-2 1,043

-169

47

-169

-167

-83

-218

-5 1,215

961

P/L with all subsidies

9,507

-287

234 1,555

-167

65

-160

-166

-73

-218

-5 1,487 1,346

According to this result, all airports except New

4

VI
SALE
type
Concession

Conclusion

Chitose are operated at a loss in case of VI if there isThis research analyzes P/L in after verticalno subsidy. On the other hand, Concession which

integration and privatization with the use of pastneeds the lowest operating cost can reduce thedata of Airside and Landside. There are these future

airports that operating at a loss to 7 airports, and 3

issues; 1) Estimate future items with the use of past

among them is airports in isolated island. New

data and analyze future privatization use thesechitose, Hakodate, Asahikawa, and Obihiro canitems. 2) Analyze about operating some airports as a

privatize conditional on subsidies in SALE type. And

group to privatize more non-profitable airports.

Kushiro and Memanbetsu may be able to privatizein Concession with subsidies. However, it is needs toReferencesnote that this table is only about airport’s economic[1]

Anne Graham: ’’managing airports’’, Butterworth -[2]

Res. Group for Air Transport Policy and Practices inbalances, and the government must pay owner’s cost

Heinemann, 2008.

in case of Concession (Table3).
Table 3, Incomes and Expenses of public

SALE type
Incomes

Profit on sales

Japan:

Concession
[3]

way

of

the

futureairportH. Murakami, K. Kato, N. Takahashi, Y. Sakakibara:‘Economics of airport’, Minerva syobou, 2008.

Subsidy for location(Subsidies)theadministration’’, 2009.

Initial payment
Concession fee
Rent for lot and building

Expenses’’About

Depreciation cost
(Subsidies)18

Geometrical Shape of Arch Formed by Collapse of Undercut Slope
Student Number: 08B20281

Name: Hidehira HIRAISupervisor: Thirapong PIPATPONGSA(1) Axial force exerting along the axis of arch with angle Ψz  x tan 

1. Introduction
A stable arch formed across a pit is beneficial to the design ofan undercut slope in open-pit mining; therefore, prediction of the
maximum stable undercut width under which the slope does notcollapse is needed. The relation between a stable width and an
inclined angle has been obtained experimentally to confirm thedeveloped theoretical relations1). The actual engineering
application under this theory was immediately tested at the Mae
Moh open-pit lignite mine in Thailand in December 2011. Fig.1
shows the stage of excavation as observed in October 2012 on
the right side of the pit. At this targeted area, the stability ofundercut slope is not only depended on arch action, but is also
contributed by some reinforcement to the arch crown using arrayof shear pins. Because, design of shear pin requires reliable
estimation of geometrical arch shape, this study reveals afundamental research in an attempt to understand criteria of archshape through a series of physical models using both 10 yen
coins and moist silica sand no.61-3).

(1)(2) Uniformly distributed load w along the gravity direction
1 w 2
(2)
z
x
2 Fh
(3) Own weight of arch with constant unit weight γ  
F (3)z  h  cosh  x   1 
 Fh  (4) Uniformly distributed load ω acting perpendicular to arch
2

Fh 

 Fh 
2
z    x  





2

(4)

Table 1 Basic geometries of arch under four particularsymmetrical loading conditions
Equilibrium
Equilibrium
GeometricalFormalong x-axis
along z-axis
shape
dFy
0
(1)
Triangledl
dFy
dx
dFx dx  0w(2)Parabola
dl
dl
dFy

(3)
Catenary
dl
dFy
dFx
dy
dx
 

(4)
Segment
dl
dl
dl
dl

Fig. 1 Actual practice of undercut slope at the Mae Moh lignite
mine in Thailand (as of October 2012)
Structural adaptation instigated by relative displacement andprincipal stress reorientation implies that both continuum anddiscrete structures can form stable stacks of arch. Formation ofarch observed from free surface or scarp after removal or
excavation of basal supports provides useful information about
load transmission. Load breaking the arch is its own weight. If
static stress adjacent to the scarp is higher than the unconfinedcompressive strength or buckling resistance of the materials,
arch will eventually collapse.

2. Theoretical background
Granular media have ability to transfer load laterally to thestationary zone by readjusting the axes of principal stresses toform the load-bearing structure. This mechanism is known asarching effects4). As shown in Fig.2, a convex upward scarp
involves passive arch action where the connecting trajectories ofmajor principal stress represent arch shape confined by thedisconnecting paths of minor principal stress. Herein, thestress-free boundary represents zero minor principal stress.
Components of contact force F are Fx in the x-axis and Fz in thez-axis. The vertical and horizontal supporting forces are Fv andFh, respectively. The differential length of curve is δl.

3 Physical model tests
Typical arch shapes based on triangle, parabola, catenary andsegment, as expressed in Eqs.(1)‒(4), respectively, are employed
for back analysis in order to approximate the unknown loading
conditions. Two different laboratory tests of undercut slope
made of closed-packing of 10 Yen coins placing on acrylic plate(Fig.3) and moist silica sand no.6 placing on acrylic sheet(Fig.4) are considered. Tables 2 and 3 reported the testing
criteria of undercut slopes using 10 yen coins and moist silica
sand1,3) (density 1395 kg/m3 and water content 10%),respectively.

Fhgravity direction

x

z
Fh

x

Differential
length x  y  l
2

Fv

z

2

2

Contact
force
F

y
l
x
Fx  Fl

Fy  F

Fh

Fv

Fig. 2 Formation of stable arch where chain of contact force isconnected to transfer load laterally
Contact chains in granular materials are outcome of structural
adaptation against the applied loads. Such force patterns can beillustrated by the arrangement of granules with no shearing force
acts on any of the contact points. According to Fig.2,
equilibrium of forces along the x-axis and the z-axis canformulate differential equations of force along the axis of anarbitrary curvature. Four structures of arch providing the
maximum stability are derived under the following symmetricalloading conditions which are also summarized to Table 1.

Fig. 3 Slope model made of 10 Yen coins showing a stable scarp
after removal of supporting racks

19

mm
300
250
200
150

Unit: mm
100

0.90
0.80
0.70
0.60
0.50
0.40
0.30
0.20
0.10

50
0
-50
-100
0

Fig. 4 Slope model made of moist silica sand no.6 with stable
undercut span after symmetrical removal of the basal support

100

200

300

400

500

600

700

mm

Fig. 5 Displacement vectors plotted from PIV results of Test 5mm

4 Results
Photo images taken in normal direction to the stable scarps of
stable slope made of 10 Yen coins were digitized. Mean squareerrors (MSE) in unit of cm2 obtained from the least square
method for interpolating each of typical four arch shapes to theactual shape of scarp were calculated and compared in Table 2
and 3. Numeric values shown in bracket represent the smallest
MSE among the choices of typical arch shape. The segment
curve obviously exhibits the best fitting geometry with theobserved shapes of scarp in physical models.

400
350
300
250
Shear strain
200

0.002360.00209
0.001820.00155
0.001280.00101
0.000740.000470.00020

150
100
50
0
0

Table 2 Testing results of undercut slope using 10 Yen coins
Test №
Inclined angle
Width
LengthShapes of
(˚)
(cm)
(cm)
scarp
Case 140
Parabola
Case 250
40
Segment
Case 360
Segment
Case 4
40
Catenary
Case 550
Segment
70
Case 660
50
Segment
Case 770
Segment
Case 840
Segment
Case 950
Segment
100
Case 1060
Segment
Case 1170
Segment
Note: A coin has a diameter 2.35 cm with mass 4.5 gram

100

200

300
mm

400

500

600

Fig. 6 Contour of shear strains based on PIV results of Test 5

Fig. 7 Collapse of undercut slope in a buckling mode of Test 5
5. Conclusion
Back analysis of loading conditions of undercut slope models
using 10 Yen coins consistently confirms that the geometricalshape of arch can be estimated to a segment. Moreover, thesegmental arch indicates that the assumption of uniformly
distributed load acting along the arch used in the undercut modelof moist silica sand is acceptable. The result of this study isbenefit to a design of reinforcement to undercut slope by giving
the reliable shape of failure zone using the segmental arch.

Table 3 Testing results of undercut slope using moist silica sand
Test
Inclined Width
Length
Wf
Bf
Mode
№
angle
(cm)
(cm)
(cm) (cm)Test 160˚55.0 50.3
1
Test 265˚47.5 39.2
2
Test 370˚40.0 32.4
2
70
40
Test 475˚35.0 28.0
1
Test 580˚
37.5 24.8
2
Test 690˚
25.0 20.9
2
Note: Thickness for all series is 5 cm
Wf = actual failure width, Bf = predicted failure width
Mode 1 = scarp failure, Mode 2 = buckling

References
1) M.H. Khosravi: Arching effect in geomaterials withapplications to retaining walls and undercut slopes,
Department of International Development Engineering,Tokyo Institute of Technology, PhD Dissertation, 2012.
2) M.H. Khosravi, T. Pipatpongsa, A. Takahashi and J.
Takemura: Arch action over an excavated pit in a stable
scarp investigated by physical model tests, Soils andFoundations, 51(4), pp.723-735, 2011.
3) M.H. Khosravi, L. Tang, T. Pipatpongsa, T. Takemura and P.
Doncommul: Performance of counterweight balance on
stability of undercut slope evaluated by physical modeling,
International Journal of Geotechnical Engineering, 6(2),pp.193-205, 2012.
4) K. Terzaghi: Stress distribution in dry and saturated sandabove a yielding trap-door, The First International
Conference on Soil Mechanics and Foundation Engineering,
Cambridge, pp.307-311, 1936.

The results of image processing by particle image velocimetry(PIV) for the slope models made of moist silica sand no.6 areshown in Fig.5 as illustration. In this figure, the displacement
vectors show the movement of the slope from the initialcondition before excavation until the time just before failure of
the slope and Fig.6 shows the contour of the shear strains which
are corresponding to the displacement vectors. The lower regionin Fig.6 showing the maximum shear strain indicates a possible
location of failure. Fig.7 shows the photo of failure at the criticalundercut width (Wf). Table 3 also summarized the predictedundercut width (Bf) and mode of failure. The value of Wf
acceptably agreed with Bf using the buckling model).

20

Numerical Analysis of Natural Convetion in an Arch Vesselby a Spectral Finite Diﬀerence Scheme
Student Number : 08-23457

Name : Keita Mitani

1. Introduction

Supervisor : Yoshihiro MOCHIMARU

Energy Equation (with neglecting dissipation):

Computional ﬂuid dynamics is adapted in many ﬁelds,
e.g. engeneering, medical engeneering and mechanicaldesign. A spectal ﬁnite diﬀerence scheme is a method
in computional ﬂuid dynamics. It has higher acuracy
and higher eﬃciency than others. In this study, a heatand ﬂuid ﬂow ﬁeld in an arch vessel placed horizontally is analyzed numerically, using a spectral ﬁnite
diﬀerence scheme.

ρcpDT= κ∇2 T,
Dt

(3)

where V : velosity, t : time, ρ : density of ﬂuid, p :pressure, µ : viscosity, T : temperature, κ : thermal
conductivity, cp : speciﬁc heat at constant pressure.
Conformal mapping from physical region (x, y), (z ≡x + i y) to the analytical region (α, β)(−∞ ≤ α ≤0, |β| ≤ π), (w ≡ eα+i β ) is given by:

2. Analysis
Fluid enclosed is assumed to be incompressible andNewtonian. The wall of lower part is heated at a uniform temperature, while the other is kept at a lower
constant temperature. No slip ﬂow on the wall is assumed.

w=

1 + if,
1 − ifwhere f is an analytic function deﬁned by)(z+1
f ≡ exp a ln
+ bi
z−1

a=

π
πθ
,b = −π−θ
π−θ

3. Basic Equations
Equation of motion :DV= −∇p + µ∇2 V + ρg
Dt

(1)

Equation of continuity (under Boussinesq approximation) :∇・V = 0(5)

(6)

Then, the upside discribes an arc, and the downside
a straight line. θ : angle of circumference from thedownside to the upside.

Fig.1 Model in analysis

ρ

(4)

(2)
1
21

Fig.2 Grid points of physical region

Boundary conditions (on r = 1 ):Then Eqs.(1) and (2) lead to the vortisity transport
equation in dimensionless form with a change of variables (r ≡eα , 0 ≤ r ≤ 1):ψJζ∂ζ∂ (ζ, ψ)
1
J
+r
= √
∂t
∂ (r, β)
Gr

(

∂
∂2
r
+
r
+
∂r 2
∂r
∂β 2
2 ∂

2

)
ζ+r

T
T

∂ (T, y)∂ (r, β)
(7)

= 0
= −

(15)
2

2 ψ (rn−1 )(1 − rn−1 )= 0 ; −π ≤ β < 0= 1 ; 0 ≤ β < π,

(16)
(17)
(18)

where rn−1 is one step before the boundary (rn ) inﬁnite diﬀerence approximation.
Relation between vorticity and a stream function:

4. Numerical Results
(
)
∂2
∂
∂2
Jζ = − r2 2 + r+ψ∂r∂r ∂β 2

Calculation was made for the case Gr = 1, P r = 0.7 .

(8)Enegy Equation:∂T∂ (T, ψ)
1
√
J
+r
=
∂t
∂ (r, β)P r Gr

(
)
2
∂
∂2
2 ∂
r
+r
+
T,
∂r2
∂r ∂β 2
(9)
Fig.3 Steady-state streamlines at
Gr = 1(∆ψ = 0.0001)

where ζ : vorticity, ψ : stream function, Gr : Grashoﬀ∂(x,y).
number, P r : Prandtl number, J ≡ ∂(α,β)
The following Fourier decomposition is used:ζ (r, β, t)=∞∑ζck (r, t) cos kβ

k=0
∞
∑

+
ψ (r, β, t)=ζsk (r, t) sin kβ

(10)

k=1
∞
∑

Fig.4 Steady-state isotherms at Gr = 1(∆T = 0.2)

ψck (r, t) cos kβ

k=0
∞
∑

+

ψsk (r, t) sin kβ

5. Conclusion

(11)

k=1

T (r, β, t)

=

∞
∑

k=0
∞
∑

+
J (r, β)=(12)References
Y. MOCHIMARU, ’Eﬀectiveness of a Spectal Finite
Diﬀerence Scheme’, Computional Fluid Dynamics Review 1998, Vol.1, 379-394, (1998).

Jck (r) cos kβ+=

Tsk (r, t) sin kβ

k=1
∞
∑
k=0
∞
∑

y(r, β)In this study, a heat and ﬂuid ﬂow ﬁeld in an arch
vessel was analyzed numerically, using a spectral ﬁnite
diﬀerence scheme. This result is symmetrical.

Tck (r, t) cos kβ

Jsk (r) sin kβ

(13)

k=1
∞
∑

Yck (r) cos kβ

k=0
∞
∑

+

Ysk (r) sin kβ

(14)

k=1

2
22

ドップラーライダーで観測された流れ場の分類と特性の把握
学籍番号：09-00985

氏名：有場次郎

1. はじめに指導教官：神田学（c） Mixed （Mixed streak and fish net）
；視線
方向風速の正負の境界が乱れているが風向は大体一定である．ストリークは不明瞭である．
（d） Fish net；視線方向風速が正の領域と負の領域が入り混じっている．風向も定まらない．
（e） Front；SNR の分布図に明瞭な境界線が確認できる．
（f） The others；上記グループのいずれにも属さない．

大気境界層とは地表面の影響（摩擦，放射など）を強く受ける大気層のことを指す．大気境
界層では乱流の働きが活発であり，大気や物質の輸送が盛んである。一方で我々は基本的にこの層内で生活しており，大気境界層は都市における汚染物質拡散や局地的豪雨のメカニズムと密接な関係を有している．したがってこれらの問題を解決するためには，大気境界層内の乱流の挙動を知る必要がある．
大気境界層の流れ場を観測する測器の 1 つとしてドップラーライダーがある．ドップラーラ
イダーとは，大気中のエアロゾルに向けてレーザー光を発射して，その反射光を受信して風速を計測するリモートセンシング装置である．ドップラーライダーの観測により都市上空には，
水平ロール渦やプリューム構造といった組織的な乱流構造が存在していることが示された[1]．このような乱流構造の発生条件を解明するために本研究では，ドップラーライダーで観測された流れ場を 6 グループに目視分類し，3 つの大気指標と比較することで乱流構造の特性を探った．

2. 流れ場の目視分類法
解析の対象は，2012 年 10 月 3 日から 10 月
13 日の間にドップラーライダーの水平スキャ
ンによって観測された計 2,447 枚の画像である．
ドップラーライダー画像には視線方向風速分
布と信号対雑音比（SNR）
分布の 2 種類がある．
本研究では流れ場を“Streak”，“No streak”，“Mixed”，“Fish net”，
“Front”，“The others”の 6 グループに目視分類した．以下に各グルー
プの分類基準を，図 1 に各グループの代表的な
視線方向風速および SNR の分布を示す．
（a） Streak；視線方向風速の正負の境界が直線状で風向が一定である．風向に沿った
ストリークが目視できる．
（b） No streak；（a）と同様に視線方向風速の正負の境界が直線状で風向が一定．スト
リークは目立たない．

図 1: 目視分類された 6 グループの代表的な視線方向風速および SNR の分布図（e）Front は SNR 分布，それ以外は視線方向
風速分布である．視線方向風速は中心から遠ざかる向きが正で，中心に向かう向きが負である．

23

3. 目視分類と大気指標との関係

目視分類と 3 つの大気指標との関係を散布図．安定時（a）（c）において，水で示す（図 2）
平風速が大きいときは Streak が，𝑧 ⁄ 𝐿の絶対値が大きいときは No streak が発生する傾向があった．また不安定時（b）（d）において，水平
風速が大きいときは安定時と同様 Streak が，𝑧 ⁄ 𝐿の絶対値が大きいときは Fish net が発生する傾向があった．また Streak と Fish net のグ
ループの中間に Mixed が分布した．
次に境界層高度の大きさに着目する．境界層
高度の大きさによって明確に違ったふるまいを見せたのは No streak および The others のグループである．両者は境界層高度が 1,000m 以下の と きに は多 く発 生した が ，境 界層 高度 が1,000m 以上になると事例数は極端に減った．一方 Streak，Mixed，Fish net のグループは，境界層高度が変化しても事例数に大きな変化はな
かった．

本研究において目視分類と比較を行ったのは水平風速，大気安定度𝑧 ⁄ 𝐿，境界層高度の 3つである．水平風速とは，流れ場の水平方向の平均的な風速である．本研究では水平スキャンの視線方向風速分布のデータに，VAD 法[2]という手法を適用して算出した．
大気安定度𝑧 ⁄ 𝐿 とはモニンオブコフ相似則[3]における説明変数であり，この値の絶対値が大きいほど浮力の影響が，絶対値が小さいほど摩擦による影響が卓越している．またこの値が正だと大気は安定であり，負だと大気は不安定である．𝑧 ⁄ 𝐿はドップラーライダーと同期観測した超音波風速計のデータから算出した．
境界層高度とは大気境界層の最上部，すなわ
ち地表面の影響が及ぶ最高高度のことである．
本研究では，ドップラーライダーを鉛直断面内でスキャンさせて得た，SNR の鉛直プロファイ
ルから境界層高度を算出している．なおここでは触れないが，境界層高度を算出する際には，従来の方法に独自の改良を加えて行った．

4. おわりに目視分類と水平風速，𝑧 ⁄ 𝐿の比較より以下のことが分かった．Streak は風速が大きく摩擦の影響が強いときに，Fish net は大気が不安定なときに，No streak は大気が安定なときに，Mixed
は摩擦と浮力の影響が均衡しているときにそれぞれ発生しやすい．また No streak および The
others の境界層高度が低かったことから，The
others に分類された流れ場にも一般的性質があり，そして No streak もその性質を共有していることが示唆された．

参考文献[1] 藤吉康志，山下和也，藤原忠誠，2005：3次元走査型コヒーレントドップラーライ
ダーによる大気境界層の流れの可視化，天気，52, p.665-666．
[2] Browing, K. A., Wexler, R., 1968: The
Determination of Kinematic Properties of aWind Field Using Doppler Radar, Journal ofApplied meteorology and climatology, 7,
p.105-108.
[3] Monin, A.S., Obkhov, A.M., 1954: Basic lawsof turbulent mixing in the ground layer of theatmosphere, Tr. Akad. Nauk SSSR Geophiz.
Inst.,24, p.163-187.

図 2: 目視分類と 3 つの大気指標との関係
𝒛 ⁄ 𝑳の符号と境界層高度𝒛𝒊 によって次の 4 ケースに分けた．（a）安定時で𝒛𝒊 が小さい，（b）不安定時で𝒛𝒊 が小さい，（c）安定時で𝒛𝒊 が大きい，（d）不安定時で𝒛𝒊 が大きい．いずれの散布図
も横軸が水平風速，縦軸が𝒛 ⁄ 𝑳の絶対値で，両軸ともに対数表記である．

24

CoOx-CeO2 mixed oxide catalysts for N2O decomposition
Student Number: 09_04003 Name: Huiting WANG Supervisor: Hirofumi HINODE, Chris SALIM

1. Introduction

2.2 Catalytic activity experiment
The catalytic activity test was carried out ina fixed-bed flow reactor under atmosphericpressure. The reactant gas containing 1000
ppm N2O, O2(0% or 10%) and He(balance gas)
was fed to the reactor at a flow rate whichcorresponded to a space velocity of 16000h-1.
Gas chromatography (GC323w; GL Science
Co., with Porapak N, Porapak Q and Molecular
Sieve 13X column) was used to analyze N2O.

Nitrous oxide (N2O) is one of the maingreenhouse gases with an average life time ofabout 150 years in the atmosphere, and a net
greenhouse effect of about 310 times and 21
times of the Global Warming Potential(GWP)
of CO2 and CH4, respectively[1]. Moreover,N2O is also identified as a contributor to thedestruction of ozone layer in the stratosphere.
The abatement of N2O emissions fromindustrial plants(adipic acid production,
combustion process, etc.) is probably the mostfeasible by implementing catalytic processes[2].
Catalysts which contain cobalt oxide, are thesubject of studies as potential N2O
decomposition catalyst[3]. Combination of
CeO2 with other metal oxides often affects themobility of oxygen on their surfaces, thus insome cases, modifies the elements’s redox
ability and catalytic performance[4].
In our laboratory, Co/TiO2 catalysts using
impregnation method and TiO2-CeO2 catalysts
using manual mixing method for N2O
decomposition has been studied[5-6]. The
addition of CeO2 improved the catalytic
activity of TiO2. Considering both of theresearches, in this work, the decomposition of
N2O over CoOx-CeO2 mixed oxide catalysts wasinvestigated.

2.3 Catalyst characterization
Characterizations of the catalysts were
performed using TG analyzer (TG8120), XRD(MultiFlex), BET(Autosorb1)
and
SEM(JSM-5310LV) -EDS (JED-2140).

3. Result and discussion

Figure 1 shows the N2O conversion by
CoOx-CeO2 ① (1:1) calcined at different
temperatures for the decomposition of N2O
without O2. For all catalysts, the N2O
conversion were almost the same. The
CoOx-CeO2①(1:1) calcined at 500℃ and 600℃
showed a slightly better catalytic activity thanothers. N2O began to decompose attemperature about 300-400℃, and achieved
100% of N2O decomposition at about
650-700℃.

2. Experimental

2.1 Catalyst preparation
CoOx-CeO2 mixed oxide catalysts wereprepared by manual mixing method. The
precursor of cobalt oxide is Co(CH3COO)2 ・
4H2O (Kanto Chemical Co., Inc.). And two
types of CeO2 (JRC-CEO-1, Santoku Co. and
JRC-CEO-2, Daiichi Kigenso Kagaku Kogyo
Co., Ltd.) were used.
CoOx was prepared by calcining cobalt
precursor at 400℃ for 5 hours under air flow.
CoOx and CeO2 were mixed manually in amortar using small amount of ethanol for
better mixing. The mixture was then calcined
at a certain temperature between 400-800℃
for another 5 hours under air flow.
In order to reduce the pressure drop, thecatalysts were pelletized, crushed and sieved
into 0.71 mm to 1.00 mm particles prior toactivity test.
The following nomenclatures for the catalyst
samples are used: CoOx-CeO2 –X (Y:Z) where X
represents the type of CeO2(1 for JRC-CEO-1,
2 for JRC-CEO-2), and (Y:Z) meanstheweight ratio of CoOx:CeO2. For example,
CoOx-CeO2 ② (1:2) means CoOx and JRC–CEO-2 were mixed in the weight ratio of 1:2.

Fig. 1 Catalytic activity of CoOx-CeO2①(1:1) for thedecomposition of N2O without O2.

Figure 2 shows the N2O conversion by
CoOx-CeO2②(1:1) for the decomposition of N2O
in the condition without O2. For all catalysts,
the N2O conversion were almost the same. The
CoOx-CeO2②(1:1) calcined at 400℃ and 500℃
showed a slightly better catalytic activity than
CoOx-CeO2 ② (1:1) calcined at others. N2O25

temperature in the presence of oxygen showed
in Figure 4. The reason for that could beconsidered that O2 competed with N2O in theprocesses of adsorption on active site[8].

began to decompose at temperature about350-400 ℃ , and achieved 100% of N2O
decompositionatabout
600-650 ℃ .
CoOx-CeO2 ① (1:1) and CoOx-CeO2 ② (1:1) did
not show obvious changes when calcined at
different
temperatures,andcatalysts
containing JRC-CEO-2 showed better catalytic
activities
thancatalysts
containing
JRC-CEO-1.

Fig.4 Effect of Oxygen on the catalytic activity of
CoOx-CeO2②(5:1) for the decomposition of N2O

4. Conclusions

Fig. 2 Catalytic activity of CoOx-CeO2②(1:1) for thedecomposition of N2O without O2.

The CoOx-CeO2①(1:1) calcined at 500℃ and
600 ℃ showed a slightly better catalytic
activity than others. The CoOx-CeO2 ② (1:1)
calcined at 400℃ and 500℃ showed a slightly
better catalytic activity than CoOx-CeO2②(1:1)
calcined at others.
Catalysts containing JRC-CEO-2 showed
better catalytic activities than catalysts
containing JRC-CEO-1.
CoOx-CeO2②(5:1) without O2 showed the best
catalytic activity with the N2O decomposition
achieving 100% at around 550℃.
N2O decomposition shifted to higher
temperature in the presence of oxygen.
References[1] S. Kannan, Appl. Clay Sci. 13 (1998) 347.
[2] P. Esteves, Y. Wu, C. Dujardin, M.K.
Dongare, P. Granger, Catal. Today 176 (2011)453.
[3] C. Ohnishi, K. Asano, S. Iwamoto, K.
Chikama, M.Inoue, Catal.Today 120(2007)145.
[4] S. Imamura, M. Shono, N. Okamoto,
R.Hamada, S. Ishida, Appl. Catal. A
142(1996)279[5] P. Gao, Master thesis, Tokyo Institute ofTechnoligy(2011)[6] Y. Fukuchi, Bachelor thesis, Tokyo Instituteof Technology(2012)[7] E. Iwanek, K. Krawczyk, J. Petryk, J.W.
Sobczak, Z. Kaszkur, Appl. Catal.B
106(2011)416[8] S. Suarez, M. Yates, A.L. Petre, J.A. Martin,
P. Avila and J. Blanco, Appl. Catal.B
64(2006)302

Fig.3 Catalytic activity of CoOx-CeO2 ② in various
weight ratios calcined at 400℃without O2

Figure 3 shows the N2O conversion by
CoOx-CeO2 ② in different weight ratios ofoxides of Co and Ce and calcined at 400℃.
CoOx-CeO2 ② (5:1) showed better catalytic
activity than other catalysts with 100% of N2O
decomposition at 550℃. Mixed oxides showed
higher catalytic activities than each bare
JRC-CEO-2 and CoOx. The higher CoOx
content, the higher the catalytic activity is. As
the cobalt content increases, the beneficial
effect of the cobalt oxide-ceria interface on theactivity of the catalyst decreases. At the sametime, the share of the activity of cobalt oxide in
the overall activity increases[7].
N2O decomposition shifted to higher

26

Synthesis of Zeolites from Lake Sludge as Ammonium Adsorbent Materials
Student number: 09_04055 Name: Wang Wenjing Supervisor: Hirofumi HINODE, Chris SALIM
1. Introduction

3M NaOH

Lakes are important sources of fresh water,
however contaminants and nutrient salts areeasily accumulated because of the closed water
systems and human activities around thelakes.

2.5M NaOH Q: Quartz (SiO2)
2M NaOH

P1: Zeolite P1(Na6Al6Si10O32・12H2O)

1.5M NaOH U: Unnamed zeolite

One of the most effective methods used torehabilitate lakes is lake dredging. But theusual disposal method of dredged lake sludge isby landfill or incineration. For resource
protection, the method of recycling should beconsidered. According to previous study, lakesludge which contains silica and alumina canbe used as raw material to synthesize zeolites
by hydrothermal treatment [1]. Zeolites are
considered as good cation exchange materials,
and they are used for ammonium ion removal
in recent years[2].

(Na6(AlSiO4)6・4H2O)

1M NaOH2θ [deg] (Cu-Kα)
Fig.1 X-ray diffraction patterns of zeolites treated withdifferent alkali concentrations (Si/Al=1.9, 120℃)

150℃

In this study, zeolites from lake sludge wereprepared using hydrothermal synthesis method
in alkali solution. Characterization andNH4+-adsorption studies were conducted todetermine which among the zeolites preparedhas the best cation-exchange capacity.

Q: Quartz (SiO2)
P1: Zeolite P1

120℃(Na6Al6Si10O32・12H2O)
A: Analcime C

100℃(Na(AlSi2)O6・H2O)

2. Experimental
2θ [deg] (Cu-Kα)

In this study, the sludge from Kasumigauralake in Yibaraki was used. The sludge wasdried, crushed and sieved (–150μm), and thendried again at about 100℃ for 24 hours. Si/Al
ratio was adjusted from 1.9 to 3.0 by addition of
Na2SiO3 solution. For zeolite synthesis, lakesludge (1.6g) was treated in different alkali
solution (1, 1.5, 2, 2.5, 3M 8ml) and put insideshaker for 24 hours in a teflon reaction vessel.
Afterwards, the teflon vessel was put into a
metal container, and then placed in the oven at
100℃, 120℃, 150℃. Afterwards, the zeolites
were washed with deionized water and dried atroom temperature.

Fig.2 X-ray diffraction patterns of zeolites synthesized at
different temperatures (Si/Al=1.9, 2M NaOH)

Si/Al=3.0

Q: Quartz (SiO2)
P1: Zeolite P1

Si/Al=2.5(Na6Al6Si10O32・12H2O)
A: Analcime C

Si/Al=1.9

After synthesis, the characterization study
was done by XRD, SEM, and ICP. Finally, thecation exchange capacity was measured by
NH4+-adsorption using back titration method.

(Na(AlSi2)O6・H2O)

2θ [deg] (Cu-Kα)
Fig.3 X-ray diffraction patterns of zeolites adjusted todifferent Si/Al ratios (120℃, 2M NaOH)

3. Results and discussion
3.1 XRD analysis

Figures 1, 2 and 3 showed that quartz was
gradually dissolved in increasing concentrationof alkali solution, temperature, and Si/Al ratio.
Figures 2 and 3, showed that both increasing
temperature and Si/Al ratio can convert Zeolite
P1 (Si/Al=1.67) to Analcime C (Si/Al=2) which
has a Si-rich structure.

The XRD analysis of the zeolites synthesized
at different alkali concentrations, temperatures,and Si/Al ratios are shown in Figures1, 2 and 3.

27

3.2 Cation-exchange capacity
Table.1 shows the Cation-exchange capacity(CEC) values of zeolite products.

Langmuir adsorption isotherm
50

Table.1 CEC values of products obtained on various alkali
concentrations, temperatures, and Si/Al ratios

Qe [mg/g]

Zeolite Preparation
Conditions/Parameter

40
30

CEC value[meq/100g]

phase

1M

131.655

Zeolite P1

10

NaOH

1.5M

159.795

Zeolite P1

0

concentration

2M

170.952

Zeolite P1(Si/Al=1.9

2.5M

127.300

Zeolite P1 &

120℃)

20

0

3M

103.850

Temperature

100℃

150.840

Zeolite P1(Si/Al=1.9

120℃

170.952

Zeolite P1

2M NaOH)

150℃

80.448Analcime-C

Si/Al Ratio

1.9

170.952

Zeolite P1

(120℃

2.5

83.800Analcime-C

2M NaOH)

3.0

97.208Analcime-C

1000

1500

2000

2500

the equilibrium concentration of NH4+ Ce [mg/L]

Unnamed zeolite
Unnamed zeolite

500

Fig.5 the Langmuir adsorption isotherm of Zeolite P1

4. Conclusions
The recycle of lake sludge is successfullycompleted in this study. Types of zeolites can bechanged by adjusting alkali concentrations,temperatures, and Si/Al ratios. A relatively
higher NH4+ adsorbent Zeolite P1 wassuccessfully synthesized from lake sludge byhydrothermal synthesis. According to Fig.5, theexperimental values fit the Langmuir isotherm
very well. By using the Langmuir adsorption
model, the maximum NH4+-adsorption capacityof Zeolite P1 is measured at 46.729 mg/g. This
result is better than those of NaA (44.3mg/g)
and faujasite (37.45mg/g), so Zeolite P1 can beconsidered as a relatively higher quality
adsorbent for ammonium ion.

The highest CEC value was measured at170.952 meq/100g (Zeolite P1) which wastreated with 2M NaOH, and synthesized at
120℃ without Si/Al adjustment. In Fig.1, thissample got the highest intensity in XRD plot
which signifies that the crystal structure wascompletely formed and this is considered as thereason of high CEC. In Fig.2 and 3, Analcime C
was generated. The more substitution of Al for
Si, the larger charge deficiency occurs, and themore cations required for electrical neutrality[3].
Because the Si/Al ratio of Analcime C is morethan the one of Zeolite P1 (Si/Al=1.67),
Analcime C showed a lower CEC.

Reference:[1] Yan Shao, Synthesis of Zeolites from Lake
Sludge,International
Development
Engineering,
Schoolof
ScienceandEngineering, Tokyo Institute of Technology[2] Ayla Arslan, Sevil Veli, Kocaeli University,
Turkey (2011)

3.3 Ammonium ion adsorption[3] Japan Association of Zeolite, Zeolite
conference home page,
http://www.jaz-online.org/, Retrieved in
January 13, 2013

In this study, the adsorption time is 3h, andthe amount of adsorption was measured byback titration method. The maximum
NH4+-adsorption capacity of Zeolite P1 isdetermined by Langmuir adsorption isotherm
at 46.729 mg/g, and the coefficient ofdetermination is 0.9985. Langmuir adsorption
isotherm is plotted in Fig.5.

28

衛星雲画像を用いた海風前線の推定および都市が海風に与える影響
学籍番号：09-04196

1

氏名：大久保

はじめに

洸平

指導教官：神田学

測高度が異なっていたため、解析に先立ち高度補正

海風は高温化した都市を冷却し、汚染物質を輸送を行っている。

し、またゲリラ豪雨の引き金にもなりうるため、そ雲画像追跡法では、輝度値変化の大きいところをの挙動を把握することは重要である。都市化と、海海風前線と推定した。図 1 において今回着目する東風侵入遅延について、蒲生 1)、大和ほか 2)は、それ

京湾、相模湾から北上する前線は赤線（図 4 では緑

ぞれ東京上空におけるブロッキング効果、都市風下

色）で、その他の前線については紫色で示してある。

域での前線の停滞を、地上観測データを用いること

次に地点法との対応を見る。図 2 から、風速、風向により指摘している。しかし、都市化と海風前線のが前線を境として海風域では風速が大きく、風向が内陸侵入の実態をより詳細に把握するためには、現

海風（南風）モードになっていることが分かる。ま状の地上観測データのみでは、時空間精度が十分でた図 3 から前線上で容積絶対湿度の上昇が見られる。

あるとは言い難い。そこで、本研究においては、時以上より、雲画像追跡法により推定した海風前線位

空間解像度の高い衛星雲画像を導入し、かつ地上気

置は、地点法と対応が取れることが分かった。また、象データも用いることで、都市が海風前線侵入に与いくつかの海風イベントでは、相模湾海風前線上でえる影響について明らかにすることを目的とした。

雲が確認されず、雲画像追跡法で前線を引くことができなかったので、地点法から推定した前線を重ね、

2

雲画像追跡法と地点法との対応について

次章でのパターン分類を行った。

この章では、衛星雲画像を用いて海風前線位置を3推定する雲画像追跡法を提案する。そして雲画像追海風前線の時空間発展のパターン分類

跡法と、従来の方法である地点法（地上観測物理量

解析対象日は 2011 年、2012 年の 7，8，9 月の１の時系列データの変化から前線到着時刻を推定する

１事例である。全日、既存研究 1)を参考に抽出した方法）とを比較し対応を見る。

海風日に含まれる（日照時間 7 時間以上、日中、降

データは、雲画像追跡法においては衛星雲画像、

水なし、日平均風速南西から南東）。これを、海風前

地点法においては、関東 7 都道府県の AMeDAS 及線の時空間発展の仕方に着目し 4 種類に分類した。

び環境省大気汚染局物質広域監視システムの地上観

以下に各パターンの特徴を記す。

測データを使用した。なお風速は観測地点ごとに観①

図 1 雲画像と、推定された前線の重ね合わせ(2011 年 9 月 14 日 14 時）赤丸はさいたまを表す。緑丸は北から熊谷、練馬、東京を表す。

2 海風収束による楔形前線の楔形が一旦解消さ𝑚 𝑠
図 2 風系場と、推定された前線の重ね合わせ（図 1 と同時刻）

29

𝑔 𝑚3 /ℎ𝑜𝑢𝑟
図 3 絶対湿度の変化量と、推定された前線の重ね合わせ（図 1 と同時刻）

分類 1、分類 3 に注目し説明していく。

れ、さいたま付近において再び楔形前線を形成する（図 4(a),(b),(c)にて例証)

図 4 に分類 1 の 2011 年 9 月 14 日(a),(b),(c)、分2 海風収束による楔形前線が見られず、さいた

類 3 の 2012 年 7 月 28 日(d),(e),(f)の海風前線と領域

ま付近において初めて楔形前線が形成される内の平均気温（図中右上[℃])からの偏差、そして風2 海風収束による楔形前線が、その形状のまま

向風速とを重ね合わせた時系列データを記す。(a)さいたま方向に移流する（図 4(d),(e),(f)にて例

で、2 海風の収束による楔形前線が確認できるが、証)(b)には一旦その楔形は解消される。(c)にさいたま付

2 海風収束による収束雲がその位置にとどまる

近で再び楔形前線が確認される。(c)の楔領域におい

ここで 2 海風収束とは、東京湾と相模湾からの海て、風速低下、及び風の収束が生じ、また高温化し風の収束のことを意味する。それぞれの事例数は、ていることが確認される（図中の黄丸)。風速低下領順に 4，1，4，2 であった。

域は，高層建物群による運動量粗度の増加領域に対

②

③

④

ここからは、海風前線侵入に対する都市影響を議応している。つまり、(c)に形成される楔形前線は、論する上で重要だと考えられる、さいたま付近での2 海風収束によるメソ影響ではなく、都市（局所）

楔形前線に着目し、その形成過程を議論するために、影響による海風遅延の結果と考えられる。一方、図(a)14:00

32.4(d)13:00

4(d),(e),(f)を見るといずれの時間にも、2 海風収束に33.5よる楔形前線が風向きに対応して（図中の黄丸）さいたま付近の高温域に侵入していくことが確認できる。さいたま付近での楔形前線の形成は両ケースに共通するが、後者の事例ではメソ影響か、都市影響（熱空気力学的作用に伴う収束）かの判断が難しい。
なぜなら、関東広域での風系場が南東モードになっており、総観場の気圧傾度により移流されたとも考(b)15:00

32.3(e)14:00

33.5

えられるためだ。メソ影響のみ、もしくはメソ影響と局所影響の 2 つの異なるスケールの現象が共に影響している可能性が考えられる。

4

おわりに

本研究では、新たに雲画像追跡法を提案した。そして、前線の時空間発展を分類することで、さいたま付(c)16:00

31.5(f)15:00

32.6

近での楔形の前線形成過程を確認した。また海風侵入に対する都市影響も確認された。海風侵入のそれぞれのパターンに対する、局所影響、メソ影響を明らかにすることが今後の課題である。

参考文献
[1]

蒲生稔，海風前線位置を読み取る方法の比較，

公害，vol. 26, p31-39, 1991
[2] 大和広明，三上武彦，高橋日出男，夏季日中における首都圏のヒートアイランド現象に海風が与え

図 4 2011 年 9 月 1４日の 14 時から 16 時(a),(b),(c)および 2012 年 7 月 28 日の 13 時から 15 時(d),(e),(f)における前線と地点データの重ね合わせ

る影響，地学雑誌，vol. 120，pp. 325-340，201130

Measurement of audible sound due to crushing progression of snail fossils
Student Number: 09_08691

Name: Wenjie WUIntroduction
Sound wave has been uniquely and sensitively detected tinychanges of micro level configuration at the grain level because
vibrations in granular material are qualitatively different fromthose found in solids and liquids1-2). Crack development andparticle crushing in brittle materials can also produce acoustic
effect; therefore, acoustic emission (AE) has emerged as in-situ
techniques to evaluate rock slope stability3). Audible sound ofcrushing progress can be observed, and might be used tomonitor crack developments inside the slope of brittle materials.
The large deposit of thickness about 10 meters of snail fossils
aging around 12-13 million years has been discovered at theMae Moh open-pit lignite mine in Thailand since 2003. Due tocrushability of snail fossils, audible sound of crushing progress
can be observed and used to monitor crack developments inside
the slope of snail fossils. This research aims to find thecharacteristic of the crushing sounds of snail fossil particles
under various conditions.

Supervisor: Thirapong PIPATPONGSA

1.

Table1 Composition of the sample
Sieve No.

Diameter (mm)

Mass (gram)

#20
#40
#60
#100
#200

0.85
0.425
0.25
0.15
0.075

30.3
19.7
9.0
4.6
2.4

Percent
retain
46 %
30 %
14 %
7%
3%

Table 2 Experiment program where σ´ is vertical stressCondition
Shearing
Sample No. σ´ (kPa)OCR
mode
B-DCL1450
1Constant
B-DCL15
200
1
load
Dry
B-DCL16400
1Constant
B-DCV11400
1volumeB-WCL150
1B-WCL250
8B-WCL350
4B-WCL4100
1Constantload
B-WCL5200
1
B-WCL6200
2B-WCL7300
1
B-WCL8400
1Wet
B-WCV150
1B-WCV250
8B-WCV350
4B-WCV4100
1Constant
volume
B-WCV5200
1
B-WCV6200
2B-WCV7300
1
B-WCV8400
12. Review of the previous study
The snail deposit at the Mae Moh open-pit mine, which isunder operation of Electricity Generating Authority of Thailand(EGAT), is regarded as the thickest shell bed in Thailand4). The
preservation area of snail fossils located a southwest pit has been
assigned to prevent damages from mining activities. The
detailed studies of characteristics of snail fossil are required inorder to preserve this geological heritage while developing theenergy resources. As being brittle inhomogeneous material, thelayer of snail fossils would lose its strength after exposing toexternal environment for long time. Therefore, non-invasive
monitoring using audio measurement was carried out in alaboratory level for measurement and analysis of crushing soundunder compressive loading using commercial microphone5).
Based on the previous study5), crushing sounds of snail fossils
under compressive loading were recorded by microphones withconnection to stethoscopes attached to the oedometer. Noisefiltering technique and sound wave analysis were conducted. At
laboratory level, linear relationship between the cumulative
counts and cumulative settlement can be confirmed.

Direct shear tests were performed by 2 steps: consolidationand shear. Consolidation test is undertaken by applying a normalpressure to a sample and measuring the vertical settlement.
Shear test is undertaken by applying shear force and measuring
the horizontal displacement. After complete the experiment, thesample was removed and followed by sieving analysis to detecta degree of crushing. Twenty sets of direct shear test in total
were conducted in dry and wet conditions under constantpressure and constant volume with different values ofover-consolidation ratio (OCR). The amount normal pressure isvaried by 50, 100, 200, 300, 400 kPa. After the end of
consolidation, the sample was sheared with a shear rate of 0.2
mm/min to reach a shear displacement of 8 mm. Moreover,crushing sound of fossils generated under the consolidation test
and shear test were recorded simultaneously by 3 sets ofstethoscope connected to miniature microphone attached to theplastic frame of shearing box and direct shear apparatus asshown in Fig. 1.

3. Experimental program
In extension to the previous studies which used snail
specimens, this research uses snail fossil particles as the
specimen, which were collected from subsurface outside thepreservation area and prepared by water sedimentation
technique in order to make the sample uniform. Each sample
was carefully prepared with the exact ratio of percent finer byweight of each particle size following the grading of the snail
fossils after removal of big particle as summarized in Table1.

4. Experimental results
The results of direct shear test are shown in Fig. 2 forcomparison between shearing modes in constant load andconstant volume. Similarly, comparison between samples
prepared in wet and dry conditions also show the same result;
therefore, indicating that no matter the shear test was conductedin constant load or volume condition, or whether the sample is
wet or dry, the obtained angle of friction is the same.
Sound wave data were filtered by low-pass filter (LPF) at 8
kHz, resampled with frequency 16 kHz down from 44 kHz,
finally filtered by high-pass filter (HPF) at 3 kHz using software

Fig. 1 Direct shear apparatus with 3 sets of stethoscope andmicrophone connection to three various locations
31

same distance. GOLDWAVE was used to observe the soundwave of B-DCL16. Figs.5 and 6 visualized the sound wave forboth the original and the processed data started from 29′02″ to29′18″. For comparison, the processed data based pm soundwave recorded by Channel 2 is in Fig 7. It can be considered thatthe sound wave caused by the direct shear apparatus might leadto the result of linear relationship shown in Fig.8.

GOLDWAVE. By setting a threshold value, a MATLAB script
applying acoustic emission algorithm counts the pulses which
are over the threshold line. Under consolidation test, a linear
relationship between cummulative count and settlement of dry
samples can be confrimed as shown in Fig. 3. The reason why
the result of B-DCV11 was totally different others although theprocess were exactly the same is possibly a run out of battary of
microphone. No particular relationship between cummulative
count and settlement of wet samples under consolidation test as
shown in Fig. 4. The crushing sound wave could not be perfectly
recorded by the stethoscope through water and plastic shear box
as sound medium. Some noice due to the sound caused byplacing the weights might be recorded time to time.

Fig. 5 Original sound data of B-DCL16 (CH1) during shearing

Fig. 6 Filtered sound data of B-DCL16 (CH1) during shearing

Fig. 7 Filtered sound data of B-DCL16 (CH2) during shearing
Fig. 2 Comparison between constant load and constant volume

Fig. 3 Cumulative count-Settlement of dry sample
Fig. 8 Cumulative count-horizontal displacement of B-DCL16
5. Conclusion
Linear relationship between the cumulative counts andcumulative settlement can be confirmed on for dry sample ofsnail fossil particles at laboratory level. Under shearing stage,crushing sound was affected by noises of apparatus. Soundfiltering technique is required in the subsequent research.
References
1) K. Tanimoto and Y. Tanaka: Yielding of soil as determinedby acoustic emission, Soils and Foundations, 26(3),pp.69-80, 1986.
2) C. Liu and S.R. Nagel: Sound in sand, Physical Review
Letters, 68(15), pp.2301-2304, 1992.
3) M. Otsu: Characteristics and Theory of Acoustic Emission,
2nd Edition, Morikita Publishing, 2005 (in Japanese).
4) W. Songtham, et al.: Middle Miocene Molluscan
Assemblages in Mae Moh Basin, Lampang Province,
Northern Thailand, ScienceAsia, 31(2), pp.183-191, 2005.
5) S. Kitakata: Crushing sound analysis of snail fossils undercompressive loading using audible frequency measurement,
Tokyo Institute of Technology, Bachelor Thesis, 2012.

Fig. 4 Cumulative count-Settlement of wet sample
Under shear stage, the number of the filtered waves over thethreshold line was counted every 100 second by using MATLAB.
After checking several threshold values, the value around 0.0003
was found suitable. Fig.6 shows the result of B-DCL16 (dry,constant load, shear under 400kPa). A linear relationship isfound between cumulative count and horizontal displacementfrom the sound data recorded by Channel 1, but was extremely
different from sound data recorded by Channel 2 although theresult should be similar because the stethoscopes of Channel 1
and Channel 2 were placed at two sides of the shear box and the32Image Coding with 1D Discrete Wavelet Transform using
Line-based Prediction
Student Number： 09_11888

1

Name： SHIYU ZHOU

3Introduction

Supervisor： Yukihiko YAMASHITA

1D wavelet transform with linebased prediction

The rapid development of information communication technology has enabled to use not only texts but
Since an image is a 2D signal, 2D DWT is used foralso still images or videos for communication. Howimage coding. An input image (I) is split into approxiever, the amount of information of digitized images andmated and detailed coefficients both vertically and horvideos is very huge. Therefore, the image coding methizontally. Fig. 1 (a) shows the coefficients of DWT. For
ods have been invented to reduce irrelevance and redunexamples, (H, L) indicates horizontally detail and vertidancy in image data and to store or transmit them in an
cally approximation coefficients of the input image (I),efficient form.
and (LL, LH) indicates horizontally approximation and
JPEG is one of the commonly used methods of lossy
vertically detail coefficients of (L, L).
compression for still images [3]. In the coding proHowever, if we introduce prediction to 2D DWT, thecess of JPEG, an image is usually split into blocks ofprecise prediction is difficult because of a distance of8 × 8 pixels, and each block is transformed by the disseveral pixels between the already decoded pixels andcrete cosine transform (DCT). Therefore, JPEG suffers
pixels to be decoded. To solve this problem, I propose
from a problem that visual degradation called block disa new usage of DWT, the 1D DWT, in which they are
tortion is observed in its decoded image. JPEG2000
adjacent. In 1D DWT for images, horizontal or vertisolved the problem by using the discrete wavelet transcal transform is performed alternately. Fig.1 (b) shows
form (DWT). Since DWT has no blockwise process,the coefficients of DWT. For examples, (H, I) indicates
DWT can solve the problem of block distortion [2, 4, 5].
horizontally detail and vertically not transformed coefPrediction is another key technology to enhance theficients of the input image (I).
performance of image coding. Since it is difficult tocombine prediction with the 2D DWT used in JPEG2000,
in this paper, I propose to use 1D DWT and combinewith prediction for image coding, develop a prediction
program, and evaluate it by experimental results.
(a) 2D DWT2(b) 1D DWTDiscrete wavelet transform (DWT)
Fig. 1: Coefficients of DWT

Daubechies proposed DWT for multiresolution signal analysis [1]. Its analysis process is as follows. An
input signal is split into a low and a high frequency components by a low and a high pass filters. Both components are downsampled by 1/2 and the approximationand the detail coefficients are obtained respectively. The
approximation coefficients are iteratively split as an input signal. Its synthesis process is almost the oppositeof the analysis process. After the approximation and thedetail coefficients are upsampled and filtered by a low
and a high pass filters, respectively, the original signal
is restored by the sum of both signals.

Fig.2 shows the block diagram of prediction process
in the proposed coding scheme. The prediction is doneline by line in a target image. In the prediction process,we only have to predict its detail coefficients, since theapproximation coefficients are given. The already decoded lines are called the reference lines and the line tobe predicted is called the target line.
The detail of prediction process is as follows. 1) Restore the low frequeny component of target line fromapproximate coefficients. 2) Predict the target line byusing the reference lines and the low frequency component of target line (Its detail is described below). 3)

33

Transform the predicted line by 1D DWT and obtain
detailed coefficients. 4) Quantize the detail coefficients.
5) Output the quantized difference between the true detailed coefficients and the predicted ones.

(a) Original image(b) Predictive image

Fig. 3: Original and predicted images
Table 1: Averaged estimation error of Lenna

Fig. 2: Block diagram of prediction process
For the line-based prediction, I used the linear prediction with clustering. First, I extracted blocks which
size is bx × by . Assume that 1D DWT is done horizontally, y is the vertical coordinate of the block. I set
the line at y = by to the target line, and the lines at
y = 1, 2, . . . , by − 1 to the reference lines. Then theblock can be expressed by a (bx by )-dimensional vector
u. Let m be the averaged value of all components in−m1
u, and v be the normalized vector ku
u−m1k , where 1 isa vector of which all components are one. Let p and q
be the value of target pixel and its normalized versionp−m.
ku−m1kFor the precise prediction, I cluster the normalized
blocks according to their features to make the categories,Ωc (c = 1, 2, . . . , C). The k-means method is applied
for the clustering. The prediction model in Ωc is givenbyq = hac , vi + bc ,(1)

5Conclusions

In this paper, I proposed an image coding with 1D
DWT using line-based prediction, and evaluated the proposed prediction method by experimental results. For
future works, I have to develop an image coding system to validate the prediction and investigate prediction
methods to improve the accuracy.

References[1] I. Daubechies, “Ten Lectures on Wavelets,” Society
for Industrial and Applied Mathematics, Philadelwhere ac and bc are constants for linear prediction. For
phia, 1992
training data v i and qi (i = 1, 2, . . . , Nc ) in Ωc , ac andbc are decided to minimize the following criterion:[2] M. Vetterli,J. Kovaevi, “Wavelets and Subband
Nc
Coding,” PrenticeHall, Englewood Cliffs, 1995.
X|qi − hac , v i i − bc |2 .
(2)[3] Sadayasu Ono, Junji Suzuki, “Achievement
i=1
method of comprehensible JPEG/MPEG2,”
Ohmsha, Tokyo, 1995 (in Japanese).

4

Experimental Results[4] C. S. Burrus, R. A. Gopinath, and H. Guo, “inFor experiments, I used the standard image Lenna
troduction to Wavelets and Wavelet Transforms: A(Fig. 4). Table 1 shows the mean squared error of thePrimer,” PrenticeHall, Englewood Cliffs, 1997.
target pixels. From Table 1, we can see the error decreases as the block size and the number of clusters C [5] Geoffrey M. Davis, Aria Nosratinia, “Waveletbased Image Coding: An Overview,” IEEE
increased. However, the errors are still large, we haveTrans. on Circuits and Systems for Video Technolto develop a more precise method.
ogy, vol . no. 6, pp.243-250, June. 2002.

34

A Pareto frontier analysis of household budget and environment for energy saving goods
省エネルギー商品に対する家計予算と環境負荷低減効果のパレート効率性分析

Student ID: 09 12244 Name: Ryoya SUEHARA Supervisor: Naoya ABE
1 Introduction
In 1970s, Japan had suﬀered from the oil crisis, and theJapanese government has made eﬀorts to secure the energy supply. At the same time, energy conservation activities, which has been mainly taken by private companies
for profit, has also become popular among the ordinary
people. The central government and local governments
now enhance the subsidy programs on energy saving device to households to promote energy conservation activities in the residential sector. The subsidized goods are
relatively expensive to conventional goods such as lighting, air-conditioners, etc. Besides, the energy crisis has
become a more serious issue after the 2011 Tōhoku earthquake and tsunami and the subsequent nuclear power
plant accident, the public attention is being focused on
the energy conservation activities. However, Whether
households can aﬀord to buy energy saving goods is depends on their budget.
This paper aims to evaluate the trade-oﬀ between economy and environment. More specifically, the economy
means the saving amount, and the environment means
the amount of CO2 emission reduction of a household. To
handle these aspect, multi-objective optimization problem are solved and the Pareto optimal solutions of eachproblem are obtained. Finally, the fairness of the subsidy
programs in Japan is discussed.
This paper focuses the four energy conservation goods;(1) residential photovoltaic systems(PV), (2) residential
cogeneration systems(ENE · FARM, EF), (3) solar water
heating systems(SWH), and (4) Eco-cars(Prius PHV, EC).
In addition, Ōta ward, one of the Tokyo’s 23 wards isfocused.

2.2 Households’ incomesHouseholds’ decision makings can be influenced by various things in practice. This study, however, is assumingthat their decision makings are influenced by their moneyin their hand; that is, the diﬀerence between their savingsand liabilities (real savings). The data of annual incomeof Ōta Ward, and the savings and liabilities of workers’
households in Japan are used to estimate their budget
constraints(shown in tables 2 and 3 respectively).
Table 2: Annual income of the inhabitants Ōta Ward
Less than 5,000 5,001 to 9,999 More than 10,000
58%
32%
10%
Gini index: 0.412
Table 3: savings and liabilities of workers
Income group
percentile
Homeownership rate

real savings≤ 5,200
40%
60%
2,640

households

5,201 to 9,980
50%
77%
5,820

≥ 9,98110%85%
18,420[Unit: thousand yen]

The half of the inhabitants have 2,640 thousand yen ofreal savings on average. Therefore, the budget constraint
on that income class should be less than 2,660 thousandyen.

3 Methodology
3.1 Multi-objective problem optimization
A multiobjective optimization task considers several conflicting objectives simultaneously. In such a case, there isusually no single optimal solution, but a set of alternatives
2 Subsidy programs and households’ incomeswith diﬀerent trade-oﬀs, called Pareto optimal solutions, or
2.1 Subsidy programs on energy conservation goods
non-dominated solutions. A multiobjective optimization
The inhabitants of Ōta Ward can enjoy the subsidies of- problem are described as follows;
fered by the central government, Tokyo Metropolitan
Government, and Ōta Ward. For instance, the subsidy
minimize { f1 (x), f2 (x), . . . , fk (x)}(eq. 1)on PV is determined by the capacity of the system represubject tox ∈ S.
sented by kW and the same can be said of SWH. EF andEC can not be divided, therefore, the subsidies are deter- involving k (≥ 2) conflicting objective function fi : Rn → R
mined by the number of systems installed. In this study to minimize simultaneously. The decision variable vectors x
to simplify the problem, the average installed capacities = (x1 , x2 , . . . , xn )T belong to the nonempty feasible region
are calculated with respect to PV and SWH. The average S ⊂ Rn . Objective vectors are images of decision veccapacity of PV installed in 2011 is about 4.3kW(J-PEC, tors and consist of objective function values z = f(x) =2011), and the average capacity of SWH is 3m2 (SSDA, ( f1 (x), f2 (x), . . . , fk (x))T . In multiobjective optimization,2012).
objective vectors are regarded as optimal if none of theother components can be improved without deterioration
to at least one of the other components. More precisely,
Table 1: Amount of subsidy for each gooda decision vector x0 ∈ S is called Pareto optimal if there
Goods
PV
EF SWH EC
does not exist another x ∈ S such that fi (x) ≤ fi (x0 ) for all
Subsidizer
i = 1, . . . , k and fi (x) < fi (x0 ) for at least one index j.
Central Gov. 129
470
0 672
Tokyo Metropolitan Gov. 431
0210
0
3.2 Application of the the modelŌta Ward 2160270This study focuses the some trade-oﬀs between two obTotal amount of subsidy
776 470
237 672
jective(See Table.4), for instance, the trade-oﬀ between[Unit: thousand yen per good]
the amount of the savings per year per ten thousand yen

35

(investment) and the amount of CO2 emission reduction. chosen. Therefore, the intersectional solution is choice 3,
the optimization problem is described as;
installing WHS only or installing all goods.
{1/ fa (x), 1/ fb (x)}minimize(eq. 2)

subject tox ∈ S1 .
fa (x) = (s · x)/(p · x)
fb (x) = r · x(eq. 3)5 Conclusions and policy implications
The result shows that the optimal choices varied based on
the household’s strategy and budget constraint. When ahousehold is interested only in the return of investment,
they would choose choices 1and 3, no choice or installing
only WHS. When a household is interested in both thereturn and the environment, they could choose choices 1,
3, 10, 15, 16. The governments seems to reduce the CO2
emission as well as minimize the expenditure. Therefore, the desirable choices are choice 1, 2, 3, 8, 16. As
a result, there are a few choices possibly being chosen.
Moreover, the households who do not have their own
house can not buy PV and SWH, and they mostly belong

Amount of savings per yen(Investment)
per year [yen/yen/year]

It is assumed that the decision variable vector can takeonly 0 or 1. Therefore, S1 is consisted of 24 = 16 decision
vectors. fa (x) represents the amount of savings per year,
fb (x) represents the amount of CO2 emission reduction.
The inverse numbers are taken to resolve the problem
into the model. s is the saving vector (s = [130 55 39 48])and provide the saving abilities of goods(the amount ofsaving per year). p is the price vector (p = [1,600 2,290
63 3,770]) and provide the subsidized prices (price persystem). r is the reduction vector (r = [1,611 1,500 580 829]) to the low income group. The governments should conand provide the reduction abilities of goods(the amount sider the fairness more in the process of policy design andimplementation.
of CO2 emission reduction kg per year).

6 Further study
In this paper, the impact of Feed-in-tariﬀ is not taken into
account to simplify the problem. The households which
install PV can sell their surplus electricity to electric powercompanies at the higher rate than the rate at which theybuy electricity from the power companies. Therefore,
the budget constraints could be changed only when thehousehold install residential photovoltaic systems into
its house. This impact would not be small and should beconsidered.
Moreover, The number of goods focused in this paper
are only four. In fact, households decide their consumption from a great number of goods in their daily lives.
When they purchase appliances, non-subsidized goods
can appear in their options. This change in the choice set
could influence the Pareto optimal choices.

#"
&$$$"

%$$$"

($$$"

#$$$"

'$$$"

!$"

!$$$"

%"
&"
!"

$"
$"

%$$"

("

'"
!$$$"

*"

!%$$"

'$$$"

+"
)"
'%$$"

!!"
!#"
!'"
#$$$"

#%$$"

!%"
!&"

!("
($$$"

(%$$"

%$$$"

References

CO2 Emission Reduction [kg/year][1] Japan Photovoltaic Expansion Center (J–PEC). 2011.

Figure 1: Pareto frontier[2] Solar System Development Association (SSDA). Solar System Data
Book. 2012.

4 Results
The number of optimal choices increases as the budget
constraint is relaxed. For households, the optimal solution is to choose choices 1, 3, 10, and 15. For governments, it is desired when choices 1, 2, 3, 8, and 16 is[3] Branke, J., K. Deb., K. Miettinen, and R. Słowiński. 2008. Multiobjective Optimization: Interactive and Evolutionary Approaches. Germany:
Springer.

Table 4: Results of frontiers
Max Sv Max Sv/Y Max CO2 /Y Max CO2 /SY Max CO2 /SY
DMVActual
&
&
and
and
and
PV EF SWH EC Cost [kY]Min cost Min Cost
min Cost
min Cost
max CO2
1
0
0
0
0
0
◦
◦
◦
2
0
0
0
1
3,770
3
0
0
1
0
63
◦
◦
◦
◦
◦
4
0
1
0
0
2,290
◦
◦
5
1
0
0
0
1,600
◦
◦
6
0
0
1
1
3,833
Never appeared in Pareto frontier
7
0
1
0
1
6,060
◦
8
0
1
1
0
2,353
◦
9
1
0
0
1
5,370
Never appeared in Pareto frontier
10
1
0
1
0
1,663
◦
◦
11
1
1
0
0
3,890
◦
◦
12
0
1
1
1
6,123
◦
13
1
0
1
1
5,433
Never appeared in Pareto frontier
14
1
1
0
1
7,660
◦
◦
15
1
1
1
0
3,953
◦
◦
◦
16
1
1
1
1
7,723
◦
◦
◦
CO2 : the amount of CO2 emission reduction, DMV: decision making vector, Sv: savings, SY: yen subsidized.
ChoiceMax CO2&
Min cost◦

36

Max Sv/Y
andmax CO2◦◦Min SYandmax CO2
◦
◦
◦

◦
◦

◦
◦

◦

Study on the Propagation of the Tidal Waves in the Mekong River
Student Number:09B13746

Name:Yuto TAKAHASHISupervisor:Hiroshi TAKAGI

is mainly determined by astronomical tides. Therefore, it can be

1. Introduction
There are many researches on floods in the lower Mekong

emphasized that an analysis on the flood in lower Mekong Delta

River basin mainly focusing on the influence of the rivershould be carried out in consideration of astronomical tides asdischarge from the upstream of the river, but there has so far

well as river discharges.

been little research that carefully considers the influence of ocean
tides. Due to its low-lying nature, the region is very likely
vulnerable to the influence of sea-level rise. The objective of thisstudy is to investigate the characteristics of tidal waves traveling
upstream.
2. Field survey
The authors carried out the field survey on March 2012 andobtained the bathymetric data and the water level of the MekongRiver in Vietnam. Figs.1 and 2 show scenes of the measurement.
These data are used in the numerical simulation.
Fig.3 The lower Mekong River basin

Amplitude (m)1DinhAn0.1

0.01

Fig.1 Scene of measurement for depth of river (Left)

0.5

1

M2
S2

O1
K1

2

1

Amplitude (m)Fig.2 Scene of measurement for water level (Right)

3. Analysis of the characteristics of the water level in the

10 15

30

days/cycle

180

360

Ssa

Sa

180

360

Ssa

Sa

CanTho0.1

0.01

lower Mekong River

0.5
M2
S2

The purpose of this study is to investigate the characteristics
of water level (WL) associated with tides in Can Tho city which

1

O1
K1

2

10 15

days/cycle

30

Fig.4 Amplitude spectrum (Top: Dinh An, Bottom: Can Tho)is the largest city in the Mekong delta region and located about
80km from the coast. Fig.3 shows the location of Can Tho.

4.

In addition to WL measured by the authors themselves, theNumerical simulation of tidal waves along the river

4.1 Conditions of the simulation

data of water level at every 15 minutes observed by MRC

Since the detailed data necessary for the numerical(Mekong River Commission) was obtained and used for a basic

simulation is scarce in the Mekong River, it is necessary toanalysis. The data analysis have been made using the data at two

postulate some of these data such as river discharge (Q),WL stations in Dinh An which is a small town nearby the coast

Manning’s coefficient of roughness (n), and the depth of river (h).

as well as Can Tho.

Using the data of annual average amount of river discharge

The amplitude spectrum is calculated by the fast Fourier

published by MRC, the estimated river discharges for threetransform for WL between 2009 and 2010. The result is shown inseasons (flood, normal, and dry seasons) are estimated to be

Fig.4. Both cities show that the amplitudes around half day, one

Q=20000, 6000, 3000 m3/s respectively. Manning’s coefficient ofday, and one year are remarkable. It appears that theseroughness is n=0.03 and 0.06, the depth of river is h=20 and 15m

amplifications are caused by astronomical tides with tidal

which are roughly estimated by the bathymetric survey by theconstituents of M2/S2, O1/K1, Sa respectively. It is natural that

authors. Fig.5 shows the computational domain of the Mekong

WL at Dinh An indicates tidal patterns because it is located at a

River and an example of WL during a tide propagates. The tidal

coast. However, it should be noted that WL at Can Tho, which ispropagation was simulated using Delft3D-FLOW originallyabout 80km upstream from the coast, also indicates similar

developed by Delft University of Technology. The results were

results. This means that the fluctuation of water level at Can Tho

compared with the observed WL by the authors on March 2012.

37

Fig.5 The computational domain of the lower Mekong River(10-Mar-2012 11:00).

4.2 Tidal damping

6. Conclusions

・

It is necessary to investigate the interaction mechanism

Computational result with Q=6,000m3/s, h=15m, n=0.03

between river flows and tidal waves in order to estimate WL in

agrees well with the observed data at Can Tho in terms of

the lower river basin precisely. It is known that the amplitude ofboth amplitude and phase.

・

tidal wave is decreased as the amount of river discharge

The numerical simulation shows that the tidal damping isincreases because the velocity of going downstream is added andaccelerated by a narrow part of the river where tidal

the run-up is suppressed [1]. In the present research, thevelocity significantly increases, whereas the influence ofnumerical simulation is performed in order to identify how river

river discharges is not remarkable.

flows and tidal waves affect with each other.

Simulation(n=0.03)
Simulation(n=0.06)
Observed WL
Field Survey

1.5

5. Result and discussion
A series of simulations with changes in Q, n, and h were1.0

2012/03/11 00:00

2012/03/10 18:00

terms of both amplitude and phase.

2012/03/10 12:00-0.5

2012/03/10 06:00

with n=0.03 agrees well with the observed data at Can Tho in

2012/03/10 00:00

0.0

too large in the case of n=0.06, whereas the computational result

2012/03/09 18:00

or n=0.06. The result shows that the tidal damping appears to be

0.5

2012/03/09 12:00

performed for the settings of Q=6,000m3/s, h=15m, and n=0.03

2012/03/09 06:00

Water Level(m)these results are shown in Fig.6 in which the simulation was

2012/03/09 00:00

performed in order to identify the sensitivity to WL. Some ofDate and Time
-1.0

Calculating the amplitude of tidal waves for 20 points atinterval of 5km from the coast to the upper boundary of theFig.6 The result of simulation and the comparison of calculated resultmodel, the damping (or amplification) ratio of the amplitudewith measured value.

(amplitude at a certain point divided by the value at the coast) is0

obtained as shown in Fig.7. Regardless of the difference in the

20

40

60

80

100

120

1.1

Q=20,000
Q=6,000
Q=3,000
Width of River

1.0

The damping ratio of tidal wave

slight increase up to about 50km and starts rapidly decreasing
after 50km. Fig.7 also shows the width of the Mekong River at
each point. It shows that the tidal damping tends to beaccelerated as the width of river become narrow. This trend canbe explained by the fact that the bottom friction that leads todecreasing the amplitude will increase as the tidal velocity getsincrease when the wave travels along a narrow part of the river.

0.9

4000

3000

0.8
0.7

2000
0.6

Width of River (m)river discharge, the tidal amplitude demonstrates a stable or

0.5
1000
0.4

On the other hand, the damping effect by the difference of0.3

the river discharge seems to be unclear. One of the possible

0

20

40

60

80

100

120

Distance from Estuary (m)reason is that the horizontal 2D simulation used in the present
research underestimates the vertical turbulence that may bringFig.7 The damping ratio of tidal wave and the width of the riverabout a significant loss of energy.
References[1] A.C. Horrevoets, H.H.G. Savenije, J.N. Schuurman and S. Graas.
The influence of river discharge on tidal damping in alluvial
estuaries, Journal of Hydrology, pp.1-3, 2004

38

Production of monomers for biodegradable plastics withlignin as raw material
Student Number:09_15277 Name: Ngoc Linh TRAN Supervisor: Ryuichi EGASHIRA, Hiroaki HABAKI
1. Introduction
Bio-base polyesters exhibiting smooth degradation havebeen widely studied as environmentally friendly polymericmaterial. They can be derived from renewable feedstock,
thereby reducing greenhouse gas emissions. For instance,
polyhydroxyalkanoates and lactic acid (raw materials forpolylactic acid) can be produced by fermentative
biotechnological processes using agricultural products andmicroorganisms. However the rate of these fermentative

temperature above 350℃[4] and the temperatures of pyrolysis
were fixed at 500℃, 600℃, 700℃, 800℃ with lignin, and at
600℃ with rice straw, respectively.

reactions is small, so production of these monomers inbulk is still difficult.
Kaneko et al. [1] reported that coumaric acid derivative
homo-polymers which are one type of bio-base polyester
showed smooth in-soil degradation. Figure 1 shows chemical
structures of various coumaric acid derivative monomers andsynthesis of their polymers [1]Figure 2. Schematic diagram of pyrolysis experimental system
Table1 Experimental condition of thermal treatment of lignin
Feed
Lignin (Rice straw)
Mass of feed
[g]
12~15
Atmosphere
Flow rate of[mL/min]50
atmosphere
Flow rate ofmoisture[mL/min]1
Temperature
T [℃]500, 600, 700, 800
Holding timet [h]1Extraction of coumaric acid from liquid products.
At extraction experiment, ether and then aqueous sodiumhydrogen carbonate solution were added in liquid products one
by one. The obtained acids in liquid products should be
converted to sodium salt forms after reacting with sodiumhydrogen carbonate, and dissolved in aqueous phase. Each
phase of this process was analyzed by UV spectroscopy. The
experimental condition of extraction is showed at table 2

Figure 1: Chemical structures of various coumaric acid derivative
monomers and synthesis of their homo-polymers

These coumaric and ferulic acids are contained as the main
phenyl propane construction monomers in lignin [2]. Lignin is
one of the major substances in plant together with cellulose andhemicellulose; however there are a great quantity of lignin
residue in the world such as rice straw and pulp industry
residue. Production of biodegradation plastics from lignin
residue is one of effective methods to resolve the environmental
problems.
In this work, the preparation of coumaric and ferulic acids
from lignin is studied. The pyrolysis of lignin under steam
atmosphere was conducted to decompose the complexity
construction of lignin structure to simple substances. Liquid
products from pyrolysis of lignin were analyzed by UV
spectroscopy and high performance liquid chromatography(HPLC) to study the existence of coumaric acid and ferulic acid.
After that, extraction experiment was conducted to recover theobtained acids. Moreover, the pyrolysis of the agriculture
residue of rice straw, which contains about 19~21% lignin by
mass, was also studied.

Table2 Experimental condition of extraction of liquid products
Feed
Liquid productsAqueous phase
Aqueous solution ofOrganic phase
Ether
Concentration aqueous solution
ofCNaHCO3 ,aqe0[kmol/m3]
0.01
Volume of feed
Vfeed [m3]3.0×10-5
Volume of aqueous phase
Vaq [m3]
3.0×10-5
Volume of organic phase
Vorg [m3]3.0×10-5
Temperature
T [℃]
298

Analysis of liquid products. Liquid products from thepyrolysis were analyzed by UV Probe 2.43 and HPLC
(Shimada. Co). At HPLC analysis, the mobile phase ismethanol. The mobile phase was delivered at a rate of 0.5
mL/min. The detection was monitored at 280nm.

2. Experimental
Material. Lignin powder was purchased from Kasei Kanto
Co. Coumaric acid and ferulic acid were purchased from Wako
pure chemical Co. Rice straw was obtained from Hiroshima.
Dried rice straw was ground and sieved to obtain particles
ranging between 1.7 and 0.35 mm.
Pyrolysis. Figure 2 shows the schematic diagram ofpyrolysis experimental system. The sample holder was filledwith lignin or rice straw, and placed at the center of thepyrolysis reactor. Nitrogen and water were supplied during thepyrolysis. All of the gases flowing out from the end of the actor
were sent to the liquid product trap to obtain the liquid productof crude pyroligneous acid (PCA). Non-condensable gases weredischarged to open air. The experimental conditions for
pyrolysis of lignin and rice straw are showed at Table 1. The
thermal decomposition of lignin should take place at

3. Results and Discussion
Pyrolysis
The products obtained from the pyrolysis were solid, liquid,and gas ones. Because of the quantity of moisture supplied
during heat treatment, mass of liquid products were greater thanmass of feed material. In Figure 3, the blue and red points show
the yield of solid products and fractional yield of liquid products
relative to material.
In pyrolysis of lignin, there was a significant reduction inyield of solid products at higher temperature. In contrast, it only
had a little change in the fractional yield of liquid products
relative to feed material. Therefore, it was considered that the39

gas products increased proportionally at higher temperature. In
the temperature range from 500℃ to 800℃, the pyrolysis
temperature had significant effects on the yields of solid and gas
products, but it had little effects on the yield of liquid products.

products and sodium hydrogen carbon solution presented thepeak of 5.5 min which did not present at liquid products sample
only. The result of HPLC analysis is showed at Figure 5
Then, the concentration of coumaric acid was quantified.
The fractional yields of coumaric acid relative to feed material
were shown in Table 3.
Table 3 Concentration of sodium salts of coumaric acid
Samples

Concentration coumaric acidin liquid products[mg/ml]

Fractional yield of coumaricacid relative to feed

Lignin 500℃

0.465

0.0010Lignin 600℃

0.488

0.0012

Lignin 700℃

0.430

0.0008Lignin 800℃

0.459

0.0009

Rice straw 600℃

0.464

0.0024

Although rice straw contained lignin with only about 0.2 bymass fraction, the yield of coumaric acid was greater than thatwith lignin. For examples, from these results, it could obtainabout 2.4 g of coumaric acid from 1kg feedstock of rice straw.
The yield of coumaric acid exceeded slightly in liquid products
heat treatment at 600℃ than other temperature. Therefore,
operating heat treatment of lignin at 600 ℃ should be
preferable. However yield of coumaric acid from lignin heat
treatment products was still small.
Extraction of coumaric acid from liquid products.
Figure 6 shows UV spectroscopy of liquid products, ether
phase, aqueous phase

Figure 3. Fractional yield of solid products and liquid products relative
to feed material

Comparing between pyrolysis of lignin and rice straw at
600 ℃ , rice straw was degraded more than lignin. Liquid
products from pyrolysis of rice straw were more than two times
to these from lignin. It might be caused by the degradation ofhemicellulose and cellulose in rice straw. These substanceswere degraded at lower temperature than lignin.
Figure 4 shows the UV spectroscopy of liquid products
from pyrolysis of lignin and rice straw at 600 ℃ . The
wavelengths of the absorptions were observed at 311nm of
coumaric acid, 322nm of ferulic acid and 274 nm of phenol,
respectively.

Figure 6. UV spectroscopy of liquid products, ether phase, water phase

The distribution coefficient, Di, was defined as,
Di =Phenol was considered as the major impurity for therecovery of coumaric acid from liquid products. The Dis of
coumaric acid and phenol were measured as 2.2 and 0.055,
respectively with liquid products of thermal treatment of lignin
at 600℃. Coumaric acid was extracted in the aqueous phase
while most of phenol remained in ether phase. Although thepurity of coumaric acid is still problem, this extraction
experiment method can remove a quantity of impurities inliquid products away coumaric acid.
4. Conclusions
It can be concluded that coumaric acid exits in liquid
products from thermal treatment of lignin. Thermal treatment
lignin at 600℃ is most adequate to obtain more quantity of
coumaric acid. Separation of sodium salt of coumaric acid fromliquid products by HPLC is analytically succeeded. Smallquantity of coumaric acid is extracted from main impurities inliquid products. It is necessary to take account of more efficient
extraction method or thermal treatment experimental conditions
to obtain more quantity and pure of coumaric acid.

Figure 4. UV-spectroscopy of liquid products from pyrolysis of lignin
and rice straw at 600℃

There were broad absorption bands with peaks around 260-280
nm in both the samples. Because the absorption at 322 nm was
relatively small, only little amount of ferulic acid could be
generated. The absorption at around 311nm was observed inliquid products obtained from pyrolysis of lignin. The
absorption of coumaric acid was also at 311nm, as shown
above.

References[1]T.Kaneko, D.Kaneko, S.Wang, High-performance lignin-mimetic
polyesters, Plant Biotechnology 27, 2010, 243-250[2] S. M. Notley and M. Norgren, The Nanoscience and Technology ofRenewable Biomaterials, ed. L. A. Lucia and O. J. Rojas, JohnWiley& Sons, Ltd, Chichester, 2009, ch. 7, pp. 173–205.
[3] Guedira, F., Pyrolise lente de la biomasse: comportement compare
destourteux d’olives, de la bagasse de canne a sucre et la sciure de bois.
Maroc, These (Docteur) – Universite Mohamed,1988, p. 122.

Figure 5. HPLC of liquid products from pyrolysis of lignin and ricestraw at 600℃

It was found that the aqueous solutions of coumaric sodium
salt showed absorption at around 286 nm. Mixture of liquid40

Tsunami Vulnerability Assessment in Tokyo Bay Using a Numerical Simulation
Student Number：09_16970 Name：Shouhei Nagase Supervisor：Hiroshi Takagi
1. Introduction
Tokyo Bay has been considered to be safe against tsunamis
because of its semi-closed shape that would effectively reducethe incoming tsunami energy. However, it is fact that a 3m
tsunami was observed at the fishing port of the Chiba
prefecture and a 1.5m tsunami was also observed in the Sumida
River after the Great East Japan Earthquake on March 11th
2011 [1]. It should be noted that the observed tsunami heightexceeds the estimated maximum level of 1.2m which wassimulated by the Tokyo disaster prevention council in 1991.
Facing unexpectedly large tsunami, it is required toreinvestigate the influence of tsunami on Tokyo Bay with theFig.2 : Comparison between observed and calculated values

new findings and review the disaster mitigation plan inFirst, the simulation was performed for the broad domain
consideration of future possible tsunamis.
with a coarse grid which encompasses the tsunami source
The aim of this research is to investigate the characteristics
generated by the earthquake on March 2011 [2]. In order toof tsunami propagation across Tokyo Bay and to identify theverify the model, the result was checked with several tidal
local amplification mechanism of tsunami heights due to agauges along the coastline in Tokyo Bay. Second, theman-made structure such as floodgates and dykes.
simulation that is limited within the area of Tokyo Bay was
performed in order to achieve the reduction of computational

2. Numerical simulation

time which enables a number of computations repeatedly. For

In this research, the tsunami propagation was simulated using
the second case, tsunamis were generated from the open
Delft3D-FLOW originally developed by Delft University ofboundary which is placed at the baymouth. The validity of theTechnology.
model was also checked in the same way.

3. Local amplification mechanism of tsunami
This research focuses on the resonance phenomenon that canincrease tsunami heights in the bay when the natural frequencyof the bay agrees with the frequency of tsunami. It is also
hypothesized that the tsunami height would be amplified in
front of a floodgate closed due to the reflection of tsunami andmay result in increasing the flood risk round the area. To
investigate these risks, several types of test tsunamis withdifferent period of sinusoidal waves were generated from theopen boundary.
Fig.1 : Observation point

41

tsunami in front of the gate closed. However, it appears to be

Table.1: Calculation Settings
Item
Calculating area
Time step
Tide level
condition

Ground levels

3-6hours(It changes according to a cycle)
Largest mesh: the Maritime Safety Agency data
Second mesh: ditto
Third mesh: uniformly T.P. -3.2m
It creates from 10m mesh data of
Geographical Survey InstituteRelative roughness

0.025 is uniformly set to a water area

Computation time
Bathymetric data

excessive because the reflection can only explain the

OutlineLargest mesh(200m square mesh),
Second mesh(40m square mesh)
Third mesh(20m square mesh)Δt=0.1sec
Mean sea-level of Tokyo Bay(T.P.+0.0m)

amplification that should be less than 2 times. One of thereasons for this excessive amplification is that the tsunami
which breaks in the shallow water turns into a bore whichcauses mass transportation towards upstream, and eventually is
accumulated in the front of the gate.

3. Conclusions
The numerical simulation is indispensable in order toestimate a future possible tsunami. In this research, thefollowing conclusions are obtained:・The tsunami simulation for the great tsunami on March2011 was performed and validated its accuracy with theobserved water elevations.

・The period of resonance in Tokyo Bay appears to bearound 50 to 90 minutes, and the resonance causes
significant increase in tsunami heights.

Fig.3: Maximum water level・The amplification of water level in front of a floodgate
closed can be reached more than two times with that of agate opened.

It is necessary to advance the research about the following
subjects.
・

Development of the tsunami risk evaluation methods
under more complex conditions in consideration of thedamaged dykes due to earthquakes and liquefactions.

・

Numerical simulation which can incorporate the gate
operation in the process of the calculation.

Fig4:Difference in maximum water level at the Kamihiraisluice between two different gate operational patterns: openig

References[1] Jun Sasaki et al., Behavior of The 2011
Tohoku Earthquake Tsunami and Resultant Damage in Tokyo
Bay, Coastal Engineering
Jounal, Vol.54, No.1, p.26, Mar.2012.
[2] Tatsuhiko Saito et al., Tsunami source of the 2011
Tohoku-Oki earthquake, Japan:Inversion analysis based on
dispaersive tsunami simulations, GEOPHYSICAL RESEARCHLETTERS, VOL.38, Isuue 7(2011)[3] 日野幹雄, 日野和子, 東京湾の進入長波に対する振動応答特性, 海岸工学論文集, 第 11 回, 1964.

and closing
It is shown that the water level deviation tend to increase withan increase in period in Fig.3. Three out of five points show
significant increases in water level around 50 to 90 minutes thatclearly demonstrate the resonance. This result agrees with theanalysis performed by Hino et al. [3] in which the range ofnatural oscillation period of Tokyo Bay is shown in between 60
and 90 minutes. Fig.4 indicates that the water levels in front of
the floodgates are much larger when closing than when opening.
This amplification can be explained by the reflection of42

CARBONATION PHENOMENA AND CORROSION BEHAVIOR OF STAINLESS STEEL BAR
IN THE MORTAR MIXED WITH SEAWATER
Student Number: 09-17052 Name: Tomohiro Nagata Supervisor: Nobuaki Otsuki

1 Introduction

Table 1 Chemical properties of bars (mass%)

Fe
C
Si
Mn
P
S
Ni
Cr
Worldwide shortage of freshwater is expected in the future,and perspective for the availability of concrete mixed with SD295A 98.33 0.2 0.24 0.74 0.22 0.27
SUS410SD
86.7
0.01
0.28
0.34
0.2
0.02
0.18
12.57
seawater is required.
The most dangerous factor for concrete mixed with SUS304SD 71.4 0.06 0.38 1.48 0.32 0.21 8.14 18.32
seawater is corrosion of reinforcing bars. Especially in land
Table 2 Mix proportions of mortar
environment reinforcing bars are under the danger of
Unit weight(Kg/m3)
W/carbonation phenomena.
Type of Cement
GypCW
OPC BFS
Sand
Carbonation decreases the pH of concrete, moreover it issumreported that internal chloride ion from seawater is
OPC
316
632
0
0
1264
50
concentrated to the inner by carbonation phenomena [1].
%
BB
315
378
247
5
1259
Therefore, reinforcing bars in the carbonated concrete mixed
with seawater are exposed to a severe environment with low
pH and large amount of chloride ion.
To overcome this problem, the usage of stainless steel asreinforcing bar is commonly practiced. However the
corrosion behavior of stainless steel reinforcing bars incarbonated concrete mixed with seawater is not investigatedenough.
Fig.1 Schematic of specimen
Thus, the objectives of this study are set below,
density in saturated surface-dry condition:2.61g/cm3 water(1)To evaluate the effect of seawater on carbonation
absorption:1.88) was used.
phenomena from carbonation rate of mortar mixed with
2.2 Specimen
seawater.
The size of mortar specimens was 40×40×160mm and the(2)To clarify chloride ion concentration in the carbonatedtypes of covering depth of specimens were 5mm and 15mm, asmortar mixed with seawater from micro analysis
shown in Fig.1. The specimens were cured for the span defined(3)To clarify corrosion behavior of stainless steel reinforcing
by concrete standard specifications (OPC:5days-curing 、bars in mortar mixed with seawater under carbonation from
BB:7days-curing) , and for only 1day (1day-curing). Table 2
electro chemical measurement.
shows the mix properties of mortar. After curing, sides of(4)To suggest the proper combination for concrete mixed
specimens are sealed with epoxy resin (Fig.1). After that they arewith seawater under carbonation.
exposed to carbonation chamber (the density of CO2:5%

2 Experiment procedures

humidity: 60%).

2.1 Materials(1)Carbonation rate was evaluated by measuring periodicallycarbonation depth of specimen, using phenolphthalein method(JIS-A-1152). Based on that, the carbonation rate coefficient was

2.3 Measurement items and method

The materials used in this study were two types of cement,
three types of bar, and two types of mixing water. Two types of
cement were Portland blast -furnace cement class B (BB) and
Ordinary Portland Cement (OPC) Three types of bar were two
types of stainless steel bars (SUS304SD and SUS410SD) andplain bar (SD295A). SUS410SD has smaller content of
Chromium (Cr) ratio than SUS304SD. Table 1 shows thechemical properties of bars. Two types of mixing water,
freshwater and artificial seawater (ASTM-D-1141-98(2003))
were used. As fine aggregate land sand (fineness modulus:2.52

derived from √𝑡 method.
(2)To clarify concentration of chloride ion by carbonation
phenomena, micro analyses were conducted on the section ofmortar specimen by EPMA method.
(3)To clarify corrosion behavior of stainless steel bars, thecorrosion current density of specimen were measured
periodically. Moreover the relationship between corrosion
current density and carbonation rate was evaluated.
43

Fig.2 Carbonation rate coefficient (natural environment)

Fig.3 The distribution map of chloride ion of BB (1day-curing)Fig.4 The relationship between remaining depth of carbonation and corrosion current density
cement. When carbonation reaches the covering depth, current
density of SD295A in mortar mixed with seawater increased
and reached 0.2µA/cm2.

3 Experimental results
3.1Carbonation phenomena
Fig.2 shows the carbonation rate coefficient of each specimen
and data obtained by Otsuki laboratory [2]. Carbonation rate
coefficients were converted into natural environment (thedensity of CO2:0.0315%) by Uomoto formula [3]From this figure and within the scope of this study carbonation
rate of mortar mixed with seawater was almost same as that
mixed with freshwater.
Besides, Fig.3 shows the chloride ion distribution map of BB(1day-curing) mixed with seawater at the age of 16weeks. From
this figure it is confirmed that chloride ion is concentrated into
the inner part due to carbonation phenomena.

4 Conclusions(1)Carbonation coefficient of specimen mixed with seawater
was almost same as that mixed with freshwater.
(2)Chloride ion is concentrated to the inner by carbonation
phenomena in the case of mortar mixed with seawater.
(3)When mortar mixed with seawater is carbonated, stainlesssteel bars isn’t corroded except for SUS410SD embedded in
Portland blast -furnace cement class B.
(4)The combination of SUS304SD and any cement, and the
combination of SUS410SD and Ordinary Portland Cement may
be proper for concrete mixed with seawater under carbonation.

3.2 Corrosion behavior of stainless steel bars
Fig.4 shows the relationship between remaining depth ofcarbonation and corrosion current density of each specimen.
Remaining depth of carbonation is the difference between
carbonation depth and covering depth.
From this figure, in the case of stainless steel bars, corrosion
current density of SUS304SD was maintained at low level
regardless of type of cement. However, in the case of
SUS410SD, before the carbonation reached the covering depth
the corrosion current density was maintained at low level, but
after the carbonation reached the covering depth, corrosion
current density of BB mixed with seawater reached
0.2µA/cm2(corrosion criterion of FIB).
On the other hand, corrosion current density of SD295A
increased as carbonation progresses, regardless of type ofReferences[1] Koichi Kishitani etc, The Relationship between Rebar
Corrosion and Neutralization in the Concrete Containing
Chloride, Proceedings of annual conference of JCI,Vol.13、No.1.
pp.601-2, 1991
[2] Takahiro Nishida etc, Resistance to transfer of materials in
the case of concrete mixed with seawater, Proceedings of annual
conference of JCI, 2013 (under reviewing)[3] Kento Uomoto etc, Influence of Concentration Carbon
Dioxide on Carbonation Ratio of Concrete, production study,
Vol.43, pp.35-38, 1991

44

Image Compression by using Vector Quantization and
Vector-Embedded Karhunen-Loève Transform
Students Number: 09 18761 Name: Kiung Park Supervisor: Yukihiko Yamasita

I

Introduction

III Representative Vectors and Corresponding Covariance Matrixes

JPEG is an international standard for still image compression. Even if it has problems such as blocking distortion and image quality degradation in edge region, it has
been used widely because its algorithm is simple and fast,
and the compression ratio is high.
Quantization used in JPEG is scalar quantization (SQ).
Shannon’s rate-distortion theory indicates that a better
performance can be achieved by coding vectors insteadof scalars [1], [2]. Thus, if we introduce vector quantization (VQ) into JPEG, we can expect better compression
performance. However, current image coding methods by
VQ have a problem that the size of representative vectors
is small. It degrades coding efficiency.
This paper proposes a new image copmression schemeby coding the difference between input vector and its
representative vetor to overcome the problem. VEKLT(Vector-Embedded Karhunen-Loève Transform) [3] isused as the transformation instead of DCT to introduce
statistically adaptive transformation.
The rest of this paper is organized as follows. In
Section 2, we will give an overview of the proposed
image compression scheme. Some key techniques of the
proposed scheme are presented in Sections 3 and 4. The
experimental results demonstrate the advantages of the
proposed scheme in Setion 5. Finally, Setion 6 concludes
this paper.

II

A. Representative Vectors
At first, normalized input vectors are obtained by using
the same method in Fig. 1 from prepared 160 training images (512×512).GMM (Gaussian Mixture Model) which
is a clustering method is implemented. The initial mean
vectors are decided by k-means. Because the mean vector
has the highest probability in each cluster, it appears very
often in images. Thus, I sample it from training images.
In my research, I sample 400 mean vectors from every
training image. Then, 64000 candidates are obtained for
representative dataset. Finally, I choose 100 vectors fromthe candidates by k-means and normalize them again. The
last 100 normalized vectors are used as the representative
vectors for VQ.
B. Covariance Matrixes for VEKLT
Because VEKLT is a kind of KLT, a covariance matrix
is needed for making a transformation matrix. In my
research, the covariance matrix is made by the difference
vectors that have a same index. Then, if the numberof representative vectors is 100, I have to prepare 100
covariance matrixes. Therefore, I made 100 covariancematrixes from 160 training images used for making VQ
dataset.

Overview of Proposed Method

IV Implementation of VEKLTThe proposed image compression scheme is illustrated
in Fig. 1. The proposed encoder is based on JPEG and thedecoder is nearly the same as the inverse of encoder.
The proposed process is the same as that of JPEG until”Subtract 128 from each pixel”. At first, DC coefficientsof input blocks are obtained by inner product with thefirst DCT basis. And then, the block is normalized andthe index of the nearest representative vector is selectedfrom dataset by using recognizer based on MSE. Afterobtaining difference data from the representative vector,
VEKLT is performed by using the covariance matrix. Next,
they are coded by arithmetic coder. Actually, Huffman
coder is used in JPEG, but the code book isn’t optimized
in proposed scheme. Thus, I use arithmetic coder into allof transformed data for fair comparison. After comparing
the code-length of AC coefficients from DCT and the sumof code-lengths of AC coefficients from VEKLT and VQ
coefficient and index, the shorter one is selected. If DCT
is selected, its index is changed to 0 and its VQ coefficientis deleted. Finally the coder outputs data structure which
consists of DC coefficients and revised indexes and VQ
coefficients and seleted AC coefficients.

Let fi be the i-th block of image f, it can be written as
fi =< fi , ΦDC > ΦDC + < fi , Φi,AC > Φi,AC .
By normalizing fi , Φi,AC is obtained. And then, thenearest vector Φ̂i,V Q to Φi,AC is selected from VQ dataset
by using recognizer. In case of 8×8 block size, thedimension of fi is 64. Because ΦDC and Φ̂i,V Q are
perpendicular, fi can be expanded asfi=+< fi , ΦDC > ΦDC + < fi , Φ̂i,V Q > Φ̂i,V Q
62∑< fi , Φij > Φij .
j=1

where Φij is perpendicular to ΦDC and Φ̂i,V Q .
The remaining basis Φij is obtained by VEKLT. Let
Ri be the covariance matrix of fi and W = ΦDC ΦT
DC +Φ̂i,V Q Φ̂Ti,V Q , the symmetric matrix Qi is defined by
Qi = Ri − Ri W − WRi + WRi W.
Φij (j = 1, . . . , 62) is given by eigenvectors of Qi corresponding nonzero eigenvalues. The transformation matrix

45

Fig. 1.

The proposed encoder

Ti is given by
Ti = [Φi1 , . . . , Φi62 ]T
After transforming fi as Ti fi , the coefficients <fi , Φij > are obtained.

V

Experimental Results

Fig. 3.

Cropped Lena reconstructed by JPEG (bpp : 0.24)

Fig. 4.

Cropped Lena reconstructed by proposed scheme (bpp : 0.24)

A. Rate-Distortion Relation

Fig. 2.

Rate-distortion curve (Lena)

I used a standard image Lena (512×512) to demonstrate
the performance of the proposed scheme. Fig. 2 andTABLE. I show the rate-distortion relation. In all of bpp(bits per pixel) range, the proposed scheme outperforms
JPEG 0.2∼1.0dB in PSNR (Peak Signal-to-Noise Ratio).
Especially, when bpp is low, it outperforms JPEG clearly.

VI Conclusion
This paper proposes a new image compression schemeby using VQ and VEKLT. Experimental results fully
demonstrate the advantages of the proposed techniques. Its
PSNR outperforms JPEG 0.2∼1.0dB and its visual quality
in edge region is better than JPEG especially at low bpp.
However, there are still many aspects to be investigatedin the future, such as indexing method, approximation ofreal covariance matrix, number of representative vectors,
their mining method, coding for transformed blocks by
VEKLT, and so on.

TABLE I

bpp
0.24
0.50
1.00

Lena
PSNR in dB
Proposed Scheme JPEG
31.69
30.67
35.47
34.73
39.27
38.90

References[1] C. E. Shannon A mathematical theory of communication, Bell
System Tech. Journal, vol. 27, pp. 379-423, 623-656, 1948.
[2] Feng Wu, Xiaoyan Sun Image Compression by Visual Pattern Vector
Quantization(VPVQ), IEEE Data Compression Conference, pp. 123131, 2008.
[3] Toshihisa TANAKA, Yukihiko YAMASHITA Vector-Embedded KLT
and Its Application in Orientation Adaptive Coding of Images, IEICE
TRANS. FUNDAMENTALS, vol. E38-A, No.6 2000

B. Visual Quality
Images in Fig. 3 and Fig. 4 are reconstructed respectively by JPEG and the proposed scheme. We can see thevisual quality of the proposed scheme is better than JPEG.

46

Development of Millimeter-wave Radio Channel
Sounder
Student Number: 09B20284

1

Name: Pham Kinh Hung

Introduction−20

Millimeter-wave with 1 mm to 10 mm wavelengthand 30 GHz to 300 GHz frequency range is an alternative frequency bands to be applied in high speed
wireless communication systems. In this study, amillimeter-wave channel sounder has been developed
to investigate the propagation characteristics at 60
GHz bands. Because of using low cost off-the-shelf
RF transceiver, the measurement accuracy is affectedby many error factors, mainly phase noise and IQ Imbalance. This study discussed about the design of thesounding signal and compensation requirement for an
accurate measurement by investigating the effect ofphase noise and I/Q imbalance in the numerical simulation.

2

Supervisor: Jun-ichi Takada, Minseok Kim

11 GHz60 GHzPSD (dBc/Hz)

−40
−60
−80
−100
−120
−140
4

6

10

10
Frequency Offset (Hz)Figure 3: Spectrum of Phase Noise

Developed hardware

the way to calibrate or reduce the influence.

Millimeter-wave channel sounder consists of transmitter, receiver and local oscillator which supply
signal to synthesize transmit signal and carrier signal. Each transmitter consists of baseband (BB) circuit and radio frequency (RF) circuit. Fig.1 is themillimeter-wave channel sounder used in the experiment. Channel sounder is used to estimate the transfer function H(f ) between the transmitter and thereceiver. The relation of H(f ) with Fourier transforms of transmit signal X(f ), receive signal Y (f )and white noise N (f ) is represented as:
Y (f ) = H(f )X(f ) + N (f ),3Phase Noise Influence and Channel
Sounding Method

Phase noise is the frequency fluctuations of the local oscillators at the transceiver. In the frequency
domain, a pure sine wave would be represented asa delta function. However, all oscillators, that generate carrier signal in RF circuit, have phase modulated noise components. The phase noise components
spread the power of signal to adjacent frequencies,
resulting in noise sidebands, thus causing interference between sub-frequencies. To reduce influence ofphase noise to measurement result, we have to choose
appropriate sounding method.
In previous research, we have already developed
channel sounder at 11 GHz. In 11 GHz system, thephase noise power is small, so it almost does not affect to the system [1]. The transmit signal is multitone which consists of 2048 single-tones arranged on(1)

In practice, as shown in Fig.2, multiplicative error factors such as the phase noise (φt (t), φr (t)) andI/Q imbalance (ε, θ) degrade the measurement accuracy. Since we can measure X(f ), Y (f ) and N (f ),estimating H(f ) accurately requires us to define theinfluence of phase noise and IQ imbalance, then findRF

RF

Baseband (BB)

BB

X(f)

Local Oscillator (LO)Radio Frequency
(RF)

8

10

H(f)

I
0
90 ° + θ

1+ εQTransmitter

I/Q imbalance

I/Q imbalance(ε ,θ)

I

N(f)

1+ ε

(ε , θ)sin ( 2πf ct + φr (t ))sin ( 2πf ct + φt (t ))Phase noise: φ (t ) − φ (t )rFigure 1: Millimeter-wave Channel
Soundert

Figure 2: Channel Sounder Model47

090 ° + θ

Q

BB

Y(f)
ReceiverTransfer function

Impulse Response

5

1

1

0
After IQI Compensation
After Back−to−back

0
−20

−5

4

2

4

5

3

5

Level (dB)−102Level (dB)

3

−15
−20

Figure 4: IQ Imbalance Compensation Result

As can be seen Fig.2, IQ imbalance is caused by aslight discrepancy between the in-phase and quadrature components in the quadrature modulator. In
this IQ imbalance model, ε, and θ denote the amplitude difference and the phase offset of the quadratureoscillator branches, respectively. Meanwhile, E represents the error matrix of ε and θ. We also denote
w(t) and x(t) as the carrier signal and baseband signal, respectively. The distorted output signal can becalculated as:

5−T

x (t) = (Ex(t)),

(2)

−100

0
100
Frequency (MHz)

200

−100
0

0.2

0.4
Delay (µs)

0.6Conclusion

References(3)[1] M. Kim, J. Takada, Y. Konishi, ”Novel Scalable
MIMO Channel Sounding Technique and Measurement
Accuracy Evaluation with Transceiver Impairments,”
IEEE Trans. Instrum. Meas., Vol.61, No. 12, pp. 31853197, Dec. 2012[2] M.Ghoraishi, M.Kim, J.Takada, ”Influence of phasenoise on the frequency division multiplexing channelsounding,” 2010 Asia Pacific Microwave Conference,
Yokohama, Japan, Dec. 2010[3]M. Kim, Y. Konishi, J. Takada, B.Gao, ” I/Q Imbalance Compensation Technique for Quadrature Modulation by Single-Tone Testing,” submitted to IEICE
Trans.Commun.

the output signal can be generated ideally. Thatmeans when we know the error parameters in model,
we can easily eliminate IQ imbalance using digitalpre-distortion on the original input signal. By using
the spectrum analyzer, we can measure the power
spectrum of image distortion which can be use to
estimate the error parameters for taking a compensation. The detailed algorithm can be seen in [3].
In this compensation method, the model parameters are optimized to suppress the values of LIRR ,that denotes the image power reduction ratio (IRR).
LIRR =−40−200

This research proposed an appropriate channelsounding method for millimeter-wave channelsounder and a calibration method to increase accuracy of measurement result.
As future works, we plan to design 4 × 4 MIMO
millimeter-wave channel sounder to measure propagation characteristic of the millimeter-wave.

T

where z(t), zei (t) are the ideal output and the image
distortion components respectively. In the equation(2), if the baseband signal is distorted as0

After IQI Compensation
After Back−to−backbe derived that NMSE corresponds to LIRR .
Fig.4 shows the result of IQ imbalance compensation in millimeter-wave channel sounder. In thisfigure, from (1) to (5) are fundamental component,
carrier leakage, image distortion, image of third harmonic and second harmonic component. By comparing left side and right side, we can see the image
distortion has already reduced, and LIRR becomes
greater than 60 dB, thus the NMSE is less than 10−6 .
Since millimeter-wave radio frequency circuit does
not allow direct current to pass through, carrier leakage can not be compensated. Similarly, the second
harmonic component and image of third harmoniccomponent also can not be compensated, but they
are small enough compared to the fundamental component.
After compensating for both transmitter and receiver, we connect them directly to estimate the instrument characteristic, then eliminate it to obtainthe unknown channel characteristic. The result isshown in the Fig.5.

IQ Compensation= (Ew(t)) x(t)= z(t) + zei (t),

−80

−35

Figure 5: The Transfer Function and Impulse Response
after Calibration

bandwidth of 400 MHz. However, in millimeter-wave
sounder, the phase noise power becomes stronger. It
is needed to increase the interval among tones (∆f ),so that the number of tones on the whole bandwidth
will decrease because of the fixed bandwidth. On theother hand, the measurement time delay calculated
as the time signal transmit from transmitter to receiver, should be smaller than the length of symbol(1/∆f ) to avoid interference among other channels.
From those conditions, the number of tones is decided
as 256 tones. Fig.3 shows the comparison of spectrumof phase noise between 11 GHz and millimeter-wave
channel sounder (at 60 GHz).

z̃(t)

−60

−25
−30

4

−40

(ε2 + ε + 2) − 2(1 + ε) cos θE[zei (t)]=, (4)E[z(t)]
4

NMSE (Normalized Mean Square Value) is used toevaluate the effect of IQ imbalance From (4), it can

48

Identification of riboflavin producing yeast in the fermentation process for Puer tea
Student Number: 09-21088 Name: Ryutaro BESSHO

1 Introduction

Supervisor: Kiyohiko NAKASAKI

2-5 Real-time PCR
The cell density of strain YPT2 in the Puer tea
samples was measured by real-time PCR method using
the DNA extracted from the Puer tea samples as atemplate. The primers specific to strain YPT2, Y2RTf,
5′-CGAGAGGTGTGAAAGGTACAC-3’, and Y2RTr,
5′-GATGGCTGTTGCTATATCCCC-3’ were designed
and used. The real-time PCR was performed using
Smart Cycler® II system (Cepheid, California, USA) andSYBR® premix EX TaqTM (TaKaRa Biomedicals, Shiga,
Japan).

Puer tea is a unique fermented tea manufactured in
the Yunnan province of China. It has been consumed by
Chinese people for centuries, and has recently become
popular in Japan as well, as a functional beverage [1].
Pure tea is produced from green tea as raw material by
the act of microorganisms whereas black tea or oolong
tea is produced by the act of enzymes contained in thetea leaves themselves. It is known that in general, Puer
tea contain higher concentration of riboflavin, a kind ofvitamin B2 compared to the green tea, which indicatesthat riboflavin in the leaves was enriched by the act of
microorganisms during the fermentation. The objectivesof this study were the identification and the isolation of a
microorganism that has the ability to produce riboflavinusing the tea leaves as substrate.

2-6 Pure cultivation of the strain YPT2 on the tea
leavesChinese green tea leaves sterilized by 1 Mrad/h of
Co-60 gamma ray irradiation for 3 h were placed in a
500 mL Erlenmeyer flask and the water content of theleaves was adjusted to around 30 % by adding sterilized
distilled water. Strain YPT2 was cultivated for 2 d inpotato dextrose medium beforehand and was inoculated
to the leaves so that the initial cell density of strain YPT2
was about 6 log CFU/g-ds. Cultivation was performed at
30 °C for 5 days with the air flow rate of 3 mL/min. The
samples were collected daily and were subjected to themeasurements of water content, pH, cell density of thestrain YPT2 and riboflavin concentration. The celldensity of strain YPT2 was measured by dilution plating
method.

2 Materials and methods
2-1 Puer tea samples and process characterization
Tea leaf samples under the fermentation process
were obtained from a company at Fuken in China. At thefactory, harvested tea leaves were heat-treated toinactivate the indigenous enzymes in the tea leaves. Then,the leaves were softened by mixing, and dried undersunlight. To promote the microbial activity at the start offermentation, water was added to the leaves. Then theleaves were piled in a windrow shape in the fermentation
room. During the fermentation, leaves were turned overroughly once a week to ensure the homogeneity of thetea leaves.
Temperature of the leaves was measured at thecenter of the pile. Tea leaf samples were collected fromthe piles on days 0, 1, 14, 21, 35, and 50 of thefermentation process and were subjected to themeasurements of pH, water content and riboflavinconcentration, and to the microbial analyses.

3 Results and discussion
3-1 Characterization of the fermentation process
The courses of temperature, pH and water content
are shown in Fig. 1. Gradual increase in temperature by
the heat of microbial activity and instant decrease intemperature due to the turning over of the pile took placealternately during the fermentation. The pH value of thestarting material was 5.5 and fell slightly to 4.3 by day
14. After that the pH value gradually increased andfinally attained at around 5.8. The water content of the
raw material was as low as 10% at the beginning of theTemp. (°C)

2-2 Determination of riboflavin concentration in the
Puer tea samples
A 1-gram of leaf sample was homogenized in 9
mL of distilled water. Extraction conditions were thesame as described by Capozzi et al [2].

Water
content (%)

pH (-)

2-3 Microbial community analysis of Puer tea
samples
The microbial community in the leaf samples wasdetermined by denaturing gradient gel electrophoresis(DGGE) analysis [3]. In the DGGE analysis, fungal 26S
rRNA gene fragment was used. Characteristic bands inthe DGGE gel were excised and base sequences of the
DNA contained in the bands were determined.
2-4 Isolation of a yeast strain
A 3-gram of leaf sample was homogenized in
27mL of distilled water and was diluted appropriately.
The diluted solutions were spread onto rose-bengal agar
medium and were cultivated at 30°C for 2 d. After that,the colonies appeared on the plates were picked up andstreaked several times in order to purify the isolated
yeast.

80
60
40
20
0
5.5
5
4.5
4
3.5
50
40
30
20
10
0

0

10

20
30
40
Fermentation period (d)50Fig. 1 The course of temperature, pH and water content
in Puer tea pile at the tea processing factory in
China

49

process and increased to around 40% on day 1 by theaddition of water, then, it maintained around 40% during
the fermentation.

fermentation. In the attempts to isolate the fungi
corresponding to the distinct bands appeared in the
DGGE gel, strain YPT2 was isolated. When the strain
YPT2 was subjected to PCR-DGGE analysis, the band
was observed at positions identical to that of the band c(Fig. 3).

3-2 Concentration of riboflavin in Puer tea samples
The course of riboflavin concentration in the Puer
tea samples is shown in Fig. 2. The initial riboflavinconcentration was 2 µg/g-ds, and it increased to 10µg/g-ds until the end of fermentation. Since enzymes
contained in the leaves were inactivated byheat-treatment, it is considered that the riboflavin was
produced by the microorganisms in the leaves.
14
12
10
8
6
4
2
0
0
10
20
30
40
50
Fermentation period (d)Log cell density(CFU/g-ds)

Conc. of riboflavin(µg/g-ds)

3-4 Concentration of strain YPT2 in Puer tea samples
The course of cell density of strain YPT2 in the
Puer tea samples determined by real-time PCR method isshown in Fig. 4. The initial concentration of strain YPT2
was 6.4 log CFU/g-ds. As early as day 1 of fermentation,
8.1 log CFU/g-ds of strain YPT2 appeared. Then, the celldensity decreased gradually until the end of fermentationand finally, it reached 5.6 log CFU/g-ds. The cell density
of YPT2 kept high at around 8.0 log CFU/g-ds from day
1 to 21 and during the same period of fermentation, theincrease in the concentration of riboflavin was observed,
suggesting that the strain YPT2 have contributed to theriboflavin production during the Puer tea fermentation.
9

Fig. 2 The course of riboflavin concentration in Puer tea
samples
3-3 PCR-DGGE analysis
DGGE image obtained from amplification of the
26S rRNA gene fragment derived from the tea leaves isshown in Fig. 3. There were five distinct bands, band a, b,
c, d and e which appeared consistently for certain periods
during the fermentation. Among the 5 distinct bands,
band a was not reproducible when the tea leaf samples
obtained from another fermentation piles were subjectedto the PCR-DGGE analysis. Therefore, the other 4 bands
except the band a were excised from the DGGE gel andDNA sequences were determined. The nucleotide blast ofthese sequences indicated the closely related species forthe fungi corresponding to those bands as follows, band
b: Blastobotrys adeninivorans, band c: Debaryomyces
hansenii var. fabryi, band d: Aspergillus sp., band e:
Talaromyces sp. Among those 4 kinds of fungal species,
Debaryomyces hansenii corresponding to the band c isfamous for the high riboflavin producing ability and it
was deduced that the fungal strain corresponding to theband c produced riboflavin during the Puer tea

8
7
6
5
4
0

10

20
30
40
Fermentation period (d)50Fig. 4 The course of cell density of strain YPT2 in Puer
tea samples
3-5 Pure cultivation of strain YPT2 on Puer tea leaves
Riboflavin producing ability of strain YPT2 wasverified by pure cultivation on green tea leaves that were
the raw materials of Puer tea. As was not shown here, theincrease in the riboflavin concentration was observed in
the raw materials along with the growth of strain YPT2.

4 Conclusions
It was confirmed that riboflavin concentration in
Puer tea leaves increased during the Puer tea
fermentation and strain YPT2, the genus of which was
famous for riboflavin production was isolated andidentified. It was also verified that riboflavinconcentration in Puer tea leaves increased when celldensity of strain YPT2 was high and that strain YPT2
had the ability to produce riboflavin in the cultivation
using the green tea leaves as substrate.

References[1] Sano, M., et al.: Effects of Pu-erh tea on lipid
metabolism in rats, Chem. Pharm. Bull., 34,
221-228, (1986)[2] Capozzi, V., et al.: Biotechnological Production ofVitamin B2-enriched Bread and Pasta, J. Agric.
Food Chem., 59, 8013-8020, (2011)[3] Abe, M., et al.: Characteristic fungi observed in thefermentation process for Puer tea, Int. J. Food
Microbiol., 124, 199-203, (2008)
Fig. 3 PCR-DGGE fingerprint for the Puer tea samples

50

大域的射影変換を 用いたパタ ーン マッ チン グに関する 研究
学籍番号： 09_23408

1

学生氏名： 苗春亭

はじ めに

指導教員： 山下幸彦

式 (3) の最大化を 直接計算する こ と は難し いため，
若原ら は次に示す解法を 提案し た [1,2,3]。 ま ず， 2 次
元ガウ ス 関数 G(x)光学的文字読取装置 (OCR) は， スキャ ナなどによ って 入力し た 画像の文字を コ ン ピ ュ ータ に よ っ て 認識する 装置であ る 。 近年では， そ の認識は統計的処理に よ っ て 行う も のが多い。 し かし な がら ， そ のた めに は学習の際に 多数の文字データ が必要と な る 。 そのよ う な 多数の文字データ が得ら れな い場合に 対して ， パタ ーン マ ッ チン グに よ っ て 文字を 認識する 手法が提案さ れて いる 。 一般に認識対象と なる 文字 (入力文字) は， 学習し た 文字 (学習文字) と 全く 同じ ではな く ， 各種の変形を 受け て いる 。 アフ ィ ン 変換は，
2 次元座標の線形変換と 平行移動から な る 変換である 。 画像全体に対する アフ ィ ン 変換に対し て 頑健な ，
入力文字と 学習文字の相関値の算出法であ る ， 大局
的アフ ィ ン 変換 (GAT) 相関法 [1,2,3] が提案さ れ， パ
タ ーン マ ッ チン グに よ る 文字認識手法と し て は， 高
い認識率を 実現し て いる 。
アフ ィ ン 変換を 拡張し たも のに， 射影変換がある 。
射影変換は， 3 次元空間の平面を 別の平面に 透視して 得ら れる 変換である 。 本論文では， GAT 相関法において 頑健と な る ア フ ィ ン 変換を 射影変換に 拡張した， 大域的射影変換 (GPT) 相関法を 提案する 。 そ して ， GPT 相関法の具体的ア ルゴリ ズム を 示し ， 計算機実験に よ っ て そ の有効性を 確認する 。

G(x) =

Dで近似する 。 こ の式は， ガウ ス 関数でぼかし た 画像の相関を 意味し て いる 。 こ の近似に よ り ， 未知成分
A, b が画像の中ではな く ， ガウ ス関数の中に含ま れる よ う に な り ， 微分な ど が容易に 計算でき る よ う にな っ た。 さ ら に ， エッ ジ方向の同一性を 導入し ，
Z Z
JGAT =
G(Ax1 + b − x2 )DDδ(∇f (x1 ), ∇g(x2 ))f (x)g(x)dx1 dx2を 定義し ， こ の最大値を GAT 相関値と し て いる 。 こ
こ で， δ(∇f (x1 ), ∇g(x2 )) は， 2 つの画像のエッ ジ方向が等し いと き に 1， そ れ以外のと き は 0 と する 関
数であ る 。 x1 , x2 の関数 (行列ま た はベク ト ルの場
合も ある )h(x1 , x2 ) に対し て ， 式 (6) の被積分関数によ る 期待値をZ Z
h(x1 , x2 ) ≡h(x1 , x2 )G(Ax1 + b − x2 )で定義する 。
GAT 相関法では， A と b に対する JGAT の微分を0 と し て 最大値を 求める 。 こ のと き ， 得ら れた 式のガウ ス 関数 G(x) の中では， A = O, b = 0 と いう
第 0 近似を 使っ て ， A と b に 対する 線形方程式を 求めて ， 次式に よ り A と b を 求めて いる 。

(1)が成立する 。 こ こ で， A と b はア フ ィ ン 変換を 表す(2, 2)-行列と 2 次元ベク ト ルである 。 定義域が D である 画像 f (x) を ， 式 (1) のアフ ィ ン 変換で変換し た
画像を fA (x) と おけ ば，

1f (A−1 (x − b))|A|Dδ(∇f (x1 ), ∇g(x2 ))f (x)g(x)dx1 dx2

2 次元座標を 表す 2 次元ベク ト ルを x と おき ， アフ ィ ン 変換し たベク ト ルを x′ と おけ ば，fA (x) =(5)

D

GAT 相関法

x′ = Ax + b(4)を 定義する 。 そ し て ， 式 (3) を ，
Z Z
G(Ax1 + b − x2 )f (x)g(x)dx1 dx2

D2√ −γkxk2
γe(2)が成立する 。 2 つの画像 f (x) と g(x) の GAT 相関値
は， 様々 な ア フ ィ ン 変換に 対する
Z
fA (x)g(x)dx
(3)

A

=

b

=



x2 x1 T
1
Ax1
x2
+
−
1
1
x2 xT1 −

x1 xT1 −x1 x1 T1

−1
(6)

ガウ ス関数に近似を 適用し たため， 一般には式 (6) の解は真の最大点を 与え な い。 従っ て ， そ のア フ ィ ン
変換で画像を 変換し ， 式 (6) を 再計算する 。 こ の操
作を 繰り 返し ， 相関値が収束し た 時に ， そ の相関値
を GAT 相関値と し て 出力する 。

Dの最大値と し て 定義さ れる 。

1
51

3

GPT 相関法

表 1: 相関法に よ る 相関値

2 次元ベク ト ルを x に 対し て ， 射影変換し た ベク
ト ルを x′ と おけ ば， 次式が成立する 。
x′ = (Ax + b)/(1 + hc, xi)

AX
BX
CX
DX
AY
BY
CY
DY
AZ
BZ
CZ
DZ

(7)

上式で、 c = 0 と すれば， 式 (7) は， アフ ィ ン 変換を表す。 こ こ で， 部分射影変換 (PPT) を ，

x′ = x/(1 + hc, xi)(8)で定義する 。 射影変換はア フ ィ ン 変換と 部分射影変換の組み合わせで表すこ と ができ る 。 部分射影変換を 求める た めの評価基準を ，


Z Z
x1
JP P T =G− x2
1 + hc, x1 i
D D
δ(∇f (x1 ), ∇g(x2 ))f (x)g(x)dx1 dx2

Org
0.8144
0.79380.73700.6746
0.52520.48450.52110.4734
0.66330.60230.65730.5528GAT
0.96310.8796
0.94530.86180.94560.88450.92340.76200.96000.87580.94360.8540GPT
0.99010.97580.96250.8861
0.95420.94690.93980.8856
0.99650.99700.99090.9111

GPT - GAT
0.02700.09620.01720.02430.0086
0.06240.01640.1236
0.03650.12120.04730.0571

図 2: 変形 BZ に よ る 結果の画像と 定義する 。 こ の c の近似値は， GAT 相関と 同様に次式で求める こ と ができ る 。

計算機実験を 行っ た。 表 4 にその結果を 示す。 変形は，A ∼ D で表す 4 種類の部分射影変換 c = (0.01, 0),(0.02, 0), (0, 0.01), (0, 0.02)， X で 表す 縦 0.9 倍横0.75 倍の縮小変換， Y で 表す時計方向20◦ 回転， Zで表す縦横を 15◦傾け る せん断を 組み合わせた 12 種類である 。 表の値は， 初期相関値， GAT 相関値， GPT
相関値， GPT 相関値から GAT 相関値を 引いた も のである 。 ま た， BZ 場合の変換画像を 図 2 に 示す。
表 4 よ り ， GPT 相関値の方が GAT 相関値の方が，射影変換に 対し て 高い値が得ら れ， 射影変換に 対して 頑健である こ と がわかる 。

c = (hx1 , x1 ix1 )−1 (hx1 , x1 ix1 − hx1 , x2 ix1 ) (9)ア フ ィ ン 変換と 部分射影変換を 交互に 収束する まで計算し ， そ の結果を 統合し て 相関を 最大に する 射影変換と する 。 そし て ， そのと き の相関値を GPT 相関と 定義する 。 図 1 に， GPT のフ ロ ーチャ ート を 示す

5

おわり に

本論文では， 射影変換に 対し て 頑健な 相関値を 与
え る GPT 相関を 提案し ， その有効性を 計算機実験によ っ て 確認し た 。 今後の課題と し て は， ア フ ィ ン 変
換と 部分射影変換の相互の関係を 解明し ， 計算精度を 向上さ せる こ と である 。

References[1] 若原 徹, 小高 和己, 大局的ア フ ィ ン 変換を 用いた形状マッ チン グ , 信学技報, PRMU96-2, May,
1996.
[2] 若原 徹, 木村 義政, 山下 喜義, GAT 反復を用いたア フ ィ ン 変換不変の濃淡文字相関法, 信学
技報, PRMU99-31, June, 1999.
図 1: GPT フ ロ ーチャ ート[3] Toru Wakahara, Yukihiko Yamashita, “Acceleration of GAT Correlation for Distortion-tolerant Im4 実験結果
age Matching”, Proc. of ICPR 2012, pp746–749,
手書き 数字を 使い， 射影変換に よ り 変形さ せた 画
Tsukuba, Dec. 11-15, 2012.
像と 原画像の GAT 相関値と GPT 相関値を 比較する
252

Effect of beam’s shape on adhesion between a multi-beam structure and rough surface
Student number：09_24320 name：Takehiko YAZAKI surpervisor：Kunio TAKAHASHI
the side surface of the beam adheres to the rigid bodyas shown in Fig 3 (b) and this contact is called area
contact.
The total energy depends on the work ofadhesion ∆𝛾, beam’s dimensions of the length L, width
W, 𝐸0 𝐼0 and n. Now the dimensionless parameter ofthe work of adhesion is defined as follows.

1. Introduction
Geckos can quickly move on rough surface. Thisability depends on adhesion force [1] which occurs
between objects and hair structure on their feet. In
order to apply this mechanism to development ofgrip devices, analyzing adhesion force is necessary.
Sekiguchi assumed contacts between a hairstructure and objects as a contact between elastic
beam and flat surface of a rigid body, and analyzed
adhesion force [2]. Using his model, effects of a hairstructure’s dimensions or elastic module on adhesion
force could be discussed, but it could discuss only
the uniform shape of a hair structure. Judging fromthe evidence that a hair structure of geckos has ashape like a spatula, analyzing the adhesion force
which is caused by hair structure of various shapes isimportant. In order to analyze the adhesion force
which is caused by hair structure of various shapes,
the purpose of this study is to advance Sekiguchi’s
adhesion model [2] as follows.
1) Modeling a beam of various shapes andanalyzing the adhesion force which occurs
between the beam and flat surface of a rigidbody.
2) Considering contacts between a multi-beam
structure and rough surface of a rigid body bysuperposition to the adhesion force of a single
beam.

2.

𝛤 = √(∆𝛾𝑊𝐿2 )/(2 𝐸0 𝐼0 )The energy is determined if the 𝛤/ tan 𝜃 and n areset, and the adhered area is equilibrium at theminimum energy condition. From the above, therelation between the displacement d and the force f isdetermined and the maximum tensile force whichcorresponds to arbitrary d equals to the grip force.

Figure 1: Distribution of the bending stiffness.

Model

2.1 The model of the contact between a single
beam and flat surface of a rigid body.
A beam of various shapes is modeled as a beam
which is changed in the distribution of bending
stiffness EI, expressed as follows.
𝑥 𝑛
𝐸𝐼 = 𝐸0 𝐼0 (1 − )(0 ≤ 𝑛 < 1)
𝐿
𝐸0 𝐼0 is EI at the base of beam. Figure 1 shows thisfunction.
Figure 2 shows various shapes of beams in response
to n. If n=0, then beam’s shape is straight, it is calledstraight beam in this study. If n>0, then the beam’s
upper surface is curved toward the edge of thebottom of the beam in response to n, it is called
tapering beam in this study. Deformation of thetapering beam is approximately represented by linear
beam theory of a straight beam.
When the beam approaches, loading process, flat
surface of a rigid body at angle θ, the beam firstly
contacts to the rigid body at the edge of one end of
the beam as shown in Fig 3 (a). This contact is called
line contact. The rigid body presses the edge of thebeam and the beam starts deforming. Consequently,

Figure 2: Various shapes of beam in response to n.

Figure 3: Schematic illustrations of contacts between anelastic beam and flat surface of a rigid body.

2.2 The model of the contact between amulti-beam structure and rough surface.
A multi-beam structure consists of beams of the53

same shape fixed to a rigid substrate as shown in
Fig 4 (a). Beams are fixed on board at even intervals
and even angles θ. The beams are arranged atintervals wide enough to avoid mutual contact.
Rough surface of a rigid body is represented as sine
wave with amplitude A as shown in Fig 4 (b).
In order to simplify states of a contact between
several beams and rigid body surface, the contact
area of rigid body surface with several beams isapproximated as flat surface is parallel to themulti-beam’s substrate as shown in Fig 5. Because ofthis approximation, the multi-beam structure’s force
is the sum of the forces between several contact
beams and flat surface of the rigid body withdifferent displacements d.

force of a multi-beam structure is larger for the caseof consisting of straight beams.
If the amplitude of surface roughness is large, thenumber of beams which adhere to rigid body surface
depends on beam’s shape, so the structure which
consisting of beams adhere easily can cause larger
grip force. When = 𝐿 sin 𝜃 /(𝑛 − 3) , this means atapering beam can adhere more easily than a straight
beam.

Figure 4: Schematic illustrations of a multi-beam
structure and rough surface of a rigid body.

Figure 6: Relation between the normalized amplitudeof surface roughness and the normalized grip force ofmulti-beam structure.

Figure 5: Simplification of contact area of the rigidbody.

3. Analysis results and discussion
A numerical simulation was done for analyzing thegrip force of a multi-beam structure under thecondition that 𝛤/ tan 𝜃 =2.0 and the number ofbeams is 100.
Figure 6 shows the relation between the amplitudeof surface roughness and the grip force of themulti-beam structure in case of the structure
consisting of straight beams (n = 0) or tapering
beams (n = 0.99). Comparing the two structures ofmulti straight beam and multi tapering beam, when
the amplitude is small the grip force of the multistraight beam structure is relativity larger, when theamplitude is large, the grip force of the multi
tapering beam structure is relativity larger. Thistendency can be discussed by analyzing contacts
between a single beam and a flat surface of rigid
body.
If the amplitude of surface roughness is small, thegrip force of a multi-beam structure depends on theadhesion force of individual beams. Figure 7 shows
the relation between the displacement and the force
of a single beam. The adhesion force of a tapering
beam is smaller than a straight beam, so the grip

Figure 7: Relation between the normalized force of asingle beam and the normalized displacement.

4. Conclusion
The model of a beam of various shapes as shown
Fig 2 was proposed. Analysis results of the grip forceof multi-beam structures against rough surface of arigid body showed the possibility that beam’s shape
can increase application performance to surface
roughness.

Reference[1] K. Autumn, et al., “Adhesive force of a single
gecko foot-hair”, Nature 405, 2000, 681-685[2] Yu Sekiguchi, et al., “Adhesion between sidesurface of an elastic beam and flat surface of a rigidbody”, Journal of Adhesion Science and Technology
Vol. 26, No. 23, December 2012, 2615–2626

54

Influences of atmospheric conditions to EHD thruster for planetary explorations
Student Number:09-26298

Name: Tianyue Li

1. Introduction
When a high electric voltage such as tens kV isapplied to a capacitor of an asymmetric structure, as
shown in Fig.1,a net force directing to a physicallysmaller electrode is generated. This physical[1]
phenomenon is called Biefeld-Brown Effect.
Athruster based on this effect is usually called
Electrohydrodynmaics (EHD) thruster or Asymmetric
Capacitor Thruster (ACT).
Because the EHD thruster has simple structure, it issuitable as a propulsion system for in-situ planetary
explorations.
AlthoughthemechanismofBiefeld-Brown Effect is not clearly understood so far,some basic feature of the EHD thruster has been
[2]
experimentally investigated.
For example, therelationship between a force on the electrode anddimensions of the electrode are one of them. In thisstudy, the EHD thruster in various atmosphericpressures and gas species are experimentallyinvestigated.

Supervisor: Yoshihiro MOCHIMARU

current data are monitored and saved in a PC during
the measurement.

Fig.2 Experimental System
3.Experimental results
Figure 3 shows the relationship between forces on
the EHD thruster and atmospheric pressures undervarious constant applied voltages. In the cases thatthe applied voltage is relatively low, the EHD thruster
obtains a stronger force under a lower atmosphericpressure condition. However, when the applied
voltage exceeds a certain threshold, the atmosphericpressure’s influence becomes negligible. If the voltage
continues to increase, the force becomes
proportional to the atmospheric pressure.

Fig.1 Structure of EHD thruster
2.Experimental method
The experimental system is shown in Fig. 2. The
main components of this system are a model of theEHD thruster , an electronic balance, and a DC highvoltage power supply. The EHD thruster is placed in aacrylic gas chamber to change the atmosphericconditions around the EHD thruster. As a voltage
between asymmetric electrodes changes, a force ismeasured by the electronic balance. The voltage andFig.3 Force on EHD thruster vs. atmospheric Pressure

55

stronger propulsion force than that in CO 2 . However,as atmospheric pressure increases, the force in CO 2
increases with the prossure. Therefore, force in CO 2

Figure 4 shows the relationship between theatmospheric pressure and the current. The current
becomes lower as the atmospheric pressure isincreased. From Fig.3 and Fig.4, we cannot find thatthe current dominates the propulsion force of theEHD thruster . As voltage changes, the force presents
different variation trend with atmospheric pressure.

becomes bigger than that in air at some pressure
values.
Figure 6 shows the relationship between the
maximum propulsion force on the EHD thruster andthe atmospheric pressure in air or CO 2 . The
maximum force is almost proportional to theatmospheric pressure, in both air and CO 2 under thepresent experimental conditions.

Fig.4 Current vs. Pressure
Figure 5 shows the relationship between the force
and the atmospheric pressure under various applied
voltages in air or CO 2 . EHD thruster in air obtains
Fig.6 Maximal force vs. Pressure in Air or CO 2

4.Conclusion
In this study, influences of atmospheric conditions to
a EHD thruster are investigated. Besides the applied
voltage and the dimensions of the EHD thruster, theatmospheric condition is indeed an importantparameter. For improvements in the EHD thruster
performance, detailed understandings of theionization and dissociation process of gas particles
around the EHD thruster are necesary.
5.References[1]Biefeld-Brown Effect
http://blazelabs.com/l-intro.asp
access 19 Feb,2013[2]Wei Liu, “Experiment on Planetary Propulsion
System using Electrostatic Force”, graduation thesis
Tokyo Institute ofTechnology,2010
Fig.5 Force vs. Pressure on Air or CO 2

56

COMPARATIVE ANALYSIS OFAVIATION PASSENGER ATTRIBUTES IN VIETNAM
Student Number: 09-26387

Name: Le Yen Lan

Supervisor: Shinya HANAOKA

1 Introduction
Vietnam’s economy is developing rapidly
during a last few years, thus income level isalso increasing. This triggers an increasing
aviation demand. Compared with 2007, thenumber of domestic passenger in 2011
increased 75%. Then, the wave of low-costcarriers (LCCs) has come to Vietnam.
Compared to other Southeast Asian countries,
Vietnam is the latest country that LCCs
appears after Malaysia, Singapore, Thailand,
Indonesia and Philippines. In 2008, after being
renamed from Pacific Airlines, Jetstar Pacific
Airlines became the first LCC in Vietnam. The
second one is VietJet Air which started
operation in 2011. There are 2 questions in thisresearch: (1) How is aviation passenger market
changed in Vietnam after LCC operation? (2)How different are Vietnam’s LCCs compared toother Asian LCCs?
In this research, passenger questionnaire
survey was conducted in two biggest airports ofVietnam – Noi Bai airport in Ha Noi (HN) and
Tan Son Nhat airport in Ho Chi Minh City
(HCM) with two purposes: (1) To clarify
Vietnam’s LCCs’ current situation. (2) To
identify and compare the passenger attributes
between LCCs and legacy carriers in Vietnam.

Table 1- Routes, Frequency and Fare
O

D

km

VNA

VietjetJetstar

Sam-

F

AF

F

AF

F

AF

ples

HN

HCM

1719

16

117

8

73

8

71

1056

HN

DN

759

9

86

2

53

1

61

205

HN

DLat

1481

1

85

1

~

0

~

8

HCM

HN

1719

17

114

8

94

8

101

632

HCM

DN

960

8

90

3

70

3

81

252

HCM

HPh

1618

4

114

1

124

2

97

191

HCM

Vinh

1426

3

107

1

81

2

88

81

HCM

P.Qu

500

5

89

1

46

0

~

37

HCM

Hue

1065

1

91

1

91

0

~

49

HCM

N.Tr

439

4

74

1

44

0

~

41

DN

HCM

960

8

~

4

~

3

~

0

Total

2552

F: frequency ( flights/ day)
AF: Average fare ( USD)

2.2 Monthly Household Income
Most of Vietnam Airlines’ passengers are morethan 1500USD in their monthly householdincome. But most passengers of LCCs are
between 400 and 600USD (see Figure 1).
Because of the entry of LCCs, middle income
earners can travel using airlines.

2 Results
The number of collected effective samples ismore than 2500 in two airports. Major results
are shown below.

2.1 Routes and Fare:
VietJet Air has more routes than Jetstar but
in some routes Jetstar has more flights than
VietJet Air as shown in Table 1. In most routes,
the average fare of LCCs was cheaper than
Vietnam Airlines (VNA). The durations ofsurvey at Ho Chi Minh City is near Vietnamese
Lunar New Year, so fare of LCCs was increased,sometimes more expensive than Vietnam
Airlines.

Figure 1: Monthly Household Income

2.3 Travel purpose
57

Passengers use Vietnam Airlines for work
purpose, but use LCCs for other private
purposes (see Figure 2). The survey duration isnear Lunar New Year in HCM, so passengers
who use for “back home” purpose are more
numerous than usually.

Table 2: Previous transportation of airtransport’s first boarding passengers
VNA
Never 60

Mekong Jetstar

VietJetAir

Pacific

Air

10

87

106

use
Rail

35

11

106

120

Bus

22

8

94

76

Car

4

0

6

5

Other

1

0

0

1

Total

122/577 29/126

293/877 308/957

2.6 Alternative choice
Most people will choose other airlines if itmeets their schedule. More users of Vietnam
Airlines will choose same airline than otherairlines. The ratio of LCCs’ users who will
choose bus or railway is higher than Vietnam
Airlines’ users. One of possible reason is thatfare of LCCs is not much more expensive thanbus or train, but fare of Vietnam Airlines ismuch more expensive.

Figure 2: Travel purpose

2.4 Reasons of choosing airlines
Figure 3 shows the result of survey question“What are the reason to choose this airline thistime for comparing with other airlines?”. Basedon the result, passengers choose LCCs becauseof low fare, suitable flight time and convenient
ticket. On the other hand passengers use legacy
carriers because of less delay, no cancel flight,
safety, well-known and comfortable seat.

Figure 4: Alternative choice

3 Conclusion
With the entry of low-cost carriers, Vietnam
aviation market has changed significantly.
More persons who are lower income can travel
by air transport. And more people who used touse bus or train changed to use air. People canuse air transport for other private purposes
than before.
Vietnam’s economy will continue to develop
and demand for aviation will potentially
increase. It creats more opportunity for LCCs
to operate in Vietnam.

Figure 3: Reasons of choosing airlines

2.5 Transport mode before using aviation
There are many new demands on both ofVietnam Airlines and LCCs. There are many
LCCs demand shift from rail and bus (see Table
2).

58

Prediction of Indoor MIMO Channel for Specific Systems using
Ray Tracing Simulation
学籍番号: 10B15328

1

名前: チョウ

カツシ

指導教員: 高田

潤一・金

T2

T3

ミンソク

概要

現在, 無線通信技術への期待は高まる一方で, 周波数資源の需要もそれと共に増えている. その結果 3 GHz 以
下の周波数帯域は隙間なく割当られ, 周波数が枯渇するという厳しい現状になってきた. 解決策として高周波数帯における通信技術の開発が強く求められている. まT1

R2

た周波数資源を有効に利用して通信の高速化を図るた

R1

R3

めに, MIMO 技術の導入が期待される [1]. 周波数が高くなるにつれて回折の効果が小さくなり，マルチパスの効果が減少するため，特に 10GHz を超える高い周波
図 1: 百年記念館

数帯における MIMO チャネルの特性は，現在の無線通表 1: 実験パラメータ周波数11 GHz帯域幅
400 MHz
アンテナ高さ1.7 m
アンテナ
偏波多重 12 素子円形アレー

信で用いられている 6GHz 以下の周波数とはかなり異なると考えられる．
このようなチャネルの特性を調べるためには, サウン
ディング実験 [2] – [5] が広く用いられているが, 実験するためのコストが膨大である. 本論文では，レイト
レースシミュレーションによるチャネル特性の予測に

3 レイトレースシミュレーション

注目し，サウンディング実験の結果と比較することでレイトレースシミュレーションには商用シミュレータ

シミュレーションの適用範囲と誤差要因を検討するこである RapLab [4] を使用した．レイトレースシミュレーとを目的としている．

ションは幾何光学近似に基づき，電波を光に見立て，パ

まずは 11 GHz において屋内サウンディング実験を行っ

スを検索する [6]．RapLab の場合には，拘束条件とした, 次に実験環境における構造物をモデル化してレイトて最大反射回数・最大回折回数を与えることで，有限レースシミュレーションを行った. シミュレーション

個のパスが探索される．具体的なシミュレーションの結果は無限の角度・遅延時間分解能を有するため，実パラメタは表 2 の通りである．

測結果と直接比較することは困難である．本研究では，

レイトレースシミュレーションをすることにより，パレイトレースシミュレーションの結果から，アンテナ

ス損失，遅延，放射方向，到来方向などが分かるが．シ

指向性と測定信号の帯域幅を考慮して，予測される実

ミュレーション結果はそのまま実験結果と比較できな

験結果を算出し，実際の実験結果と比較検討を行った．

いため，本論文ではシミュレーション結果に実験条件

比較には遅延スペクトル及び双角度スペクトルを用いの帯域制限とアンテナ応答を適用し，予測される実験

た．比較結果から，シミュレーションの妥当性と誤差結果を求めている [7] [8]．

要因を考察した．

4 実験結果とシミュレーション結果の比較

2

MIMO 伝送は次の式でモデル化される．

屋内 MIMO チャネルサウンディング実験

東京工業大学大岡山キャンパス百年記念館（図 1）にy(f ) = H(f )x(f ) + n(f )おいて、MIMO チャネルサウンディング実験を行った。

(1)

実験で得られたデータは各アンテナ間の複素電力から

x(f ) は送信信号スペクトル，y(f ) は受信信号スペク
トル，n(f ) は雑音スペクトルであり，H(f ) は MIMO

なる伝達関数である。

チャネル伝達関数である．実験で測定した結果とシミュ実験を表 1 のようなパラメータで行った [2]。

レーションの結果に対し同じ手法で解析を行う．遅延スペクトルは伝達関数のフーリエ逆変換より求めることができる．双角度電力スペクトルは，角度に対する59表 2: シミュレーションパラメータ周波数11 GHz帯域
無限大
アンテナ高さ1.7 m
最大反射回数2回
最大回折回数1回構造物の媒質モデル コンクリート

Double−Directional Spectrum in Azimuth

Double−Directional Spectrum in Azimuth

−180
−135

−80 −180

−80

−82

−82

−135

−84

−84
−90

−86
Azimuth at Rx [deg.]

Azimuth at Rx [deg.]

−90
−45

−88

0

−90
−92

45

−86
−45

−88

0

−90
−92

45

−94
90

−94
90

−96
135

−98

180
−180 −135 −90

−45
0
45
90
Azimuth at Tx [deg.]

測定結果

135

180

−96
135

−100 180−180 −135 −90

−98
−45
0
45
90
Azimuth at Tx [deg.]

135

180

シミュレーション結果

図 3: 角度スペクトル出した．実験結果とシミュレーション結果を比較する
測定結果

ことで以下のような結論を得た．

シミュレーション結果

シミュレーションはコストは小さいが，周囲の環境を

図 2: 遅延スペクトル

想定していないためシミュレータの精度が高くない．また，レイトレースによりチャネルの特性を正しく推定できることに加え，直接波を含めいくつかの強いパス

アレー応答ベクトルを用いて，ビームフォーミングにも正しく推定できる．しかし，それぞれの位置におけより次の式で求められる．

Yf (φT , φR ) =∗aH
f (φR )Hf af (φT )22||aH
f (φR )|| ||af (φT )||るチャネルの詳細な性質や弱いパスの推定には限界が
ある．

(2)

今後の研究ではシミュレーションの精度と実用性をあ

af (φR ), af (φT ) はそれぞれ受信アンテナ，送信アンテ

げる．具体的には，小さい物体を単純な形状で近似しナのアレー応答ベクトルである [2] [3] [9]．

てモデル化することでシミュレーション確度が上がる

図 2 は送信点を T1 から T2 まで直線的に移動したときかどうかを検討する．

の測定結果及びシミュレーション結果より得られた遅延

参考文献

スペクトルである．図 2 より，モデルの材質の反射率が異なるため，レベルの差が見られるが，直接波（LoS）と電力の強い一回反射波は双方で一致している．遅延の大きいパスについては，正しく推定できなかった．
図 1 の T1 と T2 の中点 (T3) における双角度スペクト
ルが図 3 である．図 2 と同様にシミュレーションでは，
直接波と一部の一回反射波のみ推定され，電力の強いパスの角度をすべて推定するのは困難だった．シミュレーションでは建物壁面のみを考慮しているのに対して，実験では，多くの椅子やテーブルなどの物体が配置されていた，従って誤差の主原因は，レイトレースシ
ミュレーションで考慮されていない物体の影響と考えられる．レイトレースでは幾何光学近似を用いるため，比較的小さい物体や表面が複雑な形状に対する散乱を正しく取り扱うことができない．したがって，物理光学近似など，より近似度の高い解法を導入して，これら小物体の散乱を考慮することで，シミュレーションの確度を上げることができると期待される [4] [6]．

5結論

本研究ではサウンディング実験結果をビームフォーミ
ングの手法で解析した．またモデル化した環境においてレイトレースシミュレーションをし，それにアンテ
ナ指向性と測定信号の帯域幅を考慮して予測結果を算60[1] 大種武雄・小川恭孝, 分かりやすい MIMO システム技術,
オーム社, 2009.
[2] 菊間信良, アダプティブアンテナ技術, オーム社, 2003.
[3] Andreas F. Mollisch, Wireless Communications, WILEY,
Now York, USA, 2005.
[4] Raplab
Kozo
KeikakuEngineeringInc.,”http://www4.kke.co.jp/network/raplab/index.html”[5] Andrea Goldsmith, Wireless Communications, Cambridge,
New York, USA, 2005.
[6] 吉敷由起子, レイトレース法を用いた電波伝搬シミュレー
ションのための非幾何光学の手法の導入, 東京工業大学
博士論文, 2011.
[7] Andreas Richter, Estimation of Radio Channel Parameters:
Models and Algorithms, Doctoral Thesis, Ilmenau University of Technology, 2005.
[8] Martin Steinbauer, The Radio Propagation Channel –A Non-Directional, Directional and Double-Directional
Point-of-View, Doctoral Thesis, Technical University of Vienna, 2001.
[9] Yuyuan CHANG, et.al. ”Analysis of Field Measurementwith Wideband MIMO Sounder at 11 GHz Frequency”, IEICE, Technical Report, AP2012-153 Jan 2013.

−100