<style> .caption { font-weight: bold; color: #555 } </style> <style> .subcaption { color: #555 } </style> # Seismic Analysis [TOC] ## 軟體：SAC ### 下載（MAC） 1. 下載圖形顯示界面 [XQuartz](https://www.xquartz.org) 1. [申請 SAC 軟體](http://ds.iris.edu/ds/nodes/dmc/forms/sac/)，得到 `sac-102.0-mac.tar.gz` <!--[由此下載](https://drive.google.com/file/d/1qkeoUllwlXPtPWH39IlWTLBLOnwMPmwz/view?usp=sharing)（請不要外流）--> 3. 使用 Terminal（或手動）解壓縮得到 `sac` 目錄，並將檔案移動到 `usr/local/`： ```Shell ## 使用終端機操作 cd Downloads/ #cd: 移動到下載項目或壓縮檔所在位置 tar -xvf sac-102.0-mac.tar.gz #tar: 解壓縮 sudo mv sac /usr/local #mv: 移動檔案 ``` ![](https://i.imgur.com/2AEJIQ7.png) :::success **Tip:** `usr/` 默認爲隱藏狀態，在 Finder 中輸入 ⌘+⇧+. 來永遠顯示隱藏的資料夾。 ::: 1. 編輯並套用 `/.zshrc` 文件： ![](https://i.imgur.com/T07nfkd.png) 手動移動到 home 目錄（在 Finder 輸入 ⌘+⇧+H）找到文件或在 Terminal 開啟文字編輯器： ```Shell nano /.zshrc #nano: 開啟文字編輯器 ``` 貼上以下文字並儲存： ``` export SACHOME=/usr/local/sac export SACAUX=${SACHOME}/aux export PATH=~/bin:${SACHOME}/bin:/usr/local/gfortran/bin:${PATH} export SAC_DISPLAY_COPYRIGHT=1 export SAC_PPK_LARGE_CROSSHAIRS=1 export SAC_USE_DATABASE=0 ``` 在 Terminal 輸入： ```Shell source /.zshrc ``` 1. 啟動 SAC： 在 Terminal 輸入： ```Shell sac ``` ![](https://i.imgur.com/rrGazaH.png) ### 檔案格式 #### header 與測站、資料相關的基本資訊紀錄在 .sac 的 "header" 中。 使用 [`listhdr` 指令](https://seisman.github.io/SAC_Docs_zh/commands/listhdr/)取得 header。 :::success 以 InSight Apr 25,2020 BHU channel 爲例： ``` SAC> lh FILE: XB.ELYSE.02.BHU.M.2020.116.000000.SAC - 1 ------------------------------------------- NPTS = 3456004 ! 資料點數 B = 0.000000e+00 ! 資料起始時間相對參考時間的秒數 E = 1.728002e+05 ! 資料結束時間相對參考時間的秒數 IFTYPE = TIME SERIES FILE LEVEN = TRUE DELTA = 5.000000e-02 ! 資料間隔（s） IDEP = VELOCITY (NM/SEC) ! 資料類型（經過去儀器響應） DEPMIN = -1.026158e+02 ! 資料最小值（自動計算） DEPMAX = 5.253109e+01 ! 資料最大值（自動計算） DEPMEN = -5.083885e-08 ! 資料平均值（自動計算） KZDATE = APR 25 (116), 2020 ! 參考時間（不一定是資料起始時間） KZTIME = 00:00:00.035 ! 參考時間（不一定是資料起始時間） KINST = VBB Velo ! 儀器類型 KSTNM = ELYSE ! 測站名稱 CMPAZ = 1.351000e+02 CMPINC = 6.060000e+01 STLA = 4.502384e+00 ! 測站緯度 STLO = 1.356234e+02 ! 測站經度 STEL = -2.613400e+03 ! 測站海拔（m） STDP = -1.000000e-01 KHOLE = 02 LOVROK = TRUE NVHDR = 6 SCALE = 8.626220e+10 ! 系統常數（經過去儀器響應，真正的值是 10^9） LPSPOL = TRUE LCALDA = TRUE KCMPNM = BHU ! 通道名稱 KNETWK = XB ! 測站網名稱 ``` ::: ### 資料處理 #### 檔案 I/O * 讀入資料（[`read` 指令](https://seisman.github.io/SAC_Docs_zh/commands/read/)） ``` r *BH*SAC ``` 會讀入所有檔名依序包含 `BH` 和 `SAC` 的 .sac 檔案。 * 寫入資料（[`write` 指令](https://seisman.github.io/SAC_Docs_zh/commands/write/)） ``` w append .test ``` 會寫入記憶體中現有的檔案，檔案後綴加上 `.test`。 #### 基本資料處理 * 去趨勢（[`rmean` 指令](https://seisman.github.io/SAC_Docs_zh/commands/rmean/)、[`rtrend` 指令](https://seisman.github.io/SAC_Docs_zh/commands/rtrend/)、[`taper` 指令](https://seisman.github.io/SAC_Docs_zh/commands/taper/)） * 濾波（[`bandpass` 指令](https://seisman.github.io/SAC_Docs_zh/commands/bandpass/)、[`highpass` 指令](https://seisman.github.io/SAC_Docs_zh/commands/highpass/)、[`lowpass` 指令](https://seisman.github.io/SAC_Docs_zh/commands/lowpass/)） #### SCALE 和去儀器響應 地震資料原則上長度單位爲公尺。但是，地動震動的尺度通常遠小於公尺，爲了方便運算，常會將資料乘以某數值加以放大，此數值紀錄在 SAC 檔的 `SCALE` header 中，有時也稱爲 sensitivity。 通常地動訊號中，以日變化最爲明顯，震幅遠高於其他訊號。因此地震資料通常會過濾低頻訊號，稱爲「儀器響應」，如圖： </br> ![](https://i.imgur.com/Jn7ujFo.png) <span class="caption"> [儀器響應示意圖。](http://ds.iris.edu/mda/XB/ELYSE/02/BHU/?starttime=2021-06-19&endtime=2021-06-19)</span> <span class="subcaption"> 可以想像成不同頻率的波有不一樣的放大倍率（SCALE），對於這個頻道而言，0.1 Hz 以下和 10.0 Hz 以上的震動都放大倍率都較低。</span> </br> 參考上圖，`SCALE` 記爲 $8.626\times10^{10}\ \mathrm{count/(m\cdot s^{-1})}$，對於頻率 $0.1-10\ \mathrm{Hz}$ 的震動，直接將數值除以 `SCALE` 即可得到真實地動速度，單位爲 $\mathrm{m\cdot s^{-1}}$。若想要得到其他頻率的真實地動速度，則需要進行「**去儀器響應**」： 1. [取得儀器響應文件](http://service.iris.edu/irisws/resp/docs/1/builder/)，將結果儲存於文件中（建議命名爲 `RESP.ALL`） 2. 啟動 SAC，讀入目標文件後，使用 [`transfer` 指令](https://seisman.github.io/SAC_Docs_zh/commands/transfer/?highlight=transfer)： ```SAC transfer from evalresp fname RESP.ALL to vel ``` 此時 SAC 文件中的單位變爲 $\mathrm{nm}$，相當於 $\mathrm{SCALE}=1\times10^9$（但 header 不會自動更改）。 ### 其他參考資料 [參考手冊（中文）](https://seisman.github.io/SAC_Docs_zh/) ## 軟體：ObsPy ObsPy 是一個 Python 函式庫 ### 下載 1. 下載 [Anaconda](https://www.anaconda.com/products/distribution#macos) 1. 啟動 Anaconda-Navigator，前往 'Environment'，點擊左下角 'Create' 建立一個新的環境 1. 點擊最右側搜尋並下載 'ObsPy'（搜尋類別設爲 'All'） 1. 使用 Jupyter Lab 啟動 ObsPy： 點擊最右側返回 'Home'，在最上方選擇有下載 ObsPy 的環境，往下找到 'Jupyter Lab' 並點擊 'Launch' ### 其他參考資料 [Documentation (English)](https://docs.obspy.org/tutorial/index.html) [Documentation (Chinese)](https://docs.obspy.org/archive/ObsPy_Tutorial_2020-04_chinese.pdf) ## 資料：InSight ### 基本資訊 InSight 在 May, 2018 降落在火星。 </br> ![](https://i.imgur.com/8Yw5mPa.jpg) <span class="caption"> [InSight 降落位置示意圖。](https://mars.nasa.gov/insight/timeline/prelaunch/landing-site-selection/)</span> <span class="subcaption"> 降落在4.5024°N 135.6234°E（Elysium Planitia，*埃律西昂平原*）。</span> </br> 相關的觀測儀器包含 SEIS（前面提到的地震儀）、TWINS（氣象資料）等。 火星日略長於地球日。1 SOL (Mars day) = 24 hr 39 min and 35.244 sec. ### 地震資料下載 * #### 震波資料 前往 [IRIS URL Builder](https://service.iris.edu/fdsnws/dataselect/docs/1/builder/)，選擇 | 欄位 | 輸入 | 說明 | | -------- | ------- | -------------------------- | | Network | XB | | | Station | ELYSE | Scientific data[^SeedConfig] | | Location | 留白 | | | Channel | BH* | 選擇所有 BH 開頭的 channel | | Format | SAC zip | 或選擇 miniSEED | 並點擊生成的連結，形如 ``` https://service.iris.edu/fdsnws/dataselect/1/query?net=XB&sta=ELYSE&cha=*&starttime=2020-04-25T00:00:00&endtime=2020-04-25T01:00:00&format=sac.zip&nodata=404 ``` :::success **Tip:** 也可以直接編輯連結內容已達到大量下載資料的目的。 ::: * #### 儀器響應文件下載 前往另一個 [IRIS URL Builder](http://service.iris.edu/irisws/resp/docs/1/builder/)，選擇對應資料。 ### 資料格式 以 BH channel 爲例[^ChannelNaming]，共有 BHU/ BHV/ BHW 三個頻道。三個頻道並不正交，如圖： </br> ![](https://i.imgur.com/fXwNd6G.jpg) <span class="caption"> [地震儀頻道方位示意圖。](https://www.sciencedirect.com/science/article/pii/S0031920120302752)</span> <span class="subcaption">BH* 通道使用的儀器爲 VBB。</span> </br> | 通道 | 方位角（azimuth, $\phi$） | 傾角（dip, $\delta$） | | -------- | --------:| --------:| | U | $135.1^\circ$ | $-29.4^\circ$ | | V | $255.0^\circ$ | $-29.7^\circ$ | | W | $345.3^\circ$ | $-29.2^\circ$ | 方位角爲從正北順時針旋轉；傾角以向下爲正。 可以將 UVW 座標系轉換爲 NZE 座標系：[^SuppText] \begin{align} \begin{bmatrix}Z\\N\\E\end{bmatrix}= \begin{bmatrix} −\sin(\delta_U) & \cos(\delta_U)\cos(\phi_U) & \cos(\delta_U)\sin(\phi_U)\\ −\sin(\delta_V) & \cos(\delta_V)\cos(\phi_V) & \cos(\delta_V)\sin(\phi_V)\\ −\sin(\delta_W) & \cos(\delta_W)\cos(\phi_W) & \cos(\delta_W)\sin(\phi_W) \end{bmatrix}^{-1} \begin{bmatrix}U\\V\\W\end{bmatrix}. \end{align} 利用 python 實現： ```Python= import numpy as np from numpy.linalg import inv azi = np.array([135.1,255.0,345.3]) dip = np.array([-29.4,-29.7,-29.2]) azi = azi/(180)*np.pi dip = dip/(180)*np.pi A = ([[ -np.sin(dip[0]), np.cos(dip[0])*np.cos(azi[0]), np.cos(dip[0])*np.sin(azi[0]) ], [ -np.sin(dip[1]), np.cos(dip[1])*np.cos(azi[1]), np.cos(dip[1])*np.sin(azi[1]) ], [ -np.sin(dip[2]), np.cos(dip[2])*np.cos(azi[2]), np.cos(dip[2])*np.sin(azi[2]) ]]) X = np.vstack((st[0],st[2],st[2])) Y = np.matmul(inv(A),X) Z = Y[0,:] N = Y[1,:] E = Y[2,:] ``` 利用 MATLAB 實現： ```Matlab function [Z,N,E] = uvw2zne(U,V,W) azi = [135.1,255.0,345.3]; dip = [-29.4,-29.7,-29.2]; azi = deg2rad(azi); dip = deg2rad(dip); A = [ -sin(dip(1)), cos(dip(1))*cos(azi(1)), cos(dip(1))*sin(azi(1)) ;... -sin(dip(2)), cos(dip(2))*cos(azi(2)), cos(dip(2))*sin(azi(2)) ;... -sin(dip(3)), cos(dip(3))*cos(azi(3)), cos(dip(3))*sin(azi(3)) ]; X = [U;V;W]; Y = A\X; Z = Y(1,:); N = Y(2,:); E = Y(3,:); end ``` ### 其他參考資料 [Official InSight mission website](https://mars.nasa.gov/insight/) [氣象資料下載](https://atmos.nmsu.edu/data_and_services/atmospheres_data/INSIGHT/insight.html#Selecting_Data) ## Paper [^SeedConfig]:https://www.seis-insight.eu/static/naming-conventions/index.html [^ChannelNaming]:https://ds.iris.edu/ds/nodes/dmc/data/formats/seed-channel-naming/ [^SuppText]:https://agupubs.onlinelibrary.wiley.com/action/downloadSupplement?doi=10.1029%2F2020EA001317&file=ess681-sup-0001-2020EA001317-SI.pdf