一種基於多重散射波幹涉的介質變化監測算法
2024-04-12 23:07:05
1.本發明涉及的是一種彈性波信號處理方法，具體地說是介質變化監測方法。
背景技術：
2.高靈敏度、高精度地監測介質變化能夠為火山預警、地震研究、材料檢測以及工程安全性評估提供技術支撐。在具有複雜幾何構型或非均質材料傳播介質中，彈性波能量會快速從直達波場轉換至多重散射波場。相較於直達波波場，多重散射波在其長程傳播過程中對介質內部的多次採樣，因此對介質物性變化敏感。通過對多重散射波對應聲信號進行幹涉處理，可以實現對介質變化的高靈敏度、高精度監測，因此被廣泛應用於地震研究以及超聲無損檢測等領域。
3.當前多重散射波幹涉分析方法的評估重心在於傳播介質彈性波速的微量變化，傳統分析方法包括windowed cross-correlation(wcc)方法以及stretching方法。其中，wcc方法算法實現簡單，但由於其理論假設與物理場景存在差異且算法本身對於周期跳躍(cycle skipping)缺乏抑制，導致對介質變化評估的精度與穩定性不足。而stretching方法雖然穩定性與評估精度有所提高，但是受限於線性時延假設以及對信號採集設備的高精度時間同步需求，導致其算法不適用於局部速度變化評估且應用成本增加。綜上，現有多重散射波幹涉算法計算精度差強人意，且在面對局部速度變化導致的非線性時延、低信噪比或時延過大導致的周期跳躍等情況時可靠性無法保證。
4.動態時間規整(dynamic time warping，dtw)方法是在語音識別領域發展起來的一種時域聲學非線性優化方法，可通過採用逐點匹配的思路實現對音頻信號的差異性補償。dtw方法具有非線性分析潛力且對周期跳躍具有一定抑制能力，但是受限於應用場景和傳播過程間的巨大差異無法直接應用於多重散射波幹涉分析。
5.綜上所述，需要針對多重散射波信號特點，圍繞波形畸變、低信噪比以及非線性時延場景下介質變化監測需求，吸收dtw方法逐點匹配的處理思路以構建新型基於多重散射波幹涉的介質變化監測算法。
技術實現要素：
6.本發明的目的在於提供實現具有病態匹配、周期跳躍抑制能力的高精度、高效能的一種基於多重散射波幹涉的介質變化監測算法。
7.本發明的目的是這樣實現的：
8.本發明一種基於多重散射波幹涉的介質變化監測算法，其特徵是：
9.(1)信號獲取；
10.(2)信號分段；
11.(3)獲取相似信息，基於互相關處理，獲得分段信號間水平偏移量αi，基於αi進行截斷處理，獲取相似信息；
12.(4)時延計算，基於歐式距離計算彎曲窗口內誤差矩陣d，由誤差矩陣d和約束條件
b計算累積誤差矩陣c，由累積誤差矩陣c及約束條件b確定最佳規整路徑p，再由最佳規整路徑p獲取平均偏移量βi，以求取時延τi；
13.(5)介質變化監測，對最佳規整路徑下兩道信號計算其去相關係數，基於步驟(4)的時延τi與流逝時間計算其速度變化，連續接收穫取介質變化前、後信號，由步驟(1)-(5)實現介質變化監測。
14.本發明還可以包括：
15.1、步驟(1)信號獲取具體為：
16.獲取介質變化前、後多重散射波時域信號，記為u0(t)和u1(t)。
17.2、步驟(2)信號分段具體為：
18.基於總時間窗口t
all
、移動時間窗的窗長t
mw
和移動步距t
step
獲得k個分段信號u
0，i
(t)和u
1，i
(t)，1≤i≤k。
19.3、步驟(3)獲取相似信息的具體方法為：
20.水平偏移量獲取：對分段信號u
0，i
(t)和u
1，i
(t)做互相關，計算互相關函數的最大值以及最大值對應的位置，此位置即為分段信號u
0，i
(t)和u
1，i
(t)間的水平偏移量αi；
21.截斷處理：由計算得到的偏移量αi，對齊分段信號，將分段信號u
0，i
(t)和u
1，i
(t)截斷為在該移動時間窗口內幾乎重疊的部分，截斷處理後分段信號u
0，i
(t)和u
1，i
(t)長度都為mi＝t
mw-|αi|。
22.4、步驟(4)時延計算具體為：
23.計算誤差矩陣：構建大小為mi×
mi的誤差矩陣d，矩陣內部元素全為正無窮大，再基於歐式距離公式重新計算彎曲窗口內的誤差矩陣d，d(j，l)＝(u
0，i
(j)-u
1，i
(l))2；
24.計算累積誤差矩陣：約束條件b將控制彎曲窗口內累積誤差矩陣c的計算，當b＝1時，c(j，l)＝d(j，l)+min{c(j-1，l)，c(j-1，l-1)，c(j，l-1)}；當b≥2時，累積誤差矩陣c使用下式計算；
[0025][0026]
尋找最佳規整路徑：在累積誤差矩陣c中，根據約束條件值確定移動路徑，找到最佳規整路徑p＝[p1，p2，...pk]，mi≤k≤2m
i-1；從pk＝(mi，mi)開始，根據c(m
i-b，m
i-b+1)、c(m
i-1，m
i-1)和c(m
i-b+1，m
i-b)三者最小值所在的行列值即為下一點p
k-a
＝(行，列)的值，a＝1或a＝b由路徑確定，再根據p
k-a
找到下一點p
k-2a
，重複此操作，直到找到p1(1，1)；
[0027]
時延估計：對規整路徑p中每一點pc＝(ic，jc)進行(i
c-jc)運算，其中1≤c≤k，對結果取平均，即為平均偏移量βi，再加上步驟(3)計算得到的偏移量αi，即為求取的該移動時間窗口下的時延τi。
[0028]
5、步驟(5)介質變化監測具體為：
[0029]
去相關係數計算：基於u
0，i
(ic)與u
0，i
(jc)，1≤c≤k計算信號間互相關係數c(i)，去相關係數di＝1-ci；
[0030]
速度變化計算：基於步驟(3)和步驟(4)計算獲取所有分段信號u
0，i
(t)和u
1，i
(t)的時延τi，1≤i≤k和移動步距計算速度變化δvi，即δvi＝(τ
i+1-τi)/t
step
，1≤i≤k-1。
[0031]
本發明的優勢在於：
[0032]
1、與stretching方法相比可以實現非線性時延以及局部介質變化的評估，且不需要信號之間嚴格時間同步；
[0033]
2、與wcc方法相比，由於採用逐點匹配，因此時間解析度大幅度提高，而且有效地抑制了周期跳躍；
[0034]
3、相比於傳統dtw方法，針對多重散射波信號的特點，對信號進行分段處理以及使用一種新的彎曲窗口，截斷處理後彎曲窗口的寬度可成倍減小，大大提高計算效率。基于波形特性計算其水平偏移量，對信號截斷處理，有效抑制病態匹配的發生，同時，合理選擇約束條件值，在不降低計算效率的同時有效提高計算精度，減弱波形畸變帶來的影響，大大提高計算精度。
附圖說明
[0035]
圖1為本發明的流程圖；
[0036]
圖2為本發明的彎曲窗口示意圖；
[0037]
圖3為本發明的最佳規整路徑的計算圖；
[0038]
圖4為本發明的仿真信號圖；
[0039]
圖5為截斷處理前的分段信號圖；
[0040]
圖6為截斷處理後的分段信號圖；
[0041]
圖7為本發明的非線性時延量級10-2
時四種方法估計結果圖
[0042]
圖8為本發明的非線性時延量級10-2
時四種方法的cdf圖。
具體實施方式
[0043]
下面結合附圖舉例對本發明做更詳細地描述：
[0044]
結合圖1-8，本發明具體包括以下步驟：
[0045]
步驟1信號獲取：本發明採用仿真信號u0(t)，對仿真信號進行線性或非線性拉伸處理，構建介質變化後的信號u1(t)；
[0046]
步驟2信號分段：對信號u0(t)和u1(t)進行分段處理，一方面是多重散射波屬於長持時信號，且若介質變化十分微弱時，常採取插值處理，將導致信號長度過長，嚴重降低計算效率；另一方面減小病態匹配的發生概率，分段處理後信號間允許的最大偏移量限制在分段信號的長度內。首先確定總時間窗口t
all
，確保總時間窗口內多重散射波信號佔主導，且信噪比較高。設置移動時間窗的窗長t
mw
，移動時間窗指具有指定長度t
mw
的時間窗口，包含在總時間窗口內，按照設定的移動步距t
step
沿著時間軸滑動該移動時間窗，實現信號分段。需要注意的是移動時間窗口內需要包含足夠多的信息量，即包含的波形要足夠複雜，但移動時間窗的窗長過大時，計算速度將會降低，因此，需要合理的設置窗長。同時，為了保證分析結果的連續性，移動步距需要較小，使移動時間窗口內信號有重疊。假定總時間窗口t
all
＝[t
1-tw，tk+tw]，移動時間窗的窗長為t
mw
＝2tw，其中滑動移動時間窗口後該窗口的中心位置為ti，1≤i≤k，其中k為分段信號個數，本發明使移動時間窗內包含五個周期的波形
變化，移動步距設置為窗長的1/15，即t
step
＝t
i-t
i-1
＝1/15
×
t
mw
，獲取分段信號u
0，i
(t)和u
1，i
(t)，1≤i≤k；
[0047]
步驟3獲取相似信息：
[0048]
3.1：水平偏移量獲取。基於互相關處理，考慮其波形特性，計算分段信號間的水平偏移量。假定時延變量ts在[-2tw，2tw]範圍內，求取第i段分段信號u
0，i
(t)和u
1，i
(t)，t∈[t
i-tw，ti+tw]的互相關係數r(ts)，公式如下所示，互相關係數最大值所對應的時延變量ts，即為分段信號u
0，i
(t)和u
1，i
(t)的水平偏移量αi，其中k指分段信號個數；
[0049][0050]
當αi≥t
mw
/5，i≥1時，考慮整段信號，令分段個數k為1，求取分段信號u
0，i
(t)和u
1，i
(t)的時延αi，1≤i≤k。
[0051]
3.2：截斷處理。基於步驟3.1獲取的分段信號u
0，i
(t)和u
1，i
(t)間的水平偏移量αi對分段信號進行截斷處理，使該時間窗口內分段信號包含的傳播介質微結構信息基本一致，減弱dtw邊界條件確定性帶來的誤差，有效抑制病態匹配的發生。若時延αi≥0，則u
0，i
(t)的第i的時間窗口調整為[t
i-tw，ti+t
w-αi]，i≥1，u
1，i
(t)的第i段的時間窗口調整為[t
i-tw+αi，ti+tw]；反之，若時延αi＜0則u
0，i
(t)的第i段的時間窗口調整為[t
i-tw+|αi|，ti+tw]，u
1，i
(t)的第i段的時間窗口調整為[t
i-t
wi
，ti+t
w-|αi|]。截斷處理後，分段信號u
0，i
(t)和u
1，i
(t)的長度都為mi＝t
mw-|αi|；
[0052]
步驟4時延計算：
[0053]
4.1：計算誤差矩陣計算。為了減小計算量以提高計算效率，對截斷處理後的第i段分段信號間只計算彎曲窗口內的誤差矩陣。本發明中提出一種新的彎曲窗口，如附圖2所示，結合sakoe-chiba band高效率與itakura-paralleogram高精度的優勢，序列對首尾使用itakura-paralleogram，中間部分使用sakoe-chiba band，既能進一步提高計算效率，又保障結果的可靠性。圖2中表示向下取整，即不超過x的最大整數。首先定義幾個參數，k2＝2k1+1，整個窗口按照縱軸方向灰色框的變化範圍，將橫坐標劃分為7個區域，分別是(1)1≤x≤k1+1；
ꢀꢀ
(2)k1+2≤x≤k2；
ꢀꢀ
(3)k2+1≤x≤k2+k3；
ꢀꢀ
(4)k2+k3+1≤x≤m
i-(k2+k3)；
ꢀꢀ
(5)m
i-(k2+k3)+1≤x≤m
i-k2；
ꢀꢀ
(6)m
i-k2+1≤x≤m
i-k
1-1；
ꢀꢀ
(7)m
i-k1≤x≤mi，縱坐標灰色框的範圍按照劃分區域根據圖中公式進行計算。舉例說明，對2000
×
2000長的序列做100次蒙特卡洛計算總耗時，再取均值，設置彎曲窗口的寬度為100，則傳統dtw算法耗時為132.3ms，sakoe-chuba band耗時為55.8ms，itakura-paralleogram耗時為69.7ms，而本方法耗時最短為54.4ms，相比於傳統dtw算法改進率為58.9％。
[0054]
首先構建mi×
mi，1≤i≤k的單位矩陣，然後單位矩陣乘以正無窮大，使矩陣內部元素全為正無窮大，記為誤差矩陣d。通常彎曲窗口的寬度設置為分段信號長度的10％，但通過截斷處理後，彎曲窗口的寬度可成倍的減小，越小則計算效率越高，但太小也會導致無法正確找到最佳的規整路徑，因此，可以根據實際需求自行調整。本發明設置彎曲窗口的寬度n為分段信號長度的5％，即n＝mi×
0.05。對彎曲窗口內的誤差矩陣d重新計算，使用歐式距離計算誤差矩陣d，d(j，l)＝(u
0，i
(j)-u
1，i
(l))2，指u
0，i
(t)第j個元素與u
1，i
(t)的第l個元素
差值的平方即為誤差矩陣d第j行第l列的值；
[0055]
4.2：計算累積誤差矩陣。約束條件(slope constraint)用於限制規整路徑的變化，約束條件值越大，則匹配路徑越平滑，可減弱波形畸變和孤立點幹擾帶來的影響，但計算複雜度也將增加。構建一個和誤差矩陣d同樣大小的矩陣，記為累積誤差矩陣c，矩陣內部元素全為正無窮大。根據設置的約束條件b值和步驟4.1計算的誤差矩陣d重新計算彎曲窗口內的累積誤差矩陣c，本發明中設置約束條件為2，計算公式如下所示；
[0056]c1，1
＝d
1，1
[0057]ci，j
＝d
i，j
+min{c
i-1，j-1
，c
i-1，j
，c
i，j-1
}i，j≤n
[0058][0059]
4.3：尋找最佳規整路徑。在累積誤差矩陣c中，找到最佳規整路徑p＝[p1，p2，...pk]，mi≤k≤2m
i-1，路徑p表示u
0，i
(t)和u
1，i
(t)中各點的匹配關係。計算規整路徑時需要滿足以下三個條件：
[0060]
①
單調性：如果pk＝(i，j)，p
k+1
＝(i
′
，j
′
)，則有i
′
≥i，j
′
≥j，目的是保障路徑上的點必須隨著時間單調進行；
[0061]
②
連續性：如果pk＝(i，j)，p
k+1
＝(i
′
，j
′
)，則有i
′
≥i+1，j
′
≥j+1，目的是限制路徑上的點只能為相鄰點，不允許有跳躍；
[0062]
③
邊界條件確定性：p1＝(1，1)，pk＝(mi，mi)；
[0063]
根據約束條件b的設定值，從pk＝(mi，mi)逆著找到p1＝(1，1)，已知pk＝(mi，mi)，確定c(m
i-b，m
i-b+1)、c(m
i-1，m
i-1)和c(m
i-b+1，m
i-b)三者最小值所在的行列位置，若位置x＝(m
i-1
，m
i-1
)，即向左下角移動一個單位，則令a＝1，否則令a＝b。於是p
k-a
＝(x，y)，x為最小值所在的行，y為最小值所在的列。例如附圖3中，帶入b＝2，則c(m
i-2，m
i-1)、c(m
i-1，m
i-1)和c(m
i-1，m
i-2)三者中最小值為c(m
i-2，m
i-1)，不是沿著左下角移動一個單位，因此a＝b＝2，即p
k-2
＝(m
i-2，m
i-1)。根據約束條件b值和p
k-2
的坐標位置可以確定p
k-1
＝(m
i-1，mi)。以此類推，直到找到p1＝(1，1)；
[0064]
4.4：時延估計。最佳規整路徑p中pc＝(ic，jc)表示的是u
0，i
(ic)與u
1，i
(jc)匹配，因此，(i
c-jc)表示第i段分段信號u
0，i
(t)和u
1，i
(t)在時間點ic時的偏移量。對最佳規整路徑中每一點pj，1≤j≤k重複上述操作，獲取分段信號u
1，i
(t)相對於分段信號u
0，i
(t)在所有時間點上的偏移量，結果取平均求取平均偏移量βi，即βi＝[(i
1-j1)+，...，+(i
k-jk)]/k，再加上步驟3計算得到的偏移量αi，即為求取的該移動時間窗口下的時延τi，τi＝βi+αi；
[0065]
步驟5介質變化監測：
[0066]
5.1：去相關係數計算；將時間點ic、jc分別代入分段信號u
0，i
(t)和u
1，i
(t)，即u
0，i
(ic)與u
0，i
(jc)，其中1≤c≤k。使用下式求取信號u
0，i
(ic)與u
0，i
(jc)的互相關係數ci，其中μ
0，i
和σ
0，i
表示信號u
0，i
(ic)的均值和方差；而μ
1，i
和σ
1，i
表示信號u
1，i
(jc)的均值和方差；根據di＝1-ci求取去相關係數di；
[0067][0068]
5.2：速度變化計算：基於步驟3和4計算獲取所有分段信號u
0，i
(t)和u
1，i
(t)的時延τi且1≤i≤k和移動步距計算速度變化δvi，即δvi＝(τ
i+1-τi)/t
step
，1≤i≤k-1；
[0069]
下面結合具體的仿真實例做進一步說明：
[0070]
圖4中黑色虛線和灰色實線分別表示介質變化前、後信號u0(t)和u1(t)，為方便計算，以時間點數表示時間，其轉換關係為以秒為單位的時間乘以時間採樣率。設置總時間窗口t
all
＝[501，2500]，移動時間窗的窗長t
mw
＝150，移動步距t
step
＝10，則分段信號個數k為184，即獲取分段信號u
0，1
(t)和u
1，1
(t)、u
0，2
(t)和u
1，2
(t)，...，u
0，184
(t)和u
1，184
(t)。
[0071]
下面以第80段分段信號為例詳細介紹步驟3和步驟4，如圖5所示，此時u
0，80
(t)＝u0(1291：1440)和u
1，80
(t)＝u1(1291：1440)，由步驟3.1計算其水平偏移量α
80
＝-2，由於α
80
＜0，u
0，80
(t)截斷為u
0，80
(t)＝u0(1291：1438)，u
1，80
(t)截斷為u
1，80
(t)＝u1(1292：1440)，截斷後分段信號如圖6所示，相較於圖5能夠發現截斷後分段信號間高度相似，分段信號間水平偏移量基本得到補償，後續的誤差矩陣和累積誤差矩陣計算中彎曲窗口的設置可成倍減小，將大大提高計算效率。截斷處理後信號長度m
80
＝t
mw-|α
80
|＝148。由於分段信號間時延很小，選擇對分段信號進行插值處理，插值間隔為0.1，插值處理後分段信號u
0，80
(t)和u
1，80
(t)的長度為m
80
＝1471。由m
80
確定彎曲窗口的長度n，向下取整，即構建m
80
×m80
大小的正無窮大誤差矩陣d，使用歐式距離公式重新計算彎曲窗口如附圖2中灰色框內的誤差矩陣，例如構建和誤差矩陣d同樣大小的內部元素值為正無窮大的矩陣，記為累積誤差矩陣c，設置坡度約束b為2，基於誤差矩陣d和坡度約束b，根據步驟3.2的公式重新計算彎曲窗口內的累積誤差矩陣c，例如由k1＝c
99，148
+2
×d100，149
+d
100，150
＝0.03447，廳2＝c
98，149
+2
×d99，150
+d
10，0150
＝0.03896以及k3＝c
99，149
+2
×d10，0150
＝0.03672三者的最小值確定c(100，150)的值，即c(100，150)＝0.03447，同理，計算彎曲窗口內的其餘累積誤差矩陣c的元素。由累積誤差矩陣c和坡度約束確定最佳規整路徑p，已知pk＝(m
80
，m
80
)和p1＝(1，1)，其中k的值不需要先知道，可以認為是一個具有固定值的未知數，當最佳規整路徑確定後，才能確定k的大小。根據c(1469，1470)＝0.008347、c(1470，1470)＝0.008332以及c(1470，1469)＝0.008261的最小值所在的行列位置確定路徑的位置，由於c(1470，1469)最小且不在對稱位置，即最小值不是c(1470，1470)，因此a＝b＝2。於是p
k-a
＝p
k-2
＝(1470，1469)，根據坡度約束值b和p
k-2
的坐標位置確定其路徑方向是先向下移動一次，再沿對稱方向移動，即確定p
k-1
＝(1471，1470)。以此類推，直至找到p1＝(1，1)。隨後，得到k＝1713，得到p＝[p1，...，pc，...，p
1713
]，其中pc＝(ic，jc)，1≤c≤k且1≤ic≤m
80
。舉例，對p
k-2
＝(1470，1469)作行坐標減去縱坐標處理，獲得偏移量y
k-2
＝(1470-1469)＝1。同理，對規整路徑p中每一點pj，1≤j≤k做同樣處理，得到偏移量y1，y2，...，yk，最後，對該y1，y2，...，yk做均值處理，求取平均偏移量β
80
＝(y1+y2+，...，+yk)/k＝-0.4776，與步驟3計算得到的水平偏移量α
80
＝-2相加，求得第80段分段信號間的時延τ
80
，τ
80
＝β
80
+α
80
＝-2.4776。同理，按照同樣的方式處理分段信號u
0，1
(t)和u
1，1
(t)、u
0，2
(t)和u
1，2
(t)，...，u
0，184
(t)和u
1，184
(t)，獲取分段信號間時延τ1，τ2，...，τ
184
。
[0072]
使用步驟5.1中公式計算最佳規整路徑下的信號sig1和sig2的互相關係數c
80
＝0.9767，其中sig1(c)＝u
0，80
(ic)，sig2(c)＝u
1，80
(jc)，即去相關係數d
80
＝0.0233。同理，獲取所有分段信號下的去相關係數d1，...，d
184
。由公式δvi＝(τ
i+1-τi)/t
step
，1≤i≤k-1，計算速度變化。
[0073]
最後，為分析本發明的基於多重散射波的介質變化監測算法的性能表現，分別在時延線性變化與非線性變化兩種場景下，考慮在不同時延量級下，本發明提及的方法與wcc方法、stretching方法以及傳統dtw算法的性能差異。圖6為非線性時延量級10-2
時四種方法估計結果圖，能夠看出只有本方法能準確可靠地估計出該量級下的非線性時延，wcc若進行擬合，能夠大致反映出非線性時延的變化趨勢，但無法準確地估計出真實值，而傳統dtw方法的結果趨於零，說明傳統dtw已經達到極限，出現解析度不足的現象，而stretching方法不能估計出非線性時延變化，其結果為線性變化趨勢，結合圖7能進一步直觀地看出本方法的誤差最小，計算精度最高。圖7為概率分布函數cdf圖，橫坐標為估計誤差，縱坐標是累積分布函數，表示在某個估計誤差下出現的概率，橫坐標越小表示估計誤差越小，其計算精度越高。同時，結合表1線性與非線性時延不同量級下四種方法的均方根誤差rmse，rmse越小，說明方法的性能越好，能夠發現本發明提及的方法對於線性時延的估計性能稍弱於stretching方法，而在非線性時延估計時，性能最好，能夠實現高精度地估計時延變化，且無論是線性時延還是非線性時延，本方法的rmse低於傳統dtw算法的一個量級，即性能優於傳統dtw方法。總的來說，本文提及的基於多重散射波的介質變化監測算法不僅適用於介質線性變化，同時也適用於介質非線性變化，能更加全面、準確地描述介質變化。
[0074]
表1線性與非線性四種方法的rmse結果
[0075][0076]