SSeennssoorrss && TTrraannssdduucceerrss Volume 6, Special Issue August 2009 www.sensorsportal.com ISSN 1726-5479 Editors-in-Chief: professor Sergey Y. Yurish, phone: +34 696067716, fax: +34 93 4011989, e-mail: editor@sensorsportal.com Guest Editors: Subhas Chandra Mukhopadhyay, Gourab Sen Gupta and Ray Yueh Min Huang Editors for Western Europe Meijer, Gerard C.M., Delft University of Technology, The Netherlands Ferrari, Vittorio, Universitá di Brescia, Italy Editor South America Costa-Felix, Rodrigo, Inmetro, Brazil Editor for Eastern Europe Sachenko, Anatoly, Ternopil State Economic University, Ukraine Editors for North America Datskos, Panos G., Oak Ridge National Laboratory, USA Fabien, J. Josse, Marquette University, USA Katz, Evgeny, Clarkson University, USA Editor for Asia Ohyama, Shinji, Tokyo Institute of Technology, Japan Editor for Asia-Pacific Mukhopadhyay, Subhas, Massey University, New Zealand Editorial Advisory Board Abdul Rahim, Ruzairi, Universiti Teknologi, Malaysia Ahmad, Mohd Noor, Nothern University of Engineering, Malaysia Annamalai, Karthigeyan, National Institute of Advanced Industrial Science and Technology, Japan Arcega, Francisco, University of Zaragoza, Spain Arguel, Philippe, CNRS, France Ahn, Jae-Pyoung, Korea Institute of Science and Technology, Korea Arndt, Michael, Robert Bosch GmbH, Germany Ascoli, Giorgio, George Mason University, USA Atalay, Selcuk, Inonu University, Turkey Atghiaee, Ahmad, University of Tehran, Iran Augutis, Vygantas, Kaunas University of Technology, Lithuania Avachit, Patil Lalchand, North Maharashtra University, India Ayesh, Aladdin, De Montfort University, UK Bahreyni, Behraad, University of Manitoba, Canada Baliga, Shankar, B., General Monitors Transnational, USA Baoxian, Ye, Zhengzhou University, China Barford, Lee, Agilent Laboratories, USA Barlingay, Ravindra, RF Arrays Systems, India Basu, Sukumar, Jadavpur University, India Beck, Stephen, University of Sheffield, UK Ben Bouzid, Sihem, Institut National de Recherche Scientifique, Tunisia Benachaiba, Chellali, Universitaire de Bechar, Algeria Binnie, T. David, Napier University, UK Bischoff, Gerlinde, Inst. Analytical Chemistry, Germany Bodas, Dhananjay, IMTEK, Germany Borges Carval, Nuno, Universidade de Aveiro, Portugal Bousbia-Salah, Mounir, University of Annaba, Algeria Bouvet, Marcel, CNRS – UPMC, France Brudzewski, Kazimierz, Warsaw University of Technology, Poland Cai, Chenxin, Nanjing Normal University, China Cai, Qingyun, Hunan University, China Campanella, Luigi, University La Sapienza, Italy Carvalho, Vitor, Minho University, Portugal Cecelja, Franjo, Brunel University, London, UK Cerda Belmonte, Judith, Imperial College London, UK Chakrabarty, Chandan Kumar, Universiti Tenaga Nasional, Malaysia Chakravorty, Dipankar, Association for the Cultivation of Science, India Changhai, Ru, Harbin Engineering University, China Chaudhari, Gajanan, Shri Shivaji Science College, India Chavali, Murthy, VIT University, Tamil Nadu, India Chen, Jiming, Zhejiang University, China Chen, Rongshun, National Tsing Hua University, Taiwan Cheng, Kuo-Sheng, National Cheng Kung University, Taiwan Chiang, Jeffrey (Cheng-Ta), Industrial Technol. Research Institute, Taiwan Chiriac, Horia, National Institute of Research and Development, Romania Chowdhuri, Arijit, University of Delhi, India Chung, Wen-Yaw, Chung Yuan Christian University, Taiwan Corres, Jesus, Universidad Publica de Navarra, Spain Cortes, Camilo A., Universidad Nacional de Colombia, Colombia Courtois, Christian, Universite de Valenciennes, France Cusano, Andrea, University of Sannio, Italy D'Amico, Arnaldo, Università di Tor Vergata, Italy De Stefano, Luca, Institute for Microelectronics and Microsystem, Italy Deshmukh, Kiran, Shri Shivaji Mahavidyalaya, Barshi, India Dickert, Franz L., Vienna University, Austria Dieguez, Angel, University of Barcelona, Spain Dimitropoulos, Panos, University of Thessaly, Greece Ding, Jianning, Jiangsu Polytechnic University, China Djordjevich, Alexandar, City University of Hong Kong, Hong Kong Donato, Nicola, University of Messina, Italy Donato, Patricio, Universidad de Mar del Plata, Argentina Dong, Feng, Tianjin University, China Drljaca, Predrag, Instersema Sensoric SA, Switzerland Dubey, Venketesh, Bournemouth University, UK Enderle, Stefan, Univ.of Ulm and KTB Mechatronics GmbH, Germany Erdem, Gursan K. Arzum, Ege University, Turkey Erkmen, Aydan M., Middle East Technical University, Turkey Estelle, Patrice, Insa Rennes, France Estrada, Horacio, University of North Carolina, USA Faiz, Adil, INSA Lyon, France Fericean, Sorin, Balluff GmbH, Germany Fernandes, Joana M., University of Porto, Portugal Francioso, Luca, CNR-IMM Institute for Microelectronics and Microsystems, Italy Francis, Laurent, University Catholique de Louvain, Belgium Fu, Weiling, South-Western Hospital, Chongqing, China Gaura, Elena, Coventry University, UK Geng, Yanfeng, China University of Petroleum, China Gole, James, Georgia Institute of Technology, USA Gong, Hao, National University of Singapore, Singapore Gonzalez de la Rosa, Juan Jose, University of Cadiz, Spain Granel, Annette, Goteborg University, Sweden Graff, Mason, The University of Texas at Arlington, USA Guan, Shan, Eastman Kodak, USA Guillet, Bruno, University of Caen, France Guo, Zhen, New Jersey Institute of Technology, USA Gupta, Narendra Kumar, Napier University, UK Hadjiloucas, Sillas, The University of Reading, UK Haider, Mohammad R., Sonoma State University, USA Hashsham, Syed, Michigan State University, USA Hasni, Abdelhafid, Bechar University, Algeria Hernandez, Alvaro, University of Alcala, Spain Hernandez, Wilmar, Universidad Politecnica de Madrid, Spain Homentcovschi, Dorel, SUNY Binghamton, USA Horstman, Tom, U.S. Automation Group, LLC, USA Hsiai, Tzung (John), University of Southern California, USA Huang, Jeng-Sheng, Chung Yuan Christian University, Taiwan Huang, Star, National Tsing Hua University, Taiwan Huang, Wei, PSG Design Center, USA Hui, David, University of New Orleans, USA Jaffrezic-Renault, Nicole, Ecole Centrale de Lyon, France Jaime Calvo-Galleg, Jaime, Universidad de Salamanca, Spain James, Daniel, Griffith University, Australia Janting, Jakob, DELTA Danish Electronics, Denmark Jiang, Liudi, University of Southampton, UK Jiang, Wei, University of Virginia, USA Jiao, Zheng, Shanghai University, China John, Joachim, IMEC, Belgium Kalach, Andrew, Voronezh Institute of Ministry of Interior, Russia Kang, Moonho, Sunmoon University, Korea South Kaniusas, Eugenijus, Vienna University of Technology, Austria Katake, Anup, Texas A&M University, USA Kausel, Wilfried, University of Music, Vienna, Austria Kavasoglu, Nese, Mugla University, Turkey Ke, Cathy, Tyndall National Institute, Ireland Khan, Asif, Aligarh Muslim University, Aligarh, India Sapozhnikova, Ksenia, D.I.Mendeleyev Institute for Metrology, Russia Kim, Min Young, Kyungpook National University, Korea South Ko, Sang Choon, Electronics and Telecommunications Research Institute, Korea South Kockar, Hakan, Balikesir University, Turkey Kotulska, Malgorzata, Wroclaw University of Technology, Poland Kratz, Henrik, Uppsala University, Sweden Kumar, Arun, University of South Florida, USA Kumar, Subodh, National Physical Laboratory, India Kung, Chih-Hsien, Chang-Jung Christian University, Taiwan Lacnjevac, Caslav, University of Belgrade, Serbia Lay-Ekuakille, Aime, University of Lecce, Italy Lee, Jang Myung, Pusan National University, Korea South Lee, Jun Su, Amkor Technology, Inc. South Korea Lei, Hua, National Starch and Chemical Company, USA Li, Genxi, Nanjing University, China Li, Hui, Shanghai Jiaotong University, China Li, Xian-Fang, Central South University, China Liang, Yuanchang, University of Washington, USA Liawruangrath, Saisunee, Chiang Mai University, Thailand Liew, Kim Meow, City University of Hong Kong, Hong Kong Lin, Hermann, National Kaohsiung University, Taiwan Lin, Paul, Cleveland State University, USA Linderholm, Pontus, EPFL - Microsystems Laboratory, Switzerland Liu, Aihua, University of Oklahoma, USA Liu Changgeng, Louisiana State University, USA Liu, Cheng-Hsien, National Tsing Hua University, Taiwan Liu, Songqin, Southeast University, China Lodeiro, Carlos, Universidade NOVA de Lisboa, Portugal Lorenzo, Maria Encarnacio, Universidad Autonoma de Madrid, Spain Lukaszewicz, Jerzy Pawel, Nicholas Copernicus University, Poland Ma, Zhanfang, Northeast Normal University, China Majstorovic, Vidosav, University of Belgrade, Serbia Marquez, Alfredo, Centro de Investigacion en Materiales Avanzados, Mexico Matay, Ladislav, Slovak Academy of Sciences, Slovakia Mathur, Prafull, National Physical Laboratory, India Maurya, D.K., Institute of Materials Research and Engineering, Singapore Mekid, Samir, University of Manchester, UK Melnyk, Ivan, Photon Control Inc., Canada Mendes, Paulo, University of Minho, Portugal Mennell, Julie, Northumbria University, UK Mi, Bin, Boston Scientific Corporation, USA Minas, Graca, University of Minho, Portugal Moghavvemi, Mahmoud, University of Malaya, Malaysia Mohammadi, Mohammad-Reza, University of Cambridge, UK Molina Flores, Esteban, Benemérita Universidad Autónoma de Puebla, Mexico Moradi, Majid, University of Kerman, Iran Morello, Rosario, University "Mediterranea" of Reggio Calabria, Italy Mounir, Ben Ali, University of Sousse, Tunisia Mulla, Imtiaz Sirajuddin, National Chemical Laboratory, Pune, India Neelamegam, Periasamy, Sastra Deemed University, India Neshkova, Milka, Bulgarian Academy of Sciences, Bulgaria Oberhammer, Joachim, Royal Institute of Technology, Sweden Ould Lahoucine, Cherif, University of Guelma, Algeria Pamidighanta, Sayanu, Bharat Electronics Limited (BEL), India Pan, Jisheng, Institute of Materials Research & Engineering, Singapore Park, Joon-Shik, Korea Electronics Technology Institute, Korea South Penza, Michele, ENEA C.R., Italy Pereira, Jose Miguel, Instituto Politecnico de Setebal, Portugal Petsev, Dimiter, University of New Mexico, USA Pogacnik, Lea, University of Ljubljana, Slovenia Post, Michael, National Research Council, Canada Prance, Robert, University of Sussex, UK Prasad, Ambika, Gulbarga University, India Prateepasen, Asa, Kingmoungut's University of Technology, Thailand Pullini, Daniele, Centro Ricerche FIAT, Italy Pumera, Martin, National Institute for Materials Science, Japan Radhakrishnan, S. National Chemical Laboratory, Pune, India Rajanna, K., Indian Institute of Science, India Ramadan, Qasem, Institute of Microelectronics, Singapore Rao, Basuthkar, Tata Inst. of Fundamental Research, India Raoof, Kosai, Joseph Fourier University of Grenoble, France Reig, Candid, University of Valencia, Spain Restivo, Maria Teresa, University of Porto, Portugal Robert, Michel, University Henri Poincare, France Rezazadeh, Ghader, Urmia University, Iran Royo, Santiago, Universitat Politecnica de Catalunya, Spain Rodriguez, Angel, Universidad Politecnica de Cataluna, Spain Rothberg, Steve, Loughborough University, UK Sadana, Ajit, University of Mississippi, USA Sadeghian Marnani, Hamed, TU Delft, The Netherlands Sandacci, Serghei, Sensor Technology Ltd., UK Saxena, Vibha, Bhbha Atomic Research Centre, Mumbai, India Schneider, John K., Ultra-Scan Corporation, USA Seif, Selemani, Alabama A & M University, USA Seifter, Achim, Los Alamos National Laboratory, USA Sengupta, Deepak, Advance Bio-Photonics, India Shearwood, Christopher, Nanyang Technological University, Singapore Shin, Kyuho, Samsung Advanced Institute of Technology, Korea Shmaliy, Yuriy, Kharkiv National Univ. of Radio Electronics, Ukraine Silva Girao, Pedro, Technical University of Lisbon, Portugal Singh, V. R., National Physical Laboratory, India Slomovitz, Daniel, UTE, Uruguay Smith, Martin, Open University, UK Soleymanpour, Ahmad, Damghan Basic Science University, Iran Somani, Prakash R., Centre for Materials for Electronics Technol., India Srinivas, Talabattula, Indian Institute of Science, Bangalore, India Srivastava, Arvind K., Northwestern University, USA Stefan-van Staden, Raluca-Ioana, University of Pretoria, South Africa Sumriddetchka, Sarun, National Electronics and Computer Technology Center, Thailand Sun, Chengliang, Polytechnic University, Hong-Kong Sun, Dongming, Jilin University, China Sun, Junhua, Beijing University of Aeronautics and Astronautics, China Sun, Zhiqiang, Central South University, China Suri, C. Raman, Institute of Microbial Technology, India Sysoev, Victor, Saratov State Technical University, Russia Szewczyk, Roman, Industrial Research Inst. for Automation and Measurement, Poland Tan, Ooi Kiang, Nanyang Technological University, Singapore, Tang, Dianping, Southwest University, China Tang, Jaw-Luen, National Chung Cheng University, Taiwan Teker, Kasif, Frostburg State University, USA Thumbavanam Pad, Kartik, Carnegie Mellon University, USA Tian, Gui Yun, University of Newcastle, UK Tsiantos, Vassilios, Technological Educational Institute of Kaval, Greece Tsigara, Anna, National Hellenic Research Foundation, Greece Twomey, Karen, University College Cork, Ireland Valente, Antonio, University, Vila Real, - U.T.A.D., Portugal Vaseashta, Ashok, Marshall University, USA Vazquez, Carmen, Carlos III University in Madrid, Spain Vieira, Manuela, Instituto Superior de Engenharia de Lisboa, Portugal Vigna, Benedetto, STMicroelectronics, Italy Vrba, Radimir, Brno University of Technology, Czech Republic Wandelt, Barbara, Technical University of Lodz, Poland Wang, Jiangping, Xi'an Shiyou University, China Wang, Kedong, Beihang University, China Wang, Liang, Advanced Micro Devices, USA Wang, Mi, University of Leeds, UK Wang, Shinn-Fwu, Ching Yun University, Taiwan Wang, Wei-Chih, University of Washington, USA Wang, Wensheng, University of Pennsylvania, USA Watson, Steven, Center for NanoSpace Technologies Inc., USA Weiping, Yan, Dalian University of Technology, China Wells, Stephen, Southern Company Services, USA Wolkenberg, Andrzej, Institute of Electron Technology, Poland Woods, R. Clive, Louisiana State University, USA Wu, DerHo, National Pingtung Univ. of Science and Technology, Taiwan Wu, Zhaoyang, Hunan University, China Xiu Tao, Ge, Chuzhou University, China Xu, Lisheng, The Chinese University of Hong Kong, Hong Kong Xu, Tao, University of California, Irvine, USA Yang, Dongfang, National Research Council, Canada Yang, Wuqiang, The University of Manchester, UK Yang, Xiaoling, University of Georgia, Athens, GA, USA Yaping Dan, Harvard University, USA Ymeti, Aurel, University of Twente, Netherland Yong Zhao, Northeastern University, China Yu, Haihu, Wuhan University of Technology, China Yuan, Yong, Massey University, New Zealand Yufera Garcia, Alberto, Seville University, Spain Zagnoni, Michele, University of Southampton, UK Zamani, Cyrus, Universitat de Barcelona, Spain Zeni, Luigi, Second University of Naples, Italy Zhang, Minglong, Shanghai University, China Zhang, Qintao, University of California at Berkeley, USA Zhang, Weiping, Shanghai Jiao Tong University, China Zhang, Wenming, Shanghai Jiao Tong University, China Zhang, Xueji, World Precision Instruments, Inc., USA Zhong, Haoxiang, Henan Normal University, China Zhu, Qing, Fujifilm Dimatix, Inc., USA Zorzano, Luis, Universidad de La Rioja, Spain Zourob, Mohammed, University of Cambridge, UK Sensors & Transducers Journal (ISSN 1726-5479) is a peer review international journal published monthly online by International Frequency Sensor Association (IFSA). Available in electronic and on CD. Copyright © 2009 by International Frequency Sensor Association. All rights reserved. SSeennssoorrss && TTrraannssdduucceerrss JJoouurrnnaall CCoonntteennttss Volume 6 Special Issue August 2009 www.sensorsportal.com ISSN 1726-5479 Research Articles Foreword: Modern Sensing Technologies - II Subhas Chandra Mukhopadhyay, Gourab Sen Gupta, Ray Yueh Min Huang .................................. 1 GPS Navigation Processing Using the IMM-Based Extended Kalman Filter Dah-Jing Jwo, Chien-Hao Tseng ....................................................................................................... 4 A Robot with Complex Facial Expressions J. Takeno, K. Mori and Y. Naito.......................................................................................................... 18 A Bayesian Approach to Solving the Non – Invasive Time Domain Reflectometry Inverse Problem Ian G. Platt and Ian M. Woodhead ..................................................................................................... 27 Estimation of Back-Surface Flaw Depth by Laminated Piezoelectric Highpolymer Film Akinobu Yamamoto, Shiro Biwa and Eiji Matsumoto ......................................................................... 43 Studies on Gas Sensing Performance of Pure and Surface Modified SrTiO3 Thick Film Resistors V. B. Gaikwad, D. D. Kajale, Y. R. Baste, S. D. Shinde, P. K. Khanna, N. K. Pawar, D. N. Chavan, M. K. Deore, G. H. Jain ...................................................................................................................... 57 Computational Aspects of Sensor Network Protocols (Distributed Sensor Network Simulator) Vasanth Iyer, S. S. Iyengar, G. Rama Murthy, M. B. Srinivas…………………………………………… 69 Sensitivity Enhancement of Wheatstone Bridge Circuit for Resistance Measurement Tarikul Islam, Shakeb A. Khan, Sheikh S. Islam, Harsh .................................................................... 92 Experimental Evaluation of the Applicability of Capacitive and Optical Measurement Methods for the Determination of Liquid Hydrogen Volume Flow Gert Holler, Daniel Hrach, Anton Fuchs ............................................................................................. 105 Effect of Organic Vapour on Porous Alumina Based Moisture Sensor in Dry Gases Saakshi Dhanekar, Tarikul Islam, S. S. Islam, Kamalendu Sengupta, Debdulal Saha...................... 117 Compensation of Imperfections for Vibratory Gyroscope Systems Using State Observers Chien-Yu Chi and Tsung-Lin Chen .................................................................................................... 128 Authors are encouraged to submit article in MS Word (doc) and Acrobat (pdf) formats by e-mail: editor@sensorsportal.com Please visit journal’s webpage with preparation instructions: http://www.sensorsportal.com/HTML/DIGEST/Submition.htm International Frequency Sensor Association (IFSA). Sensors & Transducers Journal, Vol. 6, Special Issue, August 2009, pp. 1-3 1 SSSeeennnsssooorrrsss &&& TTTrrraaannnsssddduuuccceeerrrsss ISSN 1726-5479 © 2009 by IFSA http://www.sensorsportal.com Foreword: Modern Sensing Technologies - II This special issue on Modern Sensing Technologies - II of the Sensors & Transducer journal is primarily focused on the different aspects of design, theoretical analysis, fabrication, characterization and experimentation of different sensing technologies. This special issue comprises 10 papers carefully selected from the reviewed papers which were presented in the 3rd International Conference on Sensing Technology (ICST´ 2008), November 30 to December 3, 2008, held at National Cheng-Kung University, Tainan, Taiwan and published in the conference proceedings. This special issue presents the extended version of the conference papers. This Special Issue has focused on different aspects of modern sensing technology, i.e. intelligent measurement, information processing, adaptability, recalibration, data fusion, validation, high reliability and integration of novel and high performance sensors in the areas of magnetic, ultrasonic, vision and image sensing, wireless sensors and network, microfluidic, tactile, gyro, flow, surface acoustic wave, humidity and ultra-wide band. While future interest in this field is ensured by the constant supply of emerging modalities, techniques and engineering solutions, as well as an increasing need from aging structures, many of the basic concepts and strategies have already matured and now offer opportunities to build upon. We are very happy to be able to offer the readers of the Sensors & Transducer journal such a diverse Special Issue both in terms of its topical coverage and geographic representation. We do hope that the journal readers will find it interesting, thought provoking, and useful in their research and practical engineering work. We would like to extend our wholehearted thanks to all the authors who have contributed their work to this Special Issue. Sensors & Transducers Journal, Vol. 6, Special Issue, August 2009, pp. 1-3 2 Guest Editors: Subhas Chandra Mukhopadhyay School of Engineering and Advanced Technology (SEAT), Massey University (Turitea Campus) Palmerston North, New Zealand S.C.Mukhopadhyay@massey.ac.nz Dr. Subhas Chandra Mukhopadhyay graduated from the Department of Electrical Engineering, Jadavpur University, Calcutta, India in 1987 with a Gold medal and received the Master of Electrical Engineering degree from Indian Institute of Science, Bangalore, India in 1989. He obtained the PhD (Eng.) degree from Jadavpur University, India in 1994 and Doctor of Engineering degree from Kanazawa University, Japan in 2000. During 1989-90 he worked almost 2 years in the research and development department of Crompton Greaves Ltd., India. In 1990 he joined as a Lecturer in the Electrical Engineering department, Jadavpur University, India and was promoted to Senior Lecturer of the same department in 1995. Obtaining Monbusho fellowship he went to Japan in 1995. He worked with Kanazawa University, Japan as researcher and Assistant professor till September 2000. In September 2000 he joined as Senior Lecturer in the Institute of Information Sciences and Technology, Massey University, New Zealand where he is working currently as an Associate professor. His fields of interest include Sensors and Sensing Technology, Electromagnetics, control, electrical machines and numerical field calculation etc. He has authored over 200 papers in different international journals and conferences, edited eight conference proceedings. He has also edited four special issues of international journals as guest editor and three books with Springer-Verlag. He is a Fellow of IET (UK), a senior member of IEEE (USA), an associate editor of IEEE Sensors journal and IEEE Transactions on Instrumentation and Measurements. He is in the editorial board of e-Journal on Non-Destructive Testing, Sensors and Transducers, Transactions on Systems, Signals and Devices (TSSD), Journal on the Patents on Electrical Engineering, Journal of Sensors. He is in the technical programme committee of IEEE Sensors conference, IEEE IMTC conference and IEEE DELTA conference. He was the Technical Programme Chair of ICARA 2004 and ICARA 2006. He was the General chair of ICST 2005, ICST 2007. He has organized the IEEE Sensors conference 2008 at Lecce, Italy as General Co-chair and currently organizing the IEEE Sensors conference 2009 at Christchurch, New Zealand as General Chair. Gourab Sen Gupta School of Engineering and Advanced Technology (SEAT), Massey University (Turitea Campus) Palmerston North, New Zealand G.SenGupta@massey.ac.nz Mr. Gourab Sen Gupta graduated with a Bachelor of Engineering (Electronics) from the University of Indore, India in 1982. In 1984 he got his Master of Electronics Engineering from the University of Eindhoven, The Netherlands. In 1984 he joined Philips India and worked as an Automation Engineer in the Consumer Electronics division till 1989. Sensors & Transducers Journal, Vol. 6, Special Issue, August 2009, pp. 1-3 3 Thereafter he worked as a Senior Lecturer in the School of Electrical and Electronic Engineering at Singapore Polytechnic, Singapore. He is with Massey University, Palmerston North, New Zealand, since September 2002 as a Senior Lecturer. He is currently carrying out studies for his PhD degree. He has published over 70 papers in international journals and conference proceedings, co-authored two books on programming, edited four conference proceedings, and edited two special issues of international journals (IEEE Sensors Journal and IJISTA) as guest editor. He is a senior member of IEEE. His area of interest is robotics, autonomous systems, embedded controllers and vision processing for real-time applications. He was the General Chair of ICARA 2004 and ICARA 2006 conferences. He was the Technical Programme chair of ICST 2005 and ICST 2007 conferences. Ray Yueh Min Huang Department of Engineering Science National Cheng-Kung University Tainan, Taiwan E-mail: huang@mail.ncku.edu.tw Ray Yueh-Min Huang is a Distinguished Professor and Chairman of the Department of Engineering Science, National Cheng-Kung University, Taiwan, R.O.C. His research interests include Multimedia Communications, Wireless Networks, Embedded Systems, and Artificial Intelligence. He received his MS and Ph.D. degrees in Electrical Engineering from the University of Arizona in 1988 and 1991 respectively. He has co-authored 2 books and has published about 160 refereed professional research papers. He has completed 10 Ph.D. and over 80 M.Sc. thesis students. Dr. Huang has received many research awards, such as the Best Paper Award of 2007 IEA/AIE Conference, Best Paper Award of the Computer Society of the Republic of China in 2003, the Awards of Acer Long-Term Prize in 1996, 1998, and 1999, Excellent Research Awards of National Microcomputer and Communication Contests in 2006. He also received many funded research grants from National Science Council, Ministry of Education, Industrial Technology of Research Institute, and Institute of Information Industry. Dr. Huang has been invited to give talks or served frequently in the program committee at national and international conferences. Dr. Huang is in the editorial board of the Journal of Wireless Communications and Mobile Computing, the Journal of Internet Technology, International Journal of Internet Protocol Technology, International Journal of Ad Hoc and Ubiquitous Computing, Journal of Security and Communication Networks and serves as an associate editor for Journal of Computer Systems, Networks, and Communications as well as International Journal of Communication Systems. He was the Technical Programme Chair of Symposium on Digital Life Technologies (SDLT2007). He was the General chair of VIP2007. He is organizing the SDLT2008, PCM2008 and ICST2008. Huang is a member of the IEEE as well as IEEE communication, computer, and computational intelligence societies. ___________________ 2009 Copyright ©, International Frequency Sensor Association (IFSA). All rights reserved. (http://www.sensorsportal.com) Sensors & Transducers Journal, Vol. 6, Special Issue, August 2009, pp. 4-17 4 SSSeeennnsssooorrrsss &&& TTTrrraaannnsssddduuuccceeerrrsss ISSN 1726-5479 © 2009 by IFSA http://www.sensorsportal.com GPS Navigation Processing Using the IMM-Based Extended Kalman Filter Dah-Jing JWO, Chien-Hao TSENG Department of Communications, Navigation and Control Engineering National Taiwan Ocean University 2 Pei-Ning Rd., Keelung 202-24, Taiwan Tel: +886-2-24622192 ext. 7209, fax: +886-2-24633492 E-mail: djjwo@mail.ntou.edu.tw Received: 10 July 2009 /Accepted: 31 July 2009 /Published: 10 August 2009 Abstract: This paper presents an interacting multiple model (IMM)-based extended Kalman filter (EKF) approach for the Global Positioning System (GPS) navigation processing. The “soft-switching” IMM estimator obtains its estimate as a weighted sum of the individual estimates from a number of parallel filters matched to different motion modes of the platform. The proposed method is applied to the GPS navigation to increase the navigation estimation accuracy at the high dynamic regions while preserving (without sacrificing) the precision at the low dynamic regions. In the illustrated example, the vehicle motion is designed to cover a broad range of dynamic behaviors: constant velocity, constant acceleration, variable acceleration, and circular turn. Simulation results show that the IMM-based EKF can substantially improve overall navigation accuracy as compared to that of single model EKF. The mode probability of the proposed IMM filter is also depicted. Copyright © 2009 IFSA. Keywords: GPS, interacting multiple model (IMM), extended Kalman filter (EKF) 1. Introduction The Global Positioning System (GPS) is a satellite-based navigation system that provides a user with the proper equipment access to useful and accurate positioning information anywhere on the globe. The well-known Kalman filter [1-3], which provides optimal (minimum mean square error) estimate of the system state vector, has been widely applied to the fields of navigation such as GPS receiver position/velocity determination. Sensors & Transducers Journal, Vol. 6, Special Issue, August 2009, pp. 4-17 5 The implementation of Kalman filter requires the a priori knowledge of both the process and measurement models. It is widely known that poorly designed mathematical model for the EKF may lead to the divergence. Clearly, if the plant parameters are subject to perturbations and dynamics of the system are too complex to be characterized by an explicit mathematical model, an adaptive scheme is needed. An adaptive Kalman filter (AKF) can be utilized as the noise-adaptive filter to adjust the parameters. It is well known that the process model is dependent on the dynamical characteristics of the vehicle onto which the navigation system is placed. To fulfill the requirement of achieving the filter optimality or to preventing divergence problem of Kalman filter, the so-called AKF approach has been proposed to prevent divergence problem of the EKF when precise knowledge on the system models are not available. Many efforts have been made to improve the estimation of the covariance matrices [4, 5]. The Interacting multiple model (IMM) algorithm [6,7] has the configuration that runs in parallel several model-matched state estimation filters, which exchange information (interact) at each sampling time. The IMM algorithm has been originally applied to target tracking [8-10]. A model probability evaluator calculates the current probability of the vehicle being in each of the possible modes. A global estimate of the vehicle’s state is computed using the latest mode probabilities. This algorithm carries out a soft-switching between the various modes by adjusting the probabilities of each mode, which are used as weightings in the combined global state estimate. The covariance matrix associated with this combined estimate takes into account the covariances of the mode-conditioned estimates as well as the differences between these estimates. The AKF approach is based on parametric adaptation, while the IMM approach is based on filter structural adaptation (model switching). The IMM-based method allows the possibility of using highly dynamic models just when required, diminishing unrealistic noise considerations in non-maneuvering situations of the system. The idea of defining different models according to the situation in which the vehicle is involved is not new. In this approach, the use of an IMM method allows exploiting the benefits of high dynamic models in the problem of vehicle navigation. Two models of the vehicle dynamic are defined. The non-maneuvering model, e.g., constant velocity model (CV model, also known as the PV model) reproduces properly straight trajectories of the vehicle. The maneuvering model, e.g., constant acceleration (CA model, also known as the PVA model) considers sharp turns and brusque accelerations in the vehicle state. The IMM-EKF filter developed calculates the probability of success of each model at every filter execution scan, supplying a realistic combined solution for the behavior of the vehicle. These probabilities are calculated according to a Markovian model for the transition between maneuver states. The remainder of this paper is organized as follows. In Section 2, preliminary background on GPS navigation using the extended Kalman filter method based on a single model is described. The basic structure of the IMM estimator is introduced in Section 3. In Section 4, simulation experiments on GPS navigation are carried out to evaluate the performance of the proposed approach in comparison to those by single-model EKF. Conclusions and some potential future work are drawn in Section 5. 2. GPS Navigation Using the Extended Kalman Filter The extended Kalman filter is a nonlinear version of the Kalman filter and is widely used for the position estimation in GPS receivers. The process model and measurement model for the Kalman filter are represented as kkk k wxfx +=+ ),(1 (1a) kkk k vxhz += ),( , (1b) Sensors & Transducers Journal, Vol. 6, Special Issue, August 2009, pp. 4-17 6 where the state vector n k ℜ∈x , process noise vector n k ℜ∈w , measurement vector m k ℜ∈z , and measurement noise vector m k ℜ∈v . In Equation (1), both the vectors kw and kv are zero mean Gaussian white sequences having zero crosscorrelation with each other: ikkikE δQww =][ T , ikkikE δRvv =][ T , kandiallforE ik 0vw =][ T , (2) where ][⋅E represents expectation, and superscript “T” denotes matrix transpose. kQ is the process noise covariance matrix, kR is the measurement noise covariance matrix, and kΦ is the state transition matrix. The symbol ikδ stands for the Kronecker delta function: ⎩ ⎨ ⎧ ≠ = = ki ki ik ,0 ,1 δ The discrete-time extended Kalman filter algorithm is summarized as follow: - Correction steps/measurement update equations: 1T 1| T 1| ][ − −− += kkkkkkkkk RHPHHPK (3) ]ˆ[ˆˆ 1|1|| −− −+= kkkkkkkk zzKxx , ),ˆ(ˆ 1| kkkk xhz =− (4) 1|| ][ −−= kkkkkk PHKIP (5) - Prediction steps/time update equations: )ˆ(ˆ ||1 kkkk xfx =+ (6) kkkkkkk QΦPΦP +=+ T ||1 (7) Equations (3-5) are the measurement update equations, and Equations (6-7) are the time update equations of the algorithm from k to step 1+k . These equations incorporate a measurement value into a priori estimation to obtain an improved a posteriori estimation. In the above equations, kP is the error covariance matrix defined by ])ˆ)(ˆ[( T kkkkE xxxx −− , in which kx̂ is an estimation of the system state vector kx , and the weighting matrix kK is generally referred to as the Kalman gain matrix. The Kalman filter algorithm starts with an initial condition value, 0 1|ˆ −kkx and 0 1| −kkP . When new measurement kz becomes available with the progression of time, the estimation of states and the corresponding error covariance would follow recursively ad infinity. The linear approximation equations for system and measurement matrices are obtained through the relations kk k |x̂xx fΦ =∂ ∂ ≈ ; 1|ˆ −=∂ ∂ ≈ kk k xxx hH (8) Further detailed discussion can be referred to Gelb [1] and Brown and Hwang [2]. The flow chart for the GPS navigation processing using extended Kalman filter approach is shown in Fig. 1. It is noted that the Sensors & Transducers Journal, Vol. 6, Special Issue, August 2009, pp. 4-17 7 predicted state shown in Equation (6) has been replaced by its linear version, i.e., kkkkk ||1 ˆˆ xΦx =+ , in this research work. Initialize state vector and state covariance matrix: 0 1|ˆ −kkx and 0 1| −kkP start Compute Kalman gain matrix from state covariance and estimated measurement covariance: 1T 1| T 1| ][ − −− += kkkkkkkkk RHPHHPK Multiply prediction error vector by Kalman gain matrix to get state correction vector and update state vector: ]ˆ[ˆˆ 1|1|| −− −+= kkkkkkkk zzKxx Update error covariance [ ] 1|| −−= kkkkkk PHKIP Predict new state vector and state covariance matrix kkkkk ||1 ˆˆ xΦx =+ kkkkkkk QΦPΦP T +=+ ||1 Fig. 1. Flow chart of the GPS navigation using the extended Kalman filter. 3. The Interacting Multiple Model Algorithm The IMM algorithm uses model (Markov chain state) probabilities to weight the input and output of a bank of parallel Kalman filters at each time instant. The approach which was initiated in takes into account a set of models ( kΦ and/or kQ ) to represent the system behavior patterns or system model. The overall estimates is obtained by a combination of the estimates from the filters running in parallel based on the individual models that match the system modes. In each cycle it consists of four major steps: interaction (mixing), filtering, mode probability calculation, and combination, as illustrated in Fig. 2. As a structure adaptation algorithm, the IMM estimators can substantially improve navigation accuracy during vehicle maneuvering (such as circular motion and acceleration) as well as during constant velocity straight-line motion over the conventional EKF. The IMM algorithm takes into account a set of models to represent the system behavior patterns or system model. The approach uses model (Markov chain state) probabilities to weight the input and output of a bank of parallel Kalman filters at each time instant. The overall estimates is obtained by a combination of the estimates from the filters running in parallel based on the individual models that match the system modes. In each cycle it consists of four major steps: interaction (mixing), filtering, mode probability calculation, and combination. The IMM Sensors & Transducers Journal, Vol. 6, Special Issue, August 2009, pp. 4-17 8 has the following properties: (1) it consists of a low bandwidth filter for the nearly uniform motion and a high bandwidth filter for the maneuvering situations; (2) these filters interact (exchange information) with time-varying weights; (3) the final estimate is a combination of each filter’s estimate, with the weights being the mode probabilities; (4) the weights for interaction and combination are based on which model fits better the data and other factors, such as the expected transition from one mode to another. 1 1|1 1 1|1 , ˆ −−−− kkkk Px 2 1|1 2 1|1 , ˆ −−−− kkkk Px Estimate & covariance combination Interaction (Mixing) EKF (Model 1) EKF (Model 2) Model probability update kz 1 kΛ 2 kΛ 1 | 1 | ,ˆ kkkk Px 2 | 2 | ,ˆ kkkk Px kkkk || ,ˆ Px kk |µ 1|1 −− kkµ 01 1|1 01 1|1 ,ˆ −−−− kkkk Px 02 1|1 02 1|1 ,ˆ −−−− kkkk Px Fig. 2. The IMM structure (one cycle with two models). The IMM algorithm includes the following stages: - Calculation of mixing probabilities. The probability that mode iM was in effect at 1−k given that jM is in effect at k conditioned on 1−kz is given by }|{},|{1 1111 | 1|1 zMzMMµ c −−−−−− = k i kk i k j k j ji kk PP (9) which are the mixing probabilities and can be written as µPµ c i kij j ji kk 1 | 1|1 1 −−− = , rji ,....,1, = , (10) where the normalizing constants ∑ = − = r i i kijj 1 1µPc (11) Sensors & Transducers Journal, Vol. 6, Special Issue, August 2009, pp. 4-17 9 - Model interaction (Mixing). Starting with x̂ 1 i k − , the mixed initial condition for the filter matched to j kM is computed µxx ji kk r i i kk j kk | 1|1 1 1|1 0 1|1 ˆˆ −− = −−−− ∑= , rj ,...,2,1= (12) The corresponding covariance is }][][{ ˆˆˆˆ 0 1|11|1 0 1|11|11|1 1 | 1|1 0 1|1 Tj kk i kk j kk i kk i kk r i ji kk j kk xxxxPµP −−−−−−−−−− = −−−− −⋅−+=∑ (13) - Model-matched filtering. The likelihood function corresponding to the r filters ],|[ 1−= k j kk j k p zMzΛ (14) are computed using the mixed initial Equation (12) and the associated covariance Equation (13) as ],,|[ 0 1|1 0 1|1ˆ PxΛ Mz j kk j kk j kk j k p −−−− = - Mode probability update. The probability of each mode is updated by cµ Λ C j j k j k 1 = , (15) where cΛC j r j j k∑ = = 1 (16) is the normalization constant, and ⎥⎦ ⎤ ⎢⎣ ⎡−= − j k j kk Tj kj kk j k υPυ P Λ 1 | | )( 2 1exp ||2 1 π , (17) where j kkk j k 1|ˆ −−= zzυ . - State estimate and covariance calculation through combination. Combination of the model-conditioned estimates and covariance is done according to the mixture equations given by ∑ = = r j j k j kkkk 1 || ˆˆ µxx , rj ,....,1= (18) { }T kk j kkkk j kk j kk r j j kkk ][][ ˆˆˆˆ ||||| 1 | xxxxPµP −⋅−+=∑ = (19) The filters in the bank of parallel filters could be the EKFs to describe their individual behavior in each maneuvering stage, resulting in the IMM-based extended Kalman filter (IMM-EKF). Sensors & Transducers Journal, Vol. 6, Special Issue, August 2009, pp. 4-17 10 4. GPS Navigation Performance Evaluation Simulation experiments have been carried out to evaluate the performance of the IMM-EKF approach in comparison with the conventional methods for GPS navigation processing. The computer codes were developed by the authors using the Matlab® software. The commercial software Satellite Navigation toolbox by GPSoft LLC was employed for generating the satellite positions and pseudoranges. The block diagram of the GPS navigation processing using the IMM-EKF is shown in Fig. 3. Reference GPS kx̂ kz IMM-EKF Measurement prediction )ˆ(ˆ 1|1| −− = kkkk xhz kzδ+ - kx̂δ Fig. 3. GPS navigation processing using the IMM-EKF. The dynamic process of the GPS receiver in lower dynamic environment can be represented by the PV model. In such case, the GPS navigation filter with three position states, three velocity states, and two clock states are considered, so that the state to be estimated is a 18× vector. The process model governed is given by ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ + ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ = ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ 8 7 6 4 2 8 7 6 5 4 3 2 1 8 7 6 5 4 3 2 1 0 0 0 00000000 10000000 00000000 00100000 00000000 00001000 00000000 00000010 u u u u u x x x x x x x x x x x x x x x x & & & & & & & & , (20) where 1x , 3x , 5x represent the east, north, and vertical position; 2x , 4x , 6x represent the east, north, and vertical velocity; and 7x and 8x represent the receiver clock offset and drift errors, respectively. The state transition matrix for the model can be found to be ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ ∆ ∆ ∆ ∆ = 10000000 1000000 00100000 0010000 00001000 0000100 00000010 0000001 t t t t kΦ (21) Sensors & Transducers Journal, Vol. 6, Special Issue, August 2009, pp. 4-17 11 Consider only the pseudorange observables are available, the measurement connection matrix kH are the partial derivatives of the predicted measurements with respect to each state, which can be constructed as follows: ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ = 01000 01000 01000 )()()( )2()2()2( )1()1()1( n z n y n x zyx zyx k hhh hhh hhh MMMMMMMM H (22) While considering the vehicle maneuvering or turning, the PVA model will be more suitable than the PV one. The PVA model adds three acceleration variables to the original PV model, leading to the 111× state vector [2]. The following transition probability matrices of the Markov chain were used ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ = 99.001.0 01.099.0 Pij The simulation scenario is as follows. The experiment was conducted on a simulated vehicle trajectory originating from the position of North 25.14923442 degrees and East 121.77749768 degrees at an altitude of 0 meter. This is equivalent to [ ]T2694047.794911032.253042299.42- m in the WGS-84 ECEF coordinate system. The location of the origin is defined as the (0,0,0) m location in the local tangent East-North-Up (ENU) coordinate frame. The three dimensional simulated vehicle trajectory is shown as in Fig. 4. Furthermore, the description of the vehicle motion is listed in Table 1. The vehicle motion is designed to cover a broad range of dynamic behaviors (constant velocity, constant acceleration, variable acceleration, and circular turn). The vehicle position and velocity in the east, north, and vertical components is provided in Figs. 5 and 6, respectively, for providing better insight into vehicle dynamic information in each time interval. Fig. 7 shows the ground track of the trajectory. The satellite geometry for the simulation is shown in Fig. 8. Figs. 9 and 10 provide the GPS navigation results for the IMM-EKF approach as compared to that of standard EKF. It can be seen that substantial estimation accuracy improvement is obtained by using the proposed technique, discussed as follows: Fig. 4. Three dimensional trajectory of the simulated vehicle. Starting point Sensors & Transducers Journal, Vol. 6, Special Issue, August 2009, pp. 4-17 12 Table 1. Description of the vehicle motion. Time interval (sec) Motion [0-50] Constant velocity straight line [51-150] Climb for 100 seconds [151-350] Circular motion with descending [351-400] Climb for 50 seconds [401-450] Clockwise turn [451-500] Constant velocity straight line level flight [501-700] Circular motion level flight [701-850] Constant velocity straight line level flight [851-950] Clockwise turn level flight [951-1100] Constant velocity straight line level flight [1101-1150] Counter-clockwise turn [1151-1300] Descend for 150 seconds [1301-1400] Constant velocity straight line Fig. 5. Vehicle position in the east, north, and vertical components. Fig. 6. Vehicle velocity in the east, north, and vertical components. Sensors & Transducers Journal, Vol. 6, Special Issue, August 2009, pp. 4-17 13 Fig. 7. Ground track of the trajectory. Fig. 8. The satellite geometry for the simulation. In the time intervals where the vehicle is not maneuvering and is conducting constant-velocity straight-line motion, the navigation accuracy based on the EKF and that on IMM-EKF are equivalent. In the time intervals where the vehicle is maneuvering, the mismatch of the model leads the conventional EKF to navigation performance degradation. Using the IMM-EKF, a position accuracy improvement of approximately 90% in the horizontal and over 50% in the vertical components, respectively, can be achieved; a velocity accuracy improvement of approximately 80% in the horizontal and over 10% in the vertical components, respectively, can be achieved. The position root mean squared errors are shown in Fig. 11. Table 2 shows the comparison of navigation accuracies for the two approaches. To inspect the correctness of the filter switching capability, the mode probability of the proposed IMM filter is depicted in Fig. 12. Table 2. Comparison of navigation accuracies for the two approaches (in meters). East North Altitude Position (m) 10.5293 8.7938 3.8128 EKF Velocity (m/s) 3.1817 2.9315 0.3668 Position (m) 0.8387 0.8948 1.5248 IMM-EKF Velocity (m/s) 0.5991 0.5593 0.2952 4 10 13 7 22 14 8 1 Sensors & Transducers Journal, Vol. 6, Special Issue, August 2009, pp. 4-17 14 (a) East (b) North (c) Altitude Fig. 9. Position errors for the IMM-EKF as compared to the EKF. EKF IMM-EKF EKF IMM-EKF EKF IMM-EKF Sensors & Transducers Journal, Vol. 6, Special Issue, August 2009, pp. 4-17 15 (a) East (b) North (c) Altitude Fig. 10. Velocity errors for the IMM-EKF as compared to the EKF. EKF IMM-EKF EKF IMM-EKF EKF IMM-EKF Sensors & Transducers Journal, Vol. 6, Special Issue, August 2009, pp. 4-17 16 (a) Position (b) Velocity Fig. 11. Comparison RMSE for (a) position and (b) velocity using IMM-EKF and EKF. Fig. 12. Model probability of the IMM-EKF approach. Sensors & Transducers Journal, Vol. 6, Special Issue, August 2009, pp. 4-17 17 5. Conclusions and Future Work The paper presents an alternative approach for designing an adaptive Kalman filter and provides an example in the Global Positioning System. An interacting multiple model based method has been presented to improve the results obtained by the extended Kalman filter implementation for GPS navigation processing. Two models that represent different dynamic conditions (maneuvering and non-maneuvering) have been carried out. Simulation experiments for GPS navigation have been provided to illustrate the accessibility. The navigation accuracy based on the proposed method has been compared to the EKF and has demonstrated substantial improvement in both navigational accuracy and tracking capability. Future research can be conducted to undertake the implementation of the following issues. Incorporation of the artificial intelligence such as fuzzy logic and neural network for further performance improvement is accessible. Utilization of the adaptive Kalman filter or nonlinear filters to replace the standard Kalman filter is another way for further performance improvement. Furthermore, a suggested direction for potential application is investigation of GPS RAIM (Receiver Autonomous Integrity Monitoring) in civil aviation. Acknowledgements The authors gratefully acknowledge the support of the National Science Council of the Republic of China through Grant NSC 96-2221-E-019-007 and NSC 97-2221-E-019-012. References [1]. A. Gelb, Applied optimal estimation, M.I.T. Press, MA, 1974. [2]. R. G. Brown, P. Y. C. Hwang, Introduction to random signals and applied Kalman filtering, (3rd edition), John Wiley & Sons, New York, 1997. [3]. P. Axelrad, R. G. Brown, GPS navigation algorithms, in: B. W. Parkinson, J. J. Spilker, P. Axelrad and P. Enga, (Ed.), Global Positioning System: Theory and Applications, Vol. I, AIAA, Washington DC, Chap. 9, 1996. [4]. R. K. Mehra, Approaches to adaptive filtering. IEEE Trans. Automat. Contr., Vol. AC-17, 1972, pp. 693-698. [5]. Mohamed, A. H. & Schwarz, K. P., Adaptive Kalman filtering for INS/GPS, Journal of Geodesy, 73, 1999, 4, pp. 193-203. [6]. X. R. Li, Y. Bar-Shalom, Estimation with Applications to Tracking and Navigation, John Wiley & Sons, Inc. 1993. [7]. X. R. Li, Y. Bar-Shalom, Performance Prediction of the Interacting Multiple Model Algorithm, IEEE Transactions on Aerospace and Electronic System, Vol. 29, No. 3, 1993, pp. 755-771. [8]. X. R. Li, Y. Bar-Shalom, Design of an Interacting Multiple Model Algorithm for Air Traffic Control Tracking, IEEE Transactions on Control System Technology, Vol. 1, No. 3, 1993, pp.186-194. [9]. L. A. Johnston, V. Krishnamurthy, An Improvement to the Interacting Multiple Model (IMM) Algorithm, IEEE Transactions on Signal Processing, Vol. 49, No. 12, 2001, pp.2909-2923. [10]. Y.-S. Kim, K.-S. Hong, An IMM algorithm for tracking maneuvering vehicles in an adaptive cruise control environment, International Journal of Control, Automation, and Systems, Vol. 2, No. 3, 2004, pp. 310-318. [11]. Lee, B. J., Park, J. B., Lee, H. J. and Joo, Y. H., Fuzzy-logic-based IMM algorithm for tracking a maneuvering target, IEE Proceedings-Radar Sonar Navig., 152, 1, 2005, pp. 16-22. ___________________ 2009 Copyright ©, International Frequency Sensor Association (IFSA). All rights reserved. (http://www.sensorsportal.com) Sensors & Transducers Journal, Vol. 6, Special Issue, August 2009, pp. 18-26 18 SSSeeennnsssooorrrsss &&& TTTrrraaannnsssddduuuccceeerrrsss ISSN 1726-5479 © 2009 by IFSA http://www.sensorsportal.com A Robot with Complex Facial Expressions J. Takeno, K. Mori and Y. Naito Meiji University Graduate School, Kawasaki-shi, Japan Received: 10 July 2009 /Accepted: 31 July 2009 /Published: 10 August 2009 Abstract: The authors believe that the consciousness of humans basically originates from languages and their association-like flow of consciousness, and that feelings are generated accompanying respective languages. We incorporated artificial consciousness into a robot; achieved an association flow of language like flow of consciousness; and developed a robot called Kansei that expresses its feelings according to the associations occurring in the robot. To be able to fully communicate with humans, robots must be able to display complex expressions, such as a sense of being thrilled. We therefore added to the Kansei robot a device to express complex feelings through its facial expressions. The Kansei robot is actually an artificial skull made of aluminum, with servomotors built into it. The face is made of relatively soft polyethylene, which is formed to appear like a human face. Facial expressions are generated using 19 servomotors built into the skull, which pull metal wires attached to the facial “skin” to create expressions. The robot at present is capable of making six basic expressions as well as complex expressions, such as happiness and fear combined. Copyright © 2009 IFSA. Keywords: Robot, Kansei robot, humanoid robots, facial expressions 1. Introduction Humanoid robots, such as Honda’s ASIMO, are being actively developed in recent years. These robots are compact and light, and some can walk on two legs, and even run. They can communicate with humans although not very fluently, and carry light objects in their two hands. Industrial robots are working in factories, and Healthcare robots are helping people in nursing homes. The Actroid is an increasingly popular type of receptionist robot working at corporate information counters. This trend will grow in the future, with many more robots working together with humans. Sensors & Transducers Journal, Vol. 6, Special Issue, August 2009, pp. 18-26 19 Humans and robots must have a deep mutual understanding when working together, and communication is important for establishing a good relationship between them. If robots could intuitively respond like humans by expressing their feelings when communicating with humans, the relationship between humans and robots would greatly improve and robots would be able to play more important and innovative roles in a variety of fields when working together with humans. Studies are going on at Massachusetts Institute of Technology (Kismet) and Tokyo University of Science (Saya) [1], [2]. Those studies use techniques that avoid addressing the issue of consciousness itself. In contrast, our study is backed by our idea of languages and their association-like flow of consciousness and kansei information related to individual languages. The authors have incorporated artificial consciousness into a robot, achieved linguistic associations resembling a flow of consciousness [3], and developed the Kansei robot that expresses its feelings according to the associations [4]. The authors’ objective was to develop a robot that responds to words emotionally. 2. About Kansei The Japanese word “kansei” is a concept that encompasses both sensibility and emotion in English. General ideas about kansei typically include: (1) Sensitivity of the sensory organs that generate sense and perception in response to external stimuli; (2) Experience invoked and governed by sensations, including feelings, impulses and desires accompanying sensations; (3) Sensory desires controlled by reason. Information that stimulates and affects human consciousness is defined as kansei information in the authors’ present study. The authors have already developed a system to systematically extract association information and kansei information for cognized objects from Internet websites and process the derived data. In our present study, we have defined a series of processes in this system as Kansei Information Processing. We calculate the kansei value based on this Kansei Information Processing, and have the robot make facial expressions conforming to the calculated kansei [5]. 3. About Kansei Information Processing This chapter describes our method of Kansei information processing. First, to represent Kansei information by a computer, we define concepts of "association Kansei network (abbreviated to AK network)" and "consciousness network" that model human thinking and knowledge representation of memory. The AK network is created by connecting words that can be associated from a given word. Since the network depends on personal experience and memory, different networks are formed depending on the person. For example, if the object is "grass," the associated words could be "leaf," "herb," or "flower." Then, we define a consciousness network using the AK network. For example, if you see "grass," there are associated words such as "leaf," "herb," and "flower" on the AK network. At this time, we extract a network that consists of these associated words with "grass" as the central word. We define this extracted part of the AK network as the "consciousness network." In other words, we define the consciousness network as a subgraph of the AK network that appears as a stream of consciousness and we think that the components of the network change every moment as time elapses. This means that changes of the graph state caused by the generation, growth, and decline of the consciousness network Sensors & Transducers Journal, Vol. 6, Special Issue, August 2009, pp. 18-26 20 due to external and internal stimuli such as visual sensations and vocal sounds. In the example above, when the object of interest shifts from "grass" to "flower," a consciousness network with "grass" as the center shifts to a consciousness network with "flower" as the center. We define this as the stream of consciousness (Fig. 1). Fig. 1. AK network and consciousness network. We consider that actions in daily life use the consciousness network described above. First of all, in this research we think that the center of consciousness is always one object. Generally, a human, however, performs multiple actions at the same time. Currently, we think that multiple consciousness networks are generated and active, and that the brain processes these consciousness networks in parallel. For example, we suppose that a human might perform three actions at the same time: "eating," "watching television," and "having a conversation." At this time, consciousness networks related to "eating," "watching television," and "having a conversation" exist respectively and a stream of consciousness to each action occurs. However, we think that the central object of consciousness is always only one object, the object is in one consciousness network, and the center of the consciousness shifts when the cognition circumstances or priority changes.(Fig. 2). Fig. 2. Example of shifting the center of consciousness. Sensors & Transducers Journal, Vol. 6, Special Issue, August 2009, pp. 18-26 21 In this study, we limited the consciousness network to only one network at this stage for the sake of simplification. We achieved a robot with consciousness by using a method of changing an associable word successively as time elapsed to change the consciousness network dynamically, and at the same time, by calculating the change of Kansei of each associated word. Put differently, the robot can change its facial expression according to a change in the Kansei value caused by an association change of the robot's internal language. Using this method, the robot can react with feeling. For the facial expressions of the robot, we decided on six basic expressions based on the ideas of P. Ekman: "happiness," "disgust," "fear," "anger," "sadness," and "surprise." First, we describe the definition of association value. When a human is conscious of one object, he or she associates many other things from it. Some of these things may be associated strongly at a given moment, but some things may be associated weakly at a later time. We define this concept of the numerical value of the degree of association as the "association value." Then, we describe the definition of the Kansei value. The Kansei value is a concept relating to the numerical value of the degree of various types of Kansei and mental images that a person has subjectively for an object or an event. The Kansei value is a concept of consciousness, no matter whether the object elicits happiness, sadness or disgust. These types of Kansei may be mixed, and the method for achieving such complex Kansei in a robot is also a target of this research. Then, we describe the association Kansei value (AK value). The AK value is a Kansei value owned by a given word associated from a keyword that a human is currently conscious of. Then, we describe the process of creating numerical values of the association value, Kansei value and AK value. As described above, the association value is the strength of the relationship between the nodes on a consciousness network. As a method of extracting this association value, we use text information from the Internet. We explain the method below. (1) Search a site that includes a keyword of the object of consciousness. (2) Extract a sentence that includes this keyword. (3) Consider all words in this sentence as candidates for words associated (associated words) from the keyword of the center of the current consciousness. At this time, store the number of occurrences of the keyword and of the candidates of the associated word in the database. (4) Calculate the probability of the occurrence of candidates of the associated words in a sentence that includes the keyword, and evaluate it as association value. Then, we describe the numerical values of the Kansei value. In a similar way to the calculation of the association value, we calculate the value based on text information from the Internet. ＜1＞ Steps (1) and (2) are similar. ＜2＞ In the database, store the number of occurrences of the keyword that is the center of consciousness in the sentence and the number of occurrences of the evaluation factor relating to happiness ("joyful," "interesting," and "joy"), the evaluation factor relating to disgust ("dirty," "bad- Sensors & Transducers Journal, Vol. 6, Special Issue, August 2009, pp. 18-26 22 tasting," and "rotten") and the evaluation factor relating to each type of feeling such as sadness, anger, and fear. ＜3＞ Calculate the probability of the occurrence of the evaluation factor of each type of feeling in a sentence that includes the keyword, and set it to be the Kansei value of each type of feeling. Then, we compare the Kansei values and define the largest value as the Kansei value of the robot relating to the keyword. Each Kansei value can be used for various purposes, but we will not describe them here. 4. Facial Expressions When researching facial expressions, we noted that many theories are presented in some fields based on psychology. In this research, we adopted the theory of P. Ekman and W.V. Friesen [6]. P. Ekman and others have stated that humans have many expressions, and that the basic expressions are surprise, happiness, sadness, anger, disgust, and fear - six expressions in total - and that other expressions can be represented by combining these six basic expressions. Based on this theory, we thought that if the six basic expressions could be achieved, complex expressions could also be represented by combining them. 5. Robot Structure to Represent Facial Expressions As an interface for robot communications that would approximate the capability of humans, we aimed to develop a robot that would have a physical structure similar to that of humans. Our robot has an artificial skeleton made of aluminum covered with artificial skin made of polyethylene, which is relatively flexible and can be used to represent facial expressions in parallel with the Kansei system by controlling small servomotors mounted inside the artificial skeleton that pull the artificial skin into facial gestures associated with various expressions (Fig. 3). Fig. 3. Developed robot. Sensors & Transducers Journal, Vol. 6, Special Issue, August 2009, pp. 18-26 23 6. About the System We think that human consciousness exists as a circulating stream of "cognition," "representation" and "behavior," and we attempted to achieve a robot with this stream of consciousness. Cognition: Information is input in a computer using a keyboard or a microphone. Input using a microphone is recognized with a voice board and the entered keyword is sent to the PC to pass the information to the Kansei information processing section of the system responsible for decision processing. Representation: The Kansei information for a keyword retrieved from the Internet is processed. To determine the feeling to generate and represent the proper facial expression according to the feeling on the robot hardware, the decision result is sent to the robot device section responsible for representation processing. Behavior: Facial expressions and voice output according to the Kansei information obtained in the decision process are generated. Control signals are sent to the servomotors that pull the skin of the robot's face, and text to be spoken by the robot is sent to the voice board. By repeatedly executing this sequence of three processes, we are able to achieve a series of streams of "cognition," "representation," and "behavior" that are simple but the most characteristic points of this research (Fig. 4). Fig. 4. System overview. 7. Generation of Facial Expressions As described above, we use the model of six basic expressions proposed by P. Ekman, and we achieved a robot capable of expressing the six basic facial expressions of surprise, happiness, sadness, anger, disgust, and fear. In addition, we used the concept of the "synthesis of facial expressions" proposed by P. Ekman to generate complex facial expressions. An example of a complex expression is a bitter smile. According to P. Ekman, a person expresses a bitter smile to hide true feelings of happiness or anger. We can also consider it to be a synthesis of happiness and fear. As described here, Sensors & Transducers Journal, Vol. 6, Special Issue, August 2009, pp. 18-26 24 we thought that various human facial expressions could be achieved by combining the six basic expressions and we successfully created these facial expressions on the robot. (Fig. 5) Neutral Angry & Sadness Fig. 5. Complex facial expression. Here, we describe the process of generating complex facial expressions. As described above, a complex expression can be achieved by synthesizing the six basic expressions. To generate complex expressions on a robot, we must consider the linkage of words and the stream of consciousness, that is, the stream of the conversation and the story. Then, we defined the concept of an "accumulation of feeling." An accumulation of feeling means an accumulation of the Kansei value of each type of feeling in the transition of consciousness in the stream of the robot's consciousness. Like human memory, however, the feeling about an event that has just occurred is different from the feeling about an event that occurred one hour before or one day before. In other words, we think that the strength of the feeling declines to a certain extent as time elapses. This is called the "decline of feeling." By using the accumulation and decline of feeling, we can generate complex expressions on the robot and achieve representation of feelings over the long-term. Then, we defined the "impression value" and achieved the accumulation and decline of feeling by calculating the "impression Kansei Value" (hereafter referred to as IK value). We will describe the impression value and the impression Kansei value below. 8. Impression Value and Impression Kansei Value First, we describe the impression value. Every human has a proper "impression" about an object or an event. In many cases, the human decides on a thought or action regarding the object according to an "impression." Alternatively, the selection of an action related to the object is influenced by the "impression." In other words, we can think that the "impression" owned by an object (word) itself is an important element for generating consciousness and representing feelings. The numerical value of the "degree of impression owned by the object itself" in the subjective memory and experience of a person is defined as the "impression value." Sensors & Transducers Journal, Vol. 6, Special Issue, August 2009, pp. 18-26 25 By using this impression value, feelings can be accumulated and declined, and synthesized facial expressions can be generated. For this purpose, we define the "impression Kansei value" that holds simultaneously the "Kansei" and "impression" of an object by using the impression value that determines the impression of an object, in addition to the Kansei value that shows temporary feelings. By preparing the impression Kansei value and accumulating the Kansei for certain time in the stream of consciousness, we can represent not only momentary changes in feeling and consciousness but also short-term feelings with the time axis considered. Also we can perform formalization to represent a "synthesized complex feeling" combined from the six basic types of feeling and achieve a "change of consciousness in the stream of conversation." 9. Problems and Considerations To change the consciousness and represent the feelings of the robot in this research, we use a database that consists of association word groups and the Kansei information owned by each word obtained by our analysis of text data extracted from the Internet. When the Kansei information processing is performed for a central word of the current consciousness, its association word group and the Kansei evaluation factor are calculated, and then based on them the Kansei value is calculated. Also by shifting the center of consciousness to the association word with the highest association value, the Kansei information processing can be performed repeatedly. We achieve an artificial stream of consciousness based on the occurrence of words that are associated sequentially in this way. However, we aim finally to achieve a stream of consciousness as the context of a sentence as well as a chain of words. If a robot can be conscious of a long sentence as well as the words, and can perform Kansei information processing, the robot can achieve not only natural conversation but also a natural generation of facial expressions in the stream of the conversation. We think that this will lead to the creation of higher communication between humans and robots. At present, the robot used in our research can already represent complex facial expressions, and we aim to achieve natural conversation between humans and robots in the future. 10. Summary We have developed a robot that generates facial expressions according to a stream of artificial consciousness and the stream of feelings generated by the consciousness. In this paper, we explained Kansei and the consciousness network that is the stream of artificial consciousness and facial expressions. We also described the structure and generation of the facial expressions achieved by the robot in this research. Finally, we described the process of generating complex facial expressions. References [1]. Breazeal, C., Designing Sociable Robots (Intelligent Robots and Autonomous Agent), MIT Pr., 2002. [2]. Kobayashi, H. Hara, F., Facial Interaction between Animated 3D Face Robot and Human Beings, International Conference on System, Man and Cybernetics, 1997. Sensors & Transducers Journal, Vol. 6, Special Issue, August 2009, pp. 18-26 26 [3]. Imai, Y. Takeno, J., Method For Extracting Associated Words, Association Value and Sensitivity Value From the Internet, The 8th World Multi-Conference on Systemics, cybernetics and Informatics, The International Institute of Informatics and Systemics (IIIS), Proceeding Vol. I, 2004, pp. 160 -165. [4]. Yamanaka, M. Kurokawa, S. Yamano, A. Takagi, A. Takeno, J., Change in Consciousness and Facial Expression of Humanoid Robot, The 8th World Multi-Conference on Systemics, cybernetics and Informatics, The International Institute of Informatics and Systemics(IIIS), Proceeding Vol. I, 2004, pp. 183-188. [5]. Ogiso, A. Kurokawa, S. Yamanaka, M. Imai, Y. Takeno, J., Expression of Emotion in Robots Using a Flow of Artificial Consciousness, Computational Intelligence in Robotics and Automation, CIRA, 2005, 10, pp. 421-426. [6]. Ekman, P. Friesen, W. V., Unmasking the Face, Prentice- Hall Inc., 1975. ___________________ 2009 Copyright ©, International Frequency Sensor Association (IFSA). All rights reserved. (http://www.sensorsportal.com) Sensors & Transducers Journal, Vol. 6, Special Issue, August 2009, pp. 27-42 27 SSSeeennnsssooorrrsss &&& TTTrrraaannnsssddduuuccceeerrrsss ISSN 1726-5479 © 2009 by IFSA http://www.sensorsportal.com A Bayesian Approach to Solving the Non – Invasive Time Domain Reflectometry Inverse Problem Ian G. PLATT and Ian M. WOODHEAD Lincoln Ventures, PO Box 133, Lincoln, Christchurch 7640, New Zealand Tel.: +64 3 325 3748, fax: +64 3 325 3725 E-mail: platti@lvl.co.nz Received: 10 July 2009 /Accepted: 31July 2009 /Published: 10 August 2009 Abstract: Non Invasive Time Domain Reflectometry (TDR) may be used to estimate the volumetric moisture content, θv, with depth for a variety of sample materials. The forward physical model is couched in terms of a moment method where integration is performed over a descretized sample space to estimate the measured propagation time tp down a pair of parallel transmission lines. We show that the inverse solution to this, which recovers relative permittivity and thus θv, is greatly facilitated by a simplification of the system geometry via, 1) realistically modeling the prior density of the sample, 2) using this prior with the inherent system symmetry to reduce the number of required discretization cells, and 3) determining a physically meaningful reduction operator to allow a coarse discretization mesh to be used. The observational equation is expressed in the Bayesian paradigm with the most accurate and robust solution obtained using the Conditional Mean of the posterior distribution constructed via a Monte Carlo method. Results of simulation show that the method is capable of providing accurate estimates of the moisture density profile down to a depth of 100 mm with an error of less than 4 %. Further, the reduction in the number of descretized cells required to accurately estimate these profiles means that the inversion procedure is quick enough to enable the real time application of the equipment, a fundamental requirement in the development. Copyright © 2009 IFSA. Keywords: Time domain reflectometry, Bayesian, Non invasive, Moisture content, Transmission lines, Electromagnetic sensor 1. Introduction The ability to calculate moisture content in soils has been developed over the past 50 years, primarily to maximize water use efficiency in agriculture and to understand subsurface hydrological processes. Sensors & Transducers Journal, Vol. 6, Special Issue, August 2009, pp. 27-42 28 Three basic methods have been developed over this period, namely; gravimetric, neutron scattering and electromagnetic [1], [2]. Gravimetric methods are accurate but have the disadvantage that the sample must be removed, while neutron scattering sensors must be closely monitored and cannot be left unattended, since they pose a health risk. The role of electromagnetic sensors, which suffers neither of these disadvantages, has thus become of increasing importance for soil moisture content measurement. Most of the current electromagnetic systems estimate the relative permittivity, εr, of the sample by deploying two or more metallic probes into the soil to measure the capacitance or impedance between them, or as in the case Time Domain Reflectometry (TDR), the time delay of a pulse sent down a pair of probes acting as a transmission line. The εr of the sample is then related to its volumetric moisture content, θv by an empirical relationship such as those derived by [3] and [4]. Since water has a high relative permittivity, εr 80, compared to values of εr ≤ 5 for most other components that make up typical samples the pulses velocity change is most sensitive to the amount of water present. The measured permittivity can be related to the samples dielectric properties by a real part representing the stored energy of the system and an imaginary part describing its energy dissipation: , (1) where is the stored energy of the system, is the electrical conductivity between the probes, is the molecular relaxation of the water molecules and is the frequency of the probe signal. In this work we assume a lossless media so that the imaginary part of (1) becomes negligible and the effective permittivity for the system is . In those cases where this is not justifiable [5] have developed a method whereby allowance for the conductivity part of the loss can be easily incorporated. For the invasive TDR system a pair of parallel transmission lines are inserted into the sample and a short pulse is sent down them generating a static electric field. This field causes the sample dielectric to produce an opposing polarization field that reduces the phase velocity of the pulse. The velocity of the pulse can be given in terms of the permittivity and permeability of the sample by: , (2) where is the relative permittivity, c is the speed of light and is the relative permeability (which is equal to one for most soils). By measuring the time the pulse takes to travel down the transmission line of length, d, and back again after a reflection from the end, the propagation time can be related to from (2) (with ) by: (3) Woodhead [6] extended this TDR arrangement by deploying the transmission lines external to the sample, effectively creating a non-invasive system. Being non-invasive the system can be used on a greater variety of materials including timber, concrete, road seal etc, with the added advantage that it does not alter the medium (.e.g. cause water wells around inserted probes) being tested. For this new system geometry the relationship defined by (2) and (3) no longer applies. Using a tomographic technique based upon discretisation of the x- y plane through the system geometry (Fig. 1), a new equation relating the pulse propagation time, , to the sample permittivity, , can be constructed [6] of the form: Sensors & Transducers Journal, Vol. 6, Special Issue, August 2009, pp. 27-42 29 (4) This equation represents the forward construction of the physical system where the form of is specified in section 2. Since we seek an estimate of from a measured value of it is the inverse of this equation that is required: (5) In common with many other inverse problems, the solution of (5) is somewhat problematic. The forward function is both non-linear and ill-conditioned thus an estimate of by the inverse process requires some form of regularized optimization technique. Further, since no analytic solution to the forward model is available, the method of moments technique in which the sample is discretised, is used as an approximation to the solution. The level of discretisation, or mesh size, will of course have an effect upon the accuracy of this approximation. The objective of this paper is to describe a methodology for improving the accuracy and efficiency of the inverse defined by (5) via a Monte Carlo simulation. 2. Background The steps in defining the function f in the forward model of (4) are derived by Woodhead [6] which we briefly recount here. When an incident electric field, is imposed upon a dielectric of volume , the dipole moment P, of the volume produces an electric field, . The total electric field is then: (6) Using the constitutive relationship this becomes, after some re-arrangement: , (7) where is the permittivity of free space and is the relative permittivity of the dielectric. may be calculated at a point r from the polarization element P via: (8) So that (7) becomes: (9) and this may be expressed in terms of a linear operator, L, acting upon P, [10] as: (10) Sensors & Transducers Journal, Vol. 6, Special Issue, August 2009, pp. 27-42 30 Upon discretisation into N cells, as depicted in Fig. 1, the polarization field of (8) becomes: (11) where m is the cell at which the polarization field from cell n (source cell) is to be calculated. Since, by model assumption, the z-component is constant, expanding over the gradient function gives: , (12) where and are the x and y differences between the points m and n respectively. Note that when m = n (i.e. the cells own polarization field) components become: (13) By using the method of moments with a subsectional basis function defined as non zero only for the source cell, Harrington [10] showed that (10) can be written in the discrete matrix form: (14) with the components of L given by the square bracketed part of (12) for and by (15) when , where . Since both L and P are functions of the independent variable, (14) is non-linear. Woodhead [6] extended the discretised model of (14) by introducing the concept of a transmission line to produce the incident field , so that the pulse velocity, v, due to the dielectric interaction can be given by the telegraphers equation (e.g. [11]): (16) so that , (17) Sensors & Transducers Journal, Vol. 6, Special Issue, August 2009, pp. 27-42 31 where is the total voltage between the transmission lines, is the component due to the polarization field induced by the dielectric sample and is the component due to the incident electric field directly from the transmission lines; d is the length of the transmission line; is the propagation time of the pulse from the beginning of the line to its reflection at the end and back to the beginning; ρ is the line charge density; µ is the magnetic permeability; b is the spacing between the parallel lines, and a is the diameter of the transmission lines. The only unknown on the right hand side of (17) is the polarization induced voltage, , contained in the total voltage Equation (17) may thus be succinctly represented by: (18) Following the method of [7] both the transmission lines and sample can be imbedded in the discretisation space, depicted by Fig. 1, so that: , (19) where M is the number of cells between the transmission lines and w is the width of the cell. Fig. 1. The 2D geometry of discrete cells and parallel transmission lines used to define the system matrix equation. The cells between the transmission lines are those (M of them) included in the calculation of . Each of the cells has a width w. For any particular value of permittivity in each sample cell, the polarization P within it can be calculated from (14) since so for those cells: , (20) where the delta function indicates that only components of cells between the transmission lines are selected. We represent this in the simplifying function notation as: Sensors & Transducers Journal, Vol. 6, Special Issue, August 2009, pp. 27-42 32 (21) The forward solution defined by (4) can now be formed by the composition of the two functions Ψ and Φ so that: (22) 3. Bayesian Formulation In exploring the solution to the inverse problem outlined by (4) we use the Bayesian framework as described for example by [8] and [9]. In this representation (4) may be written in the observational form as: , (23) where now T and F represent the random variables (with and as realizations) and we have added as the measurement noise, assumed to be independent of F. An inverse solution to this equation can be given in terms of the variable densities by the well know Bayes Formula: , (24) where is the posterior density, is the likelihood density function and is the prior density function of the permittivity and will act as a regularization for different physical assumptions about the likely structural disposition of throughout the sample. The term , often called the evidence factor, can be viewed as a scaling factor that ensures a unit total probability of the posterior density function. Since it plays no part in the evaluation of the required parameters of the posterior distribution it is usually ignored so that the above equation becomes: (25) To obtain an estimate of corresponding to the measurement we can use either the Maximum A Posteriori Estimator (MAP) or Conditional Mean (CM) of the posterior distribution given by (25), with associated errors given by the covariance . Both the CM and estimates require integration over the posterior density : (26) with the estimate error given by the posterior covariance (27) while the MAP estimate: (28) Sensors & Transducers Journal, Vol. 6, Special Issue, August 2009, pp. 27-42 33 is usually achieved by some sort of optimization method. The construction of the applicable form of (25) for the TDR model will be detailed in section 3.4 after some relevant properties of the system (such as parameter independence) are established. 3.1. Required Level of Discretisation In the observational model of (23), T is the measured value of the propagation time, with realization . Since no measurements are available, model validation must be carried out using simulated data which should be generated at a high enough level of resolution that it is as close as possible to the perceived real case. For this case we content ourselves with a level of discretisation that provides a model accuracy that cannot be much improved upon by further decreases in mesh size. To determine this level of discretisation, the polarization electric field, , between the transmission lines was calculated using (22) for an increasing number of individual cells over a constant sample volume. The results of this are shown in Fig. 2. Here the polarization field converges to that given by the mesh sizes of Res 4 (defined in Table 1). Table 1. Discretisation values. nx, ny, nz, are the number of cells in the directions x, y and z respectively. Resolution Number Cell Size (m) nx ny nz Res 1 0.100 4 1 1 Res 2 0.033 12 3 3 Res 3 0.020 20 5 5 Res 4 0.014 28 7 7 Res 5 0.011 36 9 9 Fig. 2. The polarization field distribution over the half plane between the transmission lines. Each of the curves represents a different Resolution, defined by Table 1, for a constant volume sample. Since Res 5 is computationally inaccessible, even for the small number of test cases required in the current development, Res 4 is chosen as being close enough to the continuous case to well approximate the operational system. Computationally it is only accessible for a limited number of Sensors & Transducers Journal, Vol. 6, Special Issue, August 2009, pp. 27-42 34 calculations and will be used in section 3.4 to construct a resolution reduction operator q and in section 4 to test the inversion procedure. 3.2. Operational Discretisation In the previous section we established that the discretisation level, labeled as Res 4 of Fig. 2 and Table 1, provided the closest accessible level of mesh size that could represent the true environment. Computation at this level of discretisation is very slow however and certainly unacceptable for the real time applications required for the TDR. Even the pre-calculation of much of the Monte Carlo sampling described in [12] would be impractical at this level. It would therefore be advantageous to perform the calculations in the coarser mesh of Res 2, and then use a map Q to project the results into those that would have resulted from the finer mesh: , (29) where the high resolution (Res 4) and low resolution (Res 2) are represented by the subscripts h and l respectively. The observational equation (23) for the high resolution, Res 4, can be written in this terminology as: , (30) where is the propagation time and is the measurement noise. Applying the operator, Q, this becomes: (31) with (32) being the error introduced by using the coarser discretisation. In this sense the operator Q can be termed as the model reduction operator [8]. Putting the total error as: ; (33) the equation (31) can be written as: , (34) where we have made the assumption throughout this formulation that the discretisation error is approximately independent of the relative permittivity, . This is discussed in more detail in section 3.4. With the propagation time measurement and the form of the likelihood function becomes: , (35) where the subscript emphasizes that the likelihood distribution is of the same form as the distribution Sensors & Transducers Journal, Vol. 6, Special Issue, August 2009, pp. 27-42 35 given by the combined error of (33). If the distribution is close to Gaussian, which is often the case, [9] shows that the above equation becomes: , (36) where and are the expectation and covariance respectively of the total error distribution, . Here the mutual dependence of and the approximation error have been ignored, a simplification that will be addressed further in the next section. The values of can be pre-calculated because of the unique system geometry and operational constraints [12], so once the operator q is known the calculation speed of (35) is acceptable for real time application. The required posterior density is thus: (37) 3.3. Prior Densities As discussed in [12] for many cases of practical interest the moisture within the sample can be considered as being stratified in the x -z plane. Apart from reducing the number of unknown permittivity values to be calculated, plane stratification also leads to symmetry in the x - z plane that can be utilized to decrease the number of dimensions of the inverse system. In keeping with this stratified geometry of the sample we define three prior distributions for , each representative of a number of real world scenarios. In the first case, the prior is given as a positive increase in relative permittivity with depth. With the discretisation level of Res 4 it can be defined as: , (38) where , is the uniform distribution and . The second prior is similarly defined but with a decrease in permittivity with depth, so that the incremental permittivity between layers is defined by . For both of these distributions the maximum change of over the sample depth is 14. A third prior is included as a subset of the first, that of a positive linear gradient over the sample. In this case the value of the upper boundary of the sample is chosen from the distribution , while the lower sample boundary value is chosen from . for the centre of each depth layer is then obtained by linear interpolation. Each of these will be used in the examples that follow. 3.4. Mapping Operator After defining the form of the next immediate concern is to find a model reduction operator Q that maximizes the decrease of in (32) and thus in (37). For such an operator to be useful it must be only weakly dependant on the actual sample relative permittivity . The Fig. 3 is a plot of the propagation times and when the sample is discretised to the level of Res 2 and Res 4 respectively, for 200 different simulations. The values for each of the layer permittivities of Res 4 are generated by a prior distribution, defined in section 3.3. Res 2 values for each case are calculated at the cell centre points of the three layers by Lagrange Sensors & Transducers Journal, Vol. 6, Special Issue, August 2009, pp. 27-42 36 interpolation over the midpoints of the seven Res 4 cells. The Fig. 3 shows the relationship between the two propagation times for the transmission line height above the sample of . A good fit to the data for each height in this figure is achieved by the polynomial: (39) where a, b and c are chosen for each of the different heights. In terms of the Bayesian paradigm of (37) then: (40) Values of the coefficients for this polynomial are easily calculated by any of the standard fitting routines and are applicable to all measurements with the same geometrical relationships and prior relative permittivity distributions. The form of (40) is also a valid approximation for the reduction operator when the other prior distributions defined in section 3.3 are used, though the coefficients at similar heights may differ. Fig. 3. The form of the projection operator Q between the low resolution and high resolution discretisation models for a transmission line heights = 0.014 m above the sample. and are the propagation times for Res 2 and Res 4 respectively. Other heights above the sample give a similar relationship. Table 2 shows the standard deviation, , of data difference from the fitted polynomial for the different priors discussed above. These values can be used for defining the distribution in (33) and together with the measurement noise , constitute the total error, from which the covariance, of (37) is extracted. The mean error in all fitted cases is zero and since we have . Further, the low dependence of the discretisation error on shown explicitly in Fig. 3, strengthens the ability of the enhanced error model of (37) to represent the real situation. Table 2. The standard deviation, , of data fitted to (40) for the indicated transmission line heights. Given in ps ( ). 30.0 13.7 7.6 27.8 10.5 5.6 linear +ve 12.5 6.7 3.2 Sensors & Transducers Journal, Vol. 6, Special Issue, August 2009, pp. 27-42 37 Using from Table 2 and from previous experimental work [12], the Gaussian distributions and can be defined for (33). The errors are dependent upon the current transmission line height only, so: (41) is a diagonal matrix containing the total error variance for each height. Re-writing (37) with and the covariance gives the form of the posterior density used in subsequent simulations as: (42) Thus the posterior distribution of can be generated with the computationally accessible Res 2. From this distribution an estimate of for each layer can made from the CM of (26) with its associated error estimate given by the covariance defined by (27). Table 3. System parameters Parameter Value (m) Sample Depth 0.1 Transmission Line Length d = 1 Diameter a = 0.01 Separation b = 0.3 Heights (0.006, 0.001, 0.014) 3.5. Some Practical Simplification for Computation The posterior distribution may be generated by the Monte Carlo approach described by [12] with an estimate for its value given by the CM. Although in this scheme much of the likelihood function is pre-calculated using (22), the large number of calculations required together with possible real time refinements results in a heavy computational burden for the required level of discretisation. To reduce the computational effort to something manageable the previous sections looked at reducing the number of dimensions of the observational form. Another purely physical reduction in complexity can be achieved by simplifying the geometry represented in Fig. 1. The geometry depicted in Fig. 1 indicates that a potentially large number of cells represent air, with a known permittivity of . These cells, apart from those between the transmission lines, provide no extra information and can be discarded. Removing these cells from the model the polarization of the sample cells may now be related to those cells between the transmission lines by a gradient function similar to that given by the bracketed terms of (12) ([12]). If G is this gradient function, then the polarization in each sample cell ( ) may be related to the field between the transmission lines ( ) by: (43) applying it to (14) after some rearrangement gives: Sensors & Transducers Journal, Vol. 6, Special Issue, August 2009, pp. 27-42 38 , (44) where and L now refer only to the sample cells. As discussed in [12] for many cases of practical interest the moisture within the sample can be considered as being stratified in the x -z plane. Apart from reducing the number of unknown permittivity values to be calculated, plane stratification can again also lead to a further simplification in the specification of the forward model. Fig. 4 represents a sample material with N cells and consequently N unknown values of permittivity. By assuming plane stratification the number of unknown permittivities is reduced to the number of stratum (5 in the example of Fig. 4). In theory if there are N unknown values of permittivity there must be N independent measurements of the propagation time, . For example, for five stratums each having an unknown value of , five measurements of are required to specify the system of equations represented by Error! Reference source not found.. To generate an independent set of components, each measurement of is made at a different height of the transmission lines above the sample. From this argument it is clear that the planar stratification also decreases the number of required measurements. Fig. 4. Physical geometry of the TDR system. Because of the form of the gradient field, it is essential that the centers of the cells between the lines are aligned with those of the sample. Another important simplification in the geometry can be made by considering the symmetry of the system. The component values of the entire samples polarization and their effect on between the transmission lines can be represented by a reflection of the shaded quarter space shown in Fig. 4. Hence, by avoiding duplication of symmetric cells, it is possible to further reduce the components of L, G and by a factor of four. 4. Results Having formulated a solution to the operational objectives outlined in the Introduction it remains to show how well it performs. Since this is a simulation, it is sensible to test the operational proc