2011/09/01〜

...................................................................................................................................

2011/08/27

RIA Novosti, Glonass-M satellite launch postponed for additional check, August 25, 2011

8/25にプレセツク宇宙センタからSoyuz-2.1Bロケットで打ち上げ予定だったGLONASS-M衛星の打ち上げ延期。前日のProress M-12M打ち上げ失敗 (参照, 参照) を受けて追加の機器点検を行うためとのこと。再打ち上げ予定は9月上旬。
Progress M-12は8/24 13:09 UTCにバイコヌール宇宙基地からSoyuz-Uロケットで打ち上げられたが、上段ステージ燃焼中に異常が発生してエンジン停止、そのまま南シベリアに落下。ロシアは8/18にプロトン-Mロケットによる通信衛星の軌道投入にも失敗 (参照) しておりトラブル続きである。10月のGalileo IOV衛星2機の打ち上げはSoyuz-STロケットを使って行われる予定だがこれには影響しないとのこと。

...................................................................................................................................

2011/08/24

M.Ge et al., Improving carrier-phase ambiguity resolution in global GPS network solutions, Journal of Geodesy, 2005

昨日書いたNAPEOSのamibuguity resolution戦略は、ほぼこのGFZ EPOSのやり方と同じ。まあほぼ確立された手法だということ。計算量という面ではまだ改良の余地があると思うけど、調べていくとこの分野もうあんまりやること残っていないのかも。ところでMGSで線形独立性をチェックするって具体的にどうやったらいいのだろう。

...................................................................................................................................

2011/08/23

C.Rocken et al., Apex - A high accuracy service for global real-time precise point positioning with GPS and GLONASS, ION GNSS 2011 abstructs

今年のIONでVeripos APEXについて発表があるのを見つけた。軌道は過去48Hのデータを基にNAPEOSで1H毎に推定し、その軌道予報を基にしてRTnetでクロックをリアルタイム推定しているとのこと。GLONASS暦は評価中のようだが地上局がまだ35局しかない様で品質面で少し厳しいのではないかと思う。しかし今やろうとしていることにとって、とても参考になる。

-------------------------------------

下で「epoch dependentパラメータとして解く」というのが多分わからないと思うので整理もかねてメモ。

大規模な最小二乗問題ではしばしば以下の様な形の正規方程式が現れる。ここでcは観測エポック間共通パラメータ、1, 2, ..., mが各観測エポックのみに依存する (epoch dependent) パラメータである。共通パラメータは例えばambiguityや対流圏パラメータ、epoch dependentパラメータは時計やキネマティック解の座標である。

N (xc',x1',x2',...xm')'=(bc',b1',b2',..,bm')'

以上より、

Nc * xc + sum (Nci' * xi) = bc
Nci * xc + Ni * xi = bi (i=1,2,...,m)

いったんxi を消去して整理すると

xc = Nc*^-1 bc*

ここで

Nc*=Nc - sum (Nci' * Ni^-1 * Nci), bc* = bc - sum (Nci' * Ni^-1 * bi)

xcが求まれば各xiは、

xi = Ni^-1 (bi - Nci * xc) (i=1,2,...,m)

Ni, Nci, biはエポック毎に計算しHDかSSDにはきだして後で使うとすれば、メモリ常駐が必要なのはNc*, bc*のみとなる。xcの解だけが必要ならHDにはきだす必要もない (例えばiterationの1回目)。

単純な基線解析の場合二重差を作って時計を観測方程式から消去するのは普通に行われるが、ネットワーク解析におけるこの方法の問題は事前にどの基線を使うかを決定しなければいけない点である。いったんゼロ差でフロート解を解いて、ambiguity決定に有利な基線を動的に検索し、ここで初めて二重差を作る方がambiguity決定性能が上がる。この方法はESOCの基線解析ソフトNAPEOS (参照) で採用されている。

NAPEOSのambiguity決定では動的に二重差amibiguityを作って信頼度順にソートし、線形独立なものだけ上から順に採用していくという手法をとっている。ただ、MW WL→NLをsquential roundingで解いているので整数最小二乗という訳ではない。ネットワーク解の全ambiguguityを整数最小二乗で解くのは計算量的に多分無理なのだが、基線毎NLだけでも整数最小二乗にしたほうがambiguity決定確率が上がるはずである。

ということで大規模ネットワーク解をどうやって解くかというのを現在検討している訳だが、メモでまとめてみると少し頭が整理されてきた気がする。

補足: 下で大規模最小二乗ではHがメモリに乗らないからQR分解は駄目と書いてしまったが、out-of-core QR factorizationという手法も開発されている様だ (参照)。外部記憶として最近のSSDを使えば、RAIDにしなくてもシーケンシャルなら200MB/sは出るから100GBのI/Oで1000sで実用性はある。ただSSDでもランダムアクセスは極端に性能が落ちるので、アルゴリズムを十分最適化する必要がある。QR分解は計算量が増える代わりに数値誤差の影響を受けにくいのでオプションで切り替えられる様にしておくのがベストだろう。資料中に「2年前まではn=10000, m=100000は大規模問題と考えられていて、クラスタが必要だった」とあるが、確かに今なら64bit OSが普通でメモリも安くなったので、これくらいの最小二乗は普通のPCでもどうってことない。ただ、n=100000になるとまだ普通のPCではかなり厳しくて実現にはいろいろな工夫が必要になる。(13:13追記)

-------------------------------------

大きな最小二乗問題を解くという必要性があって、どれくらいの大きさなら普通の計算機で実用的な速度で解けるかを検討してみた。実行時間計測用のテストルーチン (matlab) とその実行結果。

function lsqtest
% large lsq performance test

for n=[1000 2000 5000 10000 20000 50000]
    for m=n*[10,20,50,100]
        p=n/1000; % sparseness
        H=rand(m/p,2*p); y=rand(m/p,1);
        tic
        N=zeros(n,n); b=zeros(n,1);
        for i=1:n/p-1
            j=[1:p,i*p+(1:p)];
            N(j,j)=N(j,j)+H'*H; b(j)=b(j)+H'*y;
        end
        t1=toc; tic
        R=chol(N);
        x=R\(R'\b);
        fprintf('n=%6d m=%7d: time=%7.3f %7.3fs\n',n,m,t1,toc);
    end
end
>> lsqtest
n=  1000 m=  10000: time=  0.067   0.019s
n=  1000 m=  20000: time=  0.067   0.012s
n=  1000 m=  50000: time=  0.150   0.012s
n=  1000 m= 100000: time=  0.288   0.013s
n=  2000 m=  20000: time=  0.049   0.075s
n=  2000 m=  40000: time=  0.067   0.077s
n=  2000 m= 100000: time=  0.129   0.073s
n=  2000 m= 200000: time=  0.251   0.073s
n=  5000 m=  50000: time=  0.263   0.970s
n=  5000 m= 100000: time=  0.467   0.989s
n=  5000 m= 250000: time=  1.078   1.002s
n=  5000 m= 500000: time=  2.315   1.006s
n= 10000 m= 100000: time=  0.760   7.466s
n= 10000 m= 200000: time=  1.394   7.512s
n= 10000 m= 500000: time=  3.346   7.546s
n= 10000 m=1000000: time=  7.160   7.562s
n= 20000 m= 200000: time=  1.780  58.118s
n= 20000 m= 400000: time= 33.312  59.012s

(CPU=i7 2600, RAM=8GB, AVX=off, matlab 2011a, win 7 64bit SP1)

以上見てわかるように正規方程式＋コレスキー分解である。最小二乗にはMGS等を使ったQR分解もよく使われるが、計算量が約2倍になるし、観測数が多くなるとHがメモリに乗らなくなる。正規方程式ではメモリはn x (n/2+2)でよいので、大規模問題では正規方程式にせざるを得ない。なおsparse性を利用すると正規行列Nの更新に必要な計算量を大幅に減らせる。上記では計画行列のsparse性が1/1000と仮定した。

ここで解きたい最小二乗問題は、たとえばGEONET全局 (1200局) の24Hスタティック解である。解析戦略としては、

(1) 電離層フリー搬送波位相でゼロ差観測方程式を立てる。
(2) (1) から正規方程式を立てていったんフロート解を求める。
(3) MW WL→NLの標準的な手法でambiguityを解く。
(4) (3) を拘束条件として導入し再度正規方程式を解く。
(5) 衛星・受信機時計はepoch dependentパラメータとしてメインの正規方程式外で解く。

未知パラメータ数: n=1200 (局) x (3 (座標) + 3 (対流圏) x 12 (2H毎) + 1 (ambiguty) x 32 (衛星) x 2 (パス)) = 123600
観測数: m= 1200 (局) x 10 (衛星) x 2880 (エポック) = 34560000

最小二乗計算時間は、(1.8s*m/200000*n/20000+59s*(n/20000)^3)*2 (回) =31696s = 7314236s=8.80H。
ここで正規行列計算はm, nに対し線形、正規方程式求解はnに対して3乗の割合で計算量が増えると仮定している。
また、所要メモリ量は n*(n/2+2)*8 = 61.1GB。

例えば、今、Xeon (6core) x 2 way + RAM 48GBのマシンが100万位 (参照)。ただこれは4GBモジュールで8GBモジュールはとたんに高くなって7万/枚位だから12枚96GBで合計200万位。来年だったらAVX対応Xeonは出ていて速度は6 core × 2 way x AVXで上記の約5倍。ということで頑張れば2-300万位の計算機を使って実用的な時間で計算できそうである。でも、プログラムの開発は大変そうで、まず戦略や最適化を理解してプログラムを書ける人間がいないというのが一番問題になるかも。やっぱり某プログラムもカルマンやめてバッチにしようかなあ。

...................................................................................................................................

2011/08/21

宇宙開発戦略専門調査会, 宇宙開発利用の戦略的推進のための施策の重点化及び効率化の方針について, 平成23年8月8日

先日ニュースになった宇宙開発戦略本部 専門調査会の報告がアップされている。

> 2. 最重要課題としての準天頂衛星システム

となっていて、QZSSばかり強調しているの様に読めるのだけど、何でこういう結論になったのかなあ。まあ、いずれにしても来年度以降の予算がつきそうなので、あと10年くらいはこれで食べて行けるかも。せっかく「再重要課題」となったのだから、ぜひ産学官うまく連携して有用な技術開発を進めていきたい。特にこの分野、日本において「学」の寄与が今まで相対的に少なかったので。

...................................................................................................................................

2011/08/20

RTCM 10403.1, Differential GNSS (Global Navigation Satellite Systems) Services - Version 3 + Amendments 1,2,3,4, and 5 to RTCM10403.1, July 1, 2011

RTCM 3.1 SSR規格 (RTCM 10403.1 - Amendment 5) の正式版が出ていたので購入。$80也。以下Amendment 5で修正・追加されたメッセージのリストを示す。

1015 GPS ionospheric correction differences (modified)
1016 GPS geometric correction differences (modified)
1017 GPS combined geometric and ionospheric correction differences (modified)
1019 GPS ephemerides (modified)
1031 GLONASS network RTK residual message (modified)
1057 SSR GPS orbit correction
1058 SSR GPS clock correction
1059 SSR GPS code bias
1060 SSR GPS combined orbit and clock corrections
1061 SSR GPS URA
1062 SSR GPS high rate clock correction
1063 SSR GLONASS orbit correction
1064 SSR GLONASS clock correction
1065 SSR GLONASS code bias
1066 SSR GLONASS combined orbit and clock correction
1067 SSR GLONASS URA
1068 SSR GLONASS hight rate clock correction

3.5.12.11 consistency of data and processing で

> Note: Please consider two corrections in the GLONASS ICD Version 5.1. In appendix A.3.1.2 the
> sign of the 4th term in dVy/dt correctly reads -2 sigma Vx instead +2 sigma Vx and dVz/dt the
> term in parentheses correctly reads (3 - 5 z²/r²) instead of (1 - 5 z²/r²)

とちゃんとGLONASS ICD v.5.1のバグに関する注記が追加されているのが有難い。(今までこのバグに関する公式文書なかったので)

ところで、GLONASS ICD v.5.1って以前ダウンロードできたサイトからダウンロードできなくなっていて、正式版がどこから入手できるか不明。現在のところ少し怪しげなここくらいしか入手先がないようだ。できればL3 CDMAやL1P、L2Pについてちゃんと記述された改定版ICDを公式サイトで入手できる様にしてほしいのだけど。

...................................................................................................................................

2011/08/19

RTKLIB: Support Information

ずいぶん遅くなっってしまったが溜まっていたパッチをまとめてリリースした。(rtklib_2.4.1_p1)

...................................................................................................................................

2011/08/17

Jcastニュース, 夢の30兆円ビジネス「日本版GPS」、本格スタートできるか, 2011/8/15

どうやったら「30兆円」という数字が積み上がるのか「予測」のエビデンスを知りたいのだけど。

...................................................................................................................................

2011/08/16

VMware上のUbuntu 10.04でAVXがonにならず行列計算速度が上がらない件、Ubuntuをネイティブ環境にインストールし直して速度を再計測してみた。軽く100 GFLOPオーバーである。AVX素晴らしい。なお、GotoBLAS2はオリジナルの開発は終了しておりAVXには対応していない。AVXに対応したバージョンも出回っている様でインストールしてみたが乗算で65GFLOPS位しか出ていない。

...................................................................................................................................

2011/08/14

テレビ東京, ワールドビジネスサテライト: 日本版GPSの可能性, 2011/8/2

先日テレビ東京WBSで放送されたQZSS特集がアップされている。良く知った顔が沢山出ている。殆どQZSSのプロモビデオで、これはこれで良いのではないかと思うのだけど、「GPSの測位精度は10m、QZSSは3cm」という (殆ど) 誇大宣伝だけはさすがにやめて欲しい。

...................................................................................................................................

2011/08/12

この日記・備考録はもともとは自分用のメモとしてもう6年くらい続けているのだけど、過去書いていることを読み直してみるとこの6年で対して進歩していないことが分ってへこむ。
たとえば、QZSSの軌道決定精度について書いている2005/11/18の記述を読むと、今考えても的確なこと書いてあって凄いなあと思う。実際たいした根拠があって書いてたはずではないけど。あと、長基線RTKの実装を始めた2005/12/19位の記述を読むと、最初に思いついたアイデアを今でも少しずつmodifyしながら実装しているに過ぎないことも分かる。以上を逆にいえば、良いアイデアを思いつくのはたいして難しい訳ではないけど、実際にそれを実装して、評価して、改良して、その技術を使えるようにして、普及させるためには、長い時間と工数と努力と忍耐といろんな失敗が必要であるということ。
以上まとまらないが、最初のアイデアが良いかというのも大事だけど、継続するというのはもっと大事だということかなあ。(まあ6年ぽっちで偉そうなことは言えないけど)

...................................................................................................................................

2011/08/08

日経, 日本版GPS整備「内閣府主体で」宇宙本部の専門調査会, 2011/08/08

> 政府の宇宙開発戦略本部（本部長・菅直人首相）の専門調査会は８日、日本版の全地球測位システム
>（ＧＰＳ）衛星「準天頂衛星」を、内閣府が予算要求して整備を進めるべきだとする最終報告書をまと
> めた。文部科学省や総務省が個別に要求する現行体制の見直しを提言した。近く戦略本部に報告書を提
> 出し、2012年度予算の概算要求に内閣府として関連費用を盛り込みたい考えだ。
> 最終報告では、国際宇宙ステーションの日本実験棟「きぼう」は産業につながる成果が明らかでないと
> し、16年以降は「経費節減を進めるべきだ」とも指摘した。準天頂衛星は内閣府を中心に開発・運用を
> 進める体制づくりも求めた。文科省や経済産業省など複数の省にまたがる宇宙行政をいわゆる「宇宙庁」
> として内閣府の下に一元化する構想をにじませた。ただ宇宙行政を統合して無駄を省くとの具体的な記
> 述は関係省庁との調整がつかず、改めて議論して最終報告書に盛り込む。準天頂衛星は既に１基が打ち
> 上がっている。２基目以降の整備を「可及的速やかに実施すべきだ」と提言した。

とのこと。
これ自身はめでたいことではあるが、もともとは「第1段階として1機の衛星による技術実証および利用実証を行い、その結果によってシステム実証を行う第2段階へ移行するか判断する」としていたはずので、まずは第１段階の技術実証・利用実証結果の整理が先のような気がする。

...................................................................................................................................

2011/08/06

L.Wanninger, Carrier-phase inter-frequency biases of GLONASS receivers, Journal of Geodesy, 2011

GLONASS受信機のinter-frequency bias (inter-channel bias) に関して基準局網のデータを使って傾向や時間安定性を調べている。なかなか見事な論文。細かくは後から読む。たぶんこの結果そのまま使えばGLONASSの軌道・時刻決定でアンビギュイティ解けるんじゃないかなあ。これを見ると、Ashtech (old)とTPSとJPS、LeicaとNovAtelの中身が (少なくともRFフロントエンドは) 同じ設計であることが予測できる。
3月まで海洋大にいて今ENRIにいる山田君のION2010論文が"Yamanda et al., 2010"として引用されている。

-------------------------------------

昨日書いたQZSSの搭載原子時計異常の件、GPSの原子時計異常についての論文に関して問い合わせがあって少し調べたのでメモ。

GPS搭載原子時計の特性評価については過去PTTI meetingで多数発表されている。

http://www.pttimeeting.org/archivemeetings/

そのうち異常ケース解析を行ったもの。

http://www.pttimeeting.org/archivemeetings/1999papers/paper4.pdf

ただ、「内蔵水晶のフリーランでクロックを出力し続ける」を裏付ける資料は見つからず、どうも私の勝手な思い込みだった可能性が高い。実際のGPS搭載クロック異常はもっと複雑なメカニズムで発生している様。

最近のGPS搭載クロック異常は2007/10/8 PRN19で発生しており、これについては2007/10/12の備考録でメモを書いている。ただこの際のクロックドリフトは今回のQZSSの事象よりかなり小さい (10^-10s/s位)。この本当の原因は不明だが地上系ソフトウェア更新が疑われている。(→GPS Worldの記事)

最近、ここJAXAの方も読んでおられる様で、以上、私の不正確な情報で今回のQZSSの障害究明を混乱させると困るので訂正しておく。

...................................................................................................................................

2011/08/05

GPS/GNSSシンポジウム2011, 2011/10/26-28, 東京海洋大学越中島会館

最初の案内が出たようなので貼っておく。聴講は無料なのでぜひ多数の方がご参加頂きますようよろしくお願いいたします。
秋の測地学会と日程がかぶるのだけど、私はこちらへ参加の予定 (多分なにか発表も)。

-------------------------------------

OSQZSS, 衛星クロックオフセット, 2011/8/5

先日の「みちびき」の搭載ルビジウム原子時計異常 (参照) 時の衛星クロックの挙動を解析している。これを見ると異常時以降、冗長系切り替えまでのクロックドリフトは、-3.5 x 10^-8s/s。GPS衛星搭載原子時計障害時の挙動は過去解析されているが、大体内蔵水晶のフリーランでクロックを出力し続けるので、概ねOCXO並みの安定度に落ちると考えれば今回の事象は説明がつく。「みちびき」に搭載されているルビジウム原子時計 (RAFS) は米Excelitas社製 (参照) で、GPS-IIR, IIFに搭載されているものと同じものの様だ。2台しか積んでいないのでもし復旧できないとすると今後結構厳しい運用を強いられることになる。ぜひ頑張って復活させてほしい。
ところで、今回の障害で地上受信機でQZSS信号を一時的にロストしたケースが多かった様だ。-3.5 x 10^-8s/sということは、L1周波数シフトは約55kHz。急速な周波数変動に受信機PLLでの搬送波追尾が間に合わなくなったということだと思う。でも通常の受信機では捕捉時ドップラー探索幅は±10kHz位が普通だから、クロック復旧前に再捕捉は不可能な気がするんだけど、少なくともOSQZSSでは再捕捉してるなあ。これは探索幅が広いということ?

補足: 計算間違えた。周波数シフトは約55Hzでした。これでは再捕捉は問題ない。でも逆に55Hz程度のシフトでPLLロックはずすかなあ。もしこの程度ではずすとするとちょっとした移動体ですぐはずれる気もするが。(13:47追記)

-------------------------------------

JAXAが主催する「複数GNSS対応軌道時刻推定アルゴリズム研究会」の案内。既に測位航法学会のメーリングリストにも流れた様なので、問題ないと思うので、ここにも転載してしまおう。

> 複数GNSS対応軌道時刻推定アルゴリズム研究会への参加依頼
>
> 本メールは、標記、研究会への参加依頼に関する案内メールです。
> 下記メールを読んでいただき、参加検討いただければと思います。
>
> みちびきも打上げから約１年運用を重ね、先日、放送暦の精度検証を完了し、
> アラートフラグを解除させていただきました。
> これからは精密暦の高精度化に取り組むフェーズに入りました。
>
> また、JAXAとしては、今後増加する他システムの精密測位利用のためにも、
> QZSS含む複数GNSS対応の軌道・時刻推定を、次世代衛星測位基盤技術の
> 要として取り組む予定です。
>
> 現在、上記取り組みとして、
> 複数GNSS対応の受信機60台のホスト機関、局の募集を開始、
> 複数GNSS対応の観測ネットワーク構築に着手するとともに、
> 軌道・時刻推定の新しい解析プログラムの開発にも着手しました。
> この解析ソフト開発は従来にない新しいプロセスとして、アルゴリズム検討を行う
> 有識者により構成する研究会を並行して走らせながら行う予定です。
> 本研究会は、2011年8月から月１回程度のペースで、2013年3月末までの開催を予
> 定しております。
> また、会合の実施方法として、遠方からの参加も可能となるようスカイプの多点
> 接続を活用するなど進めたいと考えております。
> 上記研究会を実施していくにあたり、皆様方の知見を是非、活用させていただき
> たいと思い、連絡させていただきました。
> 是非、ご検討いただき参加いただければと思います。
>
> なお、本研究会への参加を希望される方は、参加規約(添付1)に同意いただき、
> 参加申込書(添付2)に
> 必要事項を記入の上、下記、メールアドレス、または、FAXにて申込書を「８月
> １２日」までに送付願います。
>
> 参加申込書　送り先
>
> メールアドレス ： QZSS-AP(a)jaxa.jp
> FAX番号　　　 ： 029-868-5987
> 研究会事務局 　： 小暮　聡、三吉　基之、河手　香織

参加申込書や参加規約が添付ファイルで送られてきているけど、転載できないのでご興味がある方は事務局の方に直接ご請求ください。(メールアドレス: (a)→@)

-------------------------------------

NICT, Ionospheric disturbances after the 2011 Tohoku earthquake

NICT津川さんらによる東北大地震後のTEC変動アニメーション。データとしてはGEONETのGPS観測データを使っているとのこと。

大きな地震後に地表の振動 (あるいは地震による津波の振動) が大気に伝播し、最終的に電離圏の電子密度分布に擾乱を与える現象が、稠密な地上GPS網のデータ解析により観測できることが既にいくつかの研究として発表されている。上記アニメーションでも、東北大地震による大気の振動が電離圏を伝わる衝撃波 (?) として綺麗に検出できている。どういうメカニズムで振動が伝播するかは幾つかの仮説が出されているが、定説がある訳でもないようだ。なお、先日津川さんにお聞きしたところによると、地震以前のデータの解析では電離圏擾乱は検出できていないとのこと。今年春の連合大会で北大の日置先生が地震前の電離層擾乱を検出したと発表されている (参照) が、解析手法にどういう違いがあるのか興味のあるところである。

...................................................................................................................................

2011/08/02

S.Jazari et al., Fast integer least-squares estimation for GNSS high-dimensional ambiguity resolution using latice theory, Journal of Geodesy, 2011

整数最小二乗 (ILS) に関するいくつかのアルゴリズムを実装してその実行時間を比較している。比較している対象はLAMBDA, MLAMBDA, AEVZ (Agrell, Erikssson, Vardy, Zeger), M-VB (modification of Viterbo-Boutros), M-SE (modification of Schnorr-Euchner)。結論として最速はAEVZ、次にM-SEで、条件にもよるがこれらはLAMBDAの100倍以上の速度だとしている。なおILSは通常reduction, searchの二つのステップで実行されるがこの論文の評価の対象はsearchのみのである。

Matlabだし、LAMBDAとMLAMBDAはX,W-Changの実装を使っている様で、最適化されたコードとは言えない (AEVZでもn=50のケースで0.1秒もかかっている)。その点、必ずしも実用的な比較という訳ではない。またreduction stepも結構時間がかかるのでsearchだけ速くてもしかたないという面もある。ところでreductionはやっぱりLLLしかないのかなあ。革新的なアルゴリズムでネットワーク解の1000基線、n=30000位の整数最小二乗が1時間位で解けるとうれしいのだけど。でももしそれが可能なら、CVPの計算困難性を利用している暗号が全滅だから、たぶん無理だろうなあ。

-------------------------------------

いつのまにか8月だあ。

...................................................................................................................................

〜2011/07/31