時間序列相似性度量方法
時間序列相似性度量常用方法為歐氏距離ED(Euclidean distance)和動態時間規整DTW(Dynamic Time Warping)。總體被分為兩類: 鎖步度量(lock-step measures) 和彈性度量(elastic measures) 。鎖步度量是時間序列進行 “一對一”的比 較; 彈性度量允許時間序列進行 “一對多”的比較。
歐氏距離屬於鎖步度量。
在時間序列中，我們通常需要比較兩端音訊的差異。而兩段音訊的長度大部分是不相等的。在語音處理領域上表現為不同人的語速不同。即時同一個人不同一時刻發同一個音，也不可能具有完全相同的時間長度。而且每個人對同一個單詞的不同音素的發音速度也是不同的，有的人會把"A"這個音拖得稍長，或者"i"稍短。在這種複雜的情況下，利用傳統的歐幾里得距離是無法求得有效的兩個時間序列之間的相似性的（即距離）。
如下圖，實線波形的a點會對應於虛線波形的b’點，這樣傳統的通過比較距離來計算相似性很明顯不靠譜。因為很明顯，實線的a點對應虛線的b點才是正確的。
考慮序列長短不一，或時間步不對齊的時候，特別是在峰值的時候，歐氏距離是無法有效計算兩個時間序列的距離。DTW演算法將其中一個序列進行線性放縮、進行某種“扭曲”操作，。可以存在一對多mapping的情況，在時間序列探索過程中效果更好。適用於複雜時間序列，屬於彈性度量
dtw演算法介紹
假設我們有兩個時間序列Q和C，他們的長度分別是n和m：（實際時間序列匹配運用中，一個序列為參考模板，一個序列為測試模板）
Q = q1, q2,…,qi,…, qn ;
C = c1, c2,…, cj,…, cm ;
為了對齊這兩個序列，我們需要構造一個n x m的矩陣網格，矩陣元素(i, j)表示qi和cj兩個點的距離d(qi, cj)（也就是序列Q的每一個點和C的每一個點之間的相似度，距離越小則相似度越高。這裡先不管順序），一般採用歐式距離，d(qi, cj)= (qi-cj)2（也可以理解為失真度）。每一個矩陣元素(i, j)表示點qi和cj的對齊。DP演算法可以歸結為尋找一條通過此網格中若干格點的路徑，路徑通過的格點即為兩個序列進行計算的對齊的點。
那麼這條路徑我們怎麼找到呢？那條路徑才是最好的呢？也就是剛才那個問題，怎麼樣的warping才是最好的。
首先，這條路徑不是隨意選擇的，需要滿足以下幾個約束：
1）邊界條件：w1=(1, 1)和wK=(m, n)。任何一種語音的發音快慢都有可能變化，但是其各部分的先後次序不可能改變，因此所選的路徑必定是從左下角出發，在右上角結束。
2）連續性：如果wk-1= (a’, b’)，那麼對於路徑的下一個點wk=(a, b)需要滿足 (a-a’) <=1和 (b-b’) <=1。也就是不可能跨過某個點去匹配，只能和自己相鄰的點對齊。這樣可以保證Q和C中的每個座標都在W中出現。
3）單調性：如果wk-1= (a’, b’)，那麼對於路徑的下一個點wk=(a, b)需要滿足0<=(a-a’)和0<= (b-b’)。這限制W上面的點必須是隨著時間單調進行的。以保證圖B中的虛線不會相交。
結合連續性和單調性約束，每一個格點的路徑就只有三個方向了。例如如果路徑已經通過了格點(i, j)，那麼下一個通過的格點只可能是下列三種情況之一：(i+1, j)，(i, j+1)或者(i+1, j+1)。
滿足上面這些約束條件的路徑可以有指數個，最佳路徑是使得沿路徑的積累距離達到最小值這條路徑。這條路徑可以通過動態規劃（dynamic programming）演算法得到。具體搜尋或者求解過程的直觀例子解釋可以參考：http://www.cnblogs.com/tornadomeet/archive/2012/03/23/2413363.html
普通dtw方法引入和使用
第一步：dtw包安裝，以Anaconda為例
從開始選單開啟 Anaconda prompt 管理工具
輸入dwt安裝指令碼
pip install dtw
第二步：dtw包匯入
from dtw import dtw
第三步：dtw函式使用
首先，需要定義距離函式，其中，x為參照序列資料，y為測試序列資料
distance=lambda x,y :np.abs(x-y)
接下來，呼叫dtw方法
d,cost_matrix,acc_cost_matrix,path=dtw(x,y,distance)
dtw函式輸入引數說明
x：參照序列資料
y：測試序列資料
distance：距離函式
dtw函式輸出引數說明
d：累計距離
cost_matrix：兩個序列之間的兩兩距離
acc_cost_matrix：dp矩陣
path：最佳路徑，包含兩個陣列，path[0] 參照序列下標值，path[1] 測試序列下標值
軌跡線列印出圖
plt.imshow(cost_matrix.T,origin='lower',cmap='gray')
plt.plot(path[0],path[1])
plt.show()
dtw-python包的安裝與引用
dtw庫的使用限制太多，不夠靈活，且作圖不夠方便，主要體現運算量大、首尾必須匹配、序列間對應個數無法限定等。dtw-python包是R語言中dtw實現的python版本，基本的API是對應的，它的優勢在於能夠自定義點的匹配模式，約束條件，和滑動視窗。同時提供方便的作圖和快速的計算（C語言的核心），官方文件點選這裡。
dtw-python 包安裝指令碼
pip install dtw-python
dtw-python 包引用
from dtw import dtw,dtwPlot,stepPattern,warp,warpArea,window #不建議直接 from dtw import * 方式，避免衝突
由於在dtw包選用過程中，先安裝了dtw包，後安裝了dtw-python包，在dtw-python引用過程中，出現dtw包引用失敗情況，建議兩者保留其一。可使用包解除安裝指令碼對dtw包進行解除安裝，並在Anaconda安裝目錄下 \Anaconda3\Lib 目錄下清除掉 以 ~ 開頭命名的資料夾
pip uninstall dtw
dtw函式的引數說明
ds=dtw(x, y=None, dist_method='euclidean', step_pattern='symmetric2', window_type=None, window_args={}, keep_internals=False, distance_only=False, open_end=False, open_begin=False)
輸入引數說明：
dist_method 定義你用的距離方法，預設為’euclidean’，即歐幾里得距離。
keep_internals=True儲存內部的資訊包含距離矩陣之類的，一般我們把它設為True。
step_pattern定義了點之間的匹配模式，有好幾種，具體檢視官網。
window_type表示全域性條件約束，也有幾種模式，同檢視官網。"none", "itakura", "sakoechiba", "slantedband"，其中sakoechiba方式需要在window_args 引數中傳遞window_size引數，方式如下：window_args= {"window_size": 1}
distance_only如果設定為True，會加速計算，只計算距離，忽略匹配函式。
open_end和open_begin設為True表示無拘束的匹配，即完全的部分匹配，預設是全域性匹配，就是嚴格對應頭和尾。
輸出引數說明：
ds.distance：累計距離
ds.index1：參照序列下標值
ds.index2：測試序列下標值
dtw輸出結果列印
ds.plot(type="twoway",offset=-2)
一般有twoway和threeway兩種模式，下面詳解，offset=-2可以理解為兩根線之間的分離程度，為了方便看清可以設的大一些。下圖左側為twoway方式的輸出結果，右側為threeway方式的輸出結果。
dtw函式使用注意事項
1.資料中不能出現空置和無效值，使用上一記錄替代缺失資料時，同一時段內缺失資料不能過多，否則會造成變化趨勢匹配結果偏差較大
2.最大值和最小值可能不會做到正確匹配
dtw使用例項
背景：同一條河流上兩個斷面的汙染物濃度在變化趨勢上呈現一定的相似性。但由於河流流速的不定，上游斷面濃度的波峰波谷在下游斷面上出現的時間間隔不同，可以通過dtw方法實現兩斷面間相似變化趨勢的時間段，以及上下游斷面間汙染物濃度變化的相互關係
第一步：新增python包引用
import numpy as np
import pandas as pd
import matplotlib.pylab as plt
import math
import datetime
import time #用於時間秒轉換成時間
from dtw import dtw,dtwPlot,stepPattern,warp,warpArea,window #不建議直接 from dtw import * 方式，避免衝突
import sys
import pymssql
import xlwt
from pandas._libs.tslibs.timestamps import Timestamp
第二步：通過csv方式讀取兩個斷面汙染物的濃度監測數值
beginTime='2019/1/2 12:00'
endTime='2019/2/18 12:00'
#從csv檔案中取出資料
df=pd.read_csv('wq3.csv',encoding='utf-8', index_col='datatime') #從csv檔案中讀取時間序列資料，index_col列定義為索引物件
df.index=pd.to_datetime(df.index)
df_xy_org=df['qingzhoucun'][beginTime:endTime] #指定時間序列中下游的監測資料列
df_sy_org=df['liangtiancun'] #指定時間序列中上游的監測資料列
#具體計算中需要排除掉缺數的情況，特別是長時間的缺數，會影響到最終的匹配結果。通過使用上一次監測資料方法不合適
for i in range(0,len(df_xy_org)):
if not df_xy_org[i] or math.isnan(df_xy_org[i]): #判斷數值是否為空
df_xy_org[i] = df_xy_org[i-1]
for i in range(0,len(df_sy_org)):
if not df_sy_org[i] or math.isnan(df_sy_org[i]): #判斷數值是否為空
df_sy_org[i] = df_sy_org[i-1]
第三步：為了更有效的找到波峰和波谷，對汙染物濃度做一階差分處理
#比對資料做一階差分處理
df_xytime=np.empty(len(df_xy_org)-1,dtype=datetime.datetime)
df_xyvalue=np.empty(len(df_xy_org)-1)
for i in range(1,len(df_xy_org)-1):
df_xytime[i-1]=df_xy_org.index[i].to_pydatetime()
df_xyvalue[i-1]=(float(df_xy_org.values[i])-float(df_xy_org.values[i-1]))
df_xy=pd.Series(df_xyvalue,index=df_xytime)
df_sytime=np.empty(len(df_sy_org)-1,dtype=datetime.datetime)
df_syvalue=np.empty(len(df_sy_org)-1)
for i in range(1,len(df_sy_org)-1):
df_sytime[i-1]=df_sy_org.index[i].to_pydatetime()
df_syvalue[i-1]=(float(df_sy_org.values[i])-float(df_sy_org.values[i-1]))
df_sy=pd.Series(df_syvalue,index=df_sytime)
第四步：定義距離函式
#定義距離計算方法 dist=lambda x,y:np.abs(np.power(100,x)-np.power(100,y))
第五步：下游斷面監測資料作為參照序列，對上游斷面監測資料中選取【參照序列長度正負2的長度】的序列，迴圈驗證最相似的時間序列
#dwt計算，比對dist生成的結果中，累計距離最小的項作為最終的輸出結果
lunci=20
bestBegin=datetime.datetime.strptime(beginTime, '%Y/%m/%d %H:%M')
bestEnd=datetime.datetime.strptime(endTime, '%Y/%m/%d %H:%M')
bestCost=sys.float_info.max
bestds='undefined'
#生成待比對斷面的監測資料，長度上邊可以比下游斷面的監測資料長度多或者少兩條記錄（-2，-1，0，1，2），要求首尾必須匹配方式
diff=np.array([-2,-1,0,1,2],dtype=np.int)
for lendiff in diff:
bestLunci=-1
for ii in range(lunci):
delta_begin=datetime.timedelta(hours=(-4)*(lunci-ii))
delte_end=datetime.timedelta(hours=(-4)*(lunci-ii-int(lendiff)))
try:
tempbegin=datetime.datetime.strptime(beginTime, '%Y/%m/%d %H:%M')+delta_begin
tempend=datetime.datetime.strptime(endTime, '%Y/%m/%d %H:%M')+delte_end
tempdf_sy=df_sy[tempbegin:tempend]
ds=dtw(df_xy,tempdf_sy,dist_method=dist,keep_internals=True,step_pattern='asymmetric',window_type='sakoechiba',open_begin=False,open_end=False,window_args= {"window_size": 1}) #引數傳遞方法
print('執行輪次：'+str(ii)+'時 ds.distance='+str(ds.distance))
if bestCost>ds.distance:
bestCost=ds.distance
bestds=ds
bestBegin=tempbegin
bestEnd=tempend
bestLunci=ii
bestds.plot(type="twoway")
print(bestdsT.distance)
print('執行輪次：'+str(ii)+'成功！')
except:
print('執行輪次：'+str(ii)+'出錯！')
第六步：找到原始斷面汙染物濃度值，並找出上下游不同時刻點的對應關係情況，輸出csv檔案
#把dtw執行結果存入csv檔案中
#建立檔案
file = xlwt.Workbook()
#寫入sheet名稱,2017-05-01T06_00_00 格式
name_sheet='錦江上下游水質關係'
table = file.add_sheet(name_sheet)
title=['df_xytime','df_sytime','df_xyvalue','df_syvalue','timediff']
bestdf_sy=df_sy[bestBegin:bestEnd]
#table.write(行_0開始, 列_0開始, 值)
#寫入表頭，第一行
for j in range(len(title)):
table.write(0, j, title[j])
for i in range(len(bestds.index1)):
#index_i+1 行，index_j列
for index_i,r in enumerate(bestds.index1):
df_xytime=datetime.datetime.strptime(beginTime, '%Y/%m/%d %H:%M')+datetime.timedelta(hours=4*1) +datetime.timedelta(hours=4*index_i) #beginTime先加一個小時的原因是，前期執行了一個一階差分，時間上往後延遲了一個週期，而bestBeginTime不存在這種問題
df_sytime=bestBegin+datetime.timedelta(hours=4*int(bestds.index2[index_i]))
df_xyvalue=df_xy_org[df_xytime]
df_syvalue=df_sy_org[df_sytime]
timediff=(df_xytime-df_sytime).days*24+(df_xytime-df_sytime).seconds/3600
table.write(index_i+1,0,df_xytime.strftime('%Y-%m-%d %H:%M:%S'))
table.write(index_i+1,1,df_sytime.strftime('%Y-%m-%d %H:%M:%S'))
table.write(index_i+1,2,df_xyvalue)
table.write(index_i+1,3,df_syvalue)
table.write(index_i+1,4,timediff)
#儲存檔案，2018-09-21T15_00_00.csv 名稱格式
name_file=name_sheet+'.csv'
file.save(name_file)
參考資料：
語音訊號處理之（一）動態時間規整（DTW）_zouxy09的專欄-CSDN部落格_動態時間規整
時間序列匹配之dtw的python實現（二）_輕舟已過萬重山的部落格-CSDN部落格_dtw python 用法