基于壓縮感知流型辨識的變化靈敏度場電容層析成像方法
1.本發明屬于圖像處理技術領域，具體涉及一種基于壓縮感知流型辨識的變化靈敏度場電容層析成像方法。
背景技術：
2.在工業和生活領域中，兩相流現象是普遍存在的，故過程層析成像技術在工業測量領域發展前景非常廣闊。兩相流系統如今廣泛應用于工業中以適應現代化工業的需要。電容層析成像技術(electrical capacitance tomography，ect)是一種成熟的過程層析成像技術。ect技術是利用被測管道周圍的傳感器陣列獲取被測物質的投影數據(電容值)，并用圖像重建算法重建出被測區域的介電常數分布。ect技術相對于x射線層析成像、伽馬射線層析成像、超聲波層析成像等其他過程層析成像技術，具有硬件成本低、響應速度快、硬件結構簡單、安全性能好(非輻射)等優點，其應用前景十分廣闊。
3.ect技術的基本思想是：不同的流動物質具有不同的相對介電常數，當各流相物流動不均或濃度發生改變時，都將致使流體相對介電常數的變化，從而引起檢測極板之間電容響應的改變。檢測中采用陣列式的電容傳感器，各電容傳感器之間彼此組合可以產生多個獨立測量值。最后通過成像計算機即可重建出兩相流體的介電常數分布情況。典型的ect系統結構圖如圖1所示。
4.被測管道內部的流體介電常數分布ε(x,y)與從傳感器陣列獲取的一組電容向量ci之間的關系可以描述為：
5.ci＝∫∫dsi(x,y)ε(x,y)dxdy
???
(1)
6.其中，d是被測管道區域，(x,y)∈d表示像素點的位置；ε(x,y)為介電常數分布函數，ci(i＝1,2...m)為傳感器陣列測量所得的m個電容測量值。si(x,y)是靈敏度的分布函數。為了快速計算靈敏度場，電勢分布法是一個很好的方法。根據電勢分布法，靈敏度場可以描述為以下公式：
[0007][0008]
其中v是激勵電壓，和分別是電極i和電極j為激勵電極時的電勢分布。
[0009]
ect的圖像重構是利用測量獲得的m個獨立電容值來反演管道內介電常數分布。目前已經有多種圖像重建方法，但是由于ect靈敏度的“軟場”問題，并且觀測數據極為有限，是一個欠定問題，因此很難重建出高質量的圖像，所以圖像重建算法一直是人們研究的熱點。目前在ect系統中有兩種典型的圖像重建算法：線性反投影法(lbp)和landweber迭代法。lbp算法重建的圖像質量相對較低，但lbp算法的優點是速度快；landweber迭代法算法是一種成像精度相對較高的迭代算法，它使用最速下降法來進行收斂迭代，具有更好的圖像質量。這兩種傳統算法都是利用固定靈敏度場來求解正問題的，固定靈敏度場是假定被測電容是與整個管道內部變化沒有關系的，這顯然是不符合實際的假設。因此，為了減少
ect靈敏度場的“軟場”影響，這里提出一種基于壓縮感知流型辨識的變化靈敏度場電容層析成像方法。迄今為止，提高電容層析圖像重建的質量和精度仍然是個關鍵技術問題。
技術實現要素：
[0010]
為了克服已有技術的不足，針對電容層析成像問題，提出一種基于壓縮感知流型辨識的變化靈敏度場電容層析成像方法。首先本發明方法設定k類典型流型，將其歸一化后的ect測量電容信號作為標準樣本，用過調整介質分布和疊加隨機噪聲的方式來生成訓練樣本集。其次將隨機高斯矩陣作為測量矩陣，分別求得測試樣本和標準樣本的觀測值，并用壓縮感知信號重構算法如正交匹配追蹤算法(omp)分別求出測試樣本和標準樣本在過完備訓練樣本集e投影下的稀疏系數。最后計算測試樣本與所有標準樣本的稀疏解之間的線性相關系數，并選擇最大相關系數所對應的流型為測試樣本的歸屬流型，并用該歸屬流型的靈敏度場替代原始的固定靈敏度場進行圖像重構；提高了電容層析成像的精度。
[0011]
本發明解決其技術問題所采用的技術方案是：
[0012]
一種基于壓縮感知流型辨識的變化靈敏度場電容層析成像方法，包括以下步驟：
[0013]
步驟一，初始化，歸一化數據，被測管道內部的流體介電常數分布ε(x,y)與從傳感器陣列獲取的一組電容向量ci之間的關系可以描述為：
[0014]ci
＝∫∫dsi(x,y)ε(x,y)dxdy
???
(1)
[0015]
其中，d是被測管道區域，(x,y)∈d表示像素點的位置；ε(x,y)為介電常數分布函數，ci(i＝1,2...m)為傳感器陣列測量所得的m個電容測量值，si(x,y)是靈敏度的分布函數；
[0016]
根據電勢分布法，靈敏度場描述為以下公式：
[0017][0018]
其中，v是激勵電壓，和分別是電極i和電極j為激勵電極時的電勢分布；
[0019]
由公式(1)看出，當介電常數相近時，電容向量ci與介電常數分布函數ε(x,y)之間的關系可近似為線性；為了方便計算，對其進行離散和歸一化處理，得到ect矩陣形式的數學模型如下所示：
[0020]
c＝sg
???
(3)
[0021]
其中，c是歸一化的m
×
1的電容測量值向量，g是大小為n
×
1的灰度向量，代表管道不同介電常數的分布；s是歸一化的大小為m
×
n的靈敏度矩陣，m為測量電容值的個數。n的大小代表了成像區域像素個數，從公式(3)可以看出，ect圖像重建的問題可以簡化為通過測量電容向量c求解灰度向量g的過程；
[0022]
步驟二，生成訓練樣本集。給定k類流型的測量電容數據組成訓練樣本集，記第i類流型的訓練集ei共有ni個樣本，訓練集ei表示為：
[0023][0024]
其中，vj∈rm×1(j＝1,2
…
ni)代表第i類流型里面的第j個訓練樣本，即ect系統從第i類流型獲取的測量電容向量；m是測量電容值得個數，8極板ect系統有m＝28，通過調整介質分布和疊加隨機噪聲的方式來生成訓練樣本集，每類流型的訓練樣本數為100個，則所
有k類流型組成的訓練樣本集為：
[0025]
e＝[e1,e2…ek
]∈rm×n，
???
(5)
[0026]
其中n＝n1+n2+
…
nk；
[0027]
對于流型i，若訓練集ei完備或者過完備，則此類流型的任意測量電容向量λ可由其訓練集ei線性表示：
[0028][0029]
其中，系數αi(j＝1,2
…
ni)為實數，為系數向量，在實際應用中，由于無法預知任意測量電容向量λ的歸屬類，可以考慮將其重新表示為k類流型組成的訓練樣本集e的線性組合：
[0030]
λ＝ex0.
???
(7)
[0031]
式中，為稀疏向量，至此，由ect系統獲得的測量電容信號λ就可由以訓練樣本集組成的過完備字典變換基e來稀疏表示，這樣就使得ect測量電容信號滿足了壓縮感知理論的前提條件；
[0032]
步驟三，計算測試樣本和標準樣本的觀測值，隨機高斯矩陣b∈rd×m作為測量矩陣，則根據壓縮感知思想，針對以測量電容向量λ作為原始信號的線性測量過程為：
[0033][0034]
其中，為觀測值，是d
×
1的向量，標準樣本的觀測值按照同樣方法計算；
[0035]
步驟四，計算測試樣本和標準樣本的稀疏系數，由于觀測向量的維數d遠小于稀疏信號x0的維數，直接求解公式(8)是一個欠定問題，將其轉換為如公式(9)所示的數學模型：
[0036][0037]
通過壓縮感知的信號重構算法求解此l0范數下的最優化問題，使用正交匹配追蹤(omp)算法得到測量電容向量λ和標準電容向量λk對應于過完備字典變換基e的稀疏解x0和xk；
[0038]
步驟五，計算測試樣本和標準樣本稀疏解的線性相關系數，通過x0與xi的線性相關程度來確定測試樣本的歸屬流型，達到流型辨識的目的；
[0039]
相關系數則是衡量變量之間線性相關程度的指標，其計算方法如式(10)所示：
[0040][0041]
其中，和分別為x0和xi的平均值。根據公式(10)分別計算測試樣本與所有標準樣本對應的稀疏解之間的線性相關系數；
[0042]
步驟六，確定流型，選擇相關系數最大的歸屬流型靈敏度場，過程如下：將最大相關系數所對應的流型確定為測試樣本的歸屬流型。選擇該歸屬流型的靈敏度場確定為新的靈敏度矩陣；
[0043]
步驟七，圖像重構，輸出灰度向量g。
[0044]
進一步，所述步驟七中，用變化靈敏度矩陣和重構算法來進行圖像重構。
[0045]
再進一步，所述步驟三中，使用隨機高斯矩陣b∈rd×m作為測量矩陣，為了準確的重構出原始信號，測量矩陣必須滿足rip準則，rip準則的等價條件是測量矩陣和稀疏基不相關，當以d
×
m維的隨機高斯矩陣作為測量矩陣時，只要滿足條件：d≥cklog(m/k)，其中c是一個很小的數，k為原始信號的稀疏度，此測量矩陣在極大概率下滿足rip準則。測量矩陣可以替換，例如離散fourier矩陣和伯努利隨機矩陣。
[0046]
更進一步，所述步驟六中，為了確定流型，選擇相關系數最大的歸屬流型靈敏度場。可以用其他性能指標來選擇流型，例如最小相對誤差和最小相對電容殘差等。
[0047]
本發明的有益效果主要表現在：提高了電容層析成像的精度。
附圖說明
[0048]
圖1是ect系統結構圖。
[0049]
圖2是基于壓縮感知流型辨識的變化靈敏度場電容層析成像方法流程圖。
[0050]
圖3是七種流型的兩種算法成像效果圖。
具體實施方式
[0051]
下面結合附圖對本發明作進一步描述。
[0052]
參照圖2和圖3，一種基于壓縮感知流型辨識的變化靈敏度場電容層析成像方法，包括以下步驟：
[0053]
步驟一，初始化，歸一化數據，被測管道內部的流體介電常數分布ε(x,y)與從傳感器陣列獲取的一組電容向量ci之間的關系可以描述為：
[0054]ci
＝∫∫dsi(x,y)ε(x,y)dxdy
???
(1)
[0055]
其中，d是被測管道區域，(x,y)∈d表示像素點的位置；ε(x,y)為介電常數分布函數，ci(i＝1,2...m)為傳感器陣列測量所得的m個電容測量值，si(x,y)是靈敏度的分布函數；
[0056]
根據電勢分布法，靈敏度場描述為以下公式：
[0057][0058]
其中，v是激勵電壓，和分別是電極i和電極j為激勵電極時的電勢分布；
[0059]
由公式(1)看出，當介電常數相近時，電容向量ci與介電常數分布函數ε(x,y)之間的關系可近似為線性；為了方便計算，對其進行離散和歸一化處理，得到ect矩陣形式的數學模型如下所示：
[0060]
c＝sg
???
(3)
[0061]
其中，c是歸一化的m
×
1的電容測量值向量，g是大小為n
×
1的灰度向量，代表管道不同介電常數的分布；s是歸一化的大小為m
×
n的靈敏度矩陣，m為測量電容值的個數。n的大小代表了成像區域像素個數，從公式(3)可以看出，ect圖像重建的問題可以簡化為通過測量電容向量c求解灰度向量g的過程；
[0062]
步驟二，生成訓練樣本集。給定k類流型的測量電容數據組成訓練樣本集，記第i類
流型的訓練集ei共有ni個樣本，訓練集ei表示為：
[0063][0064]
其中，vj∈rm×1(j＝1,2
…
ni)代表第i類流型里面的第j個訓練樣本，即ect系統從第i類流型獲取的測量電容向量；m是測量電容值得個數，8極板ect系統有m＝28，通過調整介質分布和疊加隨機噪聲的方式來生成訓練樣本集，每類流型的訓練樣本數為100個，則所有k類流型組成的訓練樣本集為：
[0065]
e＝[e1,e2…ek
]∈rm×n，
???
(5)
[0066]
其中n＝n1+n2+
…
nk；
[0067]
對于流型i，若訓練集ei完備或者過完備，則此類流型的任意測量電容向量λ可由其訓練集ei線性表示：
[0068][0069]
其中，系數αi(j＝1,2ani)為實數，為系數向量，在實際應用中，由于無法預知任意測量電容向量λ的歸屬類，可以考慮將其重新表示為k類流型組成的訓練樣本集e的線性組合：
[0070]
λ＝ex0.
???
(7)
[0071]
式中，為稀疏向量，至此，由ect系統獲得的測量電容信號λ就可由以訓練樣本集組成的過完備字典變換基e來稀疏表示，這樣就使得ect測量電容信號滿足了壓縮感知理論的前提條件；
[0072]
步驟三，計算測試樣本和標準樣本的觀測值，隨機高斯矩陣b∈rd×m作為測量矩陣，則根據壓縮感知思想，針對以測量電容向量λ作為原始信號的線性測量過程為：
[0073][0074]
其中，為觀測值，是d
×
1的向量。標準樣本的觀測值按照同樣方法計算；
[0075]
步驟四，計算測試樣本和標準樣本的稀疏系數，由于觀測向量的維數d遠小于稀疏信號x0的維數，直接求解公式(8)是一個欠定問題，將其轉換為如公式(9)所示的數學模型：
[0076][0077]
通過壓縮感知的信號重構算法求解此l0范數下的最優化問題，使用正交匹配追蹤(omp)算法得到測量電容向量λ和標準電容向量λk對應于過完備字典變換基e的稀疏解x0和xk；
[0078]
步驟五，計算測試樣本和標準樣本稀疏解的線性相關系數，通過x0與xi的線性相關程度來確定測試樣本的歸屬流型，達到流型辨識的目的；
[0079]
相關系數則是衡量變量之間線性相關程度的指標，其計算方法如式(10)所示：
[0080][0081]
其中，和分別為x0和xi的平均值。根據公式(10)分別計算測試樣本與所有標準
樣本對應的稀疏解之間的線性相關系數；
[0082]
步驟六，確定流型，選擇相關系數最大的歸屬流型靈敏度場，過程如下：將最大相關系數所對應的流型確定為測試樣本的歸屬流型。選擇該歸屬流型的靈敏度場確定為新的靈敏度矩陣；
[0083]
步驟七，圖像重構，輸出灰度向量g。用變化靈敏度矩陣和重構算法來進行圖像重構。
[0084]
為了驗證本發明方法的性能，進行了仿真模擬實驗驗證。利用comsol multiphysics軟件建立不同流型仿真模型并生成仿真數據m
×
1的電容向量c和m
×
n的靈敏度矩陣s并導出數據，8種不同流型包括1/3層流、1/2層流、2/3層流、滿管、環流、核心流、泡狀流、空管。通過調整介質分布和疊加隨機噪聲的方式來生成訓練樣本集，每類流型的訓練樣本數為100個。用壓縮感知求稀疏解的方法，求出測試樣本和標準樣本的稀疏系數，最后通過計算測試樣本和不同標準樣本之間的線性相關系數，出最大的相關系數的歸屬流型，用該流型的靈敏度場替換原來的固定靈敏度場。用matlab r2019a仿真工具對數據進行處理，用圖像重構算法重建ect的圖像。三維仿真模型建立如下：管道為圓形，8個測量極板組成陣列，兩相流管道的截面被分割為40
×
40個像素，較密介質的介電常數為3，較稀介質的介電常數為1。成像算法選擇lbp算法和landweber迭代算法，信噪比為50db。分別使用變化靈敏度場和固定靈敏度場進行試驗。
[0085]
測量的8種流型分別為1/2層流、1/3層流、2/3層流、核心流、環流、泡狀流、滿管流和空管流，其中空管流的變化靈敏度場就是固定靈敏度場。兩種算法分別使用變化靈敏度場和固定靈敏度場的成像結果如圖x所示。從圖像上可以看出，相比于使用固定靈敏度場成像，使用本發明方法重建的圖像會有一個清晰的輪廓邊緣，例如，層流和環狀流。對于泡狀流來說，使用固定靈敏度場不能很好的分辨出泡狀，但本發明方法中lbp算法可以呈現出泡狀的輪廓，landweber迭代算法可以清晰的顯示出泡狀。故本發明方法重建的圖像從視覺效果上優于傳統的固定靈敏度場成像方法。
[0086]
為了對本發明方法圖像重建的性能進行數值上的評價，分析重建的圖像與原圖像的圖像相對誤差和相關系數。圖像相對誤差可定義為重建的圖像灰度向量與原始圖像灰度向量的差異程度，圖像相對誤差的公式如下:
[0087][0088]
其中，g是原始灰度向量。表示g的估計，表示為重建后的灰度向量。
[0089]
圖像的相關系數可定義為重建后的灰度向量與原始灰度向量之間的線性相關度，圖像的相關系數的公式如下:
[0090][0091]
其中，是向量g的平均值，而是向量的平均值。
[0092]
表1給出了本發明方法和傳統的固定靈敏度場成像方法在七種流型下的圖像相對誤差和相關系數。
[0093][0094][0095]
表1
[0096]
由圖3可得，本發明方法重建的圖像誤差比傳統固定靈敏度場成像方法小。這說明本發明方法重建的圖像好于傳統固定靈敏度場成像方法。本發明方法效果良好，獲得了較好的圖像。
[0097]
本說明書的實施例所述的內容僅僅是對發明構思的實現形式的列舉，僅作說明用途。本發明的保護范圍不應當被視為僅限于本實施例所陳述的具體形式，本發明的保護范圍也及于本領域的普通技術人員根據本發明構思所能想到的等同技術手段。
技術特征：
1.一種基于壓縮感知流型辨識的變化靈敏度場電容層析成像方法，其特征在于，所述方法包括以下步驟：步驟一，初始化，歸一化數據，被測管道內部的流體介電常數分布ε(x,y)與從傳感器陣列獲取的一組電容向量c
i
之間的關系可以描述為：c
i
＝∫∫
d
s
i
(x,y)ε(x,y)dxdy
????????????
(1)其中，d是被測管道區域，(x,y)∈d表示像素點的位置；ε(x,y)為介電常數分布函數，c
i
(i＝1,2...m)為傳感器陣列測量所得的m個電容測量值，s
i
(x,y)是靈敏度的分布函數；根據電勢分布法，靈敏度場描述為以下公式：其中，v是激勵電壓，和分別是電極i和電極j為激勵電極時的電勢分布；由公式(1)看出，當介電常數相近時，電容向量c
i
與介電常數分布函數ε(x,y)之間的關系可近似為線性；為了方便計算，對其進行離散和歸一化處理，得到ect矩陣形式的數學模型如下所示：c＝sg
????????????????????????
(3)其中，c是歸一化的m
×
1的電容測量值向量，g是大小為n
×
1的灰度向量，代表管道不同介電常數的分布；s是歸一化的大小為m
×
n的靈敏度矩陣，m為測量電容值的個數，n的大小代表了成像區域像素個數，從公式(3)可以看出，ect圖像重建的問題可以簡化為通過測量電容向量c求解灰度向量g的過程；步驟二，生成訓練樣本集，給定k類流型的測量電容數據組成訓練樣本集，記第i類流型的訓練集e
i
共有n
i
個樣本，訓練集e
i
表示為：其中，v
j
∈r
m
×1(j＝1,2
…
n
i
)代表第i類流型里面的第j個訓練樣本，即ect系統從第i類流型獲取的測量電容向量；m是測量電容值得個數，8極板ect系統有m＝28，通過調整介質分布和疊加隨機噪聲的方式來生成訓練樣本集，每類流型的訓練樣本數為100個，則所有k類流型組成的訓練樣本集為：e＝[e1,e2…
e
k
]∈r
m
×
n
，
????????????????
(5)其中n＝n1+n2+
…
n
k
；對于流型i，若訓練集e
i
完備或者過完備，則此類流型的任意測量電容向量λ可由其訓練集e
i
線性表示：其中，系數α
i
(j＝1,2
…
n
i
)為實數，為系數向量，在實際應用中，由于無法預知任意測量電容向量λ的歸屬類，可以考慮將其重新表示為k類流型組成的訓練樣本集e的線性組合：λ＝ex0.
?????????????????????????
(7)式中，為稀疏向量，至此，由ect系統獲得的測量電容
信號λ就可由以訓練樣本集組成的過完備字典變換基e來稀疏表示，這樣就使得ect測量電容信號滿足了壓縮感知理論的前提條件；步驟三，計算測試樣本和標準樣本的觀測值，隨機高斯矩陣b∈r
d
×
m
作為測量矩陣，則根據壓縮感知思想，針對以測量電容向量λ作為原始信號的線性測量過程為：其中，為觀測值，是d
×
1的向量，標準樣本的觀測值按照同樣方法計算；步驟四，計算測試樣本和標準樣本的稀疏系數，由于觀測向量的維數d遠小于稀疏信號x0的維數，直接求解公式(8)是一個欠定問題，將其轉換為如公式(9)所示的數學模型：通過壓縮感知的信號重構算法求解此l0范數下的最優化問題，使用正交匹配追蹤(omp)算法得到測量電容向量λ和標準電容向量λ
k
對應于過完備字典變換基e的稀疏解x0和x
k
；步驟五，計算測試樣本和標準樣本稀疏解的線性相關系數，通過x0與x
i
的線性相關程度來確定測試樣本的歸屬流型，達到流型辨識的目的；相關系數則是衡量變量之間線性相關程度的指標，其計算方法如式(10)所示：其中，和分別為x0和x
i
的平均值，根據公式(10)分別計算測試樣本與所有標準樣本對應的稀疏解之間的線性相關系數；步驟六，確定流型，選擇相關系數最大的歸屬流型靈敏度場，過程如下：將最大相關系數所對應的流型確定為測試樣本的歸屬流型，選擇該歸屬流型的靈敏度場確定為新的靈敏度矩陣；步驟七，圖像重構，輸出灰度向量g。2.如權利要求1所述的基于壓縮感知流型辨識的變化靈敏度場電容層析成像方法，其特征在于，所述步驟七中，用變化靈敏度矩陣和重構算法來進行圖像重構。3.如權利要求1或2所述的基于壓縮感知流型辨識的變化靈敏度場電容層析成像方法，其特征在于，所述步驟三中，使用隨機高斯矩陣b∈r
d
×
m
作為測量矩陣，為了準確的重構出原始信號，測量矩陣必須滿足rip準則，rip準則的等價條件是測量矩陣和稀疏基不相關，當以d
×
m維的隨機高斯矩陣作為測量矩陣時，只要滿足條件：d≥ck log(m/k)，其中c是一個很小的數，k為原始信號的稀疏度，此測量矩陣在極大概率下滿足rip準則。
技術總結
一種基于壓縮感知流型辨識的變化靈敏度場電容層析成像方法，首先，本發明方法設定k類典型流型，將其歸一化后的ECT測量電容信號作為標準樣本并在此基礎上以調整介質分布和疊加隨機噪聲來生成訓練樣本集，待測流型電容值的歸一化數據作為測試樣本。然后，通過計算測試樣本與所有標準樣本的稀疏解之間的線性相關系數，并將最大相關系數所對應的流型確定為測試樣本的歸屬流型，選擇該歸屬流型的靈敏度場替代原始的固定靈敏度場進行圖像重構。通過仿真模擬實驗表明，所得的性能指標和成像效果優于傳統的固定靈敏度場成像方法，圖像邊緣清晰平滑。故本發明方法提高了電容層析成像的精度，為電容層析成像研究提供了一種新的方法和手段。手段。手段。
