2018.12.16~12.31 RMT訊號

原始圖檔在此下載

108/01/01更新: 感謝Meteor Shin提供2018-12-29 10:43:41~12:59:54缺失圖檔, 每6秒1張

2018-12-25 11:56:44台東M4.5訊號時序

未來將此列為警報訊號來進行研究: "MR突然上升>40且Mw均值>4.0以上, 持續超過30秒"
此震跟11月澎湖地震一樣, 震前1分鐘開始有強烈訊號, 時序圖如下:

原始訊號圖影片:

2018.12.01~12.15 RMT訊號

原始圖檔在此下載

2018.11.28 21:43 澎湖餘震3.3

這個地震在大陸地震網站有紀錄, 台灣沒有, 可能是判定標準問題所以沒有列出. 澎湖正下方仍有一點訊號, 可能因震度3.3太小, 無法讓該點上方岩石產生震動, 所以少了很多點

2018.11.16~11.31 RMT訊號

 原始圖檔在此下載

2018.11.27 18:29 澎湖餘震芮氏4.7

訊號模式相同, 但下方訊號提前時間縮短為半分鐘

由此可知: 同一地點產生會相同訊號模式, 是一個符合物理原理與觀察事實的合理假設

2018.11.26 20:50 21:24 澎湖西方同地點餘震

同地點的訊號模式與上午07:57芮氏6.1相同, 下方1分鐘前先產生訊號, 然後才出現在澎湖島西方, 像是從下方往上跑, 但是中間有一段可能是地質太硬所以沒看到震動, 直接穿越過去跑到上方:

20:50 MR訊號圖:

21:24 MR訊號圖:

2018.11.26 07:57 澎湖西方芮氏6.1

澎湖西方芮氏6.1, 訊號從出震點直接穿過中央山脈,傳到台灣東方海域, 之後全台灣才開始產生後續訊號:

2018.11.01~11.15 RMT訊號

沖繩板塊: 參考https://zh.wikipedia.org/wiki/%E6%B2%96%E7%B9%A9%E6%9D%BF%E5%A1%8A
訊號路徑為張裂型板塊邊界(含左方), 對應訊號出現於台灣西北至北海域, 主要出震地點為東北方之日本與俄羅斯
訊號路徑為隱沒帶(含右方), 對應訊號出現於台灣東部外海, 主要出震地點為東北方之日本俄羅斯與太平洋

11/01 8時~11/3 6時, 可能是11/2 03:17日本~11/3 2:29俄羅斯震前震後訊號, 連動台灣宜蘭外海小震 (訊號路徑為張裂型板塊邊界)
11/03 15時, 可能是16:44印尼震前訊號
11/06 5時~8時, 可能是09:00印尼震前訊號
11/11 5時~8時, 東北角垂直線訊號, 可能是09:50俄羅斯震前訊號 (訊號路徑為隱沒帶)
11/11 12時~15時, 東北角垂直線訊號, 可能是16:28日本震前訊號 (訊號路徑為隱沒帶)
11/14 15時~15日11時, 因硬碟已滿導致資料未儲存, 於11/15 12時修復

 原始圖檔在此下載

[Python] 以tensorflow分析Mw的雙高斯混合模型

根據中央極限定理: 大量相互獨立隨機變數的均值經適當標準後為常態分布。RMT不論有無發生地震均會不斷一直計算, 以目前觀察到的Mw為例, 平時多在3.5以下, 當遠方地震傳來時, 經驗上常看到Mw>4的現象. 假設每次解出的Mw是一個隨機變數Xi, 那麼:

X=X1+X2+X3+....

不論原來的Xi是甚麼分佈, X應是常態分佈. 如果無震時的均值Mw0, 有震時的均值Mw1, 那麼長時間統計下來, 應該會看到兩個以Mw0與Mw1為均值的常態分佈混和體, 分別代表有震於無震的分佈情形 (實際分佈當然是跟地震生發生的情形有關). 以今年5月分的資料為例, 出現的次數圖形如下:

看起來的確是這麼一回事. 要如何確認呢? 最近流行AI技術, 所以也來玩玩看. 假設Y1代表無震時的分佈, Y2代表有震時的分佈, 且Y1與Y2皆為常態分佈, 兩個分佈有不同的最大值,平均值與標準差, 那麼Y=Y0+Y1應該可以找到兩個常態分佈和的方程式來fit. 從fit後的誤差大小,  可以評估符的準確度.

首先是讀取檔案:

建立 tensorflow graph: y1_pred與y2_pred為兩個常態分佈, 參數a1/u1/s1與a2/u2/s2, 一般來說都會用亂數去產生, 但這種方式常常難以收斂, 因為我們用眼睛就可以看到這兩個分佈的長相, 所以這裡直接給個大概的數字, 然後再讓tensorflow的optimizer去找, 這樣可以省掉很多時間:

執行session:這裡設定loss<0.018或太久沒有降低loss就跳出

最後輸出結果與圖形:

執行過程:

執行結果:

共執行4416次, a1_final= 0.99921775 , u1_final= 2.6671188 ,s1_final= 0.25635517
a2_final= 0.761186 , u2_final= 3.5796206 ,s2_final= 0.32652023

最後跳出時的誤差為 0.017999463, 資料48筆, 平均每筆誤差 0.00037498880798618, 看起來還不錯, 可以確認Mw是由兩個常態分佈混合而成.

求出這些參數後，未來可以用於分類問題的最佳門檻值計算。假設兩個類別的分佈:

若取樣之樣本應屬於第1類而被歸類為第2類，或是應屬於第2類而被歸類為第1類，此兩種分類錯誤的機率，如下圖的綠色區域所示，求 x ̃ 使得 f1 (x ̃)=f2 (x ̃)， x ̃ 是使系統分類錯誤機率最小的門檻值。若觀察到的Mw值大於 x ̃ ，且落在第2類(有出震)分佈的2σ範圍內，則有較高的機率是屬於已發生或即將發生之地震訊號。

把datafile改成 RMT_Mw_201802, 重新跑一次, 可以看到206花蓮地震對分佈的影響:

程式跟資料在這裡下載

2018.10.16~10.31 RMT訊號

10/24日21時~25日3時, 25日19時~26日2時持續強烈訊號, 可能是02:36日本石卷M5.7前兆 (出震後2小時訊號消失)
10/26日7時~11時北方強烈訊號, 可能是11:04俄羅斯M5.6前兆 (出震後2小時訊號消失)

 原始圖檔在此下載

2018.10.23 西藏 M4.8, 波速僅4.7Km/

00:49:33 西藏日喀則市 M4.8, 30.46N  87.76E, 深度 11Km, 距離台灣3381Km, 波速=3381/12/60=4.7Km/s, 低於一般波速, 因此特別紀錄之


最大訊號 Mw=4.4 發生於:2018-10-23 01:02:11
Mw>=3.9 持續時間: 2018-10-23 01:02:11 ~ 2018-10-23 01:10:01
Mw>=4.2 持續時間: 2018-10-23 01:02:11 ~ 2018-10-23 01:02:16

2018.10.01~10.15 RMT訊號

10/03 13時, 15時, 無對應地震
10/04 10時起, 康芮颱風靠近, 東北角連續訊號
10/08 11~13時, 疑為13:45日本硫磺島M4.9臨震訊號
10/08 13:45日本硫磺島M4.9, 造成東北方強烈訊號, 可能與10/10 17:43宜蘭M4.4有關
10/15 5~7時, 可能與07:33 印尼M4.3有關
10/15 15~20時, 未出震

原始圖檔在此下載

[Python] 從中央氣象局下載地震活動彙整列表

氣象局網站的 地震活動彙整 列表, 檢視網頁內容, 真正資料網頁連結為 https://scweb.cwb.gov.tw/Page.aspx/?ItemId=20&loc=tw&adv=1, 因為是ASP網頁, 所以需先取出ASP的傳遞參數, 資料表格之 class 為 datalist4, 使用 BeautifulSoup 即可取出表格內容

# -*- coding: utf-8 -*- """ Created on Sun Oct  15 11:53:25 2018 @author: ghosty @Program: cwbweb_list.py @Prupose: download quake data from CWB web source: https://www.cwb.gov.tw/V7/earthquake/rtd_eq.htm """ import requests from bs4 import BeautifulSoup #import dateutil    def downloadCWBweb(year,month):     url = 'https://scweb.cwb.gov.tw/Page.aspx/?ItemId=20&loc=tw&adv=1'                  agent = 'Mozilla/5.0 (Windows NT 6.3; WOW64; rv:34.0) Gecko/20100101 Firefox/34.0'     headers = {'Content-type': 'application/x-www-form-urlencoded',            'Accept': 'text/html,application/xhtml+xml,application/xml;q=0.9,*/*;q=0.8',            'User-Agent': agent}     payload = {             '__VIEWSTATE':'',             '__VIEWSTATEGENERATOR':'',             '__VIEWSTATEENCRYPTED':'',             '__EVENTVALIDATION':'',             'ctl03_ddlYear':'',             'ctl03_ddlMonth':'',             'ctl03_btnSearch':''         }     response1 = requests.post(url)     if  response1.status_code != requests.codes.ok:         print("CWB requesr fail")         return         soup = BeautifulSoup(response1.text, "lxml")     payload['__VIEWSTATE']=soup.find('input',id='__VIEWSTATE')['value']         payload['__VIEWSTATEGENERATOR']=soup.find('input',id='__VIEWSTATEGENERATOR')['value']     payload['__VIEWSTATEENCRYPTE']=soup.find('input',id='__VIEWSTATEENCRYPTED')['value']     payload['__EVENTVALIDATION']=soup.find('input',id='__EVENTVALIDATION')['value']     payload['ctl03$ddlYear']="{:4d}".format(year)     payload['ctl03$ddlMonth']="{:0>2d}".format(month)     payload['ctl03$btnSearch']=''         response2 = requests.post(url,data=payload, headers=headers)     soup2 = BeautifulSoup(response2.text, "lxml")     table = soup2.find('table', attrs={'class':'datalist4'})     rows = table.find_all('tr')     quakeData  = []     for row in rows:         cols = row.find_all('td')         cols = [item.text.strip() for item in cols]         if (len(cols)>0): #skip empty row                    quakeData.append([item for item in cols if item]) # Get rid of empty values                 return quakeData quakeData = downloadCWBweb(2018,8) for data in quakeData:     print(data)

啟用地震資料自動匯入比對系統

為了避免人工失誤與主觀因素, 從USGS下載資料, 自動匯入動態圖檔的功能, 已經啟用. 這個功能是先從USGS下載所有地理位置資料, 將英文地名翻譯成中文, 建立一個地理字典來做的, 由於常發生的位置很固定, 所以字典的數量並不多, 目前只建立約1000筆資料, 未來可以慢慢加強. 產生的資訊會根據方位來顯示, 總共有12個方向,每個方位只顯示震度最大的那一個. 根據經驗, 排除距離太遠且震度太小的地震:

(1)距離>3000且震度<4.5

(2)距離>6000且震度<5

(3)距離>9000且震度<6
(4(前一個小時地震時間與小時開始的時間秒數差>距離/6的秒數 (合理的波速)

之前以訊號尋找對應出震的方式, 無法看出遠方地震的影響程度, 現在不論有無訊號皆會顯示, 可看出該地點於該時間的影響程度有多少.

[Python] 從USGS下載地震資料

為了分析地震地點的地理名稱, 使用Python語言, 從USGS下載地震資料存成文字檔. 因 USGS 使用逗點作為分隔, 因此這裡使用 $ 錢號來做輸出資料區隔, 以避免資料混淆, 亦可匯入 excel 進行整理. 輸出資料有兩個欄位, 第一個是原始地震資料來源, 第二個是擷取的地點名稱, 打算做連續地點關聯性分析, 以及加上自動翻譯為中文, 作為資料與動畫分析時自動匯入使用

# -*- coding: utf-8 -*- """ Created on Sun Oct  7 11:53:25 2018 @author: ghosty @Program: usus_location_list.py @Prupose: download USGS data and analyze the quake location name for           (1) location relation analysis           (2) translation from English to Traditional Chinese, used for AI analysis tool development """ import csv import requests import os.path import datetime #import dateutil from dateutil.relativedelta import relativedelta #init global variables locFileName = 'usgs_loc.txt' locationList=[] quakeList=[] #    def downloadUSGS(starttime, endtime, minmagnitude, maxmagnitude):  #,latitude=121.1,longitude=23.5,maxradiuskm=20000):     url = 'http://earthquake.usgs.gov/fdsnws/event/1/query'     query = {             'format': 'csv',             'starttime': starttime,             'endtime': endtime,             'minmagnitude': '5',             'eventtype':'earthquake'             }     response = requests.get(url, params=query)     if  response.status_code != requests.codes.ok:         print("USGS requesr fail")         return     #print("encoding: %s" % response.encoding)     quakeList = list(csv.reader(response.text.split('\n'), delimiter=','))     header=quakeList[0]     #time,latitude,longitude,depth,mag,magType,nst,gap,dmin,rms,net,id,updated,place,type,horizontalError,depthError,magError,magNst,status,locationSource,magSource     quakeList.pop(0) #remove header line     print('USGS download #',len(quakeList))     return quakeList, header def addLocation(quakeData):     for quake in quakeData:         #print('quake:',quake)         if (len(quake)==0): #skip null record             continue         quakeTime = quake[0]         quakeLatitude = quake[1]         quakeLongitude = quake[2]         quakeDepth = quake[3]         quakeMag = quake[4]         quakeMagType = quake[5]         quakePlace = quake[13]         #print(quakeTime, quakeLatitude, quakeLongitude, quakeDepth, quakeMag, quakeMagType, quakePlace)         locs = quakePlace.split(',')         #print(locs) #debug for 2013-02-01         for loc in locs:             try:                 ofPos = loc.find(" of ")    #do not usr 'of' only to prevent error                 if (ofPos>0):                     loc=loc[ofPos+4:]                 #loc.replace("the ","")                 loc.replace("the ","")                 #loc.strip() #remove space, strip() fail to use                 while (loc[0]==' '):                       loc=loc[1:]                 if (loc not in locationList):                     quake = '{0} {1} {2} {3} {4} {5} {6}'.format(quakeTime, quakeLatitude, quakeLongitude, quakeDepth, quakeMag, quakeMagType, quakePlace)                     print("Add [",loc, '] from [',quake,']')                     quakeList.append(quake)                     locationList.append(loc)                      #print(len(quakeList),len(locationList))             except:                 #source data error                 print('parse fail: ',quakeTime, quakeLatitude, quakeLongitude, quakeDepth, quakeMag, quakeMagType, quakePlace)                 def saveLocation(locationList, locFileName):     #print(len(locationList))        print("export ",len(locationList)," locations")     with open(locFileName, 'w') as locfile:         for i in range(len(locationList)):             locfile.write('\"'+str(quakeList[i]) + '\"$\"' + str(locationList[i]) + '\"\n')  def loadLocation(locationList, locFileName):     if (os.path.exists(locFileName)):         with open (locFileName, 'r') as locfile:            for line in locfile:                quakeData = line.strip().split('$')                quakeList.append(quakeData[0].strip('\"'))                locationList.append(quakeData[1].strip('\"'))         print("import ",len(locationList)," locations")         #------main()-------- # Load processed location data loadLocation(locationList, locFileName)  #for all date dt1 = datetime.date(1900, 1, 1) #end_dt =  datetime.date(2000, 1, 31) end_dt =  datetime.date.today() while (dt1<end_dt): # in daterange(start_dt, end_dt):     """     # if there are too many records at one download, reduce the download period     if (dt1 < datetime.date(2018, 1, 1)):         delta_dt = relativedelta(years=1)     else:         delta_dt = relativedelta(months=1)     """     delta_dt = relativedelta(years=1)     dt2 = dt1 + delta_dt     print("\nProcessing...",dt1.strftime("%Y-%m-%d"),"~",dt2.strftime("%Y-%m-%d"))     quakeDATA, header = downloadUSGS(dt1, dt2, 4.0, 10)     addLocation(quakeDATA)      dt1 = dt2          #save processed result saveLocation(locationList, locFileName)