接觸慣性測量原件 IMU 已有將近 10 年，鑒於感測器在生活中越來越常見，與自己工作又息息相關，透過本篇記錄自己對IMU的了解，並介紹IMU的種類、如何閱讀規格書與簡單資料讀取程式的撰寫。
IMU 的規格種類
IMU 依據精度規格可以分成 Strategic (戰略)、 Navigation (導航)、Tactical (戰術) 與微機電 (Micro Electro Mechanical Systems, MEMS) 四種等級，依照下表顧名思義，左側多為國防或太空產業用途的高階IMU，最右則是民用的 MEMS IMU。
不過下表年代稍久，目前IMU精度不斷提升，售價不斷下降，實際情況已經跟下表有所不同。目前有看到 STMicroelectronics 的 MEMS IMU ASM330LHH 陀螺飄移可以達到 3 度/小時。
IMU是陀螺 (Gyro) 與加速度計 (Accelerometer) 組成，分別提供角速度與加速度觀測量。加速度計一般都是 MEMS，精度規格差異不大，而陀螺相對是比較重要的觀測量 (可以想像偏差 1 度，隨移動距離而遞增的軌跡偏移量有多麼可觀)，依運作原理可分成 MEMS、光纖陀螺和雷射陀螺三種，後兩者因技術原理，體積較大。下圖是三種陀螺的精度與應用範圍，可以對應到上面IMU的等級。IMU的成本和精度，很大程度是看使用的陀螺規格。
一般手機裡面搭載的是 MEMS IMU最便宜的等級，只要幾塊美金，精度也最差。這類低成本的 MEMS IMU也多見於掃地機器人與 AR / VR 設備，如下方左圖中的黑色方塊，約 1 個指節的大小。
下方右圖左下角則是德國 iMAR 的戰略等級IMU，是國立成功大學江凱偉教授研究室做研究使用，一般用作參考系統，提供 ground truth 軌跡。因為採用雷射陀螺，使得體積接近一台小主機。下方右圖右側兩顆黑色與一顆銀色則是江凱偉教授研究室推出的戰術與戰略級 GNSS/INS 組合導航系統。
https://store.siqma.com/bmi160-6-axis-imu-sensor-module.html
IMU 的誤差與率定
至於IMU的精度規格要怎麼看呢？IMU的誤差來源分為系統誤差與隨機誤差，前者包含：
- Bias 偏差：量測值與真值的常數偏移
- Scale factor 尺度因子：量測值與真值相較的縮放比例，真值越大，尺度因子造成誤差越大
- Non-orthogonality (misalignment) 非正交：因為陀螺與加速度計一般是 x, y, z 三軸組成，如果三軸間彼此不正交，真值會在軸間產生投影造成誤差
以圖示說明系統誤差的表現行為如下：
MEMS 陀螺的偏差補償可以在靜置條件下移除非零值，但加速度計需要確定靜置與水平條件(例如使用水準氣泡)，才能將水平軸非零值視為偏差，並移除垂直軸非重力加速度的部分。因此實務應用上多只能處理陀螺的偏差，因為IMU通常不是水平狀態。
尺度因子與非正交誤差則需要準確的參考值或真值，才有辦法根據公式計算，故IMU一般是在實驗室透過多位置率定法 (multi-positions method) 配合專門的精密轉台，提供三軸角速度精確真值，以及線性加速度台提供精確的加速度真值，再採集動態數據計算各項系統誤差。若進一步搭配溫箱考量系統誤差的溫度效應，轉台價格可達數百萬台幣。如下面這台三軸轉台：
下圖則是學生時期做的陀螺偏差隨溫度變化的率定，紅藍代表兩組不同的數據集，橫軸是溫度，縱軸是偏差，可以看出陀螺偏差在不同溫度下不是固定常數，而且有一定趨勢。
隨機誤差主要包含：
- Random Walk：高斯常態分佈的雜訊，相關時間 (correlation time) 遠小於採樣時間，行為上與白噪聲相同
- Bias instability：一種低頻雜訊，因其表現行為而被視為系統誤差 bias 的波動，是應用上重要的指標參數，也就是上面規格表的陀螺飄移
其他像是 Quantization noise、Rate ramp 等較不常見於產品規格書，就不多做討論。
具時間序列的資料其隨機誤差一般用阿倫方差 (Allan Variance) 分析，概念是將資料依特定時間長度分段，各段分別取平均，再求相鄰段平均值的差值序列，最後計算差值序列的方差。公式如下：
當用較短的時間長度分段時，看的其實是資料的高頻變化，用較長的時間長度分段時，看得就是資料的低頻變化。例如把總長 1000 秒的資料，用100 秒的時間長度切分成 10 段，各段平均計算差值序列，再計算方差，其實就是以 100 秒為單位，看這 10 段資料彼此間隨時間的變異性當作 Noise，也就是每 100 秒 1 次的變化。相對於以 1 秒切割，計算每 1 秒 1 次的變化來說，前者屬於低頻雜訊，後者屬於高頻雜訊。
阿倫方差是 
接著依據圖上線段的斜率判斷隨機誤差的類型。但一般實務狀況，你的資料不見得會有或是看得出全部類型的隨機誤差，斜率也不會完美的與上圖相符，下圖是我自己的實測資料，可以看出應有3-4種隨機誤差。
當你判斷出有哪些隨機誤差後，再根據公式將阿倫標準差轉為隨機誤差值，可查閱下表。以 random walk 為例，當圖上有斜率 -0.5 的線段時，可根據定義查找該線段與 tau = 1 交會處的阿倫標準差值 (如果斜率 -0.5 的線段沒有經過 tau = 1 的位置，則需要先求出該線段的直線方程式，再用 x = 1 代入求解 y)。換言之，用 1 秒切分資料所計算的阿倫標準差，就是 random walk，而單位則依據你輸入的資料單位而定。
以陀螺 deg/s 為例，基於 random walk 的功率譜密度 (Power Spectral density, PSD) 推演與阿倫方差的關係得到下式的結論，阿倫標準差單位是deg/s， tau 的單位是秒，所以 randon walk 的單位為 deg/sqrt(s)：
Bias instability 則與截止頻率 f0 (cutoff frequency) 有關，當斜率為零的位置與阿倫曲線交會的 tau 遠大於 f0 時，可以將該交點的阿倫標準差直接除 0.664得到 bias instability。所以 bias instability 的單位有兩種表現形式。
規格書範例
以下圖 MEMS IMU BMI 160 的規格書為例，假設 Output Noise 等同是 bias instability ，當取樣頻率 (ODR) 為 200 Hz 時， noise 為 0.07 deg/s ，但在取樣頻率 10 Hz 時，從原始 PSD 推導公式， noise 為 0.007 deg/s/sqrt(Hz) = 0.007 deg/sqrt(s) 。這說明 bias instability 兩種單位的由來。
兩者的單位轉換可以直觀的乘除 sqrt(time)，參考來源：
ADIS16485 ARW = 18 deg/hour x 1/60 = 0.3 deg/sqrt(hour). https://ez.analog.com/mems/w/documents/4509/faq-gyroscope-angle-random-walk
另外，表中的 Zero-rate offset 就是系統誤差的 bias ，不同IMU的規格書可能用 bias 也可能用 offset，但原則上我認為是一樣的東西。表中也描述了 bias 的溫度效應，即溫度每上升 1 K，bias 會增加 0.05 deg/s。
此外，某些規格書可能會看到 “turn on bias” 和 “in run bias”，前者按字面是開機就存在的偏差，故我認為是指系統誤差的 bias，後者則是運行過程中的偏差，應是指隨機誤差的 bias instability。
再來 BMI 160 表中的 Sensitivity 就等於 scale factor，以 FS125為例：
(Value_max/min – Value_typical) / Value_typical
(16.9 - 16.4) / 16.4 = 3%
(15.9 - 16.4) / 16.4 = -3%
so in summary +/- 3%也就是實際資料量測值會是真值的 +/- 1.03倍。參考來源：
最後附上自駕車常用的 NovAtel PwrPak7 系列產品規格，有沒有看得懂了？
Allan Variance Note
最後值得一提的是，當用較大的時間長度分段時，代表你的資料集被切分的段數越少，所以當資料總時間長度不足時，阿倫方差圖靠右的值其實是不可靠的，所以我自己一般是不考慮最右側的隨機誤差。這可能也是一般IMU規格書多只看到 bias instability 與 random walk 的原因吧！
補充：一般阿倫變方分析都蒐集 2-8 小時的靜態數據。
IMU Serial Port
下面用 Python 做了一個從 Serial Port 讀取 MEMS IMU數據的範例：
import serial
import binascii
＃設定 IMU serial port name 與 baudrate
imu = serial.Serial('/dev/ttyUSB0', 115200, timeout=1)
# 寫入文字檔
path = 'imu.txt'
f = open(path, 'w')
# 以 Bytes 型別發送命令給 IMU (某些IMU可能需要)
# 依規格書替換 command，可為 Ascii 型別，可能包含 Checksum，換行符號\r\n
imu.write(b'command')
while True:
# 以 Bytes 型別讀取資料
imu_data = imu.readline()
print(len(imu_data), ' ', imu_data)
# 做一些必要的判讀，像是標頭為 data_header 與語句長度為 58
# data_header 可為 ascii or binary
if imu_data[0:8] == b'data_header' and len(imu_data) == 58:
try:
# int.from_bytes 把 Python 專有的 Bytes 型別轉 int
# byteorder 有 big 跟 little 要依規格書設定
# signed 也是依據規格書看該欄位是 signed 還是 unsigned binary
imu_count = int.from_bytes(imu_data[11:13], byteorder='big', signed=False) * 1
imu_status = int.from_bytes(imu_data[13:15], byteorder='big', signed=False) * 1
# imu_data[15:17] 表示第 16-17 共兩個 bytes 是陀螺 x 軸角速度資料
# MEMS IMU 原始數據單位多為 LSB (least significant bit)
# 以下面陀螺 X 軸的範例， signed 2 bytes = 1 sign bit + 15 data bits
# 15 data bits 對應十進位範圍為 -16384 ~ 16383，對照規格書陀螺的 range，是 +/- 100 deg/s
# LSB/deg/s = 32764 / 200 = 163.84 ， deg/s/LSB = 200 / 32768 = 0.006
# 故下面可將原始輸出 LSB * 0.006，就可將單位轉換為 deg/s
imu_rate_x = int.from_bytes(imu_data[15:17], byteorder='big', signed=True) * 200 / 32768 #°/s
imu_rate_y = int.from_bytes(imu_data[17:19], byteorder='big', signed=True) * 200 / 32768 #°/s
imu_rate_z = int.from_bytes(imu_data[19:21], byteorder='big', signed=True) * 200 / 32768 #°/s
imu_acc_x = int.from_bytes(imu_data[21:23], byteorder='big', signed=True) * 100 / 32768 #m/s2
imu_acc_y = int.from_bytes(imu_data[23:25], byteorder='big', signed=True) * 100 / 32768 #m/s2
imu_acc_z = int.from_bytes(imu_data[25:27], byteorder='big', signed=True) * 100 / 32768 #m/s2
# 可以透過 binascii.b2a_hex 將 Bytes 型別轉為16進制
# 可以透過 decode('ascii') 將 Bytes 型別轉為 ascii
"""
data_gx = binascii.b2a_hex(a[15:17])
data_gx = bytes(data_gx).decode('ascii')
"""
# 寫入轉換後的數據
#Full data
f.write('{} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {} {}\n'.format(imu_count, imu_status, imu_rate_x, imu_rate_y, imu_rate_z, imu_acc_x, imu_acc_y, imu_acc_z)
except:
pass
f.close()注意範例中的IMU對象是 tamagawa seiki 的產品，需要先對IMU寫入一個 request command，其他IMU並不一定需要。另外，若以 BMI 160 為例，規格書已經直接提供 LSB/deg/s，不需要照註解的步驟計算轉換參數。
測繪與 IMU 的關係
我的專長是測繪與空間資訊 (Geomatics)，也可以用測量 (Survey) 跟製圖 (Mapping) 概括，那為什麼會跟 IMU 扯上關係呢？讓我娓娓道來…
測量製圖技術從傳統使用光學經緯儀和全站儀的地面測量，演進到 GNSS RTK 單點測量，再到移動測繪系統 (Mobile Mapping System, MMS)，或簡稱測繪車，測繪的效率大幅提升。可以快速看一下我母系的介紹影片：
像路上的看到 Google 街景車，如果增配有高精度的定位設備、相機和光達 (需要遠高於自駕車用的規格)，就能夠用攝影測量 (Photogrammetry) 對蒐集的影像資料，進行量測與製圖，此時這台車就可以稱為 MMS 。台灣有將近十台 MMS 從事測繪與自駕車用的高精地圖製作，除了測繪公司的車以外，我還看過 Garmin, Apple 等公司的 MMS。下圖引用母系的初代與二代 MMS：
https://www.grb.gov.tw/search/planDetail?id=13422724
MMS 上搭載的高精度定位設備，就是由 GNSS RTK + IMU 組成。此時的 IMU 多會搭配計算單元，加入慣性導航演算法，就稱為慣性導航系統 (Inertial Navigation System, INS)。
在載體上的應用稱為捷聯式慣性導航 (Strap-down INS) ，在車用領域則又稱為航位推算 (Dead-Reckoning) ，並依據是否使用輪速計 (Speedmeter or Odometer) ，分為無線航位推算 (Untethered Dead Reckoning, UDR) 和車用航位推算 (Automotive Dead Reckoning, ADR)。
整合的目的是讓 GNSS 與 IMU / INS 做 sensor fusion 形成互補關係， GNSS 提供 INS 誤差約制的能力，而在隧道或市區等 GNSS 訊號不佳的路段，透過 INS 維持高精度定位，同時能夠估計載體姿態 (roll, pitch and yaw)，並高頻率的輸出相關結果 (20 Hz以上)。下表說明了互補關係：
為了滿足攝影測量的應用，高精度定位設備提供相機或光達精準的外方位參數，就顯得非常重要，概念可以參考下面另一篇文章。也因此當年我被指導教授領進門，從此踏入以 IMU 為主體的定位演算法研究。
參考資源：
https://docs.python.org/3/library/stdtypes.html
https://blog.csdn.net/qq_38212787/article/details/117996702
如有錯誤歡迎指正~