Galileo HAS（high accuracy service）その1
はじめに
欧州の測位衛星Galileoは、みちびきと同様に、測位信号とともにその測位精度を高めるメッセージ（補強メッセージ）も放送しています。そのメッセージ仕様は、Galileo High Accuracy Service Signal-in-Space Interface Control Document (HAS SIS ICD) にて公開されています。
幸運にも私はGalileo E6B信号を受信できる受信機を持っていますので、この補強メッセージの解読にチャレンジしています。まずは仕様書に書かれている例を実行してみました。
Galileo HASのメッセージ構造
GalileoのE6B信号上にあるHAS（ハス、high accuracy service）は、1秒間に1メッセージを伝送し、複数メッセージを組み合わせて補強を実現します。これは、みちびきのL6D信号にて伝送されるCLAS（シーラス、centimeter-level augmentation service）メッセージや、L6E信号にて伝送されるMADOCA-PPP（マドカピーピーピー、multi-GNSS advanced orbit and clock augmentation - precise point positioning）のメッセージ伝送と同様です。
また、Galileo HASのメッセージフォーマットは、みちびきのCLASやMADOCA-PPPで用いられるCSSR（compact space state representation）を考慮しているとされています。
7.8 Relationship with other formats
The Galileo HAS SIS ICD is a self-standing document defined taking into account the formats of already existing messages for providing high accuracy corrections. In particular, it takes into account QZSS Interface Specification of the Centimeter-Level Augmentation Service (IS- QZSS-L6-001) [RD4], which uses the Specification of Compact State-Space Representation (CSSR) for Satellite-Based Augmentation Messages defined under RTCM-SC-104.
しかしながら、私は、これらのメッセージ構造はかなり異なっているように思います。Galileo HASは、特に、メッセージの一部分が消失しても自己修復できることを意識しています。移動受信において起こるであろうメッセージの部分消失への対策は、Galileo I/NAV（アイナブ、integrity navigation message）航法メッセージにもあります。また、含まれる補強情報の種別は、みちびきのCSSRではメッセージを読み進めてはじめてわかるのに対して、GalileoのCSSRではヘッダにビットイメージ形式にて固定的に記されています。
Galileo E6B信号は、1秒間に984シンボル伝送されます。受信機は、これに123行8列のデインターリーブ（連続誤りを離散誤りに変換する）を行い、符号化速度1/2、拘束長7のビタビ復号を行います（Sect. 2.2 FEC Coding and Interleaving Parameters）。このようにして得た、ビタビ復号後の492ビットメッセージのうちの「HASページ」のメッセージ長は448ビットです。このHASページは、24ビット長のヘッダーと、424ビット長（=53オクテット）の「HAS Encoded Page」からなります。ヘッダーには、復号後の全ページ数k（the number of non-encoded pages, 1から32）や、復号前ページ番号N（1から255）が記されます（Sect. 3 HAS Page Header）。HASメッセージの復号は、このヘッダー情報をもとにして、HAS Encoded Pageを解釈してゆくことで行います。HASは、kページに冗長ページを加えてNページを伝送することで、ページの部分消失を自動修復します。
復号後ページ数は32以下に制限されます。Annex Cの例では15ページです。上述の通り1ページは424ビットからなり、余りメッセージ部分にはゼロと交番符号（01の交互）で埋めます（Sect. 4. HAS Message Overview）。しかし、このあまり部分は、ゼロ埋めされるのか、交番埋めされるのかは、わかりません。
padding bits with a binary sequence of zeroes and ones “010101…” are appended to the message in order to reach a message length which is a multiple of 424 bits.
また、図2のN=255の表記について、復号前の全ページ数Nは常に255なのか、それより少ないこともあるのかについては、私にはわかりませんでした。Nが常に255ならば、1周期のHASメッセージ受信に255秒かかることになり、メッセージ部分消失対策の代償が大きすぎるように思います（CLASやMADOCA-PPPにおいて、1周期の伝送時間は30秒です）。これは、実際のHASメッセージで確認してみようと思います。
メッセージ部分消失の自動修復は、リードソロモン（Reed-Solomon）符号化により実現しています。受信機は、受信したHASページ列に生成行列の逆行列を掛けることで、メッセージ復号を行います。この逆行列計算と積計算は、ガロア体上で実施します（Sect. 6.4 HAS Reception and Decoding）。k＝2、k=3、…、k=16に対応する生成行列の逆行列をあらかじめ計算しておけば、受信語とこれらのガロア体での積計算のみでHASメッセージを復号できます。
Note that only one matrix inversion is required for all vertical words and that Galois Field arithmetic needs to be taken into account for all operations.
リードソロモン符号はブロック符号なので、この誤り訂正を行わないオプションもあります（Sect. 6. HAS Message Encoding and Decoding）。しかし、Annex C. Reed-Solomon Decoding Exampleにある例を見ると、このリードソロモン復号前後で同じメッセージ列が見当たらないことから、このリードソロモン復号を行わなければならないようです。
When any subset C’ of C, of k different pages (C’1, …, C’k), is received without error, the message can be decoded.
Reed-Solomon Decoding Exampleの実行
そこで、Annex Cに書かれているリードソロモン復号例をMATLABを用いて試してみます。ガロア体での逆行列を求めるため、MATLAB本体の他に、Communications Toolboxが必要です。
この例では、符号化前ページ数kは15です。そこで、はじめにk＝15に対応する15行15列の生成行列Dを定義します。ここでは10進数にてこの要素を与えています。
D = [
31 50 155 253 213 220 84 174 239 85 87 105 214 81 160;
113 18 35 135 205 43 156 23 127 169 162 160 15 49 202;
204 239 127 208 89 187 30 192 37 152 221 214 211 49 93;
72 7 24 67 1 245 154 234 84 179 37 96 222 33 64;
253 151 182 118 101 136 118 241 195 26 152 14 225 28 193;
114 171 242 238 47 124 59 125 65 23 39 150 161 226 5;
33 32 3 8 36 151 121 17 218 26 98 82 65 146 162;
37 190 149 41 64 68 119 19 153 51 235 147 203 136 225;
58 217 47 14 1 13 117 8 167 10 105 226 96 158 229;
169 119 204 119 80 22 46 55 120 70 39 68 156 140 150;
145 19 150 65 190 97 199 178 76 115 138 198 136 18 180;
235 120 75 39 150 196 72 209 145 27 180 77 11 2 154;
143 165 24 101 222 187 133 80 114 98 164 11 16 227 43;
15 105 201 161 101 197 235 191 127 28 238 232 231 198 234;
84 157 205 255 217 251 101 194 230 208 26 232 23 201 46
];
これをgf()にてガロア体に変換し、その逆行列GFDinvを求めます。
GFD = gf(D,8);
GFDinv = inv(GFD)
GFDinv = GF(2^8) array. Primitive polynomial = D^8+D^4+D^3+D^2+1 (285 decimal)
Array elements =
200 5 121 64 146 79 105 195 15 22 150 75 227 167 168
19 123 74 36 87 170 123 139 57 66 117 19 60 150 3
71 178 51 76 147 41 146 222 254 149 82 30 109 52 69
218 100 198 13 145 216 164 119 203 199 219 224 100 198 41
220 215 176 139 91 173 91 109 135 26 12 188 209 164 238
137 243 144 94 235 114 136 177 122 41 132 110 53 186 53
251 142 229 87 112 189 132 245 226 176 21 171 39 199 231
103 219 81 164 107 111 129 192 209 102 161 141 119 23 230
204 204 167 128 199 2 33 60 17 232 124 25 79 159 233
203 113 26 162 49 31 24 40 59 114 219 100 211 120 129
98 15 14 203 1 75 205 94 78 35 41 23 178 107 76
41 152 227 188 201 126 137 5 126 121 55 226 21 41 45
69 107 172 85 239 175 198 249 164 130 179 235 162 58 196
175 7 101 251 240 158 218 225 58 249 193 64 218 17 214
174 231 25 182 215 187 49 97 157 231 127 42 126 221 196
次に15行からなる受信語Wを定義します。列数は常に53です。
W = [
132 123 199 73 235 125 113 116 36 71 136 251 69 70 145 140 0 39 42 235 193 84 146 204 110 181 90 88 128 226 97 186 227 23 26 35 221 11 229 98 252 141 111 216 142 98 41 194 158 125 140 153 223;
52 154 227 99 77 33 11 173 50 147 166 127 182 33 1 233 221 84 48 123 198 121 237 105 155 213 12 174 174 197 100 133 243 248 22 84 12 174 206 164 198 22 146 238 91 24 202 171 181 189 162 121 57;
85 1 29 145 14 230 225 85 194 242 140 77 215 250 214 40 200 226 106 5 171 215 135 151 77 226 225 111 142 246 176 156 0 215 18 228 41 8 34 151 24 174 236 105 28 5 39 243 194 63 128 181 19;
44 163 27 35 21 83 238 106 156 122 59 255 250 132 43 45 12 243 8 9 16 185 194 2 126 136 115 220 237 47 141 167 212 35 164 47 217 206 88 195 238 68 125 44 175 49 177 138 4 213 165 186 120;
55 190 96 216 35 121 141 182 26 28 152 34 238 248 75 122 213 237 99 213 34 61 152 173 145 204 133 143 64 117 119 92 224 76 187 36 160 208 177 95 127 213 58 214 134 44 121 248 82 63 169 191 75;
187 28 69 29 89 4 160 228 22 185 43 88 154 12 86 206 43 199 115 152 40 239 11 192 73 228 145 24 154 41 63 49 40 36 224 176 100 94 31 100 152 109 111 135 185 118 207 58 18 247 59 144 33;
117 25 72 154 251 194 111 69 202 191 253 159 120 178 246 68 171 41 251 163 124 202 254 239 152 25 2 5 204 223 192 231 250 120 193 179 234 80 108 166 166 167 210 195 99 135 159 118 132 143 164 128 36;
143 12 156 52 139 203 193 61 89 3 53 84 14 168 101 194 207 61 113 59 188 39 200 99 26 41 88 222 211 134 178 117 71 15 136 150 150 65 88 124 204 128 23 28 51 166 204 221 251 63 53 44 190;
203 226 36 10 145 27 54 129 243 142 43 63 242 57 243 98 229 59 74 201 41 44 96 199 124 97 197 70 118 78 134 66 106 138 68 197 64 140 187 91 201 10 138 135 16 254 109 113 144 220 128 204 93;
29 55 158 167 195 223 144 158 158 116 87 219 101 36 71 28 189 52 215 17 199 92 176 139 74 132 108 3 25 126 46 191 226 239 14 161 44 70 247 253 202 246 58 36 35 29 77 144 52 14 217 139 221;
122 57 40 21 48 65 99 21 77 50 204 30 233 166 117 3 48 3 115 250 224 78 143 108 245 144 255 199 147 114 161 38 145 41 107 172 132 82 95 202 166 152 75 83 88 143 25 25 186 202 151 159 222;
125 19 56 207 112 92 184 147 239 181 113 209 24 245 173 57 173 51 3 160 148 255 182 92 140 168 146 194 234 61 53 190 137 15 91 228 231 9 111 222 52 62 205 189 90 185 129 222 74 19 154 94 29;
161 204 117 222 253 61 201 66 207 106 21 166 117 149 224 164 249 50 45 172 71 205 29 87 112 81 177 95 215 130 214 162 83 43 182 9 188 112 183 111 5 174 231 176 103 151 117 7 232 167 19 33 234;
207 147 205 21 140 244 31 178 149 173 157 33 161 85 130 130 237 116 136 51 54 137 106 123 126 234 208 57 145 34 116 229 209 226 26 86 63 239 245 210 21 211 61 189 43 85 215 103 160 170 234 163 56;
215 200 167 19 210 166 18 96 224 77 5 145 106 148 222 103 157 196 233 132 109 61 229 187 163 152 17 62 27 210 42 67 181 2 23 108 68 206 189 76 58 39 164 43 254 9 87 41 18 228 135 212 165
];
同様にこのWもガロア体に変換して、Dに掛け、その結果を16進表記します。
GFW = gf(W,8);
GFM = GFDinv * GFW;
compose("%02x", GFM.x)
ans =
15×53 string
column 1 to 14
"00" "0c" "c0" "0b" "20" "ff" "df" "ff" "ff" "00" "81" "00" "f7" "ff"
"ff" "68" "22" "fe" "a2" "18" "07" "c1" "93" "f7" "59" "80" "35" "fd"
"47" "f9" "03" "ff" "9d" "f7" "80" "5c" "15" "ff" "9f" "dc" "ff" "80"
"43" "fd" "b0" "23" "04" "00" "7f" "e5" "03" "0f" "f1" "ac" "40" "02"
"4f" "e1" "ff" "88" "25" "fe" "8f" "fc" "ff" "00" "48" "08" "1f" "e3"
"f0" "0f" "ff" "70" "4b" "f7" "1f" "ff" "fd" "c0" "97" "fb" "40" "0c"
"40" "56" "05" "c0" "88" "04" "00" "44" "03" "01" "2f" "e2" "7f" "ef"
...
これは、Annex C.にある復号後メッセージと一致しました。
Pythonでの実装
Pythonでこの手順をどのように実装しようかと悩んでいたところ、Dr. Eric SibertからCSSRlibにHAS Decoderが実装されていることを教えていただきました。
CSSRlibでのHAS復号の実装では、Pythonのgaloisモジュールを用いてガロア体での積や逆行列を計算しています。具体的なコードがあると、とても勉強になります。私も頑張って作ってみます。
また、k=15の生成行列は与えられていますが、それ以外の生成行列は仕様書HAS SIS ICDに書かれていませんので、自分で考えなければならないようです。大学生のときに学んだような気のするリードソロモン符号を学び直しているところです。
（たぶん）その2に続きます。
関連記事
- Galileo HAS（high accuracy service）その3 27th October 2023
- Galileo HAS（high accuracy service）その2 20th February 2023