電阻抗成像探測(cè)地下溶洞的幾種方法*

吳澤宇，劉祚秋，汪利

中山大學(xué)航空航天學(xué)院，廣東 廣州 510006

電阻抗成像方法是可以用于探測(cè)地下溶洞的一種方法，其原理是根據(jù)物體內(nèi)部不同物質(zhì)的導(dǎo)電參數(shù)（如電阻率、電容率）的不同，通過(guò)對(duì)物體表面電流、電壓的施加和測(cè)量來(lái)獲知物體內(nèi)部導(dǎo)電參數(shù)的分布，進(jìn)而重建出反映物體內(nèi)部結(jié)構(gòu)的圖像，該方法的難點(diǎn)一是固有的病態(tài)性，即邊界電壓數(shù)據(jù)的微小擾動(dòng)可能引起解的巨大變化；另一個(gè)困難是它的信息量小。如何實(shí)現(xiàn)系統(tǒng)的高精度、高分辨率和算法的快速收斂是其應(yīng)用的主要難點(diǎn)［1］。

電阻抗成像的起源可以追溯到20 世紀(jì)20 年代，地球物理學(xué)研究者提出了線性電極陣列的電阻率成像（Resistivity Imaging）技術(shù)。19 世紀(jì)70年代，生物醫(yī)學(xué)研究者提出了圓形電極陣列的斷層電阻率測(cè)量技術(shù)（Tomographic Resistivity Measurement Technique）。1978 年，Henderson 和Webster 做出了第一幅電阻抗圖像。他們使用固定于胸部上的由一大電極和與之相對(duì)的若干小電極組成的電極系統(tǒng)。通過(guò)測(cè)量從各小電極流向大電極的電流所形成的等位差，獲得了可以清楚地顯示肺臟位置的阻抗圖像。但是這還不是斷層圖像，而是類似X 胸片的透視圖像。1982 年，英國(guó)Sheffield大學(xué)的Brown 和Barber 實(shí)現(xiàn)了第一個(gè)手臂的阻抗層析圖像，開辟了電阻抗層析成像（EIT）技術(shù)這一新的研究領(lǐng)域［2］。

本文采用電阻抗成像法探測(cè)了地下溶洞，并采用Tikhonov正則化和稀疏正則化方法以解決反問(wèn)題固有的非適定性問(wèn)題；在正則化的過(guò)程中運(yùn)用了靈敏度分析法，又分別采用線性近似直接求解和迭代求解兩種方法進(jìn)行求解，并進(jìn)一步分析了幾種方法的識(shí)別精度與誤差敏感度，為實(shí)際工程中的地下溶洞識(shí)別提供了一定的參考。

1 問(wèn)題的描述

1.1 正問(wèn)題

考慮一個(gè)空間區(qū)域Ω 及其邊界?Ω［3-4］，電阻抗成像的正問(wèn)題是在已知材料電阻抗率ρ(x)，x∈Ω（或其倒數(shù)，視電阻率c(x) = 1/ρ(x)）的情況下，求解區(qū)域內(nèi)的電勢(shì)u(x)（或電壓）分布。考慮邊界?Ω上輸入輸出電流，正問(wèn)題的控制方程［5］為

其中n為邊界上的外法向向量，SI為所有邊界電流點(diǎn)（或電極）的集合，考慮電流幾乎為單點(diǎn)輸入輸出，el表示的是包含電流點(diǎn)l的微小區(qū)域，Il是點(diǎn)l處的電流量。Il＞0 表示電流輸入，Il＜0 則為電流輸出。如果某個(gè)電極方案只有點(diǎn)1 輸入單位電流，而點(diǎn)2 輸出單位電流，則有I1= 1，I2= -1，Il=0，l∈SI{1，2}。此外，還可以測(cè)量一系列邊界點(diǎn)的電勢(shì)，這些測(cè)量的電勢(shì)點(diǎn)的集合記為Su；電勢(shì)點(diǎn)Su與電流點(diǎn)SI可以有部分重合。問(wèn)題（1）～（2）是典型的Neumann 邊值問(wèn)題，為了確保解的存在性，需要滿足所有輸入輸出電流總和為0，即

其中uk表示點(diǎn)k上的電勢(shì)。

基于以上描述，考慮電阻抗成像探測(cè)地下溶洞這一問(wèn)題（如圖1 所示）。問(wèn)題所屬區(qū)域簡(jiǎn)化為圖1右側(cè)的矩形區(qū)域。電流點(diǎn)與電勢(shì)點(diǎn)不一致，但數(shù)目相等。

圖1 電阻抗成像探測(cè)地下溶洞示意圖Fig.1 Diagram of electrical impedance imaging for cavity detection

通?？梢允褂糜邢拊M(jìn)行正問(wèn)題求解，得到如下方程

其中Kj為單位電阻率時(shí)的單元?jiǎng)偠染仃嚒?/p>

1.2 反問(wèn)題

2 電阻抗成像的理論與方法

2.1 初值或平均視電阻率估計(jì)

此時(shí)，可以選擇初值為c0=c01。

2.2 靈敏度分析

接下來(lái)，將根據(jù)已有的初值c0和靈敏度矩陣進(jìn)行參數(shù)識(shí)別?？紤]使用不同的正則化方法（Tikhonov 正則化［8］和稀疏正則化［9］）以及是否迭代的因素，主要可以分為4 種方法：Tikhonov 更新法、稀疏更新法、Tikhonov迭代法、稀疏迭代法。

2.3 Tikhonov更新法

這里假設(shè)初值c0與待識(shí)別結(jié)果c比較接近，此時(shí)可以對(duì)R(c)使用Taylor展開，得到

其中正則化稀疏λ＞0 可通過(guò)L-curve 方法［10-11］進(jìn)行確定。

2.4 稀疏更新法

同樣考慮如公式（12）～（13）的線性近似的過(guò)程，但使用稀疏正則化處理非適定性。實(shí)際溶洞探測(cè)問(wèn)題中，溶洞的存在相對(duì)于探測(cè)區(qū)域通常體積較小，在空間分布上具有一定稀疏性（或更新量(c-c0)是稀疏的），因此可以使用稀疏正則化進(jìn)行處理。稀疏更新法的求解問(wèn)題為

2.5 Tikhonov迭代法

2.6 稀疏迭代法

與Tikhonov迭代法類似，稀疏迭代法則表現(xiàn)為使用稀疏更新法重復(fù)、迭代求解的過(guò)程。迭代式為

3 數(shù)值算例

選取矩形區(qū)域長(zhǎng)寬均為4 m×4 m，每邊均勻劃分為8 個(gè)單元，假設(shè)初始區(qū)域的視電阻率是均勻的，取值為c0= 10?？紤]如圖2 中三種溶洞情況（溶洞處的視電阻率等效為3，相對(duì)于土體的10，有較大的折減），其中Case 1、Case 2 為單個(gè)溶洞的情況、Case 3為多個(gè)溶洞的情況，圖塊顏色與視電阻率的對(duì)應(yīng)關(guān)系如圖右側(cè)色標(biāo)所示。

圖2 溶洞探測(cè)的三種情況Fig.2 Three cases of cavity detection

考慮如圖1右側(cè)的矩形區(qū)域和電極方案，通過(guò)給定的電流與真實(shí)視電阻率分布得到電勢(shì)結(jié)果，分別利用Tikhonov 更新法、稀疏更新法、Tikhonov迭代法和稀疏迭代法識(shí)別以上三種溶洞情況。其中Tikhonov 法的正則化參數(shù)通過(guò)L-curve 方法進(jìn)行確定，而稀疏正則化參數(shù)取10-8。

另外，考慮誤差存在時(shí)的情形，此時(shí)電勢(shì)情況u=u0(1+ fr·Randn)；其中，u0為真實(shí)電勢(shì)；Randn是與電勢(shì)矩陣同階的隨機(jī)數(shù)矩陣，其每個(gè)元素都是[-1，1]上的正態(tài)分布隨機(jī)變量；fr為變量的系數(shù)，這里取fr = 0.001。

Case 1在無(wú)誤差和有誤差情況下的識(shí)別結(jié)果分別如圖3～4 所示，Case 2 在無(wú)誤差和有誤差情況下的識(shí)別結(jié)果分別如圖5～6 所示，Case 3 在無(wú)誤差和有誤差情況下的識(shí)別結(jié)果分別如圖7～8所示。從結(jié)果中可以看出：

圖3 無(wú)誤差時(shí)Case 1的識(shí)別結(jié)果Fig.3 Results of Case 1 without error

1）Tikhonov 更新法在無(wú)誤差的情況下識(shí)別單個(gè)溶洞（圖3 和圖5）時(shí)，在溶洞位置可以識(shí)別到視電阻率的折減，但與真實(shí)情況相比，視電阻率在數(shù)值上相差較大，結(jié)果較為粗糙；而在有誤差的情況下（圖4和圖6），在溶洞位置也能看出視電阻率的差別，但與無(wú)誤差情況相比，結(jié)果的差別較大；而在識(shí)別多個(gè)溶洞的情況時(shí)（圖7 和圖8），無(wú)論是否有誤差，視電阻率的分布較為混亂，多處出現(xiàn)電阻率折減，電阻率在數(shù)值上與真實(shí)情況有較大差別，無(wú)法識(shí)別出正確的溶洞位置。

圖4 有誤差時(shí)Case 1的識(shí)別結(jié)果Fig.4 Results of Case 1 with error

圖5 無(wú)誤差時(shí)Case 2的識(shí)別結(jié)果Fig.5 Results of Case 2 without error

圖6 有誤差時(shí)Case 2的識(shí)別結(jié)果Fig.6 Results of Case 2 with error

2）稀疏更新法在無(wú)誤差的情況下識(shí)別單個(gè)溶洞（圖3 到圖6）時(shí)，可以明顯地看到溶洞位置有明顯的視電阻率折減，而在其他位置的視電阻率在數(shù)值上與真實(shí)情況十分接近；識(shí)別相鄰的多個(gè)溶洞的情況時(shí)（圖7和圖8），較大溶洞處有視電阻率的折減，但不能完全識(shí)別出溶洞形狀，數(shù)值上也與真實(shí)情況有較大差別；且較小溶洞處沒(méi)有視電阻率的折減，其他位置的視電阻率也與真實(shí)情況有一定差別，但各處視電阻率之間相差不大，且無(wú)誤差和有誤差結(jié)果相似。

圖7 無(wú)誤差時(shí)Case 3的識(shí)別結(jié)果Fig.7 Results of Case 3 without error

圖8 有誤差時(shí)Case 3的識(shí)別結(jié)果Fig.8 Results of Case 3 with error

3） Tikhonov 迭代法在無(wú)誤差的情況下識(shí)別單個(gè)溶洞（圖3 和圖5）時(shí)，在溶洞位置可以識(shí)別到視電阻率的折減，其他位置的視電阻率與真實(shí)情況在數(shù)值上存在較大差別，且視電阻率分布不均勻。在有誤差的情況下（圖4和圖6），視電阻率分布混亂，有多處位置存在明顯的視電阻率折減，無(wú)法識(shí)別出溶洞位置；而在無(wú)誤差的情況下識(shí)別多個(gè)溶洞時(shí)（圖7），能準(zhǔn)確識(shí)別到溶洞位置處的視電阻率折減，其他位置的視電阻率也較為均勻。但在有誤差的情況下（圖8），整個(gè)識(shí)別區(qū)域的視電阻率分布混亂，無(wú)法識(shí)別到溶洞位置。

4） 稀疏迭代法溶洞探測(cè)效果非常好，在無(wú)誤差的情況下識(shí)別單個(gè)與多個(gè)溶洞（圖3、圖5 和圖7）時(shí)，結(jié)果與真實(shí)情況幾乎完全一致；而在有誤差的情況下（圖4、圖6 和圖8），結(jié)果也基本保持不變。

4 結(jié) 論

Tikhonov 更新法方法簡(jiǎn)單，可以大體識(shí)別單個(gè)溶洞的位置，但結(jié)果粗糙，對(duì)誤差較為敏感，且無(wú)法識(shí)別相鄰的多個(gè)溶洞；稀疏更新法也比較簡(jiǎn)單，識(shí)別單個(gè)溶洞時(shí)結(jié)果較好，對(duì)誤差敏感度也較低，但在識(shí)別相鄰的多個(gè)溶洞時(shí)存在遺漏；Tikhonov 迭代法比Tikhonov 更新法復(fù)雜，在無(wú)誤差時(shí)能識(shí)別單個(gè)或多個(gè)溶洞的位置，但對(duì)誤差的敏感程度非常大，在有誤差情況下完全無(wú)法識(shí)別出溶洞位置；稀疏迭代法計(jì)算量大，但識(shí)別結(jié)果非常好，對(duì)誤差的敏感度也很低，無(wú)論是否有誤差，無(wú)論是識(shí)別單個(gè)還是多個(gè)溶洞，都能得到非常理想的結(jié)果。

總而言之，從結(jié)果上看，幾種方法的排序（由好到壞）為：稀疏迭代法＞稀疏更新法＞Tikhonov更新法＞Tikhonov迭代法。