第10回CEReS
環境リモートセンシングシンポジウム
資料集
Proceedings of the10th Symposium on
Remote Sensing for Environment
180°

240°

300°

0°

60°

120°

180°

130

140

60°

0°

-60°

70

80

90

100

110

120

Mean PAR (W/m2) averaged during growing season in 1986 - 1995; Presented by R. Suzuki

180°

240°

300°

0°

60°

120°

180°

60°

0°

-60°

-0.1

0.0

0.1

0.2

0.3

0.4

0.5

0.6

0.8

Mean NDVI in summer 5 months in 1986 - 1995 (CEReS NDVI data) ; Presented by R. Suzuki

2008年2月28日
千葉大学けやき会館
千葉大学 環境リモートセンシング研究センター
Center for Environmental Remote Sensing (CEReS)
Chiba University

第 １０ 回CEReS
環境リモートセンシングシンポジウム
資料集
Proceedings of the 10th Symposium on
Remote Sensing for Environment

2008年2月28日
千葉大学けやき会館

千葉大学 環境リモートセンシング研究センター
Center for Environmental Remote Sensing (CEReS)
Chiba University

共同研究報告書 平成 19 年度

は じ め

に

平成１９年度の環境リモートセンシング研究センター・シンポジウムを、平成２０年２月２８日に開催しました。本シンポジウムは毎年開催し、共同研究成果の発表と評価を行い、この機会を利用して今後の研究方針の議論と策定を行う事をおもな目的としております。平成１９年度の共同利用研
究成果発表では、口頭発表１８件、ポスター発表２３件の申し込みがありました。
平成１９年度の共同利用研究は前年度に引き続いて、プロジェクト研究と一般研究から構成されています。プロジェクト研究は、（１） 「衛星データによる地球表層環境変動の実態把握とその要因解析」、（２） 「衛星データによるユーラシア大陸の植生３次元構造の変遷を中心とする表層・植生・
土地被覆変動の研究とデータ解析・処理手法・検証データ観測手法の研究」、（３）「衛星データと地上観測ネットワークによる放射収支の評価と大気パラメータの長期変動」、（４） 「地域社会に役立つリモートセンシングの実現―多様な空間情報のシナジーによる社会基盤情報の発信―」、の主要な四つの課題、及び（５）「未来を切りひらく新しいリモートセンシングの展開」の新課題を加えて、本年度も活発な共同研究を通じて成果を上げて頂きました。
本冊子は、多くの方に成果をご覧いただくため、発表会での講演・ポスターの内容を提出していた
だき、資料集として編集したものです。この資料集の内容は、また、環境リモートセンシング研究セ
ンターのホームページにも掲載しております。
本共同研究において、衛星データや幅広い環境に関するデータを活用したリモートセンシングの研究活動を行い、成果発表を通じてご協力を頂いた研究者の皆様に感謝いたします。
今後、個々の研究が大いに発展し、センターの共同研究がますます活性化することを確信しております。
さらに、学術研究の推進体制に関連して、皆様のご協力を共同利用･共同研究の推進へと活用させていただきたいと考えております。

平成２０年３月３１日

環境リモートセンシング研究センター
センター長
西尾 文彦

第 10 回 CEReS 環境リモートセンシングシンポジウム資料集（2007 年度 CEReS 共同利用研究発表会資料集）

目次[プロジェクト-１]

1

全球 NDVI の降水量，気温，放射量の経年変化に対するレスポンス ------------------

2鈴木力英 ：(独)海洋研究開発機構地球環境フロンティア研究センター

モンゴル高原とその周辺の土壌水分・植生の時空間変動 -------------------------１

２

１

３

開發一郎 、近藤昭彦 、秋山量平 、小池俊雄 、太田 哲

5３

１．広島大学大学院総合科学研究科 ２．千葉大学 CEReS ３．東京大学工学系研究科

アジアの植生・作付分類と光合成（CO2 固定）速度平面分布の推定 ------------------

6

金子大二郎 ：松江工業高等専門学校 環境・建設工学科

地球地図・土地被覆データ整備における都市域の抽出 ---------------------------１

１

１

１

１

１２

２

梶川昌三 、筒井俊洋 、山田美隆 、阿久津修 、砂川英征 、建石隆太郎
１．国土地理院地理調査部環境地理課 ２．千葉大学 CEReS

多偏波 SAR データを用いた海氷物理量測定 -----------------------------------１

２

若林裕之 、中村和樹 、西尾文彦

１８３

１．日本大学工学部情報工学科 ２．(独)産業技術総合研究所 ３．千葉大学 CEReS

西部太平洋熱帯域と東アジアにおける気候環境の関連性に関する調査 -------------１

２１

２

菅野洋光 、Prima O.D.A

１．東北農業研究センター ２．公立大学法人岩手県立大学 ソフトウエア情報学部

Study of Coral Reef Distribution in Bali island using Satellite ｄata ------------

２２

大澤高浩、 I Ketur Suwardika ：ウダヤナ大学海洋リモートセンシング研究センター

ミクスチャモデルによる MODIS データの雲量推定 --------------------------------

２３朴 鍾杰 ：東京情報大学 環境情報学科

高分解能合成開口レーダのクロス偏波画像に見られるテクスチャを利用した針葉樹林バイオマス
計測------------------------------------------------------------------１

１

２

大内和夫 、岩切宗利 、王 海鵬

１．防衛大学校 電気情報学群情報工学科 ２．復旦大学

２４

夜間光をパラメータとする二酸化炭素排出量の推定−東アジア地域における二酸化炭素の排出量の変動について―------------------------------------------------

２６

原 政直 ：(株)ビジョンテック

[プロジェクト-２]

２７

衛星データを用いた植生純一次生産量推定の精度向上のための基礎研究 II --------

２８

村松加奈子、陳路 ：奈良女子大学共生科学研究センター

NOAA/AVHRR GAC を用いた全球時系列データセットの作成 --------------------１

２

松岡真如 、本多嘉明 、梶原康司

３２２

１．高知大学農学部農学科国際支援学コース

２．千葉大学 CEReS

カラマツ林における方向別反射特性の季節変化と日変化〜光合成モデルに必要な樹冠
構造パラメータ推定に向けて〜----------------------------------------------

３４

松山 洋、長谷川宏一 ：首都大学東京 都市環境科学研究科[プロジェクト-３]

３５

衛星画像を用いたオホーツク地域のエアロゾル光学特性の導出 -------------------１

２

朝隈康司 、山﨑輝彬 、久世宏明

３６３

１．東京農業大学生物産業学部アクアバイオ学科 ２．千葉大学大学院融合科学研究科 ３．千葉大学 CEReS

多波長マイクロ波放射計データを用いた水物質リトリーバルの研究：マイクロ波放射計
データの非静力雲解像モデルへの同化法の開発-------------------------------

４０-------------------------------------

４６

青梨和正 ：気象庁気象研究所予報研究部第２研究室

船舶搭載用新型スカイラジオメータの開発
１

１

２

２

２

２

３

４

塩原匡貴 、矢吹正教 、笹本和敏 、阿部泰三 、江井和則 、河井和弘 、内山明博 、中島映至 、村山利幸５、青木一真６、小林 拓７、村治能孝８、山野牧８
１． 国立極地研究所 ２．プリード ３．気象研 ４．東大 ５．東京海洋大 ６．富山大 ７．山梨大
８．エナジシェアリング

太陽放射観測による父島のエアロゾルの光学的特性 ----------------------------１

２

２

３

４

青木一真 、塩原匡貴 、矢吹正教 、村山利幸 、兼保直樹

１．富山大学大学院理工学研究部 ２．極地研究所 ３．東京海洋大 ４．(独)産業技術総合研究所

４９

北西太平洋域における下層雲の衛星リモートセンシ ングの航空機検証観測---------

５１

浅野正二、佐藤茉莉 ：東北大学大学院理学研究科 附属大気海洋変動観測研究センター

2007 年春季黄砂と視程映像観測 --------------------------------------------１

２

３

５４

４

木下紀正 、永松哲郎 、飯野直子 、足立祐樹

１．鹿児島大学産学官連携推進機構 ２．鹿児島大 水産学部 ３．鹿児島大 工学部 ４．熊本大

ミリ波雲レーダ FALCON による観測-雲物理量導出と最近の結果 -----------------鷹野敏明

１・２

高村民雄

２・３

１

１

１

１

１

１

５８

１

、河村洋平 、阿部英二 、山口 潤 、二葉健一 、横手慎一 、鈴木陽平 、渡辺潤一 、、熊谷 博４、大野裕一４、中西裕治１・５

１．千葉大学 大学院工学研究科 ２．千葉大学 VBL ３．千葉大学 CEReS ４．情報通信研究機構 ５．サイテック

Ashra 観測の現状報告 -----------------------------------------------------

６１

佐々木真人 ：東京大学宇宙線研究所

Temporal and Spatial Variation of Cloud Measured with a Portable Automated Lidar-------------------------------------------------------------------------１

２

２

２

６２

３

椎名達雄 、竹内延夫 、眞子直弘 、久世宏明 、内藤季和

１．千葉大学大学院融合科学研究科 ２．千葉大学 CEReS ３．千葉県環境研究センター[プロジェクト-4]

６７

Analysis on the relationship between urban heat islands and urban
development in Taipei by satellite images (First Report) --------------------１

２

６８

３

白 迎玖 、近藤昭彦 、原田一平

１．東北公益文科大学公益学部 ２．千葉大学 CEReS ３．(独)国立環境研究所社会環境システム研究領域

熱環境緩和対策の立案におけるリモートセンシングデータの活用に関する検討-------

６９

ＲＳ・ＧＩＳによる沿岸生態系主要景観の広域分布動態の解析 ----------------------

７０

岡田信行：ALT 都市環境研究所

１

２

３

４

５

５

６

仲岡雅裕 、渡辺健太郎 、山北剛久 、近藤昭彦 、小倉久子 、平間幸雄 、石井光廣

１． 千葉大大学院理学研究科 ２．千葉大大学院自然科学研究科 ３．千葉大大学院理学研究科
４．千葉大学 CEReS ５．千葉県環境研究センター ６．千葉県水産総合研究センター

筑波山頂における自動気象観測の復活
１

１

１--------------------------------------２３

上野健一 、林 陽生 、辻村真貴 、浅沼 順 、寄崎哲弘

７１

１．筑波大学生命環境科学研究科 ２．筑波大 陸域環境 ３．（有）クリマテック

IT を活用した低コスト生産のための効率的てん菜栽培手法の確立(1)効率的な収穫体系確立のためのてん菜収量推定方法の検討-----------------

７２

丹羽勝久、横堀 潤 ：株式会社ズコーシャ

Measurement of NO2 and aerosol in the atmospheric pollution using differential
optical absorption spectroscopy (DOAS) with a PC projector light source
１

２

２

３

３

４

------

７３

１

原田一平 、宮﨑正志 、片岡大祐 、久世宏明 、近藤昭彦 、浜田 崇 、一ノ瀬俊明

１．(独)国立環境研究所社会環境システム研究領域 ２．千葉大学大学院融合科学研究科
３．千葉大学 CEReS ４．長野県環境保全研究所[プロジェクト-5]

７５

衛星を用いた地殻変動監視に関する研究

火山活動のモニタリング インドネシア・メラピ火山-------------------------------------------------------------------------１

１

１

２

服部克巳 、野口恭史 、大野 望 、前田 崇 、高野 忠

７６２

１．千葉大学大学院理学研究科 ２．JAXA[一般研究]

７９

合成開口レーダーによる洋上風力発電ポテンシャルの評価-----------------------

８０

香西克俊、大澤輝夫 ：神戸大学海事科学研究科

SAR offset displacement from Sumatra-Andaman Earthquake Mw 9.2 on December 26,
2004 --------------------------------------------------------------------

８４

Ashar Muda Lubis、伊勢崎修弘 ：千葉大学大学院理学研究科

Trial of Digital Filter Photography for Alteration Mineral Detection in the HachimantaiArea, NE Japan
１--------------------------------------------------------２

３

３

３

９３

４

Myint Soe 、 建石隆太郎 、石山大三 、Krit Won-In 、高島 勲 、 Punya Charusiri
１． 秋田大学大学院工学資源学研究科 ２．千葉大学 CEReS

３．秋田大学 工学資源学部付属環境資源学研究センター ４．Chulalongkorn University

Retrieval of drop size distribution using both PR/TRMM and CPR/CloudSat data
S.K. Dash、田中 佐、上田太郎 ：山口大学理工学部研究科

９８

静止衛星データの全球合成技術の開発 ---------------------------------------

９９

操野年之、太原芳彦、奥山 新、橋本 徹、加藤浩司 ：気象庁気象衛星センターデータ処理部

泥炭分布地域における生物地球化学的環境解析へのリモートセンシング技術の利用可能性の検討（予報）１------------------------------------------------------

１

１０２

２

伊豫部勉 、原口 昭 、西尾文彦

１．北九州市立大学国際環境工学部 ２．千葉大学 CEReS

干渉合成開口レーダ(InSAR)による南極の氷河・氷床マッピング―等速変動場の仮定を用いない流動位相の導出―--------------------------

１０３＝飛行場・滑走路における積雪・凍結路面の検出＝ ---------------------------

１０４

１

１

２

河田哲郎 、木村 宏 、西尾文彦

１．岐阜大学工学部電気電子工学科 ２．千葉大学 CEReS

合成開口レーダによる道路凍結のモニタリング
１

２

長 康平 、Josaphat Tetuko Sri Sumantyo

１．株式会社ウェザーニューズ 道路気象コンテンツサービスグループ ２．千葉大学 CEReS

Monitoring Kashiwa City Using Remote Sensing ---------------------------１

２

１０５

３

大山克己 、Luhur Bayuaji 、 Josaphat Tetuko Sri Sumantyo

１． 千葉大学環境健康都市園芸フィールドセンター ２．千葉大学大学院自然科学研究科
３．千葉大学 CEReS

都市域におけるアルベードの計測 --------------------------------------------

１０６

菅原広史 ：防衛大学校 地球海洋学科

ライダーによる粒子径分布計測のための多重散乱モデル -------------------------

１０７

吉田弘樹 ：岐阜大学工学部

高空間分解能衛星データを用いた森林立ち枯れ領域判定法の開発 ---------------１

１

１１０

２

松島 大 、石戸大介 、樋口篤志

１．千葉工業大学工学部建築都市環境学科 ２．千葉大学 CEReS

第１０回ＣＥＲｅＳ環境リモートセンシングシンポジウムプログラム--------------------

１１１

プロジェクト−１

全球 NDVI の降水量、気温、放射量の経年変化に対するレスポンス

Boreal forest in Alaska

Arctictundra
Arctictundrain in
Siberia
Siberia

CEReS平成19年度共同利用研究発表会 (2008年2月28日)

New

Mongoliasteppe
Mongolia
steppeCEReS平成19年度共同利用研究発表会 (2008年2月28日)

植生指数(NDVI)データ
Twenty-year global 4-minute AVHRR NDVI dataset(provided by CEReS)

POTENTIAL

このデータはPathfinder AVHRR Land (PAL) global 10-day
composite 8 km NDVI dataを編集して作成された．PALの10dailyデータより各月について最大NDVIを選んで月別値としている．

Nemani et al. (2003)は植生のデータは使用せず，気象データのみを使用，それぞれをスケーリングすることによって植生への影響力を推定している．

CEReS平成19年度共同利用研究発表会 (2008年2月28日)

CEReS平成19年度共同利用研究発表会 (2008年2月28日)

気象データ
10年間(1986〜1995年)の気温，降水量，PARの月別値を International
Satellite Land Surface Climatology Project Initiative II (ISLSCP
Initiative-II)のデータセットから得た．

気温と降水量：CRU (Climate Research Unit, University of East Anglia)によって作成された，全球1 x 1度グリッドにおける月別値．

光合成有効放射量(PAR： 0.4〜0.7μm):
ISCCP (International Satellite Cloud Project)によって作成された，全球1 x 1度における月別値．
CEReSのデータに基づく暖候季5か月における平均NDVI(1986 〜1995年)の分布
CEReS平成19年度共同利用研究発表会 (2008年2月28日)

CEReS平成19年度共同利用研究発表会 (2008年2月28日)

暖候季5か月における平均気温(1986 〜1995年)の分布

CEReS平成19年度共同利用研究発表会 (2008年2月28日)

暖候季5か月における平均降水量(1986 〜1995年)の分布

CEReS平成19年度共同利用研究発表会 (2008年2月28日)

解析の方法
1986〜1995(10年間)において，以下のパラメータ間の経年変化の相関係数を全球陸上の各2x2度グリッドで計算した．
NDVI vs. 気温
NDVI vs. 降水量
NDVI vs. 光合成有効放射量(PAR)
もともと1x1度グリッドである値を2x2度にしたのは，相関係数の計算のためのサンプル数を40個(10年x4グリッドセル)にして有意水準を上げるた
め．すべてのパラメータは暖候季の5か月間(北半球では4〜8月，南半球
では10〜2月)の平均を使用した．
暖候季5か月における平均PAR (1986 〜1995年)の分布

CEReS平成19年度共同利用研究発表会 (2008年2月28日)

CEReS平成19年度共同利用研究発表会 (2008年2月28日)

上： 気温とNDVIとの相関係数の分布．中：降水量とNDVIとの相関係数の分布．下：PARとNDVIとの相関係数の分布．

1986〜1995年における気温とNDVI(マゼンタスケール)，降水量とNDVI(シアンスケール)，PARとNDVI(黄色スケール)との経年変化の相関係数(正の部分のみ)の分布．黒い部分は計算の対象外であることを意味する．白色はすべての組み合わせに関して相関が無かったことを意味する．
CEReS平成19年度共同利用研究発表会 (2008年2月28日)

CEReS平成19年度共同利用研究発表会 (2008年2月28日)

POTENTIAL
Nemani et al. (2003)による結果との比較．Nemani et al. (2003)は植生のデータは使用せず，気象データのみを使用，それぞれをスケーリングすることによって植生への影響力を推定している．

ACTUAL

相関係数の分布とOlsonによる土地被覆分類との比較．
CEReS平成19年度共同利用研究発表会 (2008年2月28日)

CEReS平成19年度共同利用研究発表会 (2008年2月28日)

まとめ
1986〜1995年における全球のNDVIと気温，NDVIと降水量， NDVIとPAR との経年変化の共変動性を調べた．
ツンドラや亜寒帯林ではNDVIと気温との共変動性が強いことがわかった．その他の地域では，一般にNDVIと降水量との共変動性が強いことがわかった．
ユーラシアや北アメリカ中緯度のいくつかの地域，南米のサバンナ，中国南部，オーストラリアでNDVIとPARとの強い共変動性が認められた．
こういった地域は概して雲量が大きな地域にあたる．
Nemani et al．(2003)の結果ともおおよそ一致していた．
以上の結果からは，植生に対する気候因子の影響力に対する基本的な関係を理解することができる．
しかし，オーストラリアでNDVIとPARとの相関係数が大きいなど，不自然な点もあり，今後解決されることが望まれる．

CEReS平成19年度共同利用研究発表会 (2008年2月28日)

ACTUALPOTENTIAL
Map of weighted climatic controls on net primary productivity determined from water
availability, average temperature, and cloudiness. Each data point represents 3
values of the membership functions based on annual mean temperature, waterbalance coefficient, and percentage of sunshine hours per year in one 0.5 x 0.5
degree grid cell. After Churkina and Running (1998)
CEReS平成19年度共同利用研究発表会 (2008年2月28日)

左図： 10年(1986〜1995)
平均の全球1度グリッド上のNDVI，温量指数，湿潤指数との関係．
右図：差図から計算されたNDVIの温量指数に対する傾き(赤)と湿潤指数に対する傾き(青)(Suzuki et al., 2006より)

過去の研究

過去の研究

NDVIの温量指数に対する傾き(赤)と湿潤指数に対する傾き(青)の分布．赤の地域の植生は温量，青の地域の植生は降水量に依存していることを意味する．(Suzuki et al., 2006)
CEReS平成19年度共同利用研究発表会 (2008年2月28日)

Red: Temperature dominant,
Blue: Aridity dominant,Green: Both parameters dominant

1982〜2000年(19年間)における植生指数と気温(T)，植生指数と降水量(P)との間の経年変化の相関係数．(a) PALの植生指数を使用，(b)GIMMSの植生指数を使用．(Suzuki
et al. 2007より)

CEReS平成19年度共同利用研究発表会 (2008年2月28日)

モンゴル高原とその周辺の土壌水分・植生の時空間変動

アジアの植生・作付分類と光合成（CO2 固定）速度平面分布の推定
H19 ＣEReS研究発表会：

研究の背景

アジアの植生・作付分類と光合成（CO2固定）
速度平面分布の推定金子大二郎（松江高専

１．気候変動に関する政府間パネルIPCC報告による脆弱性の指摘：
気象災害と農業生産に大きな影響が及ぶ．
２． CO2の吸収量と炭素循環の解明のためにはCO2固定量を推定する必要３．穀物生産量を推定するためには，光合成速度のほかに作付面積が必要．
・Terra衛星のMODISによって計測され大略の土地利用分類

環境・建設工学科）がされた分類データをもとに，植生を細分類する方法を検討．
・分類対象は森林と作物であり，地表水分指標LSWIと植生指標NDVIの植生形態学によって作物種を細分類．
４．Asiaの植生域の光合成速度（ CO2固定量）分布を推定
CO2単位型の光合成・単位収量モデル：
研究全体流れ図：光合成速度（CO2固定速度），
単位収量モデル，
土地被覆分類，気孔開度（CWSI)

研究順序の方針

Crop Production Research Flow

穀物収量モデルの開発と同時に，

Modeling on Crop IndicesCrop Information１）光合成速度PSNの平面分布の推定
２）単位はCO2固定量 ・・・炭素循環，CO2吸収

Meteorological Cond.
Radiation
Temperature
Cloud
Precipitation

木本でなく 草本バイオマスエネルギー

Satellite Data

Crop statistics
Yield
rice,wheat
Maize
Land use
GrowthSterility
InjurySolarradiation（実蒸発散量）/（Penman蒸発位）

CWSI

s

Eac
Ep

Ulan Baa tor

Photosynthesis
rate PSN

Monitoring of Crop Production
Crop yield
Early warning of poor harvest

C hoybalsan

N anjing
C hangsha

A lla habad
C ulcutta

Phitsa nulok
Na kho n Saw an
M adrasB angkokBangalore

Food demandand security

Crisis management

Iw am izaw a

O hoga ta
Furuka w a
PyongyangSowl
Shijiazhuang
Ta ejonT aeguJinan
K w ang ju

N agpurLand use /Land cover

Monitoring

LandcoverAll monitoring sites in SE̲Asia

Shenyang

N ew D elhif PAR

Crop Production Index
Integrated photosynthesis
i NDVI
GDD

同時に，旱魃指標BMDIと水分指標NDWIによる

Main monitoring sites in SE_Asia

ASTER
NDVI, SAVI

Photosynthesis rate３）熱収支法による水ストレス（CWSI)と気孔開度の評価法．

NOAA
MODIS

Biomass
LAI
Water stress
Spectrum
WSI

U don Thani
K hon Kaen
N akhon Ra tcha thani

JAPAN

World
English Name Latitude Longitude weather
number

Iwamizawa
Shinshinotsu
Ohogata
Furukawa
Nagaoka
Kuki
Izumo
MatsueCHINA
Shenyang
Beijing
Jinan
Shijiazhuang
Zhengzhou
Xian
NanjingShanghai
Wuhan
ChangshaNanchang
Chengdu
Yichang
North Korea Anju
Pyongyang
Sariwon
Kaesong
MONGOLIA Choibalsan
UlaanBaator

43°12.6'N
43°13.2'N
39°59.8N
38°34'N
37°27'N
36°05'N
35°19.6'N
35°27.3'N
41°46'N
39°56'N
36°41'N
38°03'N
28°36'N
34°18'N
32°00'N
31°24'N
30°37'N
28°12'N
34°43'N
30°40'N30°42'N
39°62'N
39°03'N
38°52'N
37°97'N
48°07'N
47°93'N141°49.3'E
141°39'E
139°59.1'E
140°57.3'E138°50.8'E
139°38.4'E132°44'E133°41'E
123°26'E
116°17'E
116°59'E
114°41'E
115°5'E
108°56'E
118°48'E
121°28'E
114°08'E
113°05'E
113°39'E
104°01'E
111°18'E
125°64'E
125°77'E
125°76'E
126°56'E
114°49'E
106°97'E

English NameWorld
Latitude Longitude weather
number47413 KOREA
37°57'N 126°96'E
47117
Seoul
(47412)36°30'N 127°39'E
47133
Taejon(47545)35°17'N 126°87'E
47156
Kwangju(47569) THAILANDKhen Kaen
16°26'N 102°50'E
48381(47572)17°23'N 102°48'E
48354
Udon Thani(47629)16°47'N 100°16'E
48378
Phitsanulok
47790
15°40'N 100°10'E
48400
Nakhon Sawan
47741
48407
Ubon Ratchathani 15°15'N 104°52'E
54342
13°42'N 100°34'E
48454
Bangkok
54511
19°30'N 97°82'E
48300
Mae Hong Son
54823
19°13'N 99°89'E
48310
Phayao53698
18°78'N 98°97'E
48327
Chiang Mai
57083 INDIA
28°00'N 73°29'E
42165
Bikaner
57036
28°58'N 77°19'E
42182
New Delhi
58238
27°48'N 95°01'E
42314
Dibrugarh
58362
26°10'N 91°57'E
42410
Gauhati
57494
25°45'N 81°72'E
42475
Allhabad
23°07'N 72°62'E
57679
Almodabad
42647
58606
22°65'N
88°44'E
42809
Calcutta
21°10'N
79°04'E
56294
Nagpur
42867
57461
43003
Mumbai(Bombay) 19°12'N 72°84'E
47050
16°20'N 81°14'E
43185
Machilipatnam
47058
13°00'N 80°17'E
43279
Madras12°95'N 74°82'E
47065
Mangalore
43284
47070
12°97'N 77°57'E
43295
Bangalore
44259
26°30'N 73°01'E
42339
Jodhpur
44292 * (47412)Sapporo, (47545)Akita, (47569)Sendai, (47572)Niigata, (47629)Kumagaya

光合成速度と水稲の穀物生産指標CPI

CPIU

FSter (Tc )tsa pc PAR PSNmax

PSNU dt
Michaelis-Menten type ........ Wheat，maize

apc PAR PSNmax

2

f rad _ mm

4m apc PSNmax PAR

2m

ここに PSN
βs
Tc
eLAI
Fster

：不稔の温度応答関数,

a mm PARbmm PARPrioul-Chartier type.......... Paddy rice

：光合成速度,PAR ：光合成有効放射：気孔開度，
a , b ：Michaelis-Menten定数： 葉温,
Tster ：不稔限界温度：有効葉面積指数f rad _ pc

apc PAR PSNmaxa pc PAR PSNmax 2 4m a pc PSNmax PAR
2m

fSyn：光合成温度影響関数

60

Michaelis-Menten型とPrioul-Chartier型(mgCO2/m2/hr)

日射影響関数の特性比較
502

ｆrad （gCO /ｄm /day）f rad・ f Syn (T c ・) ・
s eLAI

PSN

f rad・ f Syn (T c ・
)
・
eLAI
s

PSN
f rad _ pcth光合成速度PSNと作物

Fsyn

2

40
30

MM
PC(m=0.9)(m=0.8)(m=0.6)(m=0.4)(m=0.2)(m=0.1)

20
10

Paddy
Wheat
Corn
Orchard Grass
Forest (Sun T.)

40

20

0
0

10

20

30

40

Global Solar Radiation (MJ/m2/hr)

0
0

10

20

30

Global solar radiation （MJ/m2/day）

植生指標と
植物形態学

0.0

Surface Air Temperature Re-Analysis Data Ta

NDVI
0.5

NOAA
August,2000
SPOT May, 2003

0.9

NDVI

SPOT August, 2003

Bad harvest: May 2001

Good harvest: May 2005

NCEPからの豊作年の降雨分布

NCEP再解析データによる気温分布

Seasonal photosynthesis rate （CO2 fixation)
f rad・ f Syn (T c ・
)
・
eLAI
s

PSN

Terra衛星MODISによる土地被覆のデータとディシジョントリー法の適用

光合成温度影響関数
Sigmoidal関数によるS字型

Prioul-Chartier type.......... Paddy rice

1
1 exp k (Tc Thsv )f syn (Tc )Photosynthesis rate (mgCO2 m-2day-1)

Michaelis-Menten type ........ Wheat，maize
60

Best harvest in 1994
50
1993
1994
1995
1996
1997
1998
1999
2000
Ave.

40

30

20

10

・概略の分類はMODIS土地被覆データを利用する．細分の分類にエキスパート型ディシジョントリー法を適用する．

Worst in 1993
0
91

121

151

181

211

Julian day

Land Cover Classification

ASTER

241

ただし，植生の季節的な詳しい形態学や作物の作付け期，或いは水分指標等の情報が決定木のために必要である．IF文による独自のプログラミング．

水分指標・水ストレス指標

Paddy: Land cover classification

Shortwave Infrared Water Stress Index

M. N.
College

By Fensholt (2003)(SIWSI)

Paddy
test site

MODIS
B6

Matsue City

1.628-1.652μm

B2

NOAA

MODIS

0.841-0.876

SIWSI

2

6

6

2

= new NDWI or

LSWI by Xiao

ASTER
SWIR
B4

1.60-1.70μm
Sensible band for waterB3

0.76-0.86

SIWSI

3

4

4

3

エキスパート型
ディシジョントゥリー法による分類：
if Mod̲lulc=12 thenbegin

{ Crop land 12 }if ( (NDWI1̲XYave[i̲SPOT] > NDWI̲WaterLowLimit){ for South China }and (Water̲Con1 > NDVI1̲XYave[i̲SPOT] )and ( NDVI2̲XYave[i̲SPOT]> NDVI1̲XYave[i̲SPOT] ) )or ( (NDWI2̲XYave[i̲SPOT] > NDWI̲WaterLowLimit)and ( Water̲Con2 > NDVI2̲XYave[i̲SPOT])and ( NDVI3̲XYave[i̲SPOT]> NDVI2̲XYave[i̲SPOT] ) )

0.6

-0.1
or ( (NDWI3̲XYave[i̲SPOT] > NDWI̲WaterLowLimit)and ( Water̲Con3 > NDVI3̲XYave[i̲SPOT])and ( NDVI4̲XYave[i̲SPOT]> NDVI3̲XYave[i̲SPOT] ) ){ for Norh China paddy }or ( (NDWI4̲XYave[i̲SPOT] > NDWI̲WaterLowLimit)

Distribution of LSWI in Asian countries using
SPOT VEGETATION data.

and ( Water̲Con4 > NDVI4̲XYave[i̲SPOT])and ( NDVI5̲XYave[i̲SPOT]> NDVI4̲XYave[i̲SPOT] ) )then{ Water and growth definition }begin{ Second Paddy definition }LC[i]:=13;{ paddy }エキスバート型ディシジョントゥリー法による分類：（MODIS土地利用分類）＋（作物分類）

AsiaについてNDVIとLSWIにより決定木を作成

Solar Radiation

Geostationary Meteorological Satellite Himawari
January 2000

GMS Himawari

August 2000
January

AugustNorth Hemisphere

抽出された穀物域のAsiaの光合成速度（CO2固定） PSN分布計算されたAsiaの光合成速度（CO2固定） PSN分布

2001.08.01

2001.08.01

Chinaの穀物地帯における光合成速度季節変化

日本国内における単位面積収量予測（大潟）（春小麦，冬小麦，トウモロコシ，水稲）

6

50

4

40

1993MM
PC
1994MM
PC
1996MM
PC
2000MMPC
Average32
1

PSN (gCO 2/m2 /day)

CYI (t/ha)

5

30
20

0
91

121

151

181

211

241

-1

Nanjing
Shijiazhuang
Jinan
Shenyang

10-2

Julian Day

0

-3

91

Chinaの穀物地帯と森林・草原における光合成速度の季節変化

121

151

181

211

Julian day

241

Chinaの計算された穀物生産指標CPIの季節変化

6000
Nanjing
Shijiazhuang
Jinan
Shenyang

5000

2000
1800

4000

1600

3000

1400
2

2000
1000
0
91

iPSN (gCO2/m 2/day)Mongolia

121

151

181

Julian day

211

241

1200
1000

Shenyang 1998
Shenyang 2000
Nanjing 1998
Nanjing 2000
Chengdu 1998
Chengdu 2000

800
600
400

1250

Thailand

CPI (gCO2/m )iPSN (gCO2/m2/day)

China

Mae Hong Son
Ulaan Baator

1000

200
0
152

750

182

212

242

Julian Day

500
250
0
91

121

151

181

Julian day

211

241

更に収量データの年数を追加して年平均CPIを求め，次に収量を標準化（CPIU)し，単位収量CYIを推定する．

結論：

Provincial grain production in China

25
20

25

Lioning
Shandong
Hebei
Jiangsu
Sichuan
Hunan

Production (Million ton )Production ( Million ton)

30

Rice
15
10
5

Lioning
Shandong
Hebei
Jiangsu
Sichuan
Hunan

20
15

1980

1990

2000

著者らは，これまでに光合成・不稔型の穀物収量モデルを開発してきた．このモデルによりAsiaにおける光合成速度と穀物生産量の平面分布を推定するために，穀物の作付を分類し，光合成速度を推定した．

Wheat
10
5
0
1970

0
1970

1980

1990

2000

2010

Year

2010

Year

Production ( Million ton )15

10

Production ( Million ton )1.5
20

Lioning
Shandong
Hebei
Jiangsu
Sichuan
Hunan

Maize

5

1.0

Lioning
Shandong
Hebei
Jiangsu
Sichuan
Hunan

1970
1980

1990

2000

Year

2010

１. MODISの提供する世界の概略的な土地被覆分類データを基に SPOTVEGETATIONによる植生指標NDVIと地表面水分指標LSWIによる植生形
態学によって，MODISの被覆分類の一つである穀物の土地被覆を，水田・冬小麦・春小麦・その他（トウモロコシを含む）に細分類した．
２．判定された土地被覆分類により，著者の開発してきた方法を用い，Asiaの光合成速度の平面分布を推定した．
３．この光合成速度から，穀物収量の平面分布の推定が可能となる．

Soybean
0.5

４．現在，気孔開度CWSIの広域推定のために、GMS日射量と世界気象デー
タを用いて水ストレス指標CWSIを計算している。

0.0
0
1970

アジアの植生・作付分類と光合成（CO2固定）速度平面分布の推定

1980

1990

Year

2000

2010

地球地図・土地被覆データ整備における都市域の抽出
地球地図・土地被覆データ整備
における都市域の抽出
2008.2.28
2008.2.28
CEReS共同利用研究発表会
CEReS共同利用研究発表会
建石隆太郎・
建石隆太郎・
梶川昌三・筒井俊洋・山田美隆・阿久津
梶川昌三・筒井俊洋・山田美隆・阿久津 修・研川英征
修・研川英征発表内容
１．共同利用研究について
２．地球地図プロジェクト
３．地球地図・土地被覆データ整備における都市域の抽出
４．全球土地被覆データの作成と公開

千葉大学環境リモートセンシング研究センター
千葉大学環境リモートセンシング研究センター
国土交通省
国土交通省 国土地理院
国土地理院
地球地図国際運営委員会（ISCGM）事務局
地球地図国際運営委員会（ISCGM）事務局

共同利用研究の背景

１．共同利用研究について

国土地理院に事務局を置く地球地図国際運営委員会（ISCGM）は、各国の地図作成機関と協力して地球環境を正確に表す全球陸域のデジタル地図整備（地球地図プロジェクト）を進めている。
地球地図プロジェクト第二期データ整備は、2007年度中を目処に全球での整備を目指している。
特に、土地被覆データ及び樹木被覆率データは全球一括で整備する方針としている。
本共同利用研究は、各国の地図作成機関の協力を得つ
つ、効率的かつ高品質な土地被覆データ及び樹木被覆率
データを整備・公開する上で重要であり、土地被覆データの作成技術の高度化として位置付けられ、都市域の抽出に関する所要の検討及び手法開発は地球地図プロジェク
トの進展に大きく貢献する。

これまでの研究
平成16年度

樹木被覆率データ作成手法に関する研究
平成17年度

土地被覆分類データの検証手法に関する研究平成18年度

全球土地被覆データ作成のための公開蓄
積型グランドトゥルースデータベース（GTDB)構築に関する研究平成19年度共同利用研究
研究課題
グローバル土地被覆データの作成と公開
研究目的
地球地図プロジェクトにおける、第二期整備
データのうち全球土地被覆データについて、「都市域の抽出に関する所要の検討及び手法開発」を目的とする
研究内容
DMSPデータ(夜間可視光画像) 、人口分布デー
タの地域により可変な閾値で都市域抽出の可能性を検討する地球地図プロジェクトとは地球環境問題の解明と対策立案、持続可能な開発の実現、及び大規模災害の軽減などに貢献することを目的として、世界各国の協力により全球陸域の基盤的地理情報(=地球地図）を整備するプロジェクト

２．地球地図プロジェクト

地球地図とは地球上の全陸域を対象とするデジタル地理データセット
統一規格、非営利目的であれば無償でデータ利用可能

縮尺100万分の1 (空間分解能1Km）
8種類のデータからなる
境界、交通網、水系、人口集中地区、標高、土地利用、土地被覆
植生

データは5年ごとに更新

地球地図プロジェクトの経緯
プロジェクトの立ち上げ

環境に関する国連会議

地球サミット.「アジェンダ21」を採択

1992

ベクタデータ： 4 項目（VPF形式）

地球地図の主な進展

我が国が地球地図構想を提唱

1994

地球地図データ（１）境界

地球地図に関するJICA集団研修を開始

1996

地球地図国際運営委員会（ISCGM）の設置

1997

国連環境開発特別総会.採択文書に地球地図が言及

1998

国連と連携した地球地図参加要請キャンペーン

水系

データの
整備・
提供

地球地図データの提供開始

2000
2002

海岸線・行政界

ヨハネスブルグサミット.「実施計画」に地球地図の促進を明示
アフリカでのJICA第三国研修に地球地図セミナーを追加

2007

交通網

河川・内水面・ダムなど

人口集中域

地球地図が全球陸域をカバーする目標年次

今後の取組

162ヶ国・16地域がプロジェクトに参加（全球陸域の97%に相当）
49ヶ国・3地域の地球地図データを公開中（113ヶ国・14地域のデータを公開準備中）
道路・鉄道・空港など

全球陸域の地球地図の公開・利用の普及

2008

都市の位置・名称など*2008年2月現在

地球地図データ（２）ラスタデータ： 4 項目（BIL形式）
標高

中南米
コロンビア （2001）
パナマ （2002）
メキシコ （2003）
アルゼンチン（2006）
・・・・

土地被覆

1m 単位の高さ情報

土地利用

欧州
ラトビア （2006）
マケドニア（2006）
・・・・

地表の状態を17種類に分類

樹木被覆（植生）アジア
日本 （2000）ラオス （2000）
ネパール （2000）
スリランカ （2000）
タイ （2000）
バングラデシュ（2001）
モンゴル （2001）
フィリピン （2001）
カザフスタン （2003）
キルギス （2003）
ミャンマー （2003）
イラン（2005）
ヨルダン（2007）
・・・・

アフリカ
ケニア （2002）
ブルキナファソ （2003）
ボツワナ（2003）
スワジランド（2004）
トリスタンダクーニャ（2006）
グルジア（2007）
・・・・
大洋州
オーストラリア （2001）
サモア（2004）
南極（2006）

土地利用の様子を9種類に分類

樹木被覆率を％で表示（ 2008.2.21現在）

地球地図プロジェクトの運営
地球地図国際運営委員会 (ISCGM) がプロジェク
トを運営
17ヶ国の地図作成機関の長など20人の委員
委員長： テイラー教授（カナダ・カールトン大学）
国土地理院が事務局をつとめる

各国の地図作成機関（NMO)がそれぞれの国の地
球地図を作成する

ISCGM/WGについて
ISCGM/各WGの内容
WG１：地球地図戦略
WG２：地球地図仕様の作成・改定
WG３：データポリシー
WG４：ラスタデータの整備

ISCGM/WG４の活動
CEReS、建石教授が議長

テーマ：土地被覆データ、樹木被覆率データ等のラスタデータの全球整備・更新
2008年に第二期ラスタデータの整備手法の確立・
全球整備を目指し土地被覆データ、樹木被覆率デー
タ整備を実施中

自力での整備が困難な国は我が国などが支援作成された地球地図は、インターネット（www.iscgm.org）で公開

研究概要
都市域の抽出が満たすべき条件の検討

３．地球地図・土地被覆データ整備における都市域の抽出

既存の土地被覆データ及び関連する主題地理情報の所在、内容、特徴を調査
既存の土地被覆データの内容及び関連する主要地理情報と地
球地図・土地被覆データを比較検討
都市域の抽出が満たすべき条件を検討・決定

都市域の抽出手法等の開発
都市域の抽出手法及び検証手法を設計し、工程及び方法論の確立
都市域の抽出手法検証手法の開発

都市域の抽出手法の検証
都市域の試験抽出
精度検証、手法の妥当性評価
評価に応じた改良

既存のグローバル都市域データ
データ名

作成機関

ラスター/ベクター（解像度/縮尺）

グローバル都市域データの特徴
DCW （Digital Chart of the World )情報源

都市域データwebサイトDigital Chart of theWorld (DCW)US Defense Mapping
Agency (DMA)1993年作
成Environmental Systems
Research Institute, Inc.
ESRI)配布

ベクター
1:1,000,000

1:1,000,000 Operational
Navigation Chart (ONC)

Urbanized area andsmall villagesin
Populated Place layer

http://www.maproom.ps
u.edu/dcw/

GLCC Version 2

US Geological Survey

ラスター
1km

AVHRR 1992

Urban and Built-Up
in land cover legend

http://edcsns17.cr.usgs.
gov/glcc/

UMD Land CoverClassification

University of Maryland
1998年作成

ラスター
1km

AVHRR 1992

Urban and Builtin land
cover legend

http://www.glcf.umiacs.
umd.edu/data/landcover/

MOD12Q1 V004 Land
Cover Product

Boston University
2002年作成

ラスター
30秒

MODIS 2001

Urban and Built-Upinland cover legend

http://duckwater.bu.edu/lc/mod12q1.html

GLC2000

European Commission’s
Joint research Center
2004年作成

ラスター
1km

SPOT/VEGETATION
2000

Urban (Artificial
surfaces andassociated areas)inland cover legend

http://wwwgvm.jrc.it/glc2000/

Global Rural-Urban
Mapping Project(GRUMP) alphaCenter for International
Earth Science Information
Network (CIESIN)
2004年作成

ラスター
30秒

DMSP-OLS
1994/1995,DCW,Tactical
Pilotage Charts

Urban extents

http://beta.sedac.ciesin.
columbia.edu/gpw/index.jsp
http://beta.sedac.ciesin.
columbia.edu/gpw/docs/UR_paper_webdraft1.p
df

人口集中域はDMA データソースから作成, 1950-70年代に整備された航空図で、多くの都市の現在の大きさの正確な図はもはや提供していない。
http://www.maproom.psu.edu/dcw/

GLC2000 (Global Land Cover 2000)
都市の土地被覆クラスを抽出
http://www-gvm.jrc.it/glc2000/Products/

DMSP（defense meteorological satellite program）（防衛気象衛星プログラム)

世界の夜間光データは1994年10月-1995年3月（のデータ)から作成された。粗い(2.7Kmの) 解像度を処理したこの地図は記録（レジスター）の状態が悪く、都市の境界を膨張させるブルーミング効果を示している。
http://dmsp.ngdc.noaa.gov/

MODIS 土地被覆成果-都市範囲
データ: MODIS, DMSP 及びグリッド化された人口 (解像度 5km).
http://duckwater.bu.edu/urban/global.html

都市域の過大推定
現在の都市域より狭い

GRUMPによる日本の都市域

GLC2000による日本の都市域

GRUMPの評価

GLC2000の評価
80000
60000

estimated city size(squar e meter)

60000
40000
20000
0

2 000040 00 0a c t u a l c it y s iz e (s qu a re m e te r)

60 00080 00000e sti ma ted city size(square meter)

G L C 2 0 0 0 U r b a n A re a

40000

N ig ht tim e L ig h t s D a t a20000

メリーランド大学土地被覆データによる日本の都市域

80000

ＤＣＷによる日本の都市域

02000 0400 00600 008000 0a c tu a l c it y s iz e ( s q ua r e m e t e r )ボストン大学土地被覆データによる日本の都市域

Global Map Land Cover Version1.1による日本の都市域

地域的な都市域データ

ブルーミングの問題のため都市の大きさを過大評価する

都市の大きさを多少低く評価し、多くの市がゼロになった

都市域抽出の前提条件
都市域の定義都市域データの品質
都市域データ作成のための情報源
LCCS （Land Cover Classification System）
との関係
検証手法

建造物、
商工業・交通施設
都市域内緑地CORINE land cover（既存データにおける）都市域の定義
DCW
GLCC（=IGBP-DISCover）
ボストン大学作成土地被覆データ
GLC2000土地被覆データ
GRUMP urban extents データ
CORINE土地被覆データ
それぞれは、統一された定義でない
厳密な定義がない

都市域の定義CORINE Land Cover（土地被覆）の例
地球地図（土地被覆）の場合
地表面の単位面積に対する地表面に正斜投影した人工構造物の比率： a
土、砂など裸地の比率：ｂ
植物の比率：ｃ
このとき、a＋b＋c=1
a＞0.5 都市域
b＞ 0.5 都市域
しかし、(bまたはｃの)世界共通の閾値決定は難しい。

都市域データ作成のための情報源
及び検証

都市域の抽出手法
人口密度データ
人口密度データ

MODIS(2003)データ
人口密度データ
DMSP/OLSデータ (夜間可視光画像）
LCCS（Land Cover Classification System）の利用
Landsat (比較データ）
データ作成•Landsat ETM+ データで都市域を描写• 人口密度データを処理、
MODIS 及び DMSP
•人口密度を合成、 MODIS及び DMSP

Landsat
Landsat ETM+
ETM+
閾値処理
閾値処理
都市域候補地
都市域候補地
閾値処理
閾値処理

MODIS
MODIS NDVI
NDVI

閾値処理
閾値処理

DMSP/OLSDMSP/OLS

緑地の削除
緑地の削除

途上国集落の削除
途上国集落の削除

GLCNMOGLCNMO 都市域データ
都市域データLandsat
Landsat ETM+
ETM+ および既存の都市域データとの比較による検証
および既存の都市域データとの比較による検証

都市域抽出されたGLCNMOおよび使用データと他の都市域画像（東京）

都市域抽出されたGLCNMOおよび使用データと他の都市域画像（北京）

都市域抽出されたGLCNMOおよび使用データと他の都市域画像（ロンドン）

都市域抽出されたGLCNMOおよび使用データと他の都市域画像（米国・アトランタ）

都市域の抽出試行と検証
世界30地域の都市域抽出
閾値の決定
人口密度
（人/ km2）

MODIS/NDVI

DMSP/OLS

アジア

> 800

< 0.52

> 10

ヨーロッパ

> 500

< 0.72

> 20

アフリカ

> 800

< 0.52

>5

北米

> 500

< 0.62

> 20

米国東海岸

> 500

< 0.82

> 20

南米

> 500

< 0.52

> 10

オセアニア

> 500

< 0.52> 10

都市域抽出されたGLCNMOおよび使用データと他の都市域画像（ババングラディシュ、ダッカ）

都市域抽出されたGLCNMOおよび使用データと他の都市域画像（ニューヨーク）

試行結果

抽出都市域の検証（ロンドン）

Landsat
Landsat TMカラー合成画像（Red
TMカラー合成画像（Red Green
Green Blue:
Blue: band
band 2,
2, 4,4, 7）では都市域は紫7）では都市域は紫色で明瞭に判読できる。
色で明瞭に判読できる。

GLCNMO（ロンドン）のLandsat ETM+画による検証と
人口密度、他の都市域データとの比較

北京郊外では点在する村落、ダッカ周辺では道路網に沿って、人口密集地域
北京郊外では点在する村落、ダッカ周辺では道路網に沿って、人口密集地域が現れており、都市域より広い範囲で人口密集地域が確認できる。
が現れており、都市域より広い範囲で人口密集地域が確認できる。

・GLCNMO・Landsat ETM+画像
・人口密度
・DCW
・ボストン大学土地被覆データ
・GLC2000・GRUMP・CORINE土地被覆データ

NDVI値の低い地域（NDVI画像の黒い部分）は概ね都市域と対応するが、ア
NDVI値の低い地域（NDVI画像の黒い部分）は概ね都市域と対応するが、アトランタの例では居住地周辺に樹木が多く、NDVIは高くなっている。
トランタの例では居住地周辺に樹木が多く、NDVIは高くなっている。
抽出した都市域(GLCNMO)はDCWと比較して、アジアでは広くなっているが、抽出した都市域(GLCNMO)はDCWと比較して、アジアでは広くなっているが、欧米では狭くなっている。DCWはGLCNMOと比べ約40年程度過去のデータ
欧米では狭くなっている。DCWはGLCNMOと比べ約40年程度過去のデータであるため都市域は一般に拡大している。しかし、DCWの欧米では緑の多いであるため都市域は一般に拡大している。しかし、DCWの欧米では緑の多い
住宅地も都市域に含めているため欧米では縮小しているように見える。
住宅地も都市域に含めているため欧米では縮小しているように見える。

ほぼ合致している

ボストン大学土地被覆データは概ね抽出した都市域(GLCNMO)と似ている。
ボストン大学土地被覆データは概ね抽出した都市域(GLCNMO)と似ている。
しかし、東京では過大評価している。これはDMSP/OLS（GRUMPと同じ情報
しかし、東京では過大評価している。これはDMSP/OLS（GRUMPと同じ情報源）の影響と推測される。
源）の影響と推測される。
GLC2000はアジアでは過小評価している傾向がある。
GLC2000はアジアでは過小評価している傾向がある。
GRUMPは都市域を過大評価している。
GRUMPは都市域を過大評価している。

（アジア）

（アジア）

（ヨーロッパ）（北米）

まとめ
米国Oak Ridge National Laboratory（ORNL）が作成した人口密度データの閾値処理により都市域の候補地を抽出する。
MODIS/NDVIデータの閾値処理により緑地の多い地域を都市域から除外する。
DMSP/OLSデータの閾値処理により夜間光の無いあるいは少ない部分を都市域から除外する。
抽出された都市域（パリ）（上左の画像）とDCW(上右)とその一部（矢印の部分）のGoogle Earth画像

多偏波 SAR データを用いた海氷物理量測定

P2007-1
西部太平洋熱帯域と東アジアにおける気候環境の関連性に関する調査
The relationship of climate environment between western tropical Pacific andEast Asia
菅野洋光

Prima Oky Dicky A.

東北農業研究センター やませ気象変動研究チーム〒020-0198 岩手県盛岡市下厨川字赤平４
TEL 019-643-3408 FAX 019-641-9296
kanno@affrc.go.jp

岩手県立大学ソフトウェア情報学部〒020-0193 岩手県岩手郡滝沢村滝沢字巣子152-52
TEL 019-694-2534 FAX 019-694-2595
prima@soft.iwate-pu.ac.jp

1. はじめに
アジアの気候を考える場合、当該地域における1970年代後半の気侯変化(climate shift)を観察することが非常に重要である。climate
shiftを境にしてエルニーニョの周期が変化し、また南アジアモンスーンへ影響する大気−海洋の分布パターンもそれ以前と以降とで異なったものが卓越している。これらの影響として、東・北日本に関しては、climate shift以前ではフィリピン付近の海水面温度(SST)が局在性を持たず、ロスビー波の伝播がほとんど認められないのに対して、climate shift以降ではフィリピン付近SSTの局在化とそれに伴う対流活動の活発化、ロスビー波（PJパターン）の励起が顕在化した。また、エルニーニョの周期は、1970年代後半以降、4〜5年の長周期変動が卓越しているが、それによって影響されるPJパターンの変動により、北日本の夏季天候も1982年以降明瞭な5年の周期変動を示している。その結果、
1983,88,93,98,2003年の冷夏、1984,94,99,2000,2004年の暑夏と、周期的に冷夏・暑夏に見舞われている。一方、インドネシアは、海洋大陸として対流活動の中心的位置を占めているが、近年その対流活動の変動により、しばしば大雨や干ばつに見舞われている。本研究では、北日
本夏季天候とPJパターンとの関係について解析し、その発現域であるインドネシアと、影響域である東アジアでの気候環境について、気温と降水量を用いて考察する。

2. データSST デ ー タ は 2007 年 9 月 に リ リ ー ス さ れ た NOAA の ERSST (Extended Reconstructed Sea Surface Temperatures)データver.3を、そしてグローバル気
象データはNCEP/NCAR再解析データを利用した。気温および降水量は、観測デー
タに基づいて作成されたUniversity of Delaware (U-Del)の0.5度グリッドのデータを用いた。

3. 結果
a.Kanno(2004)により、北日本夏季天候の1982年以降の5年周期が指摘されていたところであるが、図1のように、2007年夏季までその周期に則っているようにみえる。
b.南 シナ 海と フィリ ピン 東方 海域 との 海水面 温 度 (SST) 東西 差を、 C(5-15N,
110-130E)-D(5-15N, 140-160E)で計算した。図2には、C-Dと500hPa高度との相関係数分布を示す。赤道域の広い範囲に正の相関が見られるほか、日本から東シベリアにかけての波列パターンが明瞭に把握できる。PJパターンが北日本夏
季天候の周期的な年々変動と深く関わっていることが明瞭である。

図1 八戸における夏季平均気温（上）と稚内マイナス
仙台の気圧差(PDWS)の1982年を起点とした5年ごとの重ね合わせc.インドネシア、ジャワ島におけるU-Del雨量データと地上気圧との相関係数分
布を図3に示す。ジャワ島南東方の広い範囲で負の相関を示し、雨量の増加に気圧の低下を伴うことが把握できる。また、北日本に正の相関が見られ、ジャワ島の雨量が多いときには気圧が高まる、すなわち夏季気温が高温になることが示唆される。これらについては、PJパターンの伝播による影響の可能性としても考えられるが、波源の地上雨量と影響域でスポット的に有意な関係がみられるのは非常に興味深い点である。インドネシアにおいては、2006年から現地調
査ならびに雨量観測を実施しており、今後これらの現象を詳細に分析する予定である。

4. 今後の課題
以上のような気象学的な結びつきは、食料事情と関連させると、その重要性は明らかである。インドネシアと日本と同期する天候変動パターンを把握できれば、干ばつと冷夏が同時に起こりうるのか、もしも同時に起こった場合、安定的に食料を受給できるのか否かについての分析が可能になると考えられる。
今後の課題として、農研機構による作物生育モデルと天候変動パターンに対応
した広域農作物生産予測モデルを提案し、気象災害発生時の農作物生産シナリ
オの開発を試みる。

図2 SST東西差(C-D)と500hPa高度との相関係数分布。
1982-2007年の夏季(JJA)。陰影域は危険率5%以下
で有意。

参考文献
Kanno, H., 2004, Five-year Cycle of North-South Pressure Difference asan Index of Summer Weather in Northern Japan from 1982 Onwards. Journalof the Meteorological Society of Japan, 82(2), 711-724.
図3 ジャワ島における降水量と地上気圧との相関係数分布。1982-2006年の夏季(JJA)．陰影域は危険率5%以下で有意。

Study of Coral Reef Distribution in Bali island using Satellite data
Takahiro Osawa*, I Ketur Suwardika**Center for Remote Sensing and Oceanography (CReSOS)
Udayana University

1. Background

6. Results•Sea Alga are important from view point of fishery resource andenvironmental management.

Table.1 Statistics Value on Each Spectral Class on the RGB Landsat Mapping Channel.

•Mapping of sea alga are utilized for estimation of environmental
index in coastal area.

No

Spectral Class•Optical property is consist from both land (reflection from solid
ocean bottom) and ocean (attenuation by seawater).

1•Classification is depends on optical properties of sea water above the
Alga.

Channel 1

Channel 2

Channel 3

Mean

Std

Mean

Std

Mean

Std

Coral Reeves

43

17.4

53

16.6

36.5

14.9

2
3

Sands
Algae

55
26.4

24.6
9.3

69.5
41.8

28.4
8.9

93.9
35.4

35.8
8.8

4

Wave Particles

139.6

70.4

129

66.4

129.9

66.8

2. PROBLEM FORMULA• How Alga Habitat and Benthic associate with Coral Reef Can
Be Mapping By Remote Sensing.

B ubles・ What is The Kind of Maps Distribution and ecological
Condition Of Alga In Bali Island at Year 2003.

Sands

3. AIM AND OBJECTIVES

CoralsS ands

1. To know Alga Habitat And Benthic associate with Coral Reef
Can Be Mapping by Remote Sensing.

Algae

2. To understand and Provision the method Of Alga Condition.
3. To create Map Distribution and Ecological Alga Condition
Around Bali in 2003.
Figure.1 Spectral Object Structure in Channel 1 and Channel 2 Spectral Rooms

4. FRAME WORK OF RESEARCH

Table.2 Relationship among the spectral classes and the associated community
No

Research Areas

Researched Objects

Amount of Pixels%

Area (Ha)S T A R T

In S i tu D a ta

1 Canggu
C i tr a L a n d s a t(L ev el 1G )H om o gen
S u b s tr a c t
d a ta ( s a n d )80%
S a m p le s
S et ofIn S i tu D a t a )B andsG e o m e tr i c
C o r r e c ti o n

T opograph y
M ap( T i e P o i n ts )B a n d s A tm o s p h e r e
C o r r e c ti o n(6 S C ode)A n c illa ry
D a ta2 Lebih

10,392

4.05

Sands

438

0.17

Coral Reef

90

0.04

81

Algae

239

0.09

2.151

Wave Particles

2871

2.21

2583

1204

0.93

713

0.55

Sands
Coral Reef
Algae

L a n d M a s k in g
L o u d D e te r m i n a ti o n
2 0%
T es tS et ofI n S i tu D a ta

Wave Particles

3 Nusa Penida
P e m b u a ta n C i tr a
K o m p o s it W a rn a a s li
(R G B 321)

Wave Particles
Sands
Coral Reef
Algae

S u p e rv is e d R G B 3 2 1
C la s ip ik a s i

S u p e r v i s i d m u l ti p l e d e p thin v a ria n c e in d e x

4 Candi Dasa

Wave Particles
Sands

M a x im u m lik e lih o o d
m a j o r i ty fi l te r

M a x im u m li k e lih o o d
m a j o r i ty fi l te r

P e ta K l a s i fi k a s i

P e ta K l a s i fi k a s i

9,352
394

1083
641

573

0.44

515

17,246

3.83

15,521

1275

0.28

1147

755

0.17

679

538

0.12

484

21,1201.71

19,008

1335

0.11

1201

Coral Reef

49

0

44.1

Algae

87

0.01

78.3

T in g k a t a k u ra s i

P e ta S e b a r a n d a n a n a l y s t
E k o lo g i K o n d is i A lg a

5. LOCATION AND TIME RESEARCH
The location where this research was done was in the ocean of Bali that was
included on the mapping of Landsat-7 ETM Path/Row : 116/066. These areas
include seas surrounding areas of Canggu, Benoa, Lebih Beach, Serangan,
Nusa Penida, and Lembongan. Below is the location map of these areas.
This research lasted for 6 months; it started in March and ended in
September 2006. The activities were collecting material activity and relevant
literature, supporting materials and tools, input data.
Scope Research cover satellite image interpretation of LANDSAT,
for the purpose of identifying Alga distribution in Bali Island,
continued with mapping of Alga distribution in Bali Island.

Figure 2. The Algae Habitant Map in Nusa Penida Sea, Klungkung State,by Landsat Mapping in May 2003

Figure3. Algae Inhabitant Mapping Area in Lebih Water, Gianyar State, by Landsat Mapping,in May 2003

7. Conclusion
1. The Landsat Satellite mapping acting as the source of data in mapping the expansion ofalgae in Bali shows the level of accuracy as big as 78.94% without using the water column
correction and 89.47% with using the water column correction. That accuracy is present
above the limit of the standard level, 75%.
2. The result show the differentiate 4 classes of object visibility and the associated community
on algae ecosystem. The four classes are coral reef, sand, wave particles, and algae.
3. The algae expansion using the supervised classification give algae expansion community
area in east Bali is 3228.3 ha (1437 pixels).

ミクスチャモデルによるMODISデータの雲量推定
東京情報大学 朴 鍾杰、安田 嘉純

概

実

要

広域観測データを用いて環境変化モニタリングを行う際、雲の存在が障害

験○ WIとMOD35との比較

既存の方法 ： 雲の特徴温度 ↓ 、植生指数↓
Max N
Max T
Min NIR
Cloud Mask

c

a

b

問題点 ：
1. コンポジット結果が雲の影響があるがどうか判断できない
2. コンポジット期間中の中間日が晴れても選ばない可能性がある
本研究の目的
ミクスチャモデルにより、比較的簡単な計算を用いて、画素中に含まれる雲の割合を示すWhite Index（WI）を提案する。

99

95

66

曇り

0

5

10

15

30

100

a)

Mixture Model ( Gradient Method)R(λ) = a×W(λ) + b×V(λ)＋ c×S(λ)
100 = a + b + c（1）

R(λ) = a×W(λ) + (100−a)×Oth(λ)
a = R(λ) / W(λ) − (100−a)×Oth(λ) ／ W(λ)| W(λ1) −W(λ2) | >> | Oth(λ1) −Oth(λ2) |したがって、| Oth(λ1) −Oth(λ2) | / | W(λ1) −W(λ2) | ≒ 0a ≒ | R(λ1)−R(λ2) | / | W(λ1)−W(λ2)|

（2）
（3）

（4）
（5）

| Oth(λ1) −α×Oth(λ2) | / | W(λ1) −α×W(λ2) |

（6）

b)

D(α) =
Σ (| Oth(λ1) −α×Oth(λ2) | / | W(λ1) −α×W(λ2) |)2 (7)MOD09反射率○ 時系列WIコンポジットとMVCとの比較
0.61

0.78

| (B1+B3+B4)-1.5×B7|

0.61

0.59

0.76

100

0.57

300

90

258

80

250

70
反射率(%)

c)

210

60

200

50
40

a) 2003年10月9日

150

b)

20日までのWIコンポジット結果c)

２７日までのWIコンポジット結果

30
20

100

10
0
B3

B4

B1
B2
B5
MODISバンド（波長順）

B6

B7

都市
雲

50
4

10

2

sea

city

desert bare-soil

2

1

0

流氷
砂漠
森林

海
裸地（農耕地）

λ1

λ2

α

D(α)

最大誤差土地被覆

B3| B4-B1 |

2.2

107

植生域 (-7.8)

sea-icecloud

B3+B4+B1

B7

1.5

26

都市域 (4.5)

B3+B4+B1

B6

1.3

82

都市域 (8.6)

WI =(R(B1)＋R(B3)＋R(B4)−1.5×R(B7)) ／ 210 × 100

foreste) 9日から27日までのMVC結果

ま と め

（8）

D） WIコンポジットによる更新日f) MVCによる更新日

9日

20日

25日

27日

MODISデータの画素内雲の混合比率を求めるためにWIを提案した。
WIの有効性を調べるために、混合スペクトルを合成し、WIと比較した
結果、都市域で最大4.5%の誤差が生じることが分かった。
雲の分布を表すMOD35との比較を行いWIが有効であることが分かった。また、WI画像を用いることで、時系列的にコンポジットが可能であることが分かった。

Tokyo University of Information Sciences

高分解能合成開口レーダのクロス偏波画像に見られるテクスチャを利用した針葉樹林バイオマス計測

夜間光をパラメータとする二酸化炭素排出量の推定− 東アジア地域における二酸化炭素の排出量の変動について −原
政 直（株）ビジョンテック

１．背景と目的
温室効果ガスによる地球温暖化問題に関して、京都議定書後の新たな枠組みつくりが議論されている。特に、二酸化炭素ガス（CO2）の排出量の削減については、2007年6月のG8サミットで「2050年までに世界の排出量を半減させる」という目標を軸として2009年までに新たな枠組みを作ることが確認された。京都議定書では、先進国にこのCO２の排出削減責任を課しているが、途上国とされている中国をはじめ主要排出国のすべてが削減に向かう必要がある（図1）。現在、各国のCO2排出量は統計的データとして集積されているが、その統計量の集積時期やその精度にばらつきがあり、また、データが空白となっている国や地域もある。また、このCO２は経済活動に伴うエネルギー消費、特に、化石燃料の転換エネルギーである電力の消費による比率が高い（図２）。
そこで、夜間の人工光は電力消費やCO２排出量と一定の相関がある *1,*2ことを利用して、米国極軌道周回型気象衛星DMSP (Defense Meteorological Satellite Program)のOLS( Operational Linescan System ）で観測された過去12年間の夜間人工光のデータを用いて、センサ特性の正規化を行った上で、東アジア地域における光力（輝度）とCO2排出量の変動との関係、および、経済成長の指標であるGDP（Gross Domestic Products）との関係を調べ、データ空白の地域、国の排出量の推定や精度のバラつきの補正など、客観的で均一な精度でCO2の排出量の推定を行うことを目的とした。
世界のCO２排出 量（ 200 4年）
ーエネ ルギー ・経済 統計要覧200 7年版 ー（ト ンは 炭素換算）

２．使用データ
・衛星データ2 3%： DMSP/OLS １kmメッシュ年単位グローバルデータセット（1992年-2003年、米国NGDC ）か

その他

ら

米国

1%
７２億 トン

1%

ド
イン

日
本

2%
2%
2%

5%

中国
ロシ
ア

1%

北緯60度〜南緯20度、東経60度〜170度の範囲を切り出して使用した（図3）。
・ベクターデータ ： CIP (The International Potato Center (Centro International De LA PAPA)
・統計データ
： ①国際連合世界統計年鑑（2000年版、2005年版）、国際連合統計局編②2003 Energy Statistics Yearbook United Nations

18%

EU
13 %

4%

日 本 エ ネル ギ ー 経 済 研 究 所 作 成（CO 2は エ ネル ギ ー 起 源 ）
鉄 鋼

：6.0％

セメ ン ト ：3.0 %紙 パ ル プ ：0.8%

産

その 他 ：1 1.4%

業,

21.20 %
発 電 など ,44.40%

発 電な ど運 輸,
21.50%

民生

オース トラリ ア
インド ネシ ア
ブラジ ル
その 他

15ヶ国

6%

世 界 の 分 野 別 CO ２ 排 出 量（ 2004年 実 績）

産業 分野内訳

米国
中国
EU（ 1 5 カ国 ）
ロシ ア
日本
インド
カナダ
韓国
メ キシ コ

22 %

民生,

運輸

12.90%
産業

図１．世界のCO2排出量（2004）

図２．世界の分野別CO2排出量（2004）

３．データセットの正規化
12年間のDMSP/OLS夜間光データはF10、F12、F14、F15の4機で観測されたもので、センサの劣化の影響やセンサ特性の違いにより、輝度値にバラつきがある。このバラつきを補正してセンサ間の正規化を図るために、F12の1999年の特性を基準とし、各画像中の輝度値が「0」となる地点、すなわち、周辺の光の影響のない湖沼の中心点と、年間を通じてセンサゲインの調整により、輝度の消失や飽和がない地点をランダムに選択（表1）し、基準センサと各補正対象センサで得られた同一地点の輝度値を計測し（図4）、
その直線回帰により得られた回帰パラメータ（図5）を用いて正規化した。
表1．輝度計測地点一覧

70

70

60

60

F10 1992
F10 1993
F10 1994
F12 1994
F12 1995
F12 1996
F12 1997
F14 1997
F12 1998
F14 1998
F12 1999

50
40
30

F

20
10
0
1

2

3

4

5

6

7

8

9

10

地点番号

y = 0 .9601 x
2

R = 0 .8503

50

F12 1999

40

30

40

30

20

20

10

10

0

0
0

10

20

30

40

50

0

60

10

20

30

40

y = 1 .2 657x

y = 0 .9657 x

2

R = 0.9 457

2

R = 0 .8362
40

30

20

30

20

10

10

0
0

図３．解析範囲

70

50

40

図4．衛星別年次別センサ特性のバラつき

60

60

50

F14 1999
F14 2000
F14 2001
F15 2001
F14 2002
F15 2002
F 15 2003

50

F12 199 4

F 10 19 92

60

F 12 1999

60

y = 1.082 3x2
R = 0.816 3

50

F12 1 999

衛星別年別地点別光力のバラつき

F12 19 99

70

5

10

15

20

輝度計測地点25
F 14 19 99

30

35

40

45

0

50

0

10

20

30
F 15 2 002

40

50

60

図5．F12に対する各センサ特性の相関

４．国別の夜間光とCO2排出量およびGDPの変動との関係正規化された夜間光画像（図6）を用いて、国別輝度値の年次変動（図7）、世界統計による国別のCO2排出量（図8）とGDPの年次変動（図9）を示す。また、図10には、国別輝度値とGDPと関係を示し、図11には、国別輝度値とCO2排出量との関係を示した。
200000 00
180000 00
160000 00
140000 00

１９９２

１９９5

中国
インド
インドネシア
北朝鮮
韓国
モンゴル
タイ
ベトナム
日本

120000 00
100000 00
80000 00
60000 00
40000 00
20000 0001992 1993 1 994 19 95 199 6 1997 1998 1999 2000 2 001 2002 20 03

１９９8

2001

図7．国別輝度値の変動

図9．統計による国別GDP

図8．統計による国別CO2排出量

図6．年次別正規化夜間光画像
GDP vs 光 量 （1 992年〜20 03年 ）年次別の輝度値とGDPの散布図を見ると、中国、インドは輝度値が高くなる傾向にありながらGDPが低い。これに対して、日本はGDPが高い割りに輝度値が低い。これは日本の生産効率が良い（単位GDP当たりのエネルギー消費量が低い）という見方もできる。一方、輝度値とCO2排出量の散布図を見ると、日本は輝度値とCO2排出量の差が年次によって変化していないことが伺える。これに対してインドは輝度値の変化がある割りにCO2の変化が少ない。もともと水力発電への依存度が高いことがその理由に上げられる。また、中国は輝度、CO２排出量共に大きく変化をしていることが伺える。図12は2000年の輝度値を国別に示したものである。また、図13は1992年の各国における輝度値を基準として、2000年における変化の度合いをカラーで示したものである。赤が強くなるほどその変化が大きいことを示している。これによると東南アジア地域の変化が著しく、中国、インドがそれに続いていることがわかる。日本や北朝鮮はほとんど変化がな
い。これは北朝鮮のエネルギー事情や、日本における省エネルギー化が進んでいることを裏付けるものと考えられる。

CO2 vs 光 量 （ 1 992 年〜 200 3年 ）25 ,0 00, 000

25, 0 00,0 00

中国

20, 0 00,0 00
DMSP国別DN積算

20 ,0 00, 000
DMSP国別DN積算

輝度変動を見ると、絶対量、増加率ともに中国、インドが高い値を示しており、日本は緩やかな増加基調を示している。しかしながら、2002年、2003年の輝度が急激な下降基調に変化していることが認められる。また、統計量で見るCO2排出量では、中国、インドは増加基調を、また、日本が緩やかな増加基調と傾向としては輝度変動に大よそ一致している結果が得られた。しかし、中国をはじめ1997年以降、下降基調になり、2000年から上昇基調に転じている。この原因は、
1997年にタイ・バーツが変動相場制へ移行したことを発端としたアジア通貨危機が影響しているものと考えられる。この影響は、GDPについても同様のことが言える。しかし、中国のGDPが波状変化を示し、インド、日本との差が大きく出ているなど、必ずしも輝度変化とは一致していない。

15 ,0 00, 000

10 ,0 00, 000

日本
インド
5 ,0 00, 000

15, 0 00,0 00

インド

日本

5, 0 00,0 00

韓国

0
0

150 0 000

中国

10, 0 00,0 00

03 000 00 04500 0 00

6 0 000 00

0

1 0000 0 0 200 0 000 30 000 0 0 4000 0 00 5 00 0000

GDP

CO2

図10．国別輝度値とGDPとの関係

図11．国別輝度値とCO2との関係

図13．1992年を基準とした2000年の輝度値変化度

図12．2000年の輝度値画像

５．まとめ

DMSP Nighttime Lights
15000000
12500000

Sum of Lights Index

OLSによる夜間画像で得られる輝度値は、CO2排出量の傾向を良くあらわしていることが認められた。
しかしながら、長期にわたるモニタなどを行う場合、衛星センサの特性の違いを十分に考慮した補正をするなど、精度の高い情報抽出には、さらに考慮すべき点がある。図14はF10、F12、F14、F15のセンサ特性を示すものであるが、センサ間でその特性がかなり違う。過去のデータを一連のものとして使用するためには、このセンサ特性の補正が今後の課題である。

10000000

F10
F12
F14

7500000

F15

5000000
2500000

参考文献*1 Doll, C.N.H., Muller,J.P. and Elvidge,C.D., Nighttime Imagery as a Tool for Global Mapping of Socioeconomic Parameters and Greenhouse Gas Emissions,
Ambio, 29(3), 157-162, 2000.
*2 Masanao Hara, Shuhei Okada, Hiroshi Yagi, Takashi Moriyama, Koji Shigehara and Yasuhiro Sugimori, Progress for Stable Artificial Lights Distribution
Extraction Accuracy and Estimation of Electric Power Consumption by mean of DMSP/OLS Nighttime Imagery, Remote Sensing and Earth Sciences, 1(1),31-42, 2004

0
1990

1992

1994

1996

1998

2000

2002

2004

Year

図14．F10,F12,F14,F15のOLS夜間モードセンサの特性

プロジェクト−2衛星データを用いた植生純一次生産量推定の精度向上のための基礎研究Ⅱ

NOAA/AVHRR GAC を用いた全球時系列データセットの作成

カラマツ林における方向別反射特性の季節変化と日変化
P2007-2, P-14

カラマツ林における方向別反射特性の季節変化と日変化〜光合成モデルに必要な樹冠構造パラメータの推定に向けて〜 松山 洋（首都大 都市環境） 代理発表：長谷川宏一
共同利用研究者

担当教員 本多嘉明・梶原康司[ はじめに ]方向別反射特性…一つの物体を多方向から放射観測（多方向放射観測）することにより，地表面の凹凸が引き起こす不均一な散乱を捉えたもの（Bi-directional Reflectance Distribution Function）
Heiskanen,（2006）: 衛星センサ：MISRのデータを用いて，BRDFが樹高や植被率の推定に有用であることを実証。
Chen et al.（2003）: HotspotとDarkspotの反射率から計算される指標HDSを用いて，植生の日向葉・日陰葉の割合を推定し，それをパラメータとした光合成モデルを作成。
長谷川ほか（2007） : 多方向放射観測から計算される指標HDSと植生指標NDVIを用いて葉面積指数を推定する手法を提案。
光学センサを用いたBRDF観測は，次世代のリモートセンシング技術として強い期待！！

（問題点）

森林を対象にした多方向放射観測には大規模な施設が必要であり，同一地点で様々な植物季節にこれを行った例は無い。
⇒実測研究が少なく,BRDFを用いたモニタリングを行う上で基礎的な知見があまりに少ない。

ビデオカメラ 分光放射計（MS-720）

BRDF観測の条件
・快晴に近い晴天時である（雲の影響を防ぐ）
・一様な森林の上に十分な高さに出る・短時間に観測角度を変え観測を行う（観測時間中の太陽高度の変化を防ぐ）

アーム部

ポール部

デジタル
傾斜計

センサ（雲台）部

RS232Cコード

本研究の目的図1 八ヶ岳カラマツ林に構築したタワー型多方向放射観測システムの概要
八ヶ岳に構築した多方向放射観測システムを利用し,各植物季節のカラマツ林の方向別反射特性を測定する。そして以下の二点を明らかにする。
①方向別反射特性の日変化：観測時の太陽高度が与える影響（モデル上でシミュレーションされているが，実測例は無い）（一日のうち複数回観測）②方向別反射特性の季節変化（まったく知見が無い）
2007年4月〜11月 計12回の観測を実施[ 結果・考察 ]多方向放射観測を行うと前方散乱方向では，植生面の陰影を多く捉え，反射率は低くなる。一方，後方散乱方向で，特に太陽天頂角と観測角度が一致するアングルでは，日向葉のみを捉え，反射率は高くなる。本発表では，BRDFの特徴を記述するため，前方散乱方向と後方散乱方向を結ぶ面（principal plane）の反射率を示す。またBRDFの特徴を表す指標として，
Lacaze et al.（2002）で提案されたHDSの日変化・季節変化（落葉期）について示す。HDSは植生を多方向から放射観測した際，最も高い反射率（Hotspotの反射率）と，最も低い反射率（Darkspotの反射率）から下の式で計算される指標である。

(ρHS−ρDS )
ρDS

HDS

HS

結果１ 方向別反射特性の日変化

25

各バンドで計算されたHDS

Reflecatance (%)

30

20
15
10
5

2.5

DS

RED
y = 0.0874x - 1.6943
R2 = 0.8024
GREEN
y = 0.0934x - 1.8693
2
R = 0.9178

2.0
1.5
1.0

NIR
RED
GREEN
BLUE

①BRDFは観測時刻に依存して変化し，観測時の太陽天頂角が小さくなればHotspotの反射率は減少し，Darkspotの反射率は増加した。
②指標HDSの値と観測時の太陽天頂角には線形の相関関係が見られた。

0.0
0
60

50

40

30

20

10

0

10

20

-10 -20 -30 -40 -50 -60

30
太陽天頂角

40

50

⇒多時期・多地点のBRDFを比較する際は，観測時の太陽高度条件をそろえるか，もしくはモデルを用いて太陽高度条件をそろえる必要がある。

60

観測角度 (°)

図 3 HDS の日変化（観測時太陽天頂角との関係）

図 2 方向別反射特性の日変化（9 月 20 日；盛夏）

結果2 方向別反射特性の季節変化
NIR
RED
GREEN
BLUE

35
30

20
15
10
5

NIR
RED
GREEN
BLUE

35
30

60

50

40

30

4.0

0.7

3.5

0.6

3.0

0.5

2.5

0.4

2.0

0.3

1.5

0.2

1.0

0.1

20
15
10
5

NDVI

0.8

10

50

40

30

20

10

0

-10 -20 -30 -40

View zenith angle (°)

10月18日 天頂角 54.0
NIR
RED
GREEN
BLUE

40
35
30
Reflecatance (%)

25

0.9

HDS（BLUE)

15

60

NIR
RED
GREEN
BLUE

HDS（GREEN）

4.5

9月20日 天頂角 45.5
35

HDS（RED）
5.0

20

0

30 20 10 0 -10 -20 -30 -40
View zenith angle (°)

HDS（NIR）

25

50

Reflecatance (%)（観測データの中から観測時の太陽高度条件が54°のものを選出し，晩夏〜晩秋にかけHDSの季節変化を検討した）

40

Reflecatance (%)

Reflecatance (%)

25

前方散乱

BLUE
y = 0.0673x - 1.2091
R2 = 0.8822

0.5

0

後方散乱

Principal
plane

Darkspotの反射率

NIR
y = 0.0496x - 1.3089
R2 = 0.9843

3.0

49.7
45.5
39.6
36.1
35.8

35

Hotspotの反射率

センサ

25

0.5

20

0

晩夏

15

紅葉前期

紅葉後期

①夏から秋にかけてHDSの値は減少していき，BRDFには季節依存性があることが分かった。
②晩秋において葉が落ちても幹や枝が作る陰影が存在するため，植生のBRDFの特徴（強い後方散乱と弱い前方散乱）は存在した。

⇒多時期・多地点のBRDFを比較する際には，観測時期を同じ季節（植物季節）にそろえる必要がある。

落葉後

図5 HDSの季節変化

10
5

0
60

50

40

30

20

10

0

-10 -20 -30 -40

View zenith angle (°)

10月30日 天頂角 54.0

0
60

50

40

30

20

10

0

-10 -20 -30 -40

View zenith angle (°)

11月8日 天頂角 53.7

図 4 カラマツ林における各植物季節（晩夏〜落葉期）の方向別反射特性[ 今後] ・特に太陽高度条件の変化に注目したBRDFモデルとの比較（Lacaze and Roujean 2001,G-function and HOt SpoT reflectance model）
・高空間分解能衛星への応用（パンクロ画像を利用してHDSを推定する）

プロジェクト−3

衛星画像を用いたオホーツク地域のエアロゾル光学特性の導出

1. はじめに

衛星画像を用いたオホーツク地域のエアロゾル光学特性の導出
Retrieval of aerosol optical thickness over Okhotsk areafrom satellite imagery
朝隈康司*)、山崎輝彬**)、久世宏明**)
Koji ASAKUMA, Teruaki YAMASAKI and Hiroaki KUZE

大気エアロゾルは時間的・空間的に変動が激しく、衛生リモー
トセンシングを用いた大気観測においては、各画素ごとの光学的厚さ（AOT）のみならず、エアロゾルモデル（波長依存性）
が重要。
また、衛星画像からAOTを求めるには地表面アルベドが既知である必要があるが、現実には１画素中に複数の物質が混在しており、その推定が困難（ミクセル問題）。
これまで、上記２つ問題を同時に解決するために、衛星画像
を反射特性が既知の物質に分類し、その分光反射率を用いて、AOTおよびモデルを推定するアルゴリズムを開発してきたが、満足できる分類結果が得られていない。

*

東京農業大学 生物産業学部
Fucluty of Bioindustry, Tokyo Univ. of Agriculture** 千葉大学 環境リモートセンシング研究センター
CEReS, Chiba Univ.

2. 研究目的と流れ

3. 利用データ
衛星画像 ：北海道地域、MODIS 1〜7 チャンネル。
期間：2007年1、3、6、7、8、9、10、11、12月
アルベドデータ： ASTER ライブラリ

衛星画像を用いて地表面被覆分類を行い、被服分類されたクラスの分光反射率を地表面アルベドと仮定し（参照アルベド）、AOTを算出するアルゴリズムの改良を行う。

既知のアルベドを任意の割合で混合しミクセルの分光反射率を作成し、参照アルベドと実際の地表面アルベドとの差を小さくする。

4.1 前処理 その１ 簡易大気補正

1/3 τ550 = 0.05

3/22 τ550 = 0.07

6/3 τ550 = 0.22

6/5 τ550 = 0.33

7/3 τ550 = 0.21

7/30 τ550 = 0.06

8/19 τ550 = 0.11

9/2 τ550 = 0.16

10/2 τ550 = 0.11

11/5 τ550 = 0.07

12/7 τ550 = 0.14

4.2 分類クラス（標準物質）

大気上端のアルベドは、大気効果により短波長ほど地表面
アルベドとの差が大きい。
分類クラスのアルベドとの差も同様。
6Sの標準モデル (Continental)を用いて簡易的に大気補正を施してから分類。

粗分類項目として被覆を４つにグループ化
1.水、雪氷，
100.0

2.植生，
グループ1

3.土壌，
80.0

10/2

12/7

60.0

Reflectance (%)

6/5

40.0

0.0
0.40

apparent image at top of atmosphere

0.60

0.80

1.00
1.20
Wave length (μm)

100.0
80.0

simple atmospherically corrected image

60.0

1.40

Dry Grass

60.0
40.0

0.0
0.40

1.60

0.60

1.00
1.20
Wave length (μm)

1.40

1.60

グループ4

80.0

20.0

0.80

100.0

グループ3

Dark silty loam
Broun silty loam
Gray silty clay
Dark clay loam
Grayish brown loam

40.0

0.0
0.40

Grass
Decidous

20.0

20.0

Reflectance (%)

1/3

グループ2

Conifer

80.0
Pure Ice
Tap Water
Coarse Snow
Medium Snow
Fine Snow
Frost

4.人工物

100.0

Reflectance (%)

2.

衛星の各観測波長と既知のアルベド（ASTERライブラリ等を利用）をマッチングすることにより被覆分類をおこなっていたが、粗分類と詳細分類の２段階分類を取り入れ精度向上を目指す。

Reflectance (%)

1.

TarAsphalt
Asphalt shingle
Asphalt concrete
Concrete
Teracotta Tile

60.0
40.0
20.0

0.60

0.80

1.00

1.20

Wave length (μm)

1.40

1.60

0.0
0.40

0.60

0.80

1.00

1.20

Wave length (μm)

1.40

1.60

5.1 分類方法

4.3 分類クラス（ミクセル）

植生指標(NDVI:VI)、土壌指標(NDSI:SI)、水指標(NDWI:SI)
を用いて粗分類。
今回は、簡易大気補正画像と分類結果を目視で比較、年間を通して尤もらしい結果になるようにルール化。

典型的な物質2種の各波長の反射率を、10%ずつ混合し、ミクセルアルベドを作成。
作成されたミクセルアルベドはそれぞれのグループに配置。

Grass
20%
40%
60%
80%
Dark silty loam

60.0

Dark silty loam
20%
40%
60%
80%
Frost

80.0
Reflectance (%)

80.0
Reflectance (%)

WI > -0.07 OR VI < 0.08→
Group≠1 AND VI > 0.35 →Group≠1, 2 AND SI > -0.2 →Group = 3 AND VI < 0.16 →いずれにも該当しないもの →

100.0

100.0

40.0

60.0
40.0
20.0

20.0

0.0
0.40

0.0
0.4

0.6

0.8

1.0
1.2
Wave length (μm)

1.4

1.6

Mixcell of Glass and Dark silty loam
→ グループ 2 と 3 に配置

0.60

0.80

1.00
1.20
Wave length (μm)

1.40

粗分類によりグループされたグループの中で最小二乗法を用いてアルベドを抽出。ただし、グループ 0 は、全てのクラス
を用いた。

1.60

Mixcell of Dark silty loam and Frost→ グループ 1 と 3 に配置

6.2 詳細分類結果①

6.1 粗分類結果
1/3

Group 1 （水、雪氷）
Group 2 （植生）Group 3 （土壌）Group 4 （人工物）
Group 0 （その他）

6/5

10/2

12/7
1/3

3/22

6/5

simple atmospherically corrected image
simple atmospherically corrected image

composed image obtained with SI assigned to R, VI to G and WI to B.

classified image

rough grouped image. Blue indicate the group 1 (water, ice and snow), Green shows the gp. 2(vegetation), Brown the gp. 3 (soil), Gray the gp.4 (urban) and Red shows others.

6.2 詳細分類結果②
8/19

7.1 AOTの導出
10/2

11/5

1.

2.

simple atmospherically corrected image

3.
4.
5.
6.

大気上端のアルベドを、AOTを変化させながら大気補正して得られるアルベドとのテーブル（ρ-τ テーブル）を作成する。
ρ-τテーブル に、大気上端のアルベドと参照アルベドに一致するAOTを求める。
continental モデルをベースにして、AOTは0から2まで変化させた。
AOTが負になる場合は、maritime、urban の順にそれぞれのテーブルに当てはめる。
いずれのテーブル内に組み合わせが無い場合は、エラーとする。
AOTが2を超えるものは、2として表示。

classified image

粒径分布を変化させ、さらに多くのテーブルを作成する方法を昨年紹介
しましたが、今回は時間の都合で省略。

1.00
0.90
0.80
0.70
0.60
0.50
0.40
0.30
0.20
0.10
0.00

atmospherically corrected albedo

atmospherically corrected albedo

7.2 AOT導出テーブルの例

0.0

0.5

1.0

1.5

2.0

aerosol optical thickness

atmospherically corrected albedo

ρ-τ table with continental aerosol model

7.3 AOT導出結果①
1/3

1.00
0.90
0.80
0.70
0.60
0.50
0.40
0.30
0.20
0.10
0.00

τ=0.0
3/22

1.0

2.0

6/5

simple atmospherically corrected image
0.0

0.5

1.0

1.5

2.0

aerosol optical thickness

ρ-τ table with maritime aerosol model

1.00
0.90
0.80
0.70
0.60

AOT map from MODIS ch.4

0.50
0.40
0.30
0.20
0.10
0.00
0.0

0.5

1.0

1.5

2.0

aerosol optical thickness

ρ-τ table with urban aerosol model

AOT map from MODIS ch.3

7.3 AOT導出結果②
8/19

τ=0.0
10/2

1.0

2.0

11/5

7.4 参考 オホーツク地域（網走周辺）拡大①
2006/8/26 without rough grouping algorithm and with only ASER library.

simple atmospherically corrected image

2007/8/19 with rough grouping algorithm and mixed spectrum albedo.
AOT map from MODIS ch. 4

atmospherically corrected image

AOT map from MODIS ch. 3

7.4 参考 オホーツク地域（網走周辺）拡大②
2006/9/26 without rough grouping algorithm and with only ASER library.

2007/10/2 with rough grouping algorithm and mixed spectrum albedo.

atmospherically corrected image

classified image

AOT-map from MODIS ch. 3

classified image

AOT-map from MODIS ch. 3

8. 考察
雲のアルベドデータを取り扱わなかったため、雪として分類されている。
海と比較し、地上と海上で異なる傾向を示す。これは地表面
反射率の推定精度に起因し、以下の傾向と一致する。
参照アルベドの大小によって以下の傾向となる。これは、放射伝達計算に用いた6Sに起因し、長波長ほど顕著になる。
大気上端のアルベドが高い場合、大気補正後のアルベドはより高い値になり、大気上端のアルベドが低い場合、大気補正後のアルベドはより低い値になる （urban は 0.09、continental は0.21、maritime は、0.043 が境界）。
モデルによっては、AOTが複数存在してしまうため、周辺画素などからの推定が必要となる。

大気上端のアルベドと参照アルベドの差が大きい場合は、AOTが負になるため、衛星画象の観測値から主成分分析などをおこない、補正可能範囲のミクセル反射率を推定する必要がある。

1.1 参考 網走地域の被覆分類とAOT-map
2006/8/26

2006/9/26

atmospherically corrected image

classified image

aot-map from MODIS ch. 3

多波長マイクロ波放射計データを用いた水物質リトリーバルの研究
Feb. 28, 2008

２００７年度CEReS
２００７年度CEReS 共同利用研究会

本研究の目的Observed
TBs

CRM forecasts

多波長マイクロ波放射計データを用いた
水物質リトリーバルの研究：
マイクロ波放射計データの非静力雲解像モデルへの同化法の開発

青梨和正（気象研）
担当教官 高村民雄

データ同化

Pr( X |Y o )Pr(Y o | X ) Pr( X ) / Pr(Y o )CRM物理量からMWR輝度温度を計算する
CRM物理量のPDFを調べる

Other
Observations
Analysis data

Outline目標
MWR輝度温度の観測オペレータ
CRM物理量の予報誤差と共分散を調べる
MWR輝度温度の観測誤差を調べる
EnKFでMWR輝度温度を使った実験
モデルシミュレーション結果からRTM計算した
TBを観測値としてつかう実験
実際の観測値をつかう実験

はじめに
MWR輝度温度の観測オペレータ
MWR輝度温度の観測オペレータ
劉氏RTM

アンサンブル予報によるCRM
物理量のPDF
PDF
アンサンブル予報によるCRM物理量のアンサンブル予報の方法
物理量のPDF
物理量のPDF

同化スキムの開発

まとめ

地上型MWR

MWR輝度温度の観測オペレータ

劉氏RTM

観測周波数 （GHｚ）
22.235, 230.35, 23.835,
26.235, 30.0,
51.250, 52.280, 53.850,
54.940, 56.660, 57.290,
58.800時間的に高頻度、但し降水時には高温度

Infrared Imager

300

10μm

Az scan at 15 deg el

240

Rain flag

220

12 per. Mov. Avg. (Az scan
at 15 deg el)

Cloud Particles

Microwave Radiometer
Radiation from Rain
Scattering by Frozen
Particles

Cloud Top Temp.
Scattering

Frozen Precip.
Snow Aggregates

200

８５GHz

１９GHｚ

Radar

180

Melting Layer

Radiation

160

2cm

0℃ Back scattering from Precip.

Rain Drops

140
6/14/03
22:00

6/14/03
16:00

6/14/03
10:00

6/14/03
4:00

6/13/03
22:00

120
6/13/03
16:00

Tb (K)

280
260

3 mm-3cm（１００−１０GHｚ）

Satellite Observation (TRMM)

放射伝達計算モデル
dTB ( , , )
TB (1
d
0

4
where

P( , , ,

'

,

'

cos ,

0 )T ()TB ( ,
K ab',

RTM計算に用いたCRMの物理量（０４年6月8日06UTC
06UTC FT=11ｈ）
温位)')d ' dK sc dz,

0SST, Winds

'

降水強度

K sc /( K ab

K sc ),

P is phase function

Radiative Transfer Code (Liu,1998)
One-dimensional model (Plane-parallel)
Mie Scattering (Sphere)
4 stream approximation
Calculate TBs for homogeneous, convective &stratiform precip with each precip types.

RTM計算結果
３０ GHｚ５３.850 GHz

２２．２３５ GHｚ

56.660 GHz

水蒸気

水蒸気プロファイルとTB22TB22-30の比較

気温プロファイルとTB53TB53-58の比較

CRM物理量のアンサンブル予報誤差分
CRM物理量のアンサンブル予報誤差分
布と共分散用いたCRM
用いたCRM

アンサンブル予報の方法（Mitchell et al, 2002)

JMANHM（Saito et al,2001)• Resolution:
5 km• Grids: 400 x 400 x 38
• Time interval: 15 s

100 members started with perturbed initial data:

X i (t0 )X c (t0 )Xi

Geostrophically-balanced perturbation Xi
3D perturbation in stream function
XY : 2D random fields with correlation

correlation(r )(1 (r / R)) * exp( (r / R ))R=1000,300,100,30,10 km (20 members)
Z : sin (mπz*/H) m=1,3,5 (34, 33, 33 members)
to perturbation in (θ, u, v, Ps), assuming geostrophic
and hydrostatic balance.
Total energy of perturbation
SD of PT ~ 1 K

Ｗ，Q
mm/h）
）
Ｗ，Qvの頻度分布(の頻度分布(全域,
全域,降水強度>3
降水強度>3mm/h
CRM.heito.20050520 ENS12(03/1/26/21UTC FT=4h)

正規化したＷ，QRの頻度分布(全域）
CRM.heito.20050520 ENS12(03/1/26/21UTC FT=4h)

W

1200

1000

GAUSS
W0000005

600

W0000010
400
W0000015

降水強度> 3 mm/h
W

全域

800

W

1200

1200

1000

1000

800

800
GAUS S

200

GAUS S

W0000020
600

0-4 .0 0 -3.00 -2.00 -1.00

W0000005

W0000025
.00

1.00

2.00

3. 00

600

W 0000005

W0000010

W 0000010

400

4.00

400
W0000015

DEV

200
0
-4.00 -3.00 -2.00 -1.00

Qr

W 0000015
200

W0000020
W0000025
.00

1.00

2.00

3.00

4.00

DEV

Qv

3000
2500

W 0000020

0-4.00 -3.00 -2.00 -1.00

W 0000025
.00

1.00

2.00

3.00

4.00

DEV

1200
1000

GAUSS

800

QV000005

600

QV000010

400

QV000015

200

QV000020

2000

Qv

1200

1000

800
GAUS S

1500

GAUSS

600

QV000005
QV000010

400

QR000005

1000

QR000010

500

QV000025

0
-4.00

-2.00
-3.00

QR000015

0
-4.00

-2.00
-3.00

DEV

.00
-1.00

2.00
1.0 0

4.00
3.00

DEV

.00
-1.00

2.00
1.00

4.00
3.00

QV000015
200

QV000020

0QV000025-4.00 -3.00 -2.00 -1.00

DEV

.00

1.00

2.00

3.00

4.00

W対Z(CFD) @ 地点Ａ地点Ａ-D(但し高度˜0.5,1,2,3,4
但し高度˜0.5,1,2,3,4 km)

対象とする事例
寒冷前線（０３／１／２７）

寒気内低気圧（０３／１／２９）

Z

CDF
5

1
exp(x2 / 2)dx
2
5

4

A

4

3

3

W18

W18

ZCDF
2

ZCDF
2

W16

W16

ZCDF
1

B

ZCDF
1

W13

W13

ZCDF

ZCDF

0

0
W9

W9

ZCDF

-1

ZCDF

-1

W7

低気圧中心
（０３／１／２７）

梅雨前線（０４／６／１）

W7

ZCDF

-2
-3

-2

-1

0

1

2

3

ZCDF

-2
-3

4

5.0

-2

-1

0

1

2

3

4

8
7

4.0
6
3.0

W18

W18

5

ZCDF
W16

1.0

C

3

ZCDF

W13

2

W13

W9

-2.0
-3

-2

-1

0

1

2

3

W7

-1

ZCDF

-2
-3

4

.005

180

ZCDF

QS16

ZCDF
QR13

.001

120

ZCDF

100

ZCDF

QR9

ZCDF

R HW216R HW213

80

ZCDF

QS9
0.000

R HW218

QR16
.002

QS13
ZCDF

QS7

60

-2

-1

0

1

2

3

4

.0005

ZCDF
-2

-1

0

1

2

3

4

.005

.0004

ZCDF

QR18
.003

QR16
.002

ZCDF

QR13
.001

ZCDF

-1

0

1

2

3

4

ZCDF
RHW216
ZCDF
90

RHW213
ZCDF
RHW29
ZCDF

70

RHW27
60
-3

ZCDF
-2

-1

0

1

2

3

4

140

120

100

R HW218
ZCDF

90

RHW218100

ZCDF

R HW216

ZCDF

QR13
.0001

-2

ZCDF

QR16
.0002

ZCDF

110

.004

ZCDF

4

RHW218100

R HW27
20
-3

120

QR18
.0003

110

R HW29

QR7
-.001
-3

3

80

ZCDF

0.000

ZCDF

2

ZCDF

40
-.001
-3

1

120

140

.001

ZCDF
0

130

160

QR18
.003

ZCDF

W7
-1

RHW2対Z(CDF)＠地点A-D

.004

.002

ZCDF

-2

(但し高度˜0.5,1,2,3,4
但し高度˜0.5,1,2,3,4 km)

QS18
ZCDF

W9
0

ZCDF

-1.0

(但し高度˜0.5,1,2,3,4
但し高度˜0.5,1,2,3,4 km)

.003

ZCDF

1

0.0

降水物理量対Z(CDF)＠地点A-D

.004

W16

ZCDF

ZCDF

D

.005

ZCDF
4

2.0

ZCDF
80

RHW21680

ZCDF

R HW213

ZCDF

ZCDF

RHW213
60

ZCDF

70
QR9

QR9

ZCDF

0.0000

ZCDF

0.000

QR7
-.0001
-3

-1

0

1

2

3

ZCDF

60

QR7

ZCDF
-2

R HW29

-.001
-3

4

ZCDF
-2

-1

0

1

2

3

4

50
-3

ZCDF
-2

-1

0

1

2

3

4

W> 0.3 m

W 対 PWC

全体

40

70
60

38
50
40

36

30
34

20
10

W

R 2 乗 = 0.6633
0.0

.5

1.0

1.5

2.0

2.5

3.0

PWC

PRECIP

32
0-10-1.0 -.530- 1.0

R 2 乗 = 0.4910
-.5

0.0

.5

1.0

1.5

2.0

2.5

ZCDF

R HW27

RHW27
20
-3

ZCDF
-2

-1

0

1

2

3

4

W対PWCの相関空間分布＠地点B(高度〜１ｋｍ）

W対降水強度,PWCの共分散
地点B (高度〜１ｋｍ)
W 対 降水強度

RHW29
40

3.0

W

300 km

W<= 0.3 m

Ensemble Kalman Filter (EnKF(EnKF))
Kalman Filter：
where K

同化スキムの開発

Xa

Xf

K (Y oK

P f H t ( HP f H t

HX f )R) 1

EnKF: use ensemble to estimate forecast
error covariance P f H t HP f H :
N

1

Pf Ht

N 1n 1( X nf

X nf ) * ( H ( X nf ) H ( X nf ))N

HP f H t

解析実験（Gaussian Sum Filter)
CNTL解析： 全物理量が正規分布すると仮定

1( H ( X nf ) H ( X nf )) * ( H ( X nf ) H ( X nf ))N 1n 1

アンサンブル平方根フィルター (EnSRF)
Snyder & Zhang(2003)
アンサンブルで予報誤差共分散を近似する
P f Ht

実験：鉛直速度が混合正規分布
W> < 0.3 m/sで２つのクラスに分ける

X
Pf Ht
wl

wl X

l

wl (Pf

1 N( X nf
N 1n 1
N

1

HP f H t

N 1n 1

X nf ) * (H ( X nf ) H ( X nf ))( H ( X nf ) H ( X nf )) * (H ( X nf ) H ( X nf ))アンサンブル平均と各メンバーの解析
l

(X

l

Xa

X na

l

X)(H(X ) H(X))t )lexp( 0.5(Yo H(X ))**2/ 2)

各クラスの解析値をwlの重みをつけ平均する解析実験の対象事例
Upper cold low case(Jan. 28-29,2003)
Ensemble forecast(21 UTC Jan.28, 2003)(100 members, sameas the previous section)
Assimilate simulated
precipitation data
(01 UTC Jan. 29,
2003)

X na

Simulated Precip datafor 03/1/29/01UTC:
CRM control +perturbation (1/28/21UTC)＋(1/28/21UTC)＋
Noise（Noise（50%）
50%）

Xa
Xn

Xf
a

WK (Y oXn
(1f

H ( X f ))WK H X n f ,R /( HP f H tR) ) 1

鉛直速度のアンサンブル平均（lev=2110m））
鉛直速度のアンサンブル平均（lev=2110m(01UTC Jan. 29, 2003)
CNTL解析（正規）真値

解析前

実験（混合正規）鉛直速度のアンサンブルのスプレッド（lev=2110m） (01UTC Jan. 29, 2003)
lev=2110m）

CNTL解析（正規）降水域では分散の大きなクラスの重みが大きくなる効果

実験(混合正規）
5

4

3

W 18
ZCDF

2

W 16
ZCDF

1

W 13
ZCDF

0
W9
ZCDF

-1

W7
-2
-3

まとめ
ＭＷＲ輝度温度の観測オペレータを開発した。
物理量PDFは混合正規分布で近似できる
鉛直速度と、PWCとの共分散がクラス毎（w><0.3m/s）に異なる。
混合正規分布を仮定した解析アルゴリズ
ムの開発

ZCDF
-2

-1

0

1

2

3

4

船舶搭載用新型スカイラジオメータの開発
平成１９年度千葉大学CEReS共同利用研究発表会, 2008.2.28（南極観測研究計画）

船舶搭載用新型スカイラジオメータの開発

スカイラジオメータ観測による広域洋上エアロゾルのDevelopment of a new shipship-board skysky-radiometer

光学特性の長期的モニタリング
研究目的
・エアロゾルの直接効果による気候影響の定量的評価
・広域のエアロゾルの光学特性の実態把握とその空間的時間的変動特性
・衛星による広域のエアロゾル・リモートセンシングの地上検証データ

研究方法

塩原匡貴，新船舶搭載スカイラジオメータ開発チーム＊

観測： 船舶搭載スカイラジオメータによる太陽直達光天空光分光観測
解析： SKYRAD.pack 解析スキーム使用*開発チーム： （極地研）塩原匡貴・矢吹正教，（プリード）笹本和敏・阿部泰三・江井和則・河井和弘，（気象研）内山明博，（東大）中島映至，（東京海洋大）村山利幸，（富山大）青木一真，（山梨大）小林拓，（エナジシェアリング)村治能孝・山野牧（千葉大学CEReS共同利用研究）

11/14
4/13

プロジェクト研究P-3：衛星データと地上観測ネットワークによる放射収支の評価と大気パラメータの長期変動

11/21

11/28

エアロゾルの光学特性に関する観測的研究

3/21

12/23
2/15

研究目的

2000-2001 (JARE42)

2001-2002 (JARE43)

2003-2004 (JARE45)

2004-2005 (JARE46)

2002-2003 (JARE44)

エアロゾルの光学的厚さを含む光学特性の広域分布について、船上
スカイラジオメータ観測に基づく解析を中心に行い、衛星観測による
海上のエアロゾル光学特性との相互比較を通して、衛星データ解析
アルゴリズムの高精度化に資するための地上検証データを提供する

本年度の研究計画
・新たに開発した船舶搭載スカイラジオメータの性能評価
・国内訓練航海および南極航海での実観測

しらせ南極航海航路上の観測データ取得状況

Inversion results for 17-19 Nov. (top panels) and 23-24 Nov. (bottom panels)

Anthropogenic

2000.11.17-19
11/17
11/18

Maritime

11/19

0= 0.94 – 1.0

Re= 1.35 – 1.37
Im= 0.0 – 0.001

dV/dlnr (cm3/cm2)0= 0.84 – 0.96
Im= 0.002 – 0.011

10-5

2000.11.17-19

Scale height = 2.0 km

Re

1E-1 0

1.45
-0.02
1

1.40

10 -9

Re= 1.50 – 1.58
Im= 0.0002 – 0.011

11/24

Mixed
0= 0.72 – 1.0

Ship tracks of the JARE 42 cruise in 2000 and the air trajectories
Aerosol optical properties in boxes were measured by Yabuki et al. (2003)(modified after Yabuki et al., JMSJ, 2003)

dV/dlnr (cm3/cm2)0= 0.76 – 0.99

10-5

0.6

10

0.1
1
Radius ( m)1 E-1 0

-0.03

1.35
400 500 675 870 1020
Wavelength (nm)102000.11.23-24

Biomass burning

0.8

-0.01
1E-1 2

1E-1 30. 1

0.01

2000.11.17-19
1.0

1E-1 1

10 -7

effective size range

11/23

2000.11.17-19
0.0

Re

0.4

SSA
400 500 675
870 1020
Wavelength (nm)

400 500 675 870 1020
Wavelength (nm)2000.11.23-24

Scale height = 2.0 km

Im

2000.11.23-24

2000.11.23-24
1.0

0.0

1 E-1 1

-0.01
10 -7

0.8

1 E-1 2

1.45
-0.02

1 E-1 3
0 .1

1

10

0.1
1
Radius ( m)

-0.03

1.35

effective size range

0.01

0.6

1.40

10-9

10

Volume size distribution

Im

0.4

SSA

400 500 675
870 1020
Wavelength (nm)

400 500 675
870 1020
Wavelength (nm)

400 500 675 870 1020
Wavelength (nm)Real part of refractive
index

Imaginary part ofrefractive index

Single scattering albedo

新船舶搭載スカイラジオメータ POM-01 Mk III の開発（プリード・極地研 産学共同研究）

後継船搭載エアロゾル観測装置の開発（１）スカイラジオメータ

現有機に対する主な変更改良点

開発目的
太陽位置検出器
ＣＣＤ → ＰＳＤ

１） 第51次以降の船上モニタリング観測を念頭に置き、しらせ後継船を想定した次世代船舶搭載型スカイラジオメータを開発する。
２） 現有のスカイラジオメータについて問題が明らかとなった追尾性能や太陽位置検出精度の向上により、船舶動揺によるデータ不良および解析不能の問題を解決する。
３） これにより、毎年実施される西太平洋〜インド洋〜南極海にいたる南極航海において、エアロゾルの広域分布特性と光学特性およびその経
年変化の実態把握に資する。

姿勢制御用駆動装置
パルスモータ → ＤＣサーボ

目標制御用ＣＰＵ
単体共有 → 独立複合船体姿勢モニター
機械式傾斜センサ → 精密小型ジャイロ

１） しらせ級の船舶搭載を想定し、ピッチ角ロール角が５度までの船体動揺に対して太陽を完璧に追尾すること
２） 天空散乱光測定では測器高度角を制御することにより、等太陽高度ス
キャンを達成すること

しらせ動揺特性解析−JARE45-46の航海データ

新南極観測船“しらせ”の概要

10
9.5

Rolling

PERIOD (sec)

9

N.Pacific

8.5
8Pitching

7.5
7
6.5

Tropics

10N

6
2

AMPLITUDE (deg)Indian/S.Pacific
Antarctic

N.Pacific

10S

Tropics

55S

Indian/S.P.

Antarctic

Rolling
1.5

1Pitching
0.5

10N

0

N.Pacific

10S

Tropics

55S

Indian/S.P. Antarctic（ユニバーサル造船株式会社ホームページより）

しらせ船上での直達光観測データと動揺データ

太陽を光源とした立体視野角のキャリブレーション測定(POM-01 MK III, 500nm ch)

SKY-RADIOMETER POM-01 Mk III, SENDAI BAY, 2007/09/03
9

SIG. OUT VOLTAGE (a.u.)

0.3°

10

2.50E-04
2.25E-04
2.00E-04
1.75E-04
1.50E-04

1.00E-04

6

500nm

5
4
3

675nm

2

0
2.0

太陽周辺を 0.1°× 0.1°のメッシュ状にスキャン

5.00E-05
S15

2.50E-05

21

17

19

13

15

9

11

5

7

3

S8
S1

solid view angle = 2.54E-4 sr

1.5

ANGLE (deg)

7.50E-05

1

400nm

7

1

1.25E-04

0.00E+00

8

Rolling

1.0
0.5
0.0
-0.5

Pitching

-1.0
-1.5
10:30

10:35

10:40

10:45

10:50

TIME (LT, 2007/09/03)

10:55

11:00

11:05

仙台湾航行中の“しらせ”の動揺特性
4

Pitching

PITCHING ANGLE (deg)

3
2
1
0
-1
-2
-3
-4
4

ave. period = 8.6 sec
Rolling

3

ROLLING ANGLE (deg)

擬似動揺装置による太陽追尾性能試験

2
1
0
-1
-2
-3

ave. period = 10.7 sec-4
10:22:00

10:23:00

10:24:00

10:25:00

10:26:00

10:27:00

10秒周期、動揺角度±５°での追尾試験観測

TIME (LT, 2007/09/03)

スカイラジオメータの取付予定位置

まとめ1) 新南極観測船（新しらせ）による広域の洋上エアロゾルの光学特性観測のために、新型の船舶搭載スカイラジオメータ（POM-01 MK III）を開発製作した。
2) 開発にあたり、旧型スカイラジオメータ（ POM-01 MK II ）では不十分だった太陽追尾性能の向上のため、パルスモーターからＤＣサーボモーターへの変更および単一ＣＰＵによる制御から複合ＣＰＵによる制御への変更改良を施した。また、太陽周辺光の観測精度の向上のため、狭視野センサーをＣＣＤカメラ方式からＰＳＤ２次元センサー方式に変更
し、測定高度角制御のため応答の速いジャイロ方式の姿勢センサーを採用した。

取付予定位置での視界障害領域

3) これにより、本年9月の国内巡航期間中の試験観測では、ピッチ角、ロール角が±１〜２度程度の船体動揺に対しては、良好に太陽追尾を行い、安定した太陽直達光の観測が可能であることが確認された。
4) 擬似動揺装置を用いた地上試験では、同装置の不具合により、十分な試験が行われず、本機の性能評価に至らなかった。

（船首）今後の課題と目標改良した擬似動揺装置により、本スカイラジオメータの追尾性能評価を行う。
しらせ級舶船での洋上観測を可能にするため、動揺角±５度、動揺周期１０秒のもとでの太陽直達光観測および天空光観測を達成する。

大気放射観測による父島のエアロゾルの光学的特性

北西太平洋域における下層雲の衛星リモートセンシングの
航空機検証観測2007 年春季黄砂と視程映像観測
平成19年度CEReS共同利用研究報告会，2008.2.28

2007年春季黄砂と視程映像観測

定点自動映像観測目的

木下紀正 (鹿児島大学 産学官連携推進機構)
永松哲郎 (鹿児島大学 水産学部)
飯野直子 (熊本大学 教育学部)

火山噴煙・黄砂状態・大陸からの汚染気塊・赤潮

足立祐樹 (鹿児島大学 工学部)

方法
インターバル撮影（デジタルカメラ・監視ビデオ・パソコンWebCam）

可視・近赤外映像
広角−手動多方向撮影

視程の指標点 Ground Positioning Points

観測定点 B :鴨池港近傍建物８階
Height = 26 m above sea level

観測定点からの距離
Distance from the Observatory

多方向撮影による広角展望

手動, 近赤外/可視 DSC-V1, V3
Wide view by multi-directional phtoes

2007.3.28-31 3.28 九州-東北の広範囲で黄砂

インターバル撮影
3h毎 Ricoh Caplio G4wide, Horizontal 65゜

衛星画像 NOAA/AVHRR
2007.3.28̲13:57JST 鹿児島大受信

MT-Sat/VISSR AVI画像
3.30-31

VISSR
AVI

07.3.31_01

MT-Sat/VISSR AVI画像
07.4.1-2

07.4. 1

07.4.2

全国63カ所で黄砂

9 JST

視程は鹿児島で 3km

昨日に引き続き、沖縄から東北までの
広い範囲
関東でも横浜、東京、千葉で黄砂
気象人HP

下関で4km、
新潟、富山、金沢、輪島、福井、豊岡、松江、鳥取、相川、長崎、熊本で 5km など、
西日本や北陸で 特に見通しが悪くなった。

韓国・ソウルでも 黄砂が増加

AD-Net 情報 4.1-2

1日午前10時現在、韓国の広範囲に黄砂警報

気象人HP

07.5.25-27
4.1.
4.2
dusty

07.4.1-15

黄砂, 光化学スモッグ

07.5.25 黄砂に関する全般気象情報 第4号
25日〜27日にかけて、 西日本〜東北の広範囲で黄砂が予想される。
25日15時、すでに厳原で黄砂を観測している。
07.5.26 黄砂 沖縄〜東北北部までの75地点で5月下旬としては大規模（今年最大規模）、
07.5.27■黄砂：富山、松本、軽井沢、長野、山口、大島、
御前崎、尾鷲、東京、横浜■光化学スモッグ：
山口、福岡、長崎、熊本に光化学スモッグ注意報。
北九州市は08:45に市内全域に発令し、運動会を予定していた小学校85校はすべて延期に。
気象人HP

MT-SatＩｎｔｅｒｖａｌ

07.5.26 09-15 JST

Ricoh 07.5.26-31

視程観測のまとめと問題点
1. 鹿児島における視程 (km) 4 - 55以上
2. 近赤外と可視の関係： 一長一短あり
3. 標高の影響, デジタル解析など

標高

標高の影響

ミリ波雲レーダ FALCON による観測−雲物理量導出と最近の結果

Ashra観測の現状報告

Ashra Collaboration
Presented by 長南 勉 (ICRR)

Ashra (All-sky Survey High Resolution Air-shower detector)
高解像度監視型望遠鏡による超高エネルギー宇宙線観測► AshraとはAshraとは► Ashraサイト全景
Ashraサイト全景

宇宙線による２種の大気発光＋可視光による広視野「同時」監視観測

大気蛍光
空気シャワー
チェレンコフ光

Mauna Kea

Ashraサイト @ Mauna Loa (3300m)

・大気蛍光・・・全方向に広がる光 （図中 紫）・チェレンコフ光・・・入射粒子と、生成された2次粒子の軌
跡 （図中 赤）の方向に指向性を持った光 （図中 青）
Ashra集光器群

2つの発光を通して宇宙線の起源となる天体について調べ、宇宙創成の謎に迫る。
超新星爆発

活動銀河核(AGN)

ガンマ線バースト(GRB)

加速された陽子

p
γ

n
π+

p

p

γ

π0

μ+
νμ

γ
γ

e+
̲ νe
νμ

► Ashra集光器
Ashra集光器► Ashraの現在と今後
Ashraの現在と今後取得された画像例

Ashraによる観測が立ち上がりつつある。

直径42度の広視野を達成
大口径イメージインテンシファイア
入射した画像の輝度増幅と画像収縮を行う。画像をそのまま増幅する点がPMTとは決定的に異なる。
チェレンコフ光と大気蛍光を効率よく取得するため、入力面には誘電多層膜
フィルタをコーティングしている。

最終的には、トリガーによる宇宙線撮像を行い、可視光と合わせた観測に移行する。

北斗七星
full-size@Mauna Loa

まずは可視光による監視観測からスタート
し、突発天体の撮像を
目指す。

4秒露光

Ashra観測開始に伴い、大気モニターが必要になる

1台で42度径の視野を持ち、視野内で分角 Ashraにおいて最も鍵となるの精度を持った集光器。これを並べることで、デバイスの開発に成功
全天の約80%を監視する。

イメージングライダーによる大気モニター
突発的な天体現象のタイムゼロを観測可能

TEMPORAL AND SPATIAL VARIATION OF CLOUDMEASURED WITH A PORTABLE AUTOMATED LIDARTATSUO SHIINA,
Graduate School of Advanced Integration Science, Chiba University
1-33 Yayoi-cho, Inage-ku, Chiba-shi, 263-8522, JAPANNOBUO TAKEUCHI, NAOHIRO MANAGO, HIROAKI KUZE,
Center for Environmental Remote Sensing, Chiba University
1-33 Yayoi-cho, Inage-ku, Chiba-shi, 263-8522, JAPANSUEKAZU NAITO
Chiba Prefecture Environmental Research Center
1-8-8 Iwasakinishi, Ichikawa-shi, 290-0046, JAPANA portable automated lidar (PAL) system, which conducts full-time operation and all-weather observation through the window, has beendeveloped. Observations of long-term temporal and spatial dynamics of the atmosphere are described and the advantage of full-time operation is
discussed.

1. Introduction
Atmospheric convection has an effect on cloud formation, and it leads to heavy rain or lighting strike. It also affects diffusionof suspended substances. Changes in the atmospheric convection due to the climate change may influence the large and localscale transformation of particles such as the yellow sand. In this context, it is essential to understand the temporal and spatial
dynamics of the atmosphere, which cannot be monitored with conventional, fixed-point observation systems or
meteorological satellites.
Lidar is an appropriate tool for monitoring time and spatial dynamics of the atmosphere, especially aerosols andclouds. Although various kinds of lidar systems have so far been developed, observations are limited in terms of time spans.
Besides, observation directions are usually fixed both horizontally and vertically. These limitations are originated from thesystem stability as well as the complication of system maintenance including the laser device.
A micro pulse lidar (MPL), developed by Spinhirne in 1993, is a compact lidar system that provides easy operation andlong-term observation [1]. Using a laser-diode pumped laser of micro joule output energy, MPL ensures the eye-safety
features. Signal-to-noise ratio was improved by narrowing the receiver field-of-view (FOV). However, this makes it difficult
to adjust the laser beam within the receiver’s FOV. Since the same telescope is used to both transmit and receive the laser
beam, a small amount of the emitted beam back-reflected from the beam splitter often damages the detector.
In this paper, we describe a portable automated lidar (PAL) system, which we have developed to conduct full-time
operation and all-weather observation through the laboratory window [2]-[4]. The PAL system has an automated correction
mechanism for misalignment of the overlap between the transmitted laser beam and the receiver FOV. Hence the system isable to operate in a stable and stand-alone way. In addition, we have recently installed the scanning mechanism by attachinga horizontal stage to the PAL system. This improvement contributes greatly to monitoring the two-dimensional structure of
the atmosphere nearly instantaneously.

2. PAL system
The PAL system is a variation of micro-pulse lidar (MPL) system. The system configuration is shown in Fig. 1 and its
specifications are summarized in Table 1. Since the transmitted energy is 15 J, the system is nearly eye-safe at the expenseof weak signals (lidar echo). To attain enough signal-to-noise ratio, the background light due to sky radiance must be
eliminated with a narrow-bandwidth filter (0.5 nm) and a narrow FOV of 0.2 mrad. At the same time it is essential to keep
the good overlap between the laser beam and the telescope FOV. Misalignment of the overlap, however, sometimes occurs
from changes in the ambient temperature and accidental disturbances. The system has the auto alignment mechanism, inwhich the laser beam is scanned vertically and then horizontally within the receiver’s FOV and the maximum in the return
signal (a certain range near the peak of the A-scope) is sought every 15 min.
The detector is a photo-multiplier operated in the photon counting mode (Hamamatsu photonics K.K. R1924P). The
lidar echoes are accumulated by a scaler (Stanford Research Systems SR430). The spatial resolution is 24 m and the
maximum observation range is 24 km (altitude 15km). The observation is made though the vertical window of the laboratory,

leading to the capability of measurement under all weather conditions. The observation data have been accumulated since theyear 2004. The system status can be checked and the data can be downloaded though the Internet.
A built-in rotation stage for horizontal scanning has recently been installed. As the PAL system is fabricated as amonolithic structure including a laser head and a detector, the rotation stage was “inserted” under all the optical systems. The
scanning observation of a range of 25 degrees is conducted every hour, interrupting the continuous measurement for about
6 min. The PAL system is operated in Chiba Prefecture Environmental Research Center, with its beam pointed northward at
the elevation angle of 38 deg. The center is located on the east of Tokyo bay, about 10 km south of Chiba University. There isan industrial area and a busy load on the seaside (west of the center).

Auto Alignment
Mechanism

Mirror Controller

Laser Head
Telescope

Laser Driver

38 deg.
PMT

ScalerHorizontal Scanning
Mechanism

Rotation stage
Controller

Window
Rotation Stage

Fig.1 System configuration of portable automated lidar.

Table.1 Specification of PAL

Laser

LD pumped Nd:YAG Laser
Pulse Power 15 J
Wavelength 532 nm

Detector
Telescope

Photo-multiplier(photon counting mode)
Schmidt-Cassegrain
Aperture 20 cm diameter
Field of View 0.2 mrad

ScalerResolution 24 m
Range 24 km
Averaging 10 or 20 s

3. PAL Observation
The main advantage of the continuous and long-term observation is capturing the local weather change that takes place in atime scale of several hours. Especially, the system can monitor the onset and recovery of bad weather conditions and changes
in polluted airs. These features are largely dependent on the site locations and conditions (urban/rural/mountains/waters).
Two examples of characteristic results from the viewpoint of long-term cloud observation are shown in the following.

Figure 2 is the result observed during 0-12 h local time on October 7, 2006. The
weather map of Fig. 2(a) [5] shows that the low pressure has moved northward passing along
the east coat of Japan, involving stationary and cold fronts. The PAL data in Fig. 2(b) also
show that the long-lasting rain from the day before stopped and the cloud gradually gained
altitude. The PAL data shown here are all corrected by the squared distance. Relatively large
echo appeared under the cloud till 7:00 (local time) in the morning. On that day, temperatureand humidity largely changed at 7:00 (local time). Wind direction was northwest, and its
speed was 8 m/s. The 10-h cloud elevation indicates the passage of highly developed low
pressure. Figure 3 is the result observed during 0-12 on September 18, 2006. It was a windyday. Low clouds of less than 1km altitude appeared during 0-5 h. They raised the altitude up
to 2 km during 5-8 h. The lidar echo from these clouds was sparse and largely fluctuating inaltitude. Sharp downturn of the cloud altitudes during 8-10 h was due to the rainfall.
1

3
Altitude [km]Fig.2(a) Weather map over Japanon October 7, 2006.

1.5

0
0

0

2

4

6
Time [hour]

8

10

12

Fig.2(b) 12-hour cloud long-term observation result : October 7, 2006. Temp. 21 deg, Hum. 35%.
1

Altitude [km]

3

1.5

0
0

0

2

4

6
Time [hour]

8

10

12

Fig.3 12-hour cloud long-term observation result: September 18, 2006. Temp. 24 deg, Hum. 82%.

Altitude [km]Examples of long-term temporal motion of the atmosphere are shown in Fig.4. Figure 4(a) is the result observed
during 0-12 h local time on September 21, 2006. The atmospheric boundary layer and cloud were captured at the altitude of 2
km and 4 km, respectively. The structure was stable and showed little change till 9:00, while the relatively large echo
appeared and raised its altitude from the ground during 9-12 h. In accordance with the change, cloud appeared at the altitudeof 1.5-2 km. This condition continued till 16:00. Figure 4(b) is the result observed during 0-12 h on December 23, 2006. The
cloud appeared at the altitude of 6 km and lowered its altitude gradually from 0 to 6 h. Another cloud appeared on theboundary layer at the altitude of 1.5km starting from 6:00. The boundary layer reduced the altitude down to 0.3-0.5 km.
Furthermore during 8-10 h, another cloud appeared on the lowered boundary layer. Obviously those results demonstrate thebenefit of long-term observation. The change in temperature, wind, and the local-climatological influence of the site location
will also be reflected in the observation data.

1

6
3
0
0

0

2

4

6

8

10

12

Time [hour]

Altitude [km](a) September 21, 2006. Temp. 24.7deg, Hum. 59%
1

6
3
0
0

0

2

4

6

8

Time [hour](b) December 23, 2006. Temp. 11.9deg, Hum. 57%
Fig. 4 12-hour atmosphere long-term observation results.

10

12

The result of temporal and horizontal-scanning observations on July 2, 2007 is shown in Fig. 5. On the day, thecloudy weather from the preceding day gradually worsened and is started to rain in the evening. Time-height indication resultof Fig. 5(a) shows that cloud moved slowly in the altitude range of 1-1.5 km during 0-8 h. The cloud altitude lowered in 8-12h, while another thin echo appeared under the cloud layer. It rained in 15-19h (Chiba city).
3
Altitude [km]

1

1.5
0

0
0
3

1

2

5

4

3

6

7

9

8

10

11

12

Altitude [km]

1

1.5
0

0
12

13

14

16

15

17

18

19

Time [hour]

21

20

22

23

24

(a) 24-hour long-term observation result
00：39

01：39

02：39

03：39

04：39

05：39

06：39

07：39

08：39

09：39

10：39

11：39

12：39

13：39

14：39

15：39

16：39

17：39

18：39

19：39

20：39

21：39

22：39

23：39

(b) 24-hour horizontal scanning observation result
100

1020

15
0

5

10

15

20

80

Air Pressure [hPa]

Humidity [%]Temperature [deg]

90
20

1010
4
1000

990
0

2

5

10

15

20

0

Time [hour]Time [hour](c) Temperature & Humidity(d) Pressure and Wind Speed data

Fig. 5 24-hour long-term and horizontal scanning observation result : July 2, 2007.

Wind Speed [m/s]625Temperature/humidity variation shown in Fig. 5(c) and pressure/ wind speed variation in Fig. 5(d) also indicate the samechange in the atmosphere activity, particularly the change of humidity in 0-12 h and 15-19 h, and the change ofpressure/wind speed in 15-19 h. The spatial distributions of the lidar echo obtained by the horizontal scanning are shown in
Fig. 5(b). Although the scanning data are also corrected for the squared distance, it is not corrected for the elevation angle.
Thus, the graphs are plotted in the beam propagation distance. The basic features of cloud echoes agree well with thetemporal variation in Fig.5(a), while the spatial structures of cloud are clearly detected in 9-12 h and 12-15 h by virtue of thehorizontal scanning for the first time. The advantage of the horizontal scanning in understanding the local atmosphere will
be fully exploited by deducing 3-dementional spatial information. In Fig. 5(a), the horizontal scanning time periods of 6 min
are shown with blanks. Alternatively, the scanning data can also be used as part of the temporal data, filling those blanks.
4. Summary
The PAL system has continued the uninterrupted, autonomous observations for nearly 4 years. The additional inclusion of thehorizontal scanning capability enables us to apply the system to new types of targets: spread of industrial smokes and dust
distributions from busy roads are good example of such applications. The system will also be useful to elucidate yellow dust
activity and the pollen density distributions. In the near future, we are planning to install multi-wavelength and multipolarization capabilities to the PAL system.

References
1. Spinhirne, “Micro pulse lidar”, IEEE Transaction Geoscience and Remote Sensing, Vol.31, No.1, pp.48-55, 1993
2. N. Lagrosas et al., “Correlation study between suspended particulate matter and portable automated lidar data”, Aerosol
Science, 36, pp.439-454, 2005
3. G. Bagtasa, N. Takeuchi, S. Fukagawa, H. Kuze, T. Shiina, S. Naito, A. Sone, H. Kan, “Mass extinction efficiency fortropospheric aerosols from potable automated lidar and β-Ray SPM counter”, Proceedings of 23rd International Laser Radar
Conference, 3P-30, pp.499-502, 2006
4. G. Bagtasa, C. Liu, N. Takeuchi, H. Kuze, S. Naito, A. Sone, H. Kan, “Dual-site lidar observations and satellite data
analysis for regional cloud characterization”, Optical Review, 14, pp.39-47, 2007
5. http://www.jma.go.jp/jma/indexe.html

プロジェクト−4RS・GIS による沿岸生態系主要景観の広域分布動態の解析
RS・GISによる沿岸生態系主要景観の広域分布動態の解析
仲岡雅裕（千葉大・理学研究科）・渡辺健太郎（千葉大・自然科学研究科） ・山北剛久（千葉大・理学研究科）・近藤昭彦（千葉大・CEReS）・
小倉久子（千葉県環境研究センター ）・平間幸雄（千葉県環境研究センター ） ・石井光廣（千葉県水産総合研究センター ）結果１：アマモ場とカキ礁の長期変動

はじめにアマモ場の長期変動：アマモ場間で大きく変異
カキ礁の長期変動：1980〜1990年ごろに形成

アマモ場やカキ礁：内湾域の主要景観
高い一次生産量
水産動物の生息場所としての機能
沿岸環境（水質や底質）の改変作用開発や水質汚染による減少が問題

2.0E+06

富津干潟
北条海岸

面積 (m 2)

1.5E+06
1.0E+06
5.0E+05

富津干潟のアマモ場→

0.0E+00

1985

1990

8.0E+04

目的

6.0E+04

2000

2005

2010

1995

2000

2005

2010

2000

2005

2010

2

面積 (m )

東京湾に点在するアマモ場・カキ礁を対象に、
１：RS/GISによりその長期変動を明らかにする。
２：東京湾の水質の変動との関連性を検討する。

1995

波左間
富津港

沖ノ島のアマモ場のGIS解析

4.0E+04

2.0E+04

0.0E+00

1985

三番瀬・猫実川河口のカキ礁→

1990

島戸倉
丑山
沖ノ島

5.0E+03

面積 (m 2)4.0E+03

三番瀬
カキ礁

方法

3.0E+03
2.0E+03
1.0E+03

東京湾千葉県沿岸のアマモ場8ヶ所、カキ礁
3ヶ所について、航空写真より1970年代以降の面積の長期変動を解析
アマモ場、カキ礁周辺における水質（水温，
透明度、栄養塩濃度）の広域かつ長期的な変動の解析←東京湾赤潮・青潮調査および東京湾水質調査から水質データを抽出，データベース化
して解析

盤洲
アマモ場

0.0E+00

1985

1990

1995

アマモ場面積の長期変動
富津干潟
富津港
アマモ場
8000

旧江戸川河口付近
猫実川河口付近
市川航路奥部2

面積(m )

6000

沖ノ島
北条海岸
アマモ場

4000

2000

津浜
丑山
島戸倉
アマモ場

0
1965

1970

猫実川河口部のカキ礁のGIS解析

調査地点→

結果２：水質の広域長期変動

1975

1980

1985

1990

1995

2000

2005

2010

カキ礁面積の長期変動

結果３：水質とアマモ場面積の関係
空間変異：アマモ場面積は水温、透明度と負、窒素濃度と正の関係の傾向（統計的に有意ではない）
1.E+07

1.E+07

1.E+06

1.E+06

1.E+06

2

2

アマモ場面積(m )1.E+07

アマモ場面積(m )アマモ場面積(m 2 )水質の空間変異：湾奥ほど透明度は低下、窒素濃度は増加

1.E+05
1.E+04
1.E+03
1.E+02

1.E+05
1.E+04
1.E+03
1.E+02

1.E+01

1.E+01

1.E+00

1.E+00171819

平均水温 (℃）

1.E+04
1.E+03
1.E+02
1.E+01020

1.E+05

2

4

6

8

10

平均透明度 (m）12

1.E+00
0.00

14

0.20

0.40

0.60

0.80

1.00

平均窒素濃度 (mg/l）

時間変異：多くのアマモ場で無相関：富津干潟、富津港の一部に相関
アマモ場面積と環境要因の相関 (r)

18
17

1.0

全窒素濃度（mｇ/l)

水温（℃）

透明度（m)

19

5
4
3

16
15
1988

199 2

19 96

200 0

20 04

2
19 88

0.9
0.8
0.7

199 6

2000

200 4

1988

19 92

19 96

200 0

500 0

猫実川j河口カキ礁面積(m 2 )カキ礁面積(m 2)

4000

3000

猫実川j河口2.0

2. 5

透明度 (m)

2

16

17

18

平均水温（℃）

19

1.0E+06

5.0E+05

0.0E+00
0.5

0.6

0.7

0.8

窒素濃度 (mg/l)

0.9

1

3.0

3. 5

300 0

0.7

アマモ場の面積変動のパターンは、アマモ場間で大きく変異
アマモ場面積と水質の時間変化の関連性：多くの場合無相関→より局所的な環境条件が各地のアマモ場の動態に影響しているアマモ場面積と水質の空間変異の関連性：内湾浅海域に堆積物底が多い→水質が悪い内湾域ほどアマモの生育適地が広いことによる擬似相関？

400 0

水質とカキ礁面積（猫実川河口）の有意な相関→水質の改善がカキ礁の拡大維持に関連しているかも可能性

200 0

2000
1.5

2.0E+04

まとめと考察カキ礁面積と環境要因の相関 (r)0.01<p<0.05
p>0.05

4.0E+04

1.5E+06200 4

猫実川河口のカキ礁の面積は、透明度と正の相関、窒素濃度と負の相関
江戸川河口のカキ礁は無相関←データ数不足(n=4)の可能性

カキ礁
水温 透明度 全窒素量
猫実川河口 0.540 0.938-0.937
旧江戸川河口 0.325 0.304-0.076富津干潟

富津港

6.0E+04

0.0E+00

結果４：水質とカキ礁の関係

5000

2.0E+06

8.0E+04p<0.001
0.01<p<0.05
p>0.05

0.6
0.5

1992

透明度 全窒素量0.561-0.702-0.5180.809
0.158-0.1350.264
0.284-0.081-0.377-0.060-0.045
0.041
0.215

アマモ場面積(m 2 )

6

20

水温
0.779
-0.179-0.630-0.4590.465
0.612-0.025

アマモ場面積(m )

アマモ場
富津港
富津干潟
丑山
島戸倉
北条海岸
沖の島
波左間

水質の長期変動：水温・透明度は増加、窒素濃度は減少

0.8

0.9

1.0

1.1

窒素濃度 (mg/l)

1.2

1.3

1. 4

アマモ場、カキ礁の周辺景観を含めたより高解像度での解析の必要性

謝辞
本研究を行うにあたり下記の方々にお世話になりました。ここに記してお礼申し上げます。
千葉大学群集生態学研究室のメンバー；千葉県水産総合研究センターのスタッフ；金田漁業協同組合；富津漁業協同組合；天羽漁業協同組合；館山船形漁業協同組合
本研究は文部省科学研究費補助金(No. 18201043)、千葉県、環境省（自然環境保全基礎調査）の支援によって行われています。

連絡先
仲岡 雅裕 （千葉大学大学院理学研究科）
nakaoka@faculty.chiba-u.jp
http://life.s.chiba-u.jp/nakaoka/index_jp.htm

Measurement of NO2 and aerosol in the atmospheric pollution using differential
optical absorption spectroscopy (DOAS) with a PC projector light source
1Ippei Harada, 1Masashi Miyazaki, 1Daisuke Kataoka, 1Hiroaki Kuze, 1Akihiko Kondoh,
2Takashi Hamada, 3Toshiaki Ichinose
1Center for Environmental Remote Sensing, Chiba University, Japan
2Nagano Environmental Conservation Research Institute, Japan

3National Institute for Environmental Studies, Japan

DOAS measurement with a projector flashlight in Seoul city

Background70

Recently the atmospheric pollution has been improved in Japan as a whole, when compared withsituations two decades ago. Some problems, however, remain in places such as urban roadside areas,
where the environmental standards have not been achieved. Therefore, efforts are required for monitoring
anthropogenic air pollution, especially the combustion products such as nitrogen oxides and particulatematters. The main pollution species in urban areas in Japan are nitrogen dioxide (NO2) and suspended
particulate matter (SPM), usually referred to as PM10. The conventional point sampling of thesepollutants at ground stations leads to concentrations for local environments. In order to evaluate average
pollution conditions, it is valuable to obtain additional information of regional concentrations measured
over a certain distance, e.g., several hundred meters to several kilometers.
The purpose of this study is to demonstrate the capability of the PC projector as a DOAS light source.
The long-path, continuous measurements have actually been performed, and the resulting concentrationsof pollutants are compared with the data from nearby ground-based monitoring stations.

DOAS
Ground System

NO2 operation (ppb)

60
50
40
30

0
12:00

13:00

14:00

15:00

16:00

17:00

18:00

19:00

The map of study area

Seoul city (August, 2007)

For the observation in Seoul measurement was set up at about 30 m height. The projector flashlight was used
as a light source and the distance of the light source and observation was 1.4km. Therefore, the value of
DOAS and Ground system data were same.

3.0

The absorption of NO2 gas species gives a structureof the order of 0.2 in terms of the optical thickness.
In order to extract the NO2 concentration,
contributions from molecules and aerosol particles
can be subtracted by simply applying a linear fit tothe back-ground. Another important aspect fromthis figure is that, by measuring the light intensity, it
becomes feasible to evaluate the aerosol
contribution in the DOAS signal itself.

Gas (NO2)

2.5
2.0
1.5
Aerosol
1.0
0.5
Air molecule
400

425

DOAS移動平均DOAS
そらまめ(鍋屋田)
AEROS

いかほ別館

The light source ofPC projector

30

20
10
鍋屋田自排局

AEROS (Ground system)

0-10
2007/11/21 18:00

2007/11/22 18:00

2007/11/23 18:00

2007/11/24 18:00

Observed variation of NO2 concentration
4.5

信州大学

0.04

Observation point
DOAS
DOAS
そらまめ(鍋屋田
AEROS )
3

0.03

1.5

0.02

0
200 7/11/23 0: 00

The map of study area

Nagano city (November, 2007)

3

3.5

DOAS
40

3)
SPM 濃度 (mg/m (mg/m)
SPM concentration

In the conventional long-path DOAS method,a continuously emitting light source is employed,
and the source (or occasionally a retroreflector) is
placed at a certain distance from the observationsite.
A simple detection system is devised by means ofa telescope and a compact CCD spectrometer.

50NO2濃度 (ppb)
NO2 concentration(ppb)

Differential Optical Absorption Spectroscopy (DOAS)

DOAS measurement with a projector flashlight in Nagano city

aerosol
Optical光学的厚さ
thicknessτ(τaerosol )We report the results of three recent campaigns: Seoul city in South Korea (August 2007), Nagano city
(November 2007) and Chiba city (December 2007) in Japan. Seoul city has the area of about 605.33
km2 (2006) (nearly the same as the Tokyo area) with the population of about 10.35 million (2006). The
restoration of the Cheong-gye stream in the midtown has contributed to alleviate the urban pollution
situations. The environmental standards of NO2 and SMP have been mostly attained in Nagano city,while heavy traffic still causes problems in the urban Chiba city area. In these measurements, wemeasured optical thickness due to NO2 absorption and aerosol extinction in the lower troposphere(atmospheric boundary layer) using nearly horizontal optical paths in a height range of 15-100m fromthe ground level.

Optical thickness

Freaya Town
(23F)

The light source ofPC projector

DOAS data measured in every five minutes.
Ground System data measured in every one hour.

The diurnal change of NO2 concentration 30 August,2007

0.0

1400m

Observation
point

20
10

Study Area

Cheonggye
streamShinhan Bank
(9F)

The comparison was difficult in Nagano city between theNO2 concentration from DOAS and the ground (AEROS)
data, mainly due to the small concentration of pollutants,
and in addition, occasional snow fall in the mountainregion where the PC projector light source was placed.

0.01
2007/11 /24 0:00

2007/11/23 12:00

Comparison between aerosol optical thickness and SPM concentration

450

W avelength (nm)DOAS measurement with a projector flashlight in Chiba city

Optical thickness of each component detected by the DOAS spectrometer

80

DOAS measurement with an obstruction flashlight in the central Tokyo

70
60

10 Aug (DOAS)
10日
11日
11 Aug (DOAS)
12日
12 Aug (DOAS)
13日
13 Aug (DOAS)
そらまめ君
Ground System

NO2NO2濃度(ppm)Concentration (ppb)

50

40

The light source of an
obstruction flashlight30

NO2 concentration(ppb)
NO2濃度 (ppb)

60

千葉大学
千葉大学

The light source ofPC projector

50
40
30
20

DOAS

10

DOAS移動平均DOAS

02007/12/27 6:00

20

AEROS (Ground system)
千葉市役所自排局
千葉市役所自排局

そらまめ(千葉市役所)
AEROS
2007/12/27 18:00

2007/12/28 6:00

2007/12/28 18:00

千葉ポートタワー
千葉ポートタワーObservation point

Observed variation of NO2 concentration

11日
11 Aug

12日Aug
1213日
13 Aug

Diurnal change of NO2 concentration 10-13 August,2007
Since the distance between the ground observation point and light path was 5 km ,it is understandable that the difference between the DOAS and ground data wassignificant. During the night to the morning on August 10 and the night on August
11, high NO2 concentration was observed because of the rush hour on the regional
highways due to the family reunion in this particular season. Since the light froman obstruction flashlight was not available in the nighttime, the DOAS method
could not capture this high concentration.

Observation point

The map of study area

Tokyo city (August, 2007)

DOAS measurement with a projector flashlight

Power consumption of the projector is 220 W, andthe emission covers the visible wavelength.

The drawback of the use of aviation obstruction light asa DOAS source is that in accordance with the regulation,
the light intensity is diminished at dusk and dawn, andduring the nighttime, blinking red lights replace theflashlights. Thus, the DOAS measurement is limited tothe daytime. Moreover, the measurement cannot be
carried out where no obstruction flashlight is situated.
Alternatively, here we propose the use of acommercially available PC projector as a white-lightsource. This is relatively inexpensive, yet the possibilityof unattended, continuous operation is quite suitable forthe DOAS measurement. This light source is portable,
and it can be used during both daytime and night time.

0.1
DOAS DOAS

4

そらまめ(千葉市役所)
AEROS

0.08

3

0
10日
10 Augτ aerosol (τ aerosol)
Optical
光学的厚さ
thickness

5(mg/m )
SPM濃度
SPM
concentration(mg/m3)

10

3

0.06

2

0.04

1

0.02

0
2007/12/27 9:00

2007/12/27 15:00

2007/12/27 21:00

0
2007/12/28 3:00

The map of study area

Chiba city (December, 2007)
A reasonable temporal correlation was found between
the result of the long-path measurement using DOAS
method and the data of NO2 concentration from aground station near each optical path in Chiba city.
Moreover, a reasonable correlation was foundbetween the aerosol optical thickness using DOAS
method and the data of SPM concentration from aground station in these cities.

Comparison between aerosol optical thickness and SPM concentration

Conclusions
A reasonable temporal correlation was found between the result of the long-path measurement
using DOAS method and the data of NO2 concentration from a ground station near each
optical path in Seoul city and in Chiba city. Moreover, a reasonable correlation was foundbetween the aerosol optical thickness using DOAS method and the data of SPM concentrationfrom a ground station in these cities. However, the comparison was difficult in Nagano city
between the NO2 concentration from DOAS and the ground (AEROS) data, mainly due to thesmall concentration of pollutants, and in addition, occasional snow fall in the mountain region
where the PC projector light source was placed.

プロジェクト−5

衛星を用いた地殻変動監視に関する研究
火山活動のモニタリング インドネシア･メラピ火山
服部克巳、野口恭史、大野望(千葉大学)、前田崇、高野忠(JAXA)

1 研究背景と目的Rp)の値との差分値である評価関数 S を定義する。

視等の物理観測や火山ごとに過去の噴火様式や噴火(ⅳ)各位置の評価関数 S をシーンごとに並べ、S の最
大値 S max 、平均値 S を算出した。それらの値から S

間隔等の解析が必要であり、共通の監視手法はない。

の特異性を表す乖離率δを次のように定義した。[1]

従来の火山噴火の予知は火山性微動や傾斜計の監

しかし、火山の噴火に伴い地表面の温度が上昇するこSSとはよく知られている。地表面の温度異常ならば、衛星センサを用いて観測することが可能である。衛星を用いれば高頻度、広域に監視でき、インドネシアやフィリピン等の島国という環境に位置する火山を連続監視することも可能である。そこで、本研究は人工衛

また S max に対応するδを最大乖離率

max と定義

し各シーンにおける S max 、 max の分布図を作成する。

星 Aqua に 搭 載 さ れ て い る 赤 外 線 セ ン サ

4

MODIS(Moderate
Resolution
Imaging
Spectroradiometer)データを用いて火山噴火に関連

前)の二次元強度マップを図 1 に表す。カラースケー

する地表面の温度異常を検出することを目的とした。

ルは放射輝度を整数化したものであり、図の中心はメ

解析結果
前節の(ⅰ)で作成した 2006 年 5 月 10 日(噴火の 4 日

ラピ山頂である。図１の山頂付近の A の領域で高い

2 MODIS データ

値を確認できる。

人工衛星 Aqua は太陽同期準回帰軌道を使用し、同次に上記の領域内で最も大きな変化が見られた山

一地点を同一条件で観測することが可能であり、観測頂の 1km 東の輝度値の時系列変化を図２に示す。こ領域内の赤道を通過する昇交軌道時(ascend)は常時

こで、縦軸は図 1 のカラースケールと等しい値であり、昼間であり、降交軌道時(descend)は常時夜間となる。

横軸は時間である。上の矢印の期間が噴火期間を表す。

ほぼ毎日インドネシア上空を通過し夜間データを得

噴火期間の少し前から輝度値の上昇が見られ、噴火期ることができる。また、夜間データには太陽の影響が間で大きく上昇していることがわかる。噴火後に高い含まれておらず、地球内部からの熱情報のみを取得す

値を保持するのは、熱源がゆっくり冷えていったためることができる。

である。

MODIS は観測対象の異なる 36 band を持つ。本研

この値の上昇が火山活動に関連するものであるか

究では band 20 を使用する。この band は中心波長確認するために山頂付近以外の場所との輝度変化と3.750μm、空間分解能 1000 m で、比較的大気の影の比較を行った。
例として図 1 の B の位置の時系列変響を受けにくく地上温度の観測が可能である。

化を図 3 に示す。図 3 では噴火活動期間中に値の上昇は見られない。したがって、山頂付近の値の上昇は火3 データ解析手順
インドネシア･メラピ火山はジャワ島中部に位置し、山噴火活動に起因することが強く示唆される。
また、領域内のほとんど全ての時系列変化から季節

過去に 1930 年、1994 年、2006 年に大きな被害噴火

変動と思われるトレンドが見られた。季節変動や空間を発生させている。本研究では人工衛星 Aqua が飛行的な変動要素(気象条件や土壌水分)の影響を取り除き、中である 2006 年のメラピ火山の噴火を解析対象に取図 1 中の C のような普段から高い値をとるような場り上げた。

所と区別するために手順(ⅲ)(ⅳ)に従い解析を行った。

次に解析手順について簡潔に記す。

山頂位置で S max 、 max を記録したシーンの S 値と乖(ⅰ) MODIS データから山頂を中心に±0.5°の領域

離率δの２次元分布図をそれぞれ図 4、図 5 に示す。

の熱源の強さを表す二次元強度マップを作成する。

両図は 2006 年 6 月 3 日の分布図であり、溶岩ドーム(ⅱ)作成した領域内の任意の位置の時系列図を作成が崩壊する前日である。図 4 より山頂位置 D で高いし、解析期間は 2005 年 4 月~2007 年 10 月迄とする。

値が確認できる。また、山頂から離れた位置 E でも(ⅲ)季節変動や空間的な変動要素を取り除くため、各同様に高い値がある。しかし、図５の最大乖離率を見シーンの領域内の任意の位置(Focal point : Fp)の値

る E’の地点の値は D’よりも小さく、E’よりも D’の特から東西南北に 5 km 離れた位置(Reference point :異性が高いことがわかった。

図 1 : 2006 年 5 月 10 日の二次元強度マップ

図5:5max の分布図(2006 年 6 月 3 日)

まとめ
MODIS band 20 の夜間データを用いてメラピ火山の噴火活動に関連する温度異常を検知することができた。また、評価関数 S 、 S max や乖離率

、

max を

導入することで、山頂付近の噴火活動に関連する温度異常のみ検出することができた。したがって MODIS
図 2 : 山頂の 1km 東の位置の時系列図

データを用いた火山噴火活動の監視は有効である。
今後は解析期間をのばしてメラピ山に関する輝度変動の特異性についてさらに調査する必要がある。また、噴煙との違いを明確にするために、可視の画像とも比較検討する必要がある。さらに本研究にて開発した手法を他の火山に適用し、有効性を確かめる必要が
ある。

参考文献
図 3 : 山頂から離れた位置の時系列図[1] Maeda Takashi, Data Processing for Extracting
Local and Faint Changes from a Satellite-borne
Microwave Radiometer Data, IEEE Transactions
on Geoscience and Remote Sensing, 2007

図 4 : S max の分布図(2006 年 6 月 3 日)

一般研究

合成開口レーダーによる洋上風力発電ポテンシャルの評価
合成開口レーダーによる洋上風力発電
ポテンシャルの評価

陸上

沿岸域(SAR, MM5)外洋(GCM, 散乱計）

○香西克俊, 大澤輝夫
(神戸大学海事科学研究科)
平成19年度CEReS共同利用研究発表会 2008年2月28日 千葉大学 けやき会館

内 容

浮体式風力発電（日本など（構想））

１．はじめに２．データ及び方法
３．結果と考察４．まとめ

着底式風力発電（デンマークなど）
図 海上風況調査手法と海域による分類

目 的
Cバンド及びLバンド合成開口レーダーを利用した海上
風速推定の精度比較を行うことにより、日本沿岸域に適した洋上風力資源評価の新しい方法を開発すること
Fig.2. First offshore? wind farmin Japan at Setana.(600KW x 2)
（http://www.mlit.go.jp/kowan/kaihatuka/
wind_hp/jirei-world）

ASAR geometrically-corrected
NRCS image（Sep.24, 2003, 01h
09m（UT), 500m resolution)

Fig.3. Horns Rev （Denmark）
２ＭＷ × ８０
（http://www.hornsrev.dk）

Specifications of ENVISAT/ASAR

Specifications of ALOS/PALSAR

Mode
Image mode (IM)
Product
Precision product
Beam
IS2
Incidence angle 18.7~26.2deg
Swath
107.7kmpolarization VV
Pixel spacing12.5m
Orbit
Descending
PeriodFeb.,03-Mar.,06
No. of scene 35

Mode
Fine beam
Incidence angle 23.8~36.0deg
Swath40~70kmpolarization HH
Pixel spacing 6.25m/12.5m
Orbit
A/D
PeriodOct.,10-Dec.,13
No. of scene 12

図 幾何補正済みPALSAR画像（Fine beammode, 2007年12月13日13h27m（UT)、円内
は白浜海洋観測塔の位置を示す。白浜風向
333度、風速11.8m/s、推定風速5.0m/s)

Table 4. Available observational equipments
Observation
instrument
ADCP ×2

Profiling
range

Manufac
ture

0~20
m

RD
Instrume
nts
ALEC
ElectronicsElectromagn
etic current
meter×2

2005 May

モニン・オブコフ相似則から導かれる風速鉛直プロファイルは次式で表される．
ここでΨ は次式で表される無次元関数である．

Notes
ADCP frequency 1200Hz

Current speed/Wave height
Temperature/Salinity/Depth/Conductivity/Dissolved Oxygen/Turbidity/Chlorophyll-a

0~200
m

ALEC
ElectronicsTemperature/Salinity/Depth

Ultrasonicwave gauge

FS
15m

KAIJO
SONIC

Wave speed/Wave height

Ultrasonicanemometer

0~30
m/s

KAIJO
SONIC

Wind speed/wind direction

YAMAH
A

15m long from stem to stern

Observation
ship

ln

z
z0

(1)
2

1

2

1

2

2 tan 1

2

ただし，

2

(1 19.3 )1 / 4

0

(2)

0

Sea-Bird
Electroni
cs

STD

u*

ln

Current
speed/Temperature/Salinity

0~200
m

CTD

u ( z)ADCP frequency 1200Hz
Wave ADCP

IOTechnic23m

高度補正の方法 (Monin-Obkov similarity)

無次数ζ(=z/L)は，バルクリチャードソン数Ribの関数として次式で求める．

Rib
4.5

10 Rib 1
11.547 Rib
1 5 Rib

1

( Rib

0)

( Rib

0)

(3)

粗度長を2.0×10 -4mと仮定すれば，次式より高度10mでの風速値が得られる．
ln

10
z0

ln

23
z0

U10

Fig.8. Marine tower at Shirahama.

(4)

U 23

Wind speed estimation using CMOD4 and LMOD
ASAR/PALSAR Raw image
Amplitude to PowerPower image
Speckle noise removal

方 法

BEST (Basic Envisat SAR
Toolbox, V4.0.5,2007, ESA)ERDAS

Noise-removed image
Backscattering
image generation

NRCS (sigma nought) image
CMOD4 (Stoffelen and Anderson, 1997)
LMOD (Isoguchi and Shimada, 2007)+In situ winddirection

Wind speed estimation at 10m height

CMOD4 σ0=f (wind speed、 relative wind direction、incidence angle)
σ0lin=b0(1+b1cosχ+b3tanhb2cos2χ)1.6(Stoffelen and Anderson, 1997)

0

LMOD σ0=f (風速W、 相対風向φ、入射角θ)σ0=A0(1+A1cosφ+A2cos2φ)（A0,A1,A2は風速と入射角の関数）-2(Isoguchi and Shimada、2007)

0

-6

-10

-8

Sigma nought (dB)

NRCS (dB)-4-9.36dB

-10
-12
-14

-20

-30

-16
0

4.3m/s

5

10

15

20

25

30

W in d s p e ed (m / s )Fig.10. Relationship between Normalized Radar Cross Section(σ0) and wind speed（May,7,2004、Incidence angle=25.8deg、Relative wind direction=344.3deg）

-40
0

5

10

15

20

25

W in d sp e e d (m / s)

図 後方散乱係数（NRCS）と風速の関係（入射角36度、相対風
向253度、2007年12月13日の例）

Azimuth direction
Radar antenna

結果と考察
Crosswind direction （χ=270°)
Downwind direction（χ=180°)Ran
ge direc tion

Upwind direction（χ=0°)
Crosswind direction （χ=９0°)

Fig.5. Geometry of SAR viewing and relative wind directionχ.

18

RMSE=1.93m/s
BIAS=0.13m/s

8
Estimated wind speed (m/s)

Estimated wind speed (m/s)

16
14
12
10
8
6
4

RMSE=2.61m/s
BIAS=-1.47m/s

6
4
2
0

2

0

0
0

2

4

6

8

10

12

14

2

16

4

6

8

10

12

Insituwind speed (m/s)

In situ wind speed (m/s)

Fig.13. Results of comparison between in situ and the estimatedbased on CMOD4 algorithm

図 LMODに基づく推定風速の白浜現場風速に対する検証結果（円内は2007年12月13日）

Image date Sep.24, 2003
CMOD-4+in situ
(Sep.24, 2003)

LMOD+in situ
(Dec.13, 2007)

Fig.12. Examples of estimated offshore wind speed distribution based on CMOD4and LMOD with in situ wind directions.(500m grid)
Shirahama observed wind speed & direction=11.0m/s, 59.2deg
Shirahama estimated wind speed=5.4m/s (CMOD4+in situ)

NRCS (dB)

CMOD4-upwind

CMOD4-crosswind

まとめ

4
2
0
-2
-4
-6
-8
-10
-12

(1) Cバンド及びLバンド合成開口レーダーを利用した海上風速
推定の精度検証を行った結果、RMSEはそれぞれ1.93m/s,
2.61m/sとなり、両者とも推定風速は実測風速に対し過少評価であった。

0

5

10

15

20

25

Wind speed (m/s)
Fig.15. Normalized radar cross section (NRCS) as a function of windspeed for the fixed radar incidence angle 21 degrees.

謝 辞
本研究は欧州宇宙機関（ESA）との共同研究プロジェクト（Offshore wind resource assessments using SAR and MM5 overJapanese coastal waters, C1P4068）に基づき、ENVISAT/ ASAR
データを取得した。また、本研究は宇宙航空研究開発機構（JAXA)との共同研究に基づき、ALOS/PALSARデータを取得した。本研
究は京都大学防災研究所との共同研究による成果である。本研
究は平成19年度文部科学省科学研究費補助金（基盤研究（B)（一般）「合成開口レーダー及びメソ気象モデルを用いた洋上風力資源調査手法の開発」課題番号19360406） 及び若手研究（A） 「洋上風力発電に関わる気象研究拠点としての風況モニタリングステー
ションの開設」課題番号19686052の支援を受けた。

(2) 実測風向が北西から東の場合のCMOD4及びLMODによる海上風速推定は大きな誤差を示した。これは北西と東方向から吹く風の吹送距離が短く、さらに海岸地形により変化を受けているためと考えられる。

SAR offset displacements from Sumatra-Andaman Earthquake
Mw 9.2 on December 26, 2004
Ashar Muda Lubis and Nobuhiro Isezaki
Geophysics Laboratory Faculty of Science Chiba University, Chiba, Japan
Tel: +81-43-290-2854, Fax: +81-43-290-2859
E-mail: asharml@graduate.chiba-u.jp

Abstract— The Sumatra-Andaman eathquake, the epiceter about 163 km western part of
Sumatra island, occured on Decemeber 26, 2004 at 7:58 a.m. local time. The magnitude of theearthquake was 9.1 on the moment magnitude scale and became one of the largest magnitudeof the earthquake in world. We use the benefit tools of Amplitude Maching Pixeal of
Synthetic Aperture Radar (SAR) data to make a rough estimation of crustal displacement
associated with this earthquake. Fours of C- band synthetic aperture radar (SAR) images
acquired by ERS-2 synthetic aperture radar instrument are analysed. In order to produce SLC(Single Look Complex) and to obtain offset displacement, we use a comertial GAMMA
software depeloved by GAMMA Remote Sensing Research and Consulting AG, Switzerland.
The results of offset SAR displacement show that maximum displacement in northren part
Sumatra Island is 4-6 meters. Generally our result agrees well with long span GPS
observation.
Key words: Sumatra, offset, SAR, displacement

1. Introduction
The Sumatra-Andaman earthquake, a devastating earthquake, was occurred in Indian
ocean at western of Sumatra Island on 07:58:53 LT December 26, 2004 which has magnitude
originally recorded as Mw 9.2, lately upgraded to 9.3 (Stein and Okal, 2005). The
hypocentre of the main earthquake was at 3°19′N 95°51.24′E, approximately 160 km west of
Sumatra, at a depth of 30 km below mean sea level (see Figure 1). Unfortunattely theearthquake was followed by tsunami. It was repoted that the tsunami made landfall within
10-15 minutes, run-up heights1 of over 30 m in Banda Aceh (Guy, 2005). Run-up heights
of more than 10 m were also documented in Thailand (12.14 m) and Sri Lanka (11 m).
Documented fatalities eventually exceeded 283,000 with the heaviest losses concentrated
along the west coast of Sumatra, but more than 40,000 were also accounted for on distantshores around the Indian Ocean. More than 1 million people were displaced (USGS, 2007).

Figure 1. History earthquake in Sumatra-Andaman region, Green star is the epicenterof 2004 and 2005 event
Many Scientifics have been analyzing the size of the earthquake, slip, duration, speed,
and spatial rupture extent and tsunami generation. Such as Tobita et al., (2005) measured
displacement generated by the earthquake using Synthetic Aperture Radar (SAR) images, but
this method measures just vertical displacement and can not yield quantitative value ofdisplacement. In this paper we try to calculate horizontal displacement associated with theerathquake using amplitude offset SAR Analyzes (SAR matching image) (Tobita et al,. 2001;
Pathier et al., 2006), and map vector dispalacement around the epicenter of the earthquake.

2. Data dan method
In order to estimate crustal displacement related to the Sumatra-Andaman earthquake,
we use radar imagery data (see Table 1) collected by the C-band (5.2 GHz) SAR instruments
aboard the ERS-2 satellite. The SAR data which are acquired on 2002 and 2005 covers areaaround north of Sumatra Island (Figure 2). The raw signal data are processed using acommercial software GAMMA SAR processor whose output consists of a signal that is ameasure of the complex backscatter of a patch on the ground delayed by the travel time of thesignal from sensor to target and back.
Table 1. SAR data used for monitoring crustal displacement in in the north of Sumatra
Islands.
Location

MasterSlave

Perpendicular Baseline (m)Period (days)

Loaction 1

2002/11/12

2005/01/25

57

1109

Location 2

2002/11/28

2005/01/06

310

1074

Figure 2. SAR images location in the north of Sumatra islands

The displacement field can be measured via incoherent or coherent cross-correlationof small image chips. This method has the advantages that it does not require phaseunwrapping. In cases where there is some degree of interferometric coherence between dataacquisitions, the single-look complex (SLC) image speckles themselves give features that canbe accurately tracked. This method is complimentary to the phase-based approach since itworks well with the large displacements. Furthermore, image cross-correlation measurements
yield the two-dimensional displacement field while measurements of the phase yield
deformation only along the line-of-sight. Accuracy of the cross-correlation method isdependent on the scene content, and image chip size.This method has been successfully
applied by Michel et al,. (1999a and 199b); Tobita et al,. (2001) and Pathier et al., (2006).
Basic Idea of measuring pixel offset from two SAR images is illustrated in Figure 3. In thismethod, using two SAR images, only use amplitude SAR image from SLC image (not phase,
not interferometry), we try to find better window size of cross-correlation between master andslave image to improve SNR.

Figure 3. Schematic of basic idea measuring offset displacement from SAR image. The
residual offset in K letter means ground displacement related to the earthquake.
In general this method is very easy, but distortion due to topography on the radar
images must be intensive in research area, therefore slave SLC must be precisely transformed
to a master SLC geometry. Otherwise, the influence due to the difference of the observationgeometry including terrain height difference must be quite big, and the "offset" that I want toknow will be very small, which may be less than 1/10 of the offset due to observationgeometry difference between a slave and a master, so the new slave image which is registeredto master image should be generated then the geometry both master and slave registeredalmost same. The procedures to product the new slave image can be seen in Figure 4. In thisprocedure the Digital Elevation Model (DEM) data are needed. I mosaic the SRTM withresolution 90 m to produce terrain height data. The DEM data which used in this research are
shown in appendix 2.

Figure 4. The procedure for generating new slave image registered to master image.
Since the new slave image is obtained the calculation for measuring offset
displacement can be done as follow:
The first method to estimate the range and azimuth registration offsets filed of two
SAR images is the intensity tracking, known as cross correlation optimization procedure(Helmut et al., 1998; Gary and Mattar, 2000). The offsets are generated with a normalized
cross correlation of image patches of detected real-value SAR intensity image. The success of
this method for the local image offset depends on the presence of nearly identical feature in
two SAR image at the scale of the employed patches. When coherence is retained, the speckle
pattern of the two images is correlated, and intensity tracking with small image patches can be
performed to remarkable accuracy. Incoherent intensity tracking is also feasible but requires
large images patches. In order to increase the estimation accuracy, oversampling rates areapplied to image patches, and a two-dimensional (2-D) regression fit to model the correlation
function around the peak is deteminated with a four-point interpolation. The location of thepeak of the 2-D cross-correlation function yields the image offset. Confidence in the offset
estimate is measured by comparing the height of the correlation peak relative to the averagelevel of the correlation function to supply a correlation signal–to-noise ratio (SNR). Coarse
information on slant-range and azimuth offsets is used to guide the search of thecross-correlation maximum. For the present study with ERS-2 data, I systematically used

patch size 128 x 128 single look pixels.
The second, coherence tracking, a method of measuring the offset between SAR
images complementary to the intensity tracking is the coherence tracking, also known as thefringe visibility algorithm or coherence optimization procedure Derauw, (1999). Throughout
the SLC SAR images, small data patches are selected; a series of small interferogram withchanging offset is constructed; and coherence is estimated. The location of the coherence
maximum is determinated by subpixel accuracy by oversampling the single-look complex
SAR image patches with zero padding and by using 2-D regression function to model thecoherence peak with four-point interpolator. The magnitude of the coherence maximum
relative to the average level is used as a quality factor to help reject unsuitable patches. Coarse
information on the slant-range and azimuth offsets is used to guide the search for coherence
maximum.
3. Discussion and conlusion

Figure 5. Horizontal displacement in the north of Sumatra islands (SAR image location 1)
associated with megahurst Sumatra-Andaman earthquake December 26, 2004

The result of offset horizontal displacements from in location 1 is shown in Figure 5.
The maximum displacement is about 3-4 meters and the amplitude of offset displacements
decreases gradually with increasing distance from epicenter. In northern part of SAR image
on location 1 shows us that displacement less than error of offset measurement. For validation,
We compare this result with displacement obtained from long span GPS observation(Subarya et al., 2006) GPS displacements (see Figure 6) were obtained from measurements
collected one month or more after the earthquake, compared to measurements collectedduring various surveys between 1991 and 2001. They obtained that the maximum horizontal
GPS displacement in location SAR image 1 is 3.63 - 3.71 meters.

Figure 6. Crustal displacement related to Sumatra-Andaman observed by long term GPS
observation, GPS data from Subarya, et al., (2006)
The maximum offset displacement from pair SAR image location 2 is 5-6 meters(Figure 7). The average of offsets horizontal displacement is about 3 meter. Unfortunately
there is no GPS data available in this research area. The closest GPS observation to SAR
image in location 2 shows that horizontal displacement is about 3.20 meter (Subarya et al.,2006). In the southern part of SAR image location (Simalue Island) which is very close toepicenter of the earthquake, the displacement from GPS observation reaches to 5.8 meter. In
general the horizontal displacements derived from SAR Amplitude offset both in SAR image
location 1 and location 2 are good agreement with GPS observed by Subarya et al., (2006)

Figure 7. Horizontal displacement in the north of Sumatra islands (SAR image location 2)
associated with megahurst Sumatra-Andaman earthquake December 26, 2004
We conclude that the Amplitude Offset Analyzes (SAR matching image) can asses
crustal displacement assciated with Sumatra-Andaman earthquake. This method two offsets
amplitude image of master and slave image were succefully correlated with signal noise ratio(SNR) more than 15. Generally we obtained that the maximum displacement in horizontal
direction derived from this method is 4-6 m and this results are good agrement withdisplacement which obtained from long span GPS measurement (Subarya et al., 2006).
4. Acknowlegment
We are thankful to Mr. Takumi Onuma at JGI Japan for introducing and training
GAMMA software. Some of figure was made by using GMT software (Wessel and Smith,1998).

5. References
Derauw.D (1999), DInSAR and coherence tracking applied to glaciology: The example of
Shirase Glacier, in Proc. FRINGE’99, Liège, Belgium.
Gray, L. K. Mattar, and G. Sofko (2000), Influence of ionospheric electron density
fluctuations on satellite radar interferometry, Geophysical Research Letters,
27(10):1451–1454.
Helmut Rott, M. Stuefer, A. Siegel, P. Skvarca, and A. Eckstaller (1998), Mass fluxes anddynamics of Moreno Glacier, Southern Patagonia Icefield, Geophysical Research
Letters, 25 (9): 1407–1410.
Guy Carpenter (2005), Tsunami:Indian Ocean event and investigation into potential global
risks, http://www.guycarp.com/portal/extranet/pdf/Tsunami_Report_Final.pdf.
Michel, R., J.P. Avouac, and J. Taboury (1999a), Measuring ground displacements from SAR
amplitude images: Application to the Landers earthquake, Geophysical Research
Letters, 26: 875-878.
Michel, R., J.P. Avouac, and J. Taboury (1999b), Measuring near field coseismic
displacements from SAR images: Application to the Landers earthquake, GeophysicalResearch Letters., 26: 3017-3020.
Pathier, E.; Fielding, E. J.; Wright, T. J.; Walker, R.; Parsons, B. E.; Hensley, S (2006),
Displacement field and slip distribution of the 2005 Kashmir earthquake from SAR
imagery, Geophysical Research Letters, 33(20), L20310, 10.1029/2006GL027193.
Tobita M, Hisashi Suito, Tetsuro Imakiire, Min Kato, Satoshi Fujiwara, and Makoto
Murakami (2006), Outline of vertical displacement of the 2004 and 2005 Sumatra
earthquakes revealed by satellite radar imagery, Earth Planets Space, 58, e1–e4.
Stein. S and Okal E.M (2005), Speed and size of the Sumatra earthquake, Nature, 434:
581-582, doi:10.1038/434581a.
USGS (2007), Magnitude 9.1 off The West Coast of Northern Sumatra
http://earthquake.usgs.gov/eqcenter/eqinthenews/2004/usslav/#summary
Tobita, M., M. Murakami, H. Nakagawa, H. Yarai, S. Fujiwara, P. A. Rosen (2001), 3-D
surface deformation of the 2000 Usu eruption measured by matching of SAR images,
Geophysical Research Letters, 28(22): 4291-4294.
Wessel, P., and W. H. F. Smith, 1998, New, improved version of the Generic Mapping Tools
Released, EOS Trans. AGU, 79, 579.

Trial of Digital Filter Photography for Alteration Mineral Detectionin the Hachimantai Area, NE JAPANIntroduction
Trial of Digital Filter Photography for Alteration Mineral
Detection in the Hachimantai Area, NE JAPAN• This study focuses on hydrothermally altered materials using satellite
image (ASTER data) and trial new digital filter photography remotesensing method in the Hachimantai area, Northern JAPAN.
• Most satellite images are good quality and georeferenced so they can beloaded directly into GIS software.
• Unfortunately, most satellite systems have limited resolution, limited
orbital periods. Cloud cover adversely affects them at the time of image
acquisition.
• Alteration zones can guide exploration geologists to hidden systems or toancient spring activity and important in geothermal resource exploration
over the Hachimantai area.

Myint Soea, Tateishi Ryutarob , Ishiyama Daizoc, Krit Won-Inc,
Isao Takashimac, Punya Charusirid

a
b
c
d

Graduate School of Engineering and Resource Science, Akita University, Akita 010-8502
Japan.
Center for Environment Remote Sensing, Chiba University, 1-33, Yayoi, Inage, Chiba 2638522, Japan.
Center for Geo-Environmental Science, Akita University, 1-1 Tegatagakuen-cho, Akita 0108502 Japan
Department of Geology, Faculty of Science, Chulalongkorn University, Bangkok 10330,
Thailand28 February 2008

Objective

The Study Area
141º

• To develop a flexible, low cost remote sensing system thatcan be applied in the detection of alteration minerals.
• The aim will be met through the following specific objectives:(1) Develop a lightweight digital imaging system capableobtain high-resolution images.
(2) Demonstrate the usefulness of the filter camera system
for alteration detection.
(3) Demonstrate the utility of the filter camera system forpre-scouting fields.

•

The acidic hydrothermal alteration zones elongated along
ENE striking fractures (Sumi, 1968).
The acidic stage has been divided into three alteration
subzones based on the distribution of kaolinite, alunite andpyrophyllite (Nakamura and Sumi 1981)(1) Silicic subzone – siliceous rocks, alunite and sulphur.
(2) Silicification subzones – silicified rock, clays, sericite,alunite, gypsum, calcite, rutile, diaspore and andalusite.
(3) Argillization subzone – clay, montmorillonite, kaolin
and alunite.
Red = suitable spectral signature in Remote Sensing

14 0° 50 '0 "E

14 0° 5 5'0 "E

1 41 °0 '0 "E

4 0°0 '0" N

Hachimantai

Hachimantai

Study Area

3 9° 55 '0 "N

3 9°5 5'0 " N

40º
Mt. Iwate
Morioka

Mt. Iwate

3 9° 50 '0 "N

3 9°5 0'0 " N
14 0° 40 '0 "E

1 40 °4 5' 0"E

1 40 ° 50 '0" E

14 0° 55 '0 "E
5

•

•

•

140 °45 '0 "E

41º

•
•

Hydrothermal alteration

140 °40 '0"E

4 0° 0' 0"N

2.5

0

1 41 °0 '0 "E
5 Kil om ete rs

Hachimantai volcanic region is one of the largest geothermal provinces inJAPAN, is located 50 km northwest of Morioka city Prefecture.
Hydrothermally altered rocks are exposed by landslide on the hill.
The Geological Survey of Japan considers mapping hydrothermal alteration
zones as an extremely important element in geothermal exploration.
The white polygons are alteration zones mapped by geological survey ofJAPAN.

Methodology(1) Processing of ASTER data Image analysis- DN image to Radiance Image calibration- Radiance Image to Reflectance Image calibration- Reflectance band ratio method- Principal Component Analysis (PCA) method- Band Math Method(2) New digital filter photography method.
(3) X-ray diffraction (XRD) analysis

•This study used the Advanced Spaceborne
ThermalEmissionandReflection
Radiometer (ASTER) data level AST3A01.
It has been acquired on September 16, 2004
with a soloar incidence angle of 50.0º andazimuth angle of 156.10º .
Another one is digital camera filter
photography image. filtered through;- visible wavelength filters - 400 nm, 450
nm, 500 nm, 550 nm and 600nm.
- infrared wavelength filters – 750 nm and
950 nm.

IMAGE DATA
•
1 2

ASTER

3

4

56 7 8 9•0
0.4

0.6 0.8

1.0

1.2 1.4 1.6 1.8
Wavelength ( um )2.0 2.2

2.4

Sub-System Number of bands
Filter band
Sensorwavelength
Product
400 nmAster L1B

Spectral range(μm)

Spatialresolution (m)3 (+ 1 backward)

0.52 up to 0.86

15x15

SWIR

6

1.60 up to 2.43

30x30

TIR

5

8.125 up to 11.65

90x90Band Ratio Method• Band ratio method used this reflectance andabsorption bands characteristic.
ASTER
1 2

3

4

56 7 8 9

80
R (%)
60

Calcite40

Correct slope andorientation of the surface• horizontal radiance (LH)
• sloped terrain radiance (LT)• the sun incidence angle (i)

0.8

1.0

1.2

1.4

1.6

1.8

Radiance Image to Reflectance Image• the relative Sun-Earth distance (d)• the exoatmospheric solar irradiance (Exo λ)• the solar zenith angle (θzenith)• The most appropriate index used to extract thelaterite area is the ratio of band 2 to band 1 in
ASTER data.
Band 2 / Band 1 = New band (iron oxide image)Principal Component Analysis Method (PCA)• The principal component transformation is a multivariate statistical technique.
• This technique indicates whether the materials are represented bright or dark
pixels in the principal components according with the magnitude and sign of theeigenvectors loading.
• As we know that the iron oxide give high reflectance values in ASTER band 2
and low in band1, we look for the principal component in which the difference ofreflectance is large at table.
PC Band 4
Band 5
Band 6
Band 7 Band 8 Band 9
1 0.774759 0.301401 0.372717 0.315822 0.229668 0.132276
2 0.610152 -0.363007 -0.299115 -0.361144 -0.412267 -0.325701
3 0.053647 -0.650366 -0.311138 0.638889 0.247169 0.089823
4 -0.077688 -0.573204 0.792854 -0.188634 0.006616 -0.034026
5 -0.133375 0.151050 0.213947 0.565574 -0.583161 -0.503657

Kaolinite

0.6(ERSDAC, 2001)(Teillet et al., 1982)

Goethite

20
0
0.4

DN image to Radiance Image• Sensor gain/offset• Calibration coefficient (δ)(Mather, 1999)

VNIR• The mineral or rock unit may have high reflectancein some spectral portion, however, it may absorb in
another spectral region (Bannari et.al., 1995).

100

Pre-Processing of ASTER data

2.0

2.2

Output PCA Images

2.4

6 -0.027586 0.055036 -0.020372 -0.078707 0.613250 -0.783279

Wavelength ( um )PC2Spectral bands of the ASTER (after Hunt, 1979)
0.60
0.55
0.50
0.45

Reflectance

0.40
0.35
0.30
0.25

Hem atite

0.20
0.15
0.10
0.05

Band 1

0.00
2250

500

750

1000

1250

1500

1750

2000

2250

Band 2

Band 2 / Band1

Wavelength (nm)Vegetation Indices

Thick

•

vegetation

•
•

Bare

•

soil

The most commonly used NDVI
for estimating green vegetation
cover in Remote Sensing.
Normalized Difference Vegetation
Index (NDVI)ASTER bands 3N (NIR) and 2 (R)were then converted into apparent
reflectance values.
Using apparent reflectance images
in Red and NIR bands, the NDVI
index was computed by thestandard formula:-1<NDVI<1ρred = Reflectance in red channel
ρNIR = Reflectance in NIR channel(Rouse et al., 1974)

Alteration Sample
Collection

MSR7000 Multispectro Radiometer each mineral

X-ray diffraction (XRD) analysis

powder sample measurement•

The reflectance spectra were measured using MSR7000 Multispectro
Radiometer (covers the 280 nm to 2500 nm wavelength).
Laboratory reflectance spectroscopy, MSR 7000 can be a definitive test of
the presence of hematite and kaolinite, if the absorptions appear strong.
Mixtures of minerals with overlapping absorption bands can be difficult tointerpret with spectroscopy .
1.00

0.950
0.925

0.95
0.900

Alu nite
0.90

Reflectance

Reflectance

0.875
0.85

0.80

Kao lin e

0.850

•

•
•

Spectral measurement methods are sensitive to different abundances ofmaterials especially clay minerals.
However, the minerals quartz and low iron feldspars have no diagnostic
absorption in Visible-NIR wavelength range but XRD is very sensitive tothem.
X-ray diffraction analyses confirmed that much of the silica is thedominant mineral in the alteration area.
According to XRD analyses of samples contains iron oxide (goethite,hematite), a variety of clay minerals including kaolinite, montmorillonite,
illite and siliceous minerals.

0.825
0.800

0.75
0.775
0.70

2024

0.750
500

750

1000

1250

1500

1750

2000

2250

500

750

1000

Wavelength (nm)

1250

1500

1750

2000

Kao

2250

Wavelength (nm)75 % Kaol.
25 % Dickit< 1 % andere

1724

Sericite & clay
0.9

0.60

1324

0.55
0.50

1283

1024

0.45
0.40

1698

Phyllic alteration

Reflect ance

Reflectance

959
0.8

1807 2005
1867

0.7

724

0.35

Dickit0.30
0.25

He m atite

424

0.20
0.15

1407

124

0.10

2144

0.05
0.00

0.6
400

600

800

1000

1200

1400 1600

1800

2000 2200

2250

500

750

Wavelength (nm)

1000

1250

1500

1750

2000

2250

4

12

20

28

36

44

52

60

68

Wavelength (nm)Trial of Digital Filter Photography

Optical Filter

Equipment
Required••••

The basic elements of digital filter photography include the charge
coupled devices (CCD) digital camera, filter, filter holder and ball head
tripod in this method.
The camera can store images as uncompressed TIFF format or RAW files.
One of the highlights of the Dimage 7 is that F 2.8, 7X optical zoom
Minolta GT lens. The focal range is 7.2 – 50.8 mm.
A conventional 5 mega pixel camera actually may output 2560×1920
pixel images (4915200 pixels) because some of the pixels in the camera
are used for various measurements in image processing.

Shooting filter photograph••

Shooting filter photographs is simple, byplacing various filters in front of the mainlens of digital camera.
But there are a few things to get good photos
in field work.

- First, make sure no light can leakthrough the filter attachment mechanism
and turn off automatic flash.
- Second, set the camera to the highest
ISO rating and remember to use thetripod for image registration.
- Focusing was a problem because theinfrared is so dark, it was going to takesome work before.

• CCD camera filtered through visible wavelength filters – 400 nm,
450 nm, 500 nm, 600 nm (visible wavelength) and Infrared filter 750
nm and 950 nm.
• CORION [(Holliston, MA) S25-F0470-4M229] filter with an optical
bandwidth of ± 25 nm.

Digital filter photo station

Recognizing hydrothermal alteration on Filter Photography

Filter Photography Results and Interpretation

Outcrop-3
silicification and argillization

Outcrop-2
Iron oxide (hematite, goethite)This study examined filter image processing based on alteration outcrops andnon-alteration outcrops in fieldwork.
Shoot filter photos and create to image cube or one packet image file.
The purest pixels selected from filter photographs using the 2D scatter plot
method.
Scatter plots provides a good way to show the relationship between spectral andimage space.

•
•
•

•

•

•
•

200

200

150
100
50

150
iron oxide

100
50

500

600
700
800
Wavelength (nm)

0
400

900

500

600
700
800
Wavelength (nm)900

The purest pixels showed 2D scatter plot of 750 nm and 450 nm sensitive for
Hydroxyl minerals and 600 nm and 500 nm for iron oxide minerals.
Blue pixels represent hydroxyl minerals and red color pixels represent ironoxide minerals.

Filter Photography and Principal Component Analysis

Filter Photography False Color Composite
The color composite of filter image 600:500:450 assigned to RGB
channels respectively, with histogram equalization, clearly displayed the
distribution of iron oxide pixels (pink color pixels).
750:450:400 assigned to RGB channels respectively, displayed the
distribution of highly altered pixels )light color pixels.
This result was coincided and recognized with other alteration outcrops
filter photos image processing results of study area.

Outcrop-2
Iron oxide (hematite, goethite)

Spectral Profile
250

0
400

•

•

Outcrop-4
silicification and argillization

Spectral Profile
250
Value (Sensitivity)•

Outcrop-2
Hydroxyl mineral with
sub silicification

Value (Sensitivity)

Outcrop-1
Hydroxyl mineral with
sub silicification•

The examination of PCA eigenvector loadings decided which of theprincipal component extracted information directly related to the target.
Iron oxide image was correlation with the PC image of 500 nm and 600
nm eigenvector loading and alteration image coincided with PC 1 image.
400 nm 450 nm
500 nm 550 nm
600 nm
750 nm 900 nm
PC Band 1 Band 2Band 3
Band 4
Band 5
Band 6 Band 7
1 0.303378 0.306451 0.602092 0.521392 0.415650 0.064133 0.052970
2 0.169661 0.047326 0.548430 -0.463789 -0.247157 -0.530721 -0.332187
3 -0.212583 -0.046964 0.464213 -0.460632 -0.048967 0.516112 0.506120
4 0.612482 0.575002 -0.326311 -0.386089 0.048430 0.165799 0.094130
5 0.193584 -0.021586 0.121346 0.224672 -0.611067 0.560183 -0.457874
6 0.112707 0.138909 0.003044 0.316755 -0.614961 -0.326242 0.618909
7 0.639275 -0.742465 -0.001563 -0.031498 0.099120 0.011076 0.170678

RGB = 600:500:450

RGB = 750:450:400

Iron oxide image

Alteration image

Pink color pixels

Light color pixels

Rechecking Field work

PC 5 for iron oxide image

PC 1 for alteration image

Estimate weight percent from outcrop

Iron oxide weight image

ROI too and Spectral Angle Mapper
MethodAlteration (silica, kaolinite) weight
imageSummary

Recognizing hydrothermal alteration on
ASTER Images

140 ° 40 '0" E

140° 45 '0 "E

14 0°50 '0 "E

14 0°5 5'0 " E

1 41 °0 '0 "E

4 0°0' 0" N

Hachimantai

Iron oxide image
VNIR band 2 / band 1 (15m)

Iron oxide

Study Area

3 9°55 '0 "NKaolinite, alunite,montmorillonite

Mt. Iwate

Silica mineral
3 9°50 '0 "N
14 0°40 '0 "E

1 40 °4 5' 0" E

1 40 °50 '0"E

14 0°55 '0 "E
5Hydroxyl minerals (30 m)
Principal Conponet Analysis

•
•
•
•

The best results for alteration mapping obtained little
vegetation cover outcrops and landslide.
The new filter photo method can be used in this hilly
and landslide outcrop as field spectrometer and canbe detected up to data of outcrop.

Silica minerals (90 m)
B 13 / B 12

Summary•

New digital filter photography remote sensing method is agood trial tool for detecting signs of alteration. Because- can be detected by high resolution (miga pixel).
- could be used ground truth field checking like mobile
field spectrometer.

•This method could be used to discriminate mainly the ironoxide and among hydroxyl mineral sub silicification sub
zone.

•

Iron oxide sensitivity at 500 nm and 600 nm filter
photography and hydroxyl mineral sub silicificication image
coincided with the high albedo principal component
analysis image (PC1).

•

Can be estimated the weight of alteration mineral percent
of each outcrop.

•
•

2.5

0

1 41 °0 '0 "E
5 Kil omete rs

The most appropriate index used to extract the iron oxide areas is the ratio ofband 2 to band 1 of ASTER data.
The vegetation mapping based on NDVI has shown no or sparse vegetation
cover in the areas.
Principal Component Analysis selection is based on the examination for PCA
eigenvector loading to the for hydroxyl mineral, kaolinite, alunite and illite etc..
Silica rich areas were mapped with TIR ASTER emissivity band 13 / band 12.

The 10th Center for Environmental Remote Sensing (CEReS) Remote Sensing Symposium on February 28, 2008, Chiba University, Japan

Sisir Kumar Dash, Tasuku Tanaka and Taro Ueda
Instrument and Information Engineering Laboratory, Department of Mechanical Engineering, Faculty of Engineering, Yamaguchi University
2-16-1, Tokiwadai, Ube-shi, Yamaguchi-ken, Japan, Postal code: 755-8611, Email: skdash@yamaguchi-u.ac.jp

Abstract
Dual frequency precipitation radar (DPR) is widely used to estimate drop size distribution (DSD). The output of a radar or observed physicalsubstance is the back-scattered cross section per unit volume ( ). The rainfall rate is derived from using empirical equation. If we have two
observations of different frequencies, we can retrieve, assuming the DSD, both the averaged drop size and number of rain drop in unit volume
from . The simultaneous observation of rainfall from both Precipitation Radar (PR) of TRMM and Cloud profiling radar (CPR) of CloudSat is
obtained in the month of January 2007. Using the two radar backscattering parameter, we tried to retrieve the drop size and number of particles.
We also discuss a new method to calculate the attenuation without solving the Hitschfeld-Bordan equation.
.
Introduction

Satellite data

PR/TRMM has observed the rainfall in the tropic and subtropic regions for more than a decade. In TRMM retrieval algorithm,
Surface Reference Technique (SRT) techniques is utilized. The radar reflectivity factors are corrected for attenuation, and theinformation on the DSD is inferred from the non-Rayleigh backscattering characteristics of the hydrometeors. The assessment of
the reliability of SRT is difficult because the statistical properties of the surface echoes depend on surface type, angle of incidence,and surface winds. So we assumed the new method to calculate the attenuation without solving Hitschfeld-Bordan equation and
SRT. To estimate the raindrop size distribution and number of particles, we must think of the scattering and the back-scatteredcross section per unit volume. We calculated the Mie and Rayleigh scattering in order to get the scattering character and the backscattered cross section per unit volume from radar equation. The back-scattered cross section per unit volume is corrected fora t t e n ua t i on . W e e v a l ua t e t he op ti c a l t hi c k n e s s o f bi n da ta a nd a t te n ua t e t he bi n da ta us i ng r e c ur s i on
equation.
.

T h e s i m ul t a n e ou s o b s e r v a t i on o f P R / TR M M a n d
CPR/CloudSat and their position on 25 January 2007 isshown on Table.1. Both the data were collected using orbit
c a lc ul at ion of N ASA’s Tw o-li ne - el em e nt da t a s e ts .
Satellite
CloudSattime [UTC]10:40:25

Latitude［°］-16.988

Longitude[°]48.40302

TRMM

10:40:27-16.987

48.40153

Table. 1. Observational time and location of satellites

26000

Methodology
TRMM
CloudSat

22000

A

Calculate attenuation i using
recursion equations

The observed back-scattered cross section per unit volume is a function of number of particle(D )N and backscattering cross section bshown in equation (1) and the scatter diagram isshown in Fig.1.
.
e
b

20000
18000

N

16000

b

14000

(D )

(1)

10000
8000
6000
4000
2000
0
1E-6

1E-5

1E-4

1E-3

0

.

0

0.1

1

1

10

100

1000
2

2

5
3

The backscatter cross section per unit volume(reflectivity) [mm /m ]

b

(D )

Compute

( 2)

'

b

( D)

n

( 1) n ( 2n 1)( a n bn )4'( ) n (m ) m n ( ) n ( )''( ) n (m ) m n (m ) n ( )

bn

n

an

Calculateexp( 2ei)n(m ) n ' ( )m n ( ) n ' (m )n(m ) n ' ( )We assumed that water particle exist at 4560 [m]of the distance from surface, because of the brightband appearance as shown in Fig.1. Fig 3a and 3b
shows the measured reflectivity of TRMM andCloudSat data. The simultaneous observation
cross point is found at bin number 13 on TRMM
and bin number 516 on the CloudSat data. These
points observed the precipitation. These data wereused to calculate the new back-scattered crosssection per unit volume. Fig 4 shows thecomparison of the observedback-scattered crossei
section per unit volume
and the apparent
backmiscattered cross-section per unit volume
withreference to the distance from the surface without
attenuation. In TRMM case, the observed backscattered cross-section perei unit volume nearly
unaffected by attenuation
shown in Fig.4a. In
other hand, the observedback-scattered crossmi
section per unit volume
from CloudSat comeunder the influence of attenuation at 2880 [m] asshown in Fig.4b. This result indicates that there isan existing of the precipitation from 3200 [m] due
to attenuation influence which appeared in the next
range of the layers.
.

( 4)

N e (cp 2)

e

(5)

e

Where
is the extinction
cp 2 cross-section that computes from the first layer using N and D and theradar range resolution
shown in equation (5).
Using the attenuation
Ai , the observed back-scattered cross-section permiunit volume ei iscompute from the apparent back-scattered cross-section per unit volume. The compute isshown in equation (6).
mi

(r )

ei

(r ) Ai (r ), N,Results and Discussions

(3)

39.0
[a]

using then 1

m n ( ) n ' (m )i 1

OutputandFig. 2 Flowchart of the recursion equation method.

D

Ai

and
using attenuatedei

is theWhere m is the complex
index of refraction andis the size parameter.
n( )
n( )
wavelength, and,are the Riccati-Bessel functions.
Using all these above
equations,wecanestimate
N
and
D.
In
the
first
layer,
we
don’t
need toAi
consider attenuation
But.Aifi we computeas second layer, we must consider the attenuation of
first layer. The attenuationisrepresentedinequation
(4).
.
n 1

-88.0

fromei cross point

2

2

m2 1
D6
m2 2

4

Fig.1. Scattering diagram shows the variation of backscatteredreflectivity with reference to the distance from surface.

eiCalculate attenuated(D )We estimate
from the scattering ,characteristic where D is the particle size. In this research,
we assume the raindrop size distribution range is from 0.02 [mm] to 5.0 [mm]. Due to the inverse
relation between size parameter ( ) and wavelength ( ), we also assume that TRMM is Rayleighscattering and CloudSat is Mie scattering. The wavelength of TRMM and CloudSat radars are
21.3[mm] and 3.19[mm] respectively. The backscattering cross section for Rayleigh and Mie
scattering is shown in equation (2) and (3) respectively.
.

12000

(6)

167.0

[b]

No attenuation data of TRMM
Attenuation data of TRMM

4500

Distance from surface [m]Distance from surface [m]

5000

3a.

4000
3500
3000
2500

No attenuation data of CloudSat
Attenuation data of CloudSat

4000

3000

2000

1000

2000
1500

0
0.10

0.15

0.20

0.25

0.30

back scatterign cross section per unit volume-327.67-126.5

Fig.4 Comparison of

75.0

0.35
2

3

[mm /m ]

0

10

20

30

40

50

60

back scatterign cross section per unit volume

70
2

80
3

[mm /m ]during attenuation influence (a: TRMM; b: CloudSat)

Fig. 3. The measured Reflectivity (dBZ) of PR/TRMM [a] and CPR/CloudSat [b]
25000

1E12

TRMM
CloudSat

5000

4500

Distance from surface [m]

1E10
1E9
1E8
1E7
1000000
100000
10000
1000

4500

4000

4000

Distance from surface [m]1E11Distance from surface [m]3

Nunber of particles N[number/m ]Distance from surface [m]

24000

3500
3000
2500
2000

3500

3000

2500

ConclusionssConclusion

TRMM
CloudSat20000

We have developed a retrieval methodof DSD for dual-frequency radar using
two satellites. Also we have discussed
a method to calculate the attenuation,b y n o t s o l v i n g H it c h f e l d - Bo r d a nequation and SRT but using recursion
equations.
.

15000

10000

5000

2000

100
1500

10

0
1.0

1

1.5

2.0

2.5

3.0

3.5

4.0

4.5

Average raindrop size distribution [m m]

0.1

5.0

1500
0

10

2

0

30

40

5

0

3

Number of particles per unit volum e [number/m ]

1E-8

1E-7

1E-6

1E-5

1E-4

1E-3

0.01

0.1

1

Optical thickness e

0.01
0

1

2

3

4

5

Fig. 6 Distribution of D and N in the atmosphere window.
Fig 7. Optical thickness observed by TRMM
and CloudSat data

Particle size D[mm]Fig 5 the relationship of the parameter D and N withreference to the distance from the surface.

Attenuation computed by recursion equation shown inequation (4). Using the sum oflayers of optical
thickness, the attenuation oflayer is estimated.
Substituting the total attenuation and back-scatteredcr os s- se c tion pe r u n it v o lum e
in to eq ( 6),subsequently incorporating
in to equation (1), theraindrop size distribution D and number of the particle
per unit volume N were computed. Fig 5 shows therelationship between raindrop size distribution D andnumber of the particle per unit volume N in 3200 [m].
The cross-section point determines the value of D and N.

However, in the total column of atmosphere window in therange from 1700 to 4560 m, there is wide distribution of D
and N is computed and shown in Fig. 6.
.
.
The optical thickness is used for the attenuation. Fig. 7
shows the optical thickness with reference to the distancefrom the surface. The pattern of optical thickness at
4560m is identical for both data. The increase of optical
thickness lies between 3200m to 4560m which is having
an absolute value of 0.001 and 0.45 for TRMM andCloudSat respectively. This indicates that using the highfrequency radar for observation, we must consider theattenuation by precipitation.
.

References
H.C.van de Hulst “LIGHT SCATTERING by small particles,’” Dover
Publications, Inc. New York, pp 121-130, 1957.
Redy Mardiana, Toshio Iguchi, and Nobuhiro Takahashi “A DualFreq uenc y Rain Profiling Method Without the Use of a Surface
Reference Technique,” IEEE Trans. Geosci. Remote Sensing, vol. 35,
pp2214-2225, 2004.
R. Meneghini, T.Kozu, H.Kumagai, and W.C.Boncyk, “A study of raine st im ati on m eth ods f rom s pa c e u si ng d ua l-wav el en gth ra da r
measurements at near nadir incidence over ocean,” J.Atmos. Ocean.
Technol. Vol. 9, pp. 364-382, 1992.
.

静止衛星データの全球合成技術の開発
静止衛星データの全球合成技術の開発

•
•
•
•
•
•

本共同研究の背景と目的
千葉大学環境リモートセンシング研究センター（CEReS）
全球静止気象衛星データセット作成• ああああ （4大学連携プロジェクト）
気象衛星センター
全球衛星搭載センサー相対較正（GSICS）

気象庁気象衛星センターデータ処理部システム管理課
操野 年之
太原 芳彦
奥山 新
橋本 徹
加藤 浩司 ※（ 発表者 ）全球合成技術の開発と精度評価
気候分野へのデータ利用促進

GSICS• Global Space-Based Inter-Calibration System– 全球衛星搭載センサー相互較正システム– WMOが推進する静止、極軌道気象衛星センサーの相互較正を現業的に行う国際的な枠組み• 参加機関– NOAA, EUMETSAT, JMA, KMA, CMA,CNES– 気候分野、数値予報同化ためのバイアス除去など幅広い応用

高スペクトル分解能サウンダの利用（AIRS/IASI）• 1000以上のチャネルを持ち非常に高い波長分解能を持つ赤外探査計
Name
AIRS
IASI

CH数
2378
8400

FOV[km] バンド[μm]13.5 3.74 - 4.61, 6.20 - 8.22, 8.80 - 15.4
12 3.62-15.5

ギャップがある• 第一段階： 静止衛星の赤外チャネルの校正– 高スペクトル分解能サウンダ（hyper sounder）を利用した校正– 各静止気象衛星運用センターが比較校正を行い、結果を公開す
る– MTSAT/GOES/METEOSAT/FY vs. AIRS/IASI

AIRS

IASI

GSICSコロケーション条件• 観測時間差
15分以内• 衛星天頂角θ Abs( cosθgeo/cosθleo)-1 < 0.03
• FOVが一様
MTSAT Pixel– FOV内放射量の標準偏差が小さい– Z検定

AIRSとMTSAT-1R(IR2)の応答関数

Hyper Sounderから静止衛星放射量の推定
１.
１.
ハイパーサウンダーを用いて静止ハイパーサウンダーを用いて静止衛星の応答関数を作成
衛星の応答関数を作成（Super（SuperChannel）
Channel）

MTSAT-1R

２.不良チャネルとGap
２.不良チャネルとGap Channels
Channelsの放射量を計算済み放射伝達計算の放射量を計算済み放射伝達計算
を用いて補完(AIRS)を用いて補完(AIRS)

Hyper Sounder

Check Range

FOV Range

Sounder FOV

Super Channelの作成
MTSAT-1R
MTSAT-1R 6.8-um
6.8-umAIRS 不良チャネル

Gap channelsの応答関数

AIRS応答関数

Super Channelの最適化計算• ωiに関する最適化問題を解きSuper Channelを作成
Radiance observed by a broadband channel is

Radiance of a super channel (linear combination of sounder radiances) is
MTSAT応答関数
Gap Channels とAIRSの各チャネルで応答関数をシミュレート

They should be approximately equal for any I (v), then

To obtain optimized wi, solve

各チャネルの重み

where

計算済み放射伝達計算を用いた補完
MTSAT-1R
MTSAT-1R6.8-um
6.8-umregion
region比較結果と考察

計算済みAIRS放射量• 衛星 センサー– MTSAT-1R搭載 JAMI 赤外センサー– Aqua搭載 AIRS
– Metop-A搭載 IASI

• 期間 2007年8月 （1ヶ月）
AIRS観測放射量

AIRS/IASI
August 200709 – 10 JST
IASI (mW/m .sr.cm )2

vs.
21 – 22 JST

AIRS/IASI

IASI (mW/m .sr.cm )

-1

2

-1

MTSAT

MTSAT-1R 6.8um

MTSAT

vs.

MTSAT – IASIMTSAT-1R 10.8um

MTSAT – IASI

放射計算によるAIRS放射量09 – 10 JST

21 – 22 JST

IASI (mW/m .sr.cm )IASI (mW/m2.sr.cm-1)

2

-1

AIRS
AIRS 補完あり
補完あり

12 – 13 JST
AIRS (mW/m2.sr.cm-1)

00 – 01 JST
AIRS (mW/m2.sr.cm-1)

AIRS
AIRS 補完あり
補完あり

MTSAT

IASI
IASI 補完なし
補完なし

MTSAT

MTSAT – AIRS

IASI
IASI 補完なし
補完なし

MTSAT – AIRS

August 2007

12 – 13 JST

00 – 01 JST

AIRS (mW/m2.sr.cm-1)

AIRS (mW/m2.sr.cm-1)

August 200709 – 10 JST
IASI (mW/m2.sr.cm-1)

vs.
21 – 22 JST

MTSAT

MTSAT-1R 3.8um

21 – 22 JST

AIRS/IASIMTSAT – IASI

AIRS/IASIMTSAT – IASIvs.

MTSAT – IASIMTSAT-1R 6.8um

IASI (mW/m2.sr.cm-1)

IASI (mW/m2.sr.cm-1)

21 – 22 JST
IASI (mW/m2.sr.cm-1)

12 – 13 JST
AIRS(mW/m 2.sr.cm-1)

00 – 01 JST
AIRS (mW/m .sr.cm )
2

-1

考察

AIRS
AIRS 補完あり
補完あり

MTSAT

IASI
IASI 補完なし
補完なし

00 – 01 JST
AIRS(mW/m2.sr.cm -1)

MTSAT – AIRS

AIRS
AIRS 補完あり
補完あり

MTSAT – AIRS

IASI
IASI 補完なし
補完なし

MTSAT – AIRS

August 2007

00 – 01 JST
AIRS (mW/m2.sr.cm-1)

今後の計画• Super Channelは高い精度で応答関数を再現する。

１． 来年度から本手法を用いた
MTSAT-1Rの赤外リキャリブレーションを開始• AIRSとIASIで殆ど結果に変わりないことから、放射計算による補完は妥当であると推定される。

２． 赤外リキャリブレーションデータの評価手法を開発
３． MTSAT, FY-2, GOESの赤外画像合成データ作成
４． 可視データ較正技術の確立

干渉合成開口レーダ(InSAR)による南極の氷河・氷床マッピング−等速変動場の仮定を用いない流動位相の導出−河田 哲郎(岐阜大),木村 宏(岐阜大),西尾文彦(千葉大)目的[2008年度の目的]

干渉処理における仮定

masterInSARの原理

slave

干渉位相：

・様々な仮定における流動位相の導出

B(r2-r1)/

・衛星軌道・基線の外挿による広域の変動位相の導出

r1

f(r1,h,B,d)r2

変位dの寄与・ALOS PALSARデータと比較検証 のみを検出する為の準備衛星InSAR南極モニタリング

d
h

露岩のある
シーン

隣のシーン

変動位相

上記の等速変動の仮定は厳密には成立しない→ 干渉ペアによる基線長に関係無く，別個に流動位相の導出，初期干渉位相からの除去を行い，基線推定を行なう隣のシーン

テストエリア(ヤマト周辺)

衛星軌道・基線推定（軌道データの誤差修正）

衛星軌道・基線の外挿

南極の氷河・氷床の流動は急激な変化はないので，期間が短い干渉位相同士の流動位相成分は等しい→ 基線の長い干渉ペアの基線推定において流動
位相が，基線の短い干渉ペア96年4月と6月の流動位相と同じになるように基線推定をチューニング

・非等速変動の仮定

南極InSAR処理のフローチャート[期待]・InSARでは変動計測が可能
・連続した観測データがあれば、変動の時間変化を計測可能（JERS-1 SARの場合４４日間隔）[問題点]
・JERS-1SARでは軌道情報精度が不十分 → 基線推定
基線は，初期干渉位相から流動位相を除いた成分より推定・等速変動の仮定200ｍDEM

衛星軌道・基線
の外挿

変動位相

変動位相

仮定による変動位相の違い
等速変動の場合(1kmDEM使用)
非等速変動の場合(200mDEM使用)

197-421

+6 cm

+6 cm

0 cm

0 cm

-6 cm

-6 cm

197-422
レーダ照射
方向

飛行方向

197-423

coherence

96.03~04

96.04~06

96.06~07

96.03~04

96.04~06

96.06~07

まとめ

今後の予定

・非等速変動の仮定における流動位相は様子が全く異なる
誤った結果である → 基線の長い干渉ペアにおける地形位相，軌道位相の除去
精度の問題？
正しい結果である → 時間経過に則した氷河・氷床流動の変化？

・ALOS PALSAR データとの比較検証
・非等速変動の仮定における流動位相結果の正誤判断

200８年2月26日 作成

合成開口レーダによる道路凍結モニタリング＝飛行場・滑走路における積雪・凍結路面の検出＝
長 康平 （株式会社ウェザーニューズ）

Josaphat ʻJoshʼ Tetuko Sri Sumantyo （千葉大学）

１．概要

４．画像特性解析

本研究は、合成開口レーダによる積雪、凍結等の路上・路面状態モニタリングの道路維持管理・雪氷対策業務等への応用における技術的課題を明らかにすることを目的とする。
今期は、高速道路等の調査に先立ち、合成開口レーダ観測データを用いた、路面域の検出、路面の雪氷の検出等の可能性や課題を検討するため、比較的平坦でマルチパス等の原因となる施設・地形の少ないであろう、飛行場の滑走路等を対象に調査を行った。
また、観測当日の天候状況や滑走路の路面状態等の情報も別途入手を試みた。
本報告は、これらの調査のうち、新千歳空港（N42°46.5’ / E141°41.5’）、および、千歳飛行場（N42°47.7’ / E141°40.0’）を対象に実施した調査事例をまとめたものである。

合成開口レーダ再生処理を行ったデータの空港周辺、滑走路周辺、各滑走路の切り出
し領域毎に、画素値出現頻度（領域全体の画素数に対する、各画素値毎の画素数の割合）のヒストグラムを作成し、受信電圧相当の画素値(DN)、後方散乱係数相当のdB 値(DNdB)の分布特性を調査した。なお、dB 値(DNdB) については、画素値(DN) から規格化後方散乱係数σ0 への換算式を参考に、次式にて計算した。

( http://net.jmc.or.jp/saishiki/index.asp )( http://gc.kls2.com/airport/CTS )( http://www.stv.ne.jp/webcam/shinchitose/index.html )

千歳飛行場

空港周辺
新千歳空港

滑走路周辺

空港ターミナル屋上カメラ映像
2008年1月14日 00:59UTC

滑走路 (RJCC / 01L)

DNdB = 20 log10(DN)[dB]

それぞれの事例について、画素値(DNdB) [dB] を横軸にとったヒストグラムを図−４に示す。これらより、空港周辺領域に比べ、滑走路周辺領域、各滑走路領域では、画素値(DNdB)の小さい画素が、より多く分布していることが分かる。

滑走路

滑走路

滑走路間

滑走路間

滑走路周辺

滑走路周辺

空港周辺

空港周辺

Fig-1 空港周辺地図および空港ライブカメラ映像

２．観測データ

１偏波・高解像度観測
2008年1月14日 01:07UTC（新千歳空港）

本調査では、陸域観測技術衛星「だいち」搭載の合成開口レーダ(ALOS / PALSAR) による、高分解能モード観測の観測生データ( Level 1.0 ) を用いた。
観測日時
シーンＩＤ
観測モード
オフナディア角2008 年1 月14 日01:07UTC ALPSRP105022760 １偏波(HH)・高分解モード
41.5
2006 年10 月6 日01:01UTC ALPSRP037192760 ２偏波(HH/HV)・高分解モード
34.4

レーダ観測当日の天候は、気象庁の地域気象観測によれば、千歳において、2008 年1
月14 日には、当日午前３時から６時にかけて、時間2cm 前後の降雪があり、レーダ観測
時刻には、17cmの積雪が観測されている。一方、2006 年10 月6 日には、前日、前々日も含め、当日は降水は観測されていない。また、2008 年1 月14 日の新千歳空港の雪氷状
況は、航空局の航空情報（雪氷状況）によれば、レーダ観測時刻前後において、新千歳空
港の両滑走路とも、ところどころ乾いた雪に覆われ、積雪は1mm 程度、摩擦係数は概ね0.4 以上の状態にあった。
N 42.884
E141.377

N 42.783
E142.170

N 43.131
E141.104２偏波・高解像度観測 （ＨＨ）
2006年10月6日 01:01UTC（千歳飛行場）Fig-4 画素値出現頻度ヒストグラム

５．路面域抽出の試行
画素値分布特性を参考に閾値を設定し、２値化処理を行って、路面域の抽出を試みた。
滑走路周辺の画素値分布を参考にDN = 60 ( = 35.6 [dB]) を、また、滑走路の画素値分
布を参考にDN = 30 ( = 29.5 [dB]) を閾値に設定した。次に、これらの閾値を用い、0≦DN ＜30、30 ≦DN ＜60 の画素を抽出する２値化画像を作成した。さらに、各々の画像に以下の配色を設定し、合成画像を作成した。これにより得られた合成画像を図−５に示す。
画素値(DN)
0 ≦DN ＜3030 ≦DN ＜60

N 43.041
E141.910

配色
Red(R)
Green(G) + Blue(B)1 月14 日の事例においては、合成開口レーダ観測データの２値化画像処理により、個々の滑走路、誘導路、駐機場等の領域を抽出可能なことが確認できた。特に、除雪がなされ、数ミリ程度の乾質の雪に覆われた滑走路と、十数センチの積雪の残る滑走路沿いの区画における画素値分布特性の差異は、ヒストグラム、２値化処理画像ともに明瞭に識別可能であることが確認できた。

ＨＨ
N 42.382
E141.263

50km

N 42.280
E142.04950km

N 42.567
E140.992

１偏波・高解像度観測
2008年1月14日 01:07UTC

ＨＶ
50km

N 42.477
E141.791

２偏波・高解像度観測
2006年10月6日 01:01UTC

Fig-2 観測生データ( Level 1.0 ) の合成開口レーダ再生処理画像

３．データ処理・解析
本調査では、以下、合成開口レーダ再生処理等のデータ解析を行い、衛星搭載合成開口レーダ観測による路面検出手法、および、路面の雪氷検出に関する調査を行った。
１．観測生データ( Level 1.0 ) の合成開口レーダ再生処理（図−２）２．解析対象領域データの切り出し
３． ２偏波合成画像の作成（図−３）４．画像特性の調査
５． ２偏波画像特性の調査
６．路面域の抽出

3km

3km

１偏波・高解像度観測
2008年1月14日 01:07UTC

ＨＨ

3km

ＨＶ

２偏波・高解像度観測
2006年10月6日 01:01UTC

Fig-5 ２値化処理合成画像

２偏波・高解像度観測
2006年10月6日 01:01UTC

3km

偏波
HH
HV

配色
Red(R)
Green(G) + Blue(B)50km

Fig-3 ２偏波合成画像（ＨＨ／ＨＶ）

６．まとめ
以上の調査結果より、衛星搭載合成開口レーダ観測データの２値化処理を用いた路面域抽出手法、同データによる路面の積雪状態検出可能性の確認がなされた。
今後、抽出された路面域の画素分布特性、後方散乱係数分布特性を調査し、積雪状態の定量的把握の可能性を検討する。また、多偏波観測についても、積雪時の事例を蓄え、定量的把握のための散乱モデル化を検討する。さらに、これら手法の高速道路等の路面状態の把握への応用を検討する。

Monitoring Kashiwa City Using Remote Sensing
Katsumi Ohyama1, Luhur Bayuaji2 and Josaphat Tetuko Sri Sumantyo2
1Center for Environment, Health, and Field Sciences, Chiba University
2 Center for Environmental Remote Sensing, Chiba University
E-mail: k_ohyama@faculty.chiba-u.jp

Introduction
Regional sustainable development requires a balance among diverse social goals, for instance environmental conservation
and/or improvement vs. economic development (Zhang and Guindon, 2006). The fulfillment of single factor is not sufficient toimprove regional sustainability, and thus, comprehensive investigation with respect to the social goals is essential. A
comprehensive investigation is more important in a developing region than in a developed region because of possible
alterations in its development plans. For further improvement in the regional sustainability, transitions in the land use andenvironmental indicators (e.g., vegetation index and surface temperature) between the past and present need to be compared
quantitatively in order to influence the results of the future development plans.

Objective

Tone River

The objective of this study is to demonstrate that satellite remote sensing can be used for
determining land use, and thus, for evaluating regional sustainability at a low cost.

Target Area

Edo River

Images of an area neighboring Tsukuba Express (railway) in Kashiwa City (northwest
region of Chiba Prefecture, Japan) acquired by the Advanced Spaceborne Thermal
Emission and Reflection Radiometer (ASTER) on the Earth observation satellite TERRA(EOS-AM1) were used. In Kashiwa City, especially in the area neighboring Tsukuba
Express line, land use has been changing dramatically. The introduction of the railway
increased accessibility to Tokyo, resulting in rapid construction of buildings (e.g., shopping
malls and houses), increase in the population, and decrease in the cultivation area and/or
green area. Hence, Kashiwa-no-ha Campus Station (35°53′34″N, 139°57′9″E)
and its surrounding area (radius: 5 km) was intensively analyzed.

Kashiwa city
Nishichiba Campus
Chiba University

Tokyo Bay

Results and Discussion
Normalized Differential Vegetation Index (NDVI)
05/16/2000

04/25/2004

04/01/2001

09/03/2005

06/04/2001

05/01/2006

09/24/2001

10/13/2002

09/30/2003

Around the Kashiwa-no-ha Campus Station(center of the images), the NDVI was higher in2000 than in 2006. While the campus station
previously contained a golf course and a forest,
drastic changes in land use occurred after the
introduction of the railway. Thus far, most of
the vegetation has been destroyed. Further,buildings are continually being constructed.

11/09/2006

Land Surface Temperature (LST)
05/16/2000

04/01/2001

06/04/2001

09/24/2001

04/25/2004

09/03/2005

05/01/2006

11/09/2006

10/13/2002

09/30/2003

Although obvious differences among the LST
images were not observed, the surface and airtemperatures around the Kashiwa-no-ha Campus
Station can increase in the near future due to thedecrease in vegetation, which facilitates sensible
and latent heat exchange. To prevent theincrease in the air temperature in the area, i.e.,the so called “heat island phenomenon,” theamount of vegetation should be maintained orincreased while simultaneously striking a balance
between vegetation and city development.

Conclusion Remarks
The NDVI and LST around the Kashiwa-no-ha Campus Station (Tsukuba Express line) were obtained fromthe images acquired by ASTER on the Earth observation satellite TERRA. As we mentioned above, in Kashiwa
City, especially in the area neighboring Tsukuba Express line, land use has been changing dramatically.
Hence, some diagnostic tool for evaluating the land use was required for the sustainable development of thecity. Continuous monitoring of the target area using remote sensing will contribute to sustainable developmentof the city as well as increase the quality of life of the people living in this area.

都市域におけるアルベードの計測
菅原広史（防衛大地球海洋）
１） 研究の背景
都市域のアルベードは，ビルキャノピー内の多重散乱（ある しかし都市域では複雑な地表面形状のために反射光が非等方的
いは日陰の効果）によって低下することが模型実験や数値計に分布しており，FOVの小さいセンサーはアルベードの計測には算 に よ り示さ れ てい る（ Aida, 1982; Aida and Gotoh, 1982;

不向きである．Yoshida et al. (1995)は都市を模擬した模型においKondo, 2001）．他の土地被覆と比較して都市域のアルベード て絶対値の50%程度に及ぶbidirectional reflectanceを実測していが小さいことは衛星計測からも指摘されている（Nakagawa and る． そこで本研究では都市域を主な対象として航空機観測を行Nakayama, 1995; Brest, 1987）．

い，非等方性反射の影響を受けないアルベードを計測した．

２） 計測
図１ フライトコース

ヘリコプターに搭載した分光放射計（英弘精機 MS200）により
475, 525, 675, 750 nmでの上向放射量を計測した．計測高度は地上610mである．計測した上向き放射量から地表面の分光反射率を求めるために，エアロゾルの鉛直分布（500〜3000m）と地上で↓図２ フライトコースに沿った建物の形状パラメータの分布．グ
ロス建蔽率は敷地面積に対す
る屋根面積の割合．天空率は道路から空が見える割合（小さくなるほど狭い街路）．

の光学的厚さの計測も行っている．今回解析するデータは1994年
8月5日14:00前後（太陽天頂角36.9°）に札幌において計測したものである．フライトコースを図１に示す．

３） 解析方法
計測されたエアロゾルの粒径分布を用い，エアロゾルの複素屈折率を1.5-0.00i （Takamura et al., 1983）として，各高度での消グロス建蔽率

形状パラメータ

1.0

散係数を算出した．これにより放射量の計測高度（610m）での大気の光学的厚さが算出される．計測高度（610m）と地表面間での

0.6
0.4
0.2
0.0

放射伝達をtwo stream近似により解くことで，計測された上向き510

放射量から地表面での分光反射率を求めた．

反射率

南端（図右端）で750nmの反射率が大きくなっているのはそれぞれ北海道大学のキャンパスと真駒内付近の森林である．土地利

15
10
0
0.8

NDVI

い）のスペクトルが入っている．図５には中高層ビル街でのスペク
トルと都市を構成する建築材料のスペクトルを比較した．建築材
料も多種多様ではあるが，今回の実測値は総じてそれらより5〜7

0.6
0.4
0.2

割低い値となった．これは建物キャノピーによる多重散乱（日陰の0.05効果）にあたる．非等方性反射の影響を受けない反射率においての面積が強く影響していることがわかる．そこで，土地被覆を

めざす．

10

15
20
海岸からの距離[km]2530

図３ フライトコースに沿った分光反射率とNDVIの分布

も，キャノピー構造によるアルベード低下の効果が確認された．
図４から森林を含めたこのエリアのアルベードには，植生

デルをさらに精緻にし，石狩平野でのアルベード分布推定を

475
525
675
750

20

1.0

いう２つのグループに分かれ，その中間に農地や北大（植生が多計測された分光反射率を再現したところ，
RMSE1.8%(reflectance)で実測値と一致した．今後はこのモ

30

5用ごとの反射率スペクトルを図４に示す．大きくは森林と市街地との反射率スペクトルを推定するものである．このモデルにより，

25

25

750nmと675nmでの反射率から求めた．コース中央付近とコース

中高層ビル群
低層ビル群
住宅地
農地
森林
北大

25
20
15
10
5

reflectance [%]30

reflectance [%]成した．このモデルは植生面割合が既知として，そのエリア

20

30

得られた分光反射率と NDVIの分布とを図３に示す．NDVIは

射率を植生面積で重み付け平均する分光反射率モデルを作

15

海岸からの距離[km]４) 結果

中高層ビル街と植生からなる混在地表面とし，両者の分光反

天空率

0.8

0
0.4

0.5

0.6
0.7
wa v e length [micro meter ]0.8

図４ 土地利用ごとに平均した反射率スペクトル

アスファルト面

45
40
35
30
25
20
15
10
5
0
0.3

ビルキャノピー実測

コンクリート面

0.5
0. 7
wave len gt h [m ic ro m eter]

図５ 中高層ビル群と建材の反射率スペクトル

0. 9

ライダーによる粒子径分布計測のための多重散乱モデル
はじめにライダーによる粒子径分布計測のための多重散乱モデル

本研究ではライダーを用いて雲の粒子径分布を計測することを目指している。レーザーを雲に照射すると、
散乱光の２次元画像が得られる。その画像の偏光および強度分布は、粒子径分布を反映していると考えられる。
↓
その２次元画像から粒子径分布を得るためのモデルを構築することを目的とし、以下の研究を実施した。
・標準粒子を用いた室内実験で、モデルを検証する。
・雲にレーザーを照射し、散乱光の２次元画像とモデルによるシミュレーション結果との比較検討を行う。

千葉大学環境リモートセンシング研究センター
平成１９年度共同利用研究発表 （2008.2.28）
P偏光イメージ

S偏光イメージ

岐阜大学工学部
吉田 弘樹
1

Yos20080205-2

本研究の多重散乱モデル

異なる粒子径間の散乱を無視したモデル

Photon receiving position x（mm）（a）Small particles→

The number of the returning photons N（photon）

Photon＋

The number of the returning photons N（photon）

The number of the returning photons N（photon）

Light path

Photon receiving position x（mm）（b）Large particles

ミー散乱

PhotonPhoton receiving position x（mm）（c）Mixed particles

43

多重散乱モデルによる
モンテカルロ法シミュレーション
ＳＴＡＲＴ

Light path

光子入射
z=0

後方散乱光測定実験
Cloud
thickness

散乱(距離,角度,偏光)

ミー散乱計算

散乱断面積と粒子数密度分布に従って、光子が次に入射する粒子のサイズを決定。
散乱角の決定。ストークス・ベクトルの計算。

光子数ｎ＝ＮＹｅｓ

Beam splitter

z <0

l
P

T0

4mm

Cover

4ml

Cover glass Suspension

No

Yes

A (mm 2)

Ｙｅｓ
ＥＮＤ

Cover

R0

He‐Ne laser

λ=543nmP 0=500μW

Ｎｏ
0<z<cloud thickness
Ｎｏ

4mm

A (mm 2) Power meter

条件設定

粒子半径、屈折率、粒子径分布密度分布 、波長、総光子数Ｎ

高速化のための、微分散乱
断面積テープルの作成

Cover glass

Cap

z

d

l

A
2 l2

ｎ＝ｎ＋１
5

6

単一粒子径サンプル

粒子

数密度0.01，0.018，0.1，0.2
0.3，0.4，0.75，1.0(μm-3)

1.0×10 -4，3.15×10 -4
φ2μmポリスチレン標準粒子 6.5×10 -4，13.0×10 -4(μm-3)

1.5×10 -5，3.0×10-5
φ4μmポリスチレン標準粒子 6.0×10 -5，12.0×10 -5(μm-3)

φ0.44μmアクリル粒子

単一粒子径サンプル

φ0.44μmアクリル粒子のＳＥＭ写真
7

単一粒子径サンプルの後方散乱光（φ0.44μmアクリル粒子懸濁液）

ポリスチレン標準粒子(Std. Dev
1.0%)
Yos20080205-8
Duke Scientific Corporation

異なる粒子径のサンプル（混合懸濁液）MixtureCondition
φ0.44μmアクリル懸濁液(1.0×10-2(μm-3)) 1ml

A

φ2μmポリスチレン標準粒子懸濁液(1.3×10-3(μm-3)) 1ml
φ4μmポリスチレン標準粒子懸濁液(1.2×10-4(μm-3)) 1ml

φ0.44μmアクリル懸濁液(1.0×10-2(μm-3)) 2ml
B

φ2μmポリスチレン標準粒子懸濁液(1.3×10-3(μm-3)) 1ml
φ4μmポリスチレン標準粒子懸濁液(1.2×10-4(μm-3)) 1ml

910

Lidar屋上設備の構成

混合懸濁液（Ａ，Ｂ）の後方散乱光

View
Port
Telescope
280mm
PC

Nd:YAG(532nm, 10ns, 180mJ)

SpectroPMT
PMTmeter1m1211

多重散乱モデルシミュレーションで得られる散乱光２次元画像

散乱光２次元画像（Ｐ偏光）

平均粒子径： 0.1 m

平均粒子径： 1 m

Camera; Nikon D70S, 1s
Lens;
Tamron AF 200-400mm
Laser;
Surelite, Nd:YAG, 532nm, 10ns, 50mJ

平均粒子径： 10 m

14
13

結論
• 多重散乱モデルを構築し、モンテカルロ法による計算機シミュレーションで散乱光の偏光および強度を２次元で得られるようにした。
• 標準粒子のサンプルを用いた室内実験とシミュレーションとの後方散乱光強度を比較し、誤差が １０％であることを検証した。
• ライダーを用いて雲の散乱光２次元画像を測定した。
シミュレーション結果と比較すると、平均粒子径〜数mと推定される。
今後の課題：計測画像とシミュレーションとの誤差が少なくなるようにイタレーションすることで、粒子径分
Yos20080205-15
布を得る。

高空間分解能衛星データを用いた森林立ち枯れ領域判定法の開発
一般研究

高空間分解能衛星データを用いた森林立ち枯れ領域判定法の開発
松島 大・石戸大介 千葉工業大学工学部
樋口篤志
千葉大学環境リモートセンシングセンター

１．背景・目的
研究の背景
・モンゴル国では北部・山岳地域を中心にカラマツ林を主とする森林が存在する。
・森林面積は必ずしも広くなく国土全体のわずか８％である(UNEP, 2002)。
・近年、森林火災や蛾の幼虫による葉の食害による森林減少(UNEP, 2002)。
・「立ち枯れ」の様相を呈する。特に2003年に大規模に生じた。
・立ち枯れ域の植林をはじめている地区もある（フブスグル湖周辺）。
・立ち枯れ領域の面積を地上だけから把握するのは困難。
・緑の森林と立ち枯れた森林では視覚的にかなり異なる様相を呈するので、可視域のリモートセンシングにより判別できる可能性がある。
本研究の目的・モンゴル国北部の森林立ち枯れ領域を衛星可視光域データから判別する方法の開発
２．方法・概要（１）航空機観測によるビデオ画像から典型的な土地被覆（立ち枯れ・緑の森林・草原・裸地・水面等）を選び、ビデオ画像のRGB値とそれぞれの撮影地点近傍における衛星画像のRGB値を比較する。
（２）土地被覆ごとのRGB値の傾向によって、衛星画像のRGB値による土地被覆の判別
アルゴリズムを作成し、実際の衛星画像に適用する。
・データ
１．航空機からのビデオ画像データの解析
実施場所：モンゴル国トブ県・ヘンティ県（北緯46.6-48.4度,東経106.8-110.7度, 図３）
実施期間：2003年7月19,20,23日,8月21,22,23日
使用機材：航空機(小型単発複葉プロペラ機・AN-2, 図２), デジタルビデオカメラ(松下電器, NV-GS5), GPS受信機(古野電気, GP-31)
データ取得・処理方法：(1) 航空機底部にほぼ真下の地面を撮影するようにデジタルビデオカメラを設置(2) 飛行中はビデオカメラは地表面の連続撮影を行う。
(3) 飛行中にGPS受信機の出力をPCに出力させ記録する(4) 飛行中にビデオカメラとGPS受信機の時計のずれを線形補間する。
(5) 各土地被覆の静止画を作成、画像解析ソフトで処理し、RGB値を解析する。
２．Landsat7-ETM+ の放射輝度の解析
・画像の時期・領域：2003年7月29日,8月28・30日の日中の画像を用いた。8月28日の領域以外は航空機が飛行した空域が含まれている（図１参照）。
・フォーマット・画素：L1Gフォーマット・１ピクセル当たり25ｍ。
・波長：Ch.1(450-520nm), Ch.2(520-600nm), Ch.3(630-690nm)
・デジタルナンバーからRGB値への変換（大田,2003）(1) “Landsat7 Science Data Users Handbook (NASA発行)”によって、デジタルナン
バーから反射率に変換（この際、レーリー散乱除去の大気補正を施した）。
(2) 以下の式によりCIE-XYZ表色系の値に変換

X

k R( ) E ( ) x ( )dY

k R ( ) E ( ) y ( )dZ

k R( ) E ( ) z ( )dk

100 / E ( ) y ( )d・判別対象の土地被覆（図５）
本研究では、立ち枯れ、緑の森林、草原、裸地の４種を判別する。

図５ 本研究の解析に用いた判別対象となる画像の例。左から立ち枯れ、緑の森林、草原、裸地。

３．結果
・航空機観測値と衛星輝度値によるRGB値の比較（図６）
両者はほぼ線形に対応している。

図６ 左から立ち枯れ、緑の森林、草原の場合のRGB値。

・衛星輝度によるR値とG値の比較（図７）
森林：R＜60, R＜G。草原：R＞60, R＞G。立ち枯れ：R＜65, R＞G
裸地・水面：R＞80, R＞G+10

図７ 衛星輝度によるR値とG値の比較。

・土地被覆判別アルゴリズム＊本研究では基準値は60

衛星RGB値
N

R＞基準値＊

R,G,B＜90
R＜G

Y

R＜G

R＜B & G＜B
判別不能

立ち枯れ

緑の森林

草原

裸地

図８ 土地被覆判定アルゴリズム。

・土地被覆分布（図９）xy z は等色関数（図４）。
ただし、R(λ)は分光反射率, E(λ)は分光輝度,(3) CIE-XYZ表色系からRGB値に変換（ディジタル画像処理編集委員会,2004）
1 / 2.2
R

255(0.035064 X

0.0174YG

255( 0.01069 X

0.019777Y

0.005441Z )0.000352Z )1 / 2.2

B

255( 0.000563X

0.00197Y

0.010511Z )1 / 2.2

2003.08. 30

2003.08.28

図９ 土地被覆判別分布図。黒：緑の森林、赤：立ち枯れ、緑：草原、青；
裸地、黄：判別不能（雲など）。
2003.0
7. 29

４．今後の課題
・草原域と立ち枯れ域におけるR値のオーバーラップ(60〜65)の解決

ウランバートル

図１ 対象地域の地形図及び衛星撮像領域。
薄緑が森林域。

図２ 観測に用いた航空機(上)、及び機内におけるビデオカメラ設置状況（下・赤丸内）。

・各土地被覆の経年変化の推移の推定
参考文献
ディジタル画像処理編集委員会：ディジタル画像処理, 財団法人画像情報教育振興協会,
p.384, 2004.
大田 登：色彩工学第２版, 東京電機大学出版局, p.310, 2003.
UNEP：State of Environment, Mongolia 2002, p.113, 2002.

図３ 航空機観測の経路例(8月22日)

図４ CIE-XYZ表色系の等色関数

平成１９年度共同利用研究発表(口頭） 2008年2月28日
9:25

センター長挨拶(千葉大学けやき会館３F レセプションホール）
Proj
番号

共同利用研究者（所属）時間
am 9：30−12：15

番号発表タイトル

対応教員

口頭発表は発表時間 11 分、質疑 4 分です鈴木力英（(独)海洋研究開発機構地球環境フロンティア研
全球NDVIの降水量，気温，放射量の経年変化に対する
9:30-9:45 究センター）
P2007-1 O-1
レスポンス

近藤

9:45-10:00 開發一郎（広島大学大学院総合科学研究科）

モンゴル高原とその周辺の土壌水分・植生の時空間変P2007-1 O-2 動

朝隈康司（東京農業大学 生物産業学部アクアバイオ学
10:00-10:15 科）

P2007-3 O-3

衛星画像を用いたオホーツク地域のエアロゾル光学特性の導出

久世

金子大二郎（松江工業高等専門学校 環境・建設工学
10:15-10:30 科）

P2007-1,3 O-4

アジアの植生・作付分類と光合成（CO2固定）速度平面
分布の推定

建石

梶川昌三・筒井俊洋（国土地理院地理調査部環境地理
10:30-10:45 課）

P2007-1 O-5 地球地図・土地被覆データ整備における都市域の抽出

建石

休

10:45-11:15

近藤

憩
雲解像モデルとマイクロ波放射計の物理量のアンサンブ
ル予報誤差分布と共分散

11:15-11:30 青梨和正（気象庁気象研究所予報研究部第２研究室）

P2004-3 O-6

11:30-11:45 塩原匡貴（情報・システム研究機構 国立極地研究所）

P2004-3 O-7 船舶搭載用新型スカイラジオメータの開発

久世

11:45-12:00 青木一真（富山大学大学院理工学研究部）

P2004-3 O-8 太陽放射観測による父島のエアロゾルの光学的特性

久世

12:00-12:15 村松加奈子（奈良女子大学共生科学研究センター）

P2007-2 O-9
昼

12:15-13:00

人工衛星データを用いた植生純一次生産量推定の精度向上のための基礎研究 II

高村

梶原・本多

食

pm 13：00−16：30
合成開口レーダーによる洋上風力発電ポテンシャルの評
久世・高村価

13:00-13:15 香西克俊（神戸大学海事科学研究科）

一般研究 O-10

13:15-13:30 伊勢崎修弘（千葉大学大学院理学研究科）

SAR offset displacement from Sumatra-Andaman
一般研究 O-11 Earthquake Mw 9.2 on December 26, 2004(発表者：
Ashar Muda Lubis）

西尾

13:30-13:45 若林裕之（日本大学工学部情報工学科）

P2007-1 O-12 多偏波SARデータを用いた海氷物理量測定

西尾

Trial of digital filter photography for alteration mineral
高島 勲（秋田大学工学資源学部 付属環境資源学研究
13:45-14:00 センター）
一般研究 O-13 detection in the Hachimantai area, NE Japan (発表者：
ミン スー）

建石

14:00-15:00 POSTER（※2枚目のプログラム参照)
休

15:00-15:15

憩 （ポスターの片付けをお願いします）

浅野正二（東北大学大学院理学研究科 附属大気海洋変
北西太平洋域における下層雲の衛星リモートセンシ ング
15:15-15:30 動観測研究センター）
P2007-3 O-14の航空機検証観測（発表者： 佐藤 茉莉)

高村

15:30-15:45 木下紀正（鹿児島大学産学官連携推進機構）

久世P2007-3 O-15 2007年春季黄砂と視程映像観測

操野年之（気象庁 気象衛星センターデータ処理部システ
15:45-16:00 ム管理）
一般研究 O-16 静止衛星データの全球合成技術の開発

高村

16:00-16:15 服部克巳（千葉大学大学院理学研究科）

P2007-5 O-17 衛星を用いた地殻変動監視に関する研究高村

16:15-16:30 鷹野敏明（千葉大学 大学院自然科学研究科）

P2007-3 O-18

O・・・ORAL P・・・Poster「ミリ波雲レーダ FALCON-I による観測-雲物理量導出と最近の結果」※ 16:30-17:30の間、CEReS見学やプロジェクトの打ち合わせなどにお使いください。

高村平成１９年度共同利用研究発表（ポスター）時間

共同利用研究者（所属）

Proj
番号

番号

発表タイトル対応教員

白 迎玖（東北公益文科大学公益学部）

一般研究 P-1

衛星画像を用いた台北の都市発展と都市ヒートアイランドとの関係分析近藤岡田信行（ALT都市環境研究所）

P2007-4P-2（仮）熱環境緩和対策の立案におけるリモートセンシング
データの活用方策の検討

近藤

仲岡雅裕（千葉大学大学院理学研究科）

P2007-4

P-3

ＲＳ・ＧＩＳによる沿岸生態系主要景観の広域分布動態の解析

近藤Prima O.D.A.（公立大学法人岩手県立大学 ソフトウエア
情報学部）

P2007-1

P-4

西部太平洋熱帯域と東アジアにおける気候環境の関連性に関する調査近藤

原口 昭（北九州市立大学国際環境工学部）

一般研究 P-5泥炭分布地域における生物地球化学的環境解析へのリ
モートセンシング技術の利用可能性の検討（予報）西尾Study of Coral Reef Distribution in Bali and Sumbawa
Island Using Satellite Data西尾

大澤高浩（ウダヤナ大学海洋リモートセンシング研究セン
P2007-1
ター）

P-6

木村 宏（岐阜大学工学部電気電子工学科）

干渉合成開口レーダ(InSAR)による南極氷河・氷床マッピ
一般研究 P-7 ング ―等速変動場の仮定を用いない流動位相の導出―

長 康平（株式会社ウェザーニューズ 道路気象コンテンツサービスグループ）

一般研究 P-8 合成開口レーダによる道路凍結のモニタリング

ヨサファット大山克己（千葉大学環境健康都市園芸フィールドセン
ター）

一般研究 P-9 Monitoring Kashiwa City Using Remote Sensing

ヨサファット朴 鍾杰（東京情報大学 環境情報学科）

P2007-1 P-10 ミクスチャモデルによるMODISデータの雲量推定

建石

S.K. Dash・田中 佐（山口大学理工学部研究科）

P2004-1 P-11

14:00-15:00 上野健一（筑波大学生命環境科学研究科）

Retrieval of drop size distribution using PR/TRMM andCPR/Cloudsat data（発表者：Sisir Kumar Dash)

P2007-4 P-12 筑波山頂における自動気象観測の復活

西尾

建石

樋口

松岡真如（高知大学農学部農学科国際支援学コース）

P2007-2 P-13

「NOAA/AVHRR GACを用いた全球時系列データセット
の 作成」

本多・梶原

松山 洋（首都大学東京 都市環境科学研究科）

カラマツ林における方向別反射特性の季節変化と日変
P2007-2 P-14 化〜光合成モデルに必要な樹冠構造パラメータ推定にむけて〜(発表者：長谷川宏一)

本多・梶原菅原広史（防衛大学校 地球海洋学科）

一般研究 P-15 都市域におけるアルベードの計測

高村佐々木真人（東京大学宇宙線研究所）

P2007-3 P-16 Ashra観測の現状報告

久世内藤季和（千葉県環境研究センター）

P2007-3 P-17

可搬型自動ライダー（PAL)によるエアロゾルの時空間計
測久世吉田弘樹（岐阜大学工学部）

一般研究 P-18

ライダーによる粒子径分布計測のための多重散乱モデル

久世

丹羽勝久（株式会社ズコーシャ）

P2007-4 P-19

ITを活用した低コスト生産のための効率的てん菜栽培手
法の確立

本郷

一ノ瀬俊明（(独)国立環境研究所社会環境システム研究
領域）

P2007-4 P-20

プロジェクタ光源を利用した長光路差分吸収分光（DOAS）法による大気汚染計測

大内和夫（防衛大学校 電気情報学群情報工学科）

原 政直（(株)ビジョンテック）

松島 大（千葉工業大学工学部建築都市環境学科）

O・・・ORAL
P・・・PＯＳＴＥＲ（＊）ポスター用パネルサイズは（縦）120cm×（横）90ｃｍです。
（＊＊）パネルを会議室4に準備いたします。適宜貼付願います。

（発表者：長南 勉）

ポスタータイトル：高分解能合成開口レーダのクロス偏波
P2007-1 P-21 画像 に見られるテクスチャを利用した針葉樹林バイオマ
ス計測
夜間光をパラメータとする二酸化炭素排出量の推定
P2007-1 P-22 −東アジア地域における二酸化炭素排出量の変動について―一般研究 P-23

高空間分解能衛星データを用いた森林立ち枯れ領域判
定法の開発

近藤・久世

西尾

西尾

樋口

第 10 回 CEReS 環境リモートセンシングシンポジウム
資料集2008 年 5 月発行
編集 千葉大学環境リモートセンシング研究センター
共同利用研究推進委員会〒263-8522 千葉市稲毛区弥生町 1-33
千葉大学 環境リモートセンシング研究センター
電話 043-290-3832 FAX 043-290-3857
URL http://www.cr.chiba-u.jp/