萬 鵬，呂毅斌，王櫻子，唐勝男

（1.昆明理工大學(xué) 理學(xué)院；2.昆明理工大學(xué) 計算中心，云南 昆明 650500）

0 引言

保角變換又稱共形映射，是復(fù)變函數(shù)的重要內(nèi)容，在很多領(lǐng)域扮演著重要角色，如流體力學(xué)、醫(yī)學(xué)圖像處理、電磁理論等［1-4］。然而，除了在少數(shù)特殊情況下可以解出保角變換的解析解外，大多數(shù)工程與科學(xué)問題中求解保角變換的解析非常復(fù)雜甚至無法求解，所以對數(shù)值保角變換研究成為科學(xué)工程領(lǐng)域必不可少的課題。目前比較著名的方法有Symm［5-6］提出的積分方程法，Dai［7-10］提出的基于模擬電荷法的數(shù)值保角變換計算法，Sangawi 等［11-12］提出的帶有廣義Neumann kernel 的邊界積分法等。

共軛調(diào)和函數(shù)具有嚴謹?shù)睦碚摶A(chǔ)而且容易計算。模擬電荷法是Steinberger 計算電場時提出的方法，Amano 利用模擬電荷法將數(shù)值保角變換的近似問題轉(zhuǎn)換為共軛調(diào)和函數(shù)的逼近問題。相較于Symm 提出的積分方程法和Nass?er 提出的帶有廣義Neumann kernel 的邊界積分法，基于模擬電荷法的數(shù)值保角變換計算法避免了很多數(shù)值積分微分等計算操作，可降低數(shù)值保角變換計算的復(fù)雜度和時間。

模擬電荷點的位置及數(shù)量選取問題至今沒有統(tǒng)一方法，很多情況下數(shù)值保角變換的最終結(jié)果對模擬電荷點位置及數(shù)量選取極為敏感，對于簡單的較為規(guī)則的邊界問題而言，如圓形或橢圓外部區(qū)域，Amano 提出的基于模擬電荷法的數(shù)值保角變換計算法具有較高精度，但是對于一些較為復(fù)雜的邊界如橙形、心形等，基于模擬電荷法的數(shù)值保角變換誤差會明顯上升，所以基于模擬電荷法的數(shù)值保角變換計算法精度還有很大的改善空間。

本文對單連通橙形外部區(qū)域、單連通橢圓外部區(qū)域給出電荷點和約束點位置，使用基于（k，j）-Padé［13-14］迭代法改進的SSOR 方法［15-16］求解電荷量，以降低模擬電荷點數(shù)量選取對最終結(jié)果的影響，最后利用求得的電荷量計算出數(shù)值保角變換函數(shù)。在數(shù)值實驗中針對不同電荷點數(shù)量進行多組實驗，通過誤差比較證明改進算法的精確性，并實現(xiàn)數(shù)字圖像的保角變換。

1 相關(guān)知識

將z-plane上Jordan 閉曲線C1,C2,…,Cn外側(cè)的無窮遠點的非有界多連通區(qū)域D映射到ω-plane（含有豎直狹縫的平面）。如果z=0 在區(qū)域D內(nèi)，在正則化條件f(∞)=∞，f'(∞)=1，a0=0 下，豎直角度的保角變換在無窮遠點的Laurent 級數(shù)可唯一表示為如下形式：

保角變換可將曲線C1,C2,…,Cn映射成豎直狹縫ω=fu(z)，fu(z)可以表示為：

其中g(shù)v和hv是一對共軛調(diào)和函數(shù)。

（1）邊界條件。

根據(jù)模擬電荷法，gu(z)+ihu(z)可逼近Gu(z)+iHu(z)，它可以表示為：

（3）約束條件。

（4）正則化條件。

可得：

（5）柯西條件。復(fù)變函數(shù)（4）是區(qū)域D內(nèi)的解析函數(shù)，D內(nèi)的任意閉曲線設(shè)為，設(shè)為只包圍Cl的封閉曲線：

故有：

可得：

其中Nl為Cl內(nèi)電荷點數(shù)量，ζli(i=1,2,…,Nl)是分布在曲線Cl外部區(qū)域D中的電荷點，zmj(j=1,2,…,Nm)為在邊界Cm上的約束點，電荷量Qli由公式（10）和邊界條件（3）以及電荷點和約束點決定，Um是狹縫位置，為um的近似值。式（6）計算出的電荷量和狹縫位置用來逼近fu，即：

其中ζl0為Cl內(nèi)的一點。

將式（12）進行變形得：

根據(jù)式（6）、式（10）和式（13），聯(lián)立一次方程求電荷點以及狹縫位置U1,U2,…,Un：

2 本文方法

上述約束方程將其寫為標準線性方程：

改進的SSOR 迭代法［11-12］可以很好地改善其病態(tài)性。首先對式（1）兩邊分別乘以ωi AT：

其中，AT A為對稱正定矩陣，ωi∈(0,2),i=1,2 為松弛因子。根據(jù)參考文獻［11］、［12］對ωi AT A進行分裂，ωi AT A=Mi-Ni,i=1,2，于是有：

其中M1=D-ω1L，N1=(1-ω1)D+ω1U，M2=D-ω2U，N2=(1-ω2)D+ω2L，D=diag(a1,1,a2,2,…,aN,N)，L為負的嚴格下三角矩陣，U為負的嚴格上三角矩陣。

對式（17）兩邊乘以μiK-1，得：其中，K 為預(yù)條件算子，一般取K 為D=diag(a1,1,a2,2,…,aN,N)或單位矩陣，i=1 時，μi取，i=2 時，μi取，x值更接近準確值，λ1，λN+1分別為Mi-Ni的最大特征值、最小特征值。

由式（18）構(gòu)造兩個半步迭代公式：

根據(jù)式（18）定義矩陣多項式r(-μi(Mi-Ni))=E-μiK-1(Mi-Ni)=ri，t(-μi(Mi-Ni))=ti，q(-μi(Mi-Ni))=qi，h(-μi(Mi-Ni))=hi，得到迭代的矩陣多項式：

改進SSOR 方法收斂性如下：

證明：根據(jù)參考文獻［3］，對任意的n≥m，x＜0，μi＞0，ωi＞0，rnm(x)＜1，由于AT A是正定的，故ωi AT A=Mi-Ni也是正定的，故Mi-Ni的特征值都大于0，則-μi(Mi-Ni)的特征值都小于0，則有：

故改進的SSOR 迭代法基于（n，m）-Padé 迭代法收斂。

基于上述方法給出基于(k,j) -Pade'迭代法的SSOR相關(guān)算法。

算法流程如圖1 所示。

Fig.1 Based on the（n，m）-Padé SSOR algorithm of iterative process圖1 基于（n，m）-Padé 迭代法的SSOR 算法流程

通過以上改進的SSOR 方法得到基于（n，m）-Padé 迭代法的數(shù)值保角變換步驟：①給出每條曲線C1,C2,…,Cn上電荷點的數(shù)量Nl(l=1,2…,n)，以及電荷點的坐標和約束點的坐標(m=1,2,…,n)；②根據(jù)公式（13）構(gòu)造約束方程組；③利用改進的SSOR 法和迭代公式求解電荷量和狹縫位置；④由模擬電荷法求出保角變換的近似函數(shù)；⑤利用改進的數(shù)值保角變換計算法實現(xiàn)圖像保角變換。

3 實驗結(jié)果

在Windows10 系統(tǒng)上用Matlab 2016b 對以橙形為邊界的單連通區(qū)域進行數(shù)值實驗并分析實驗結(jié)果，然后通過該方法實現(xiàn)圖像的數(shù)值保角變換。根據(jù)復(fù)變函數(shù)中的最大模原理［18-19］得到誤差和豎直狹縫位置的誤差計算公式。

例1：當l=1 時的橙形外部區(qū)域

邊界：

邊界參數(shù)方程：

約束點：

電荷點：

其中：z0=zN,zN+1=z1

圖2 構(gòu)造40×40 網(wǎng)格，圖3 是橙形外部區(qū)域。約束點分布在邊界上，電荷點在區(qū)域的內(nèi)部，其中i是虛數(shù)單位，0＜q＜1，是調(diào)整電荷點位置的參數(shù)。

圖3 是圖2 用改進后的算法進行數(shù)值保角變換結(jié)果。從圖中可以看出變換前后的角度關(guān)系，顯示該算法的優(yōu)越性。圖中把橙形變成一條豎直的線，由圖2 和圖3 可以看到網(wǎng)格變換前后的變化。

Fig.2 Problem areas（a）圖2 問題區(qū)域

Fig.3 Fig.2 transformation results under the improved algorithm圖3 圖2 在改進算法下的變換結(jié)果

用式（20）作為誤差指標，圖4 和圖5 為分別用Amano法、LQ 分解法、雙共軛梯度法和本文方法求解電荷量，并計算出數(shù)值變換函數(shù)，最后畫出變換函數(shù)的誤差曲線。從圖4 可以看出，隨著電荷數(shù)量的增加，本文方法的誤差一直在減小，而Amano 方法在數(shù)量較小的情況下誤差明顯增加，在超過一個特定范圍之后又明顯減小，但本文方法是一直在減小，而且和Amano 方法有很大差距。從實驗數(shù)據(jù)計算可知，本文方法比原有方法精度平均提高0.101 9，最大提升率為11.2%，LQ 分解法和雙共軛梯度法雖然在一定程度上比Amano 方法好，但是整體效果還是比改進算法要差一些。從圖5 可以看到豎直狹縫位置之間的誤差。Amano 方法的誤差剛開始很大，然后突然減小，說明Ama?no 方法不穩(wěn)定。雖然本文方法和Amano 方法從誤差曲線上看差距不是很大，但本文方法誤差剛開始時在減小，當電荷數(shù)量達到某個臨界點時誤差曲線會趨于平衡狀態(tài)，說明本文方法更加穩(wěn)定。從實驗數(shù)據(jù)計算可知，本文方法比原有方法精度平均提高0.009 0，故通過圖3 與圖4 可以看出本文方法較Amano 方法更加穩(wěn)定，且本文方法可改善方程組的病態(tài)性，對于求解病態(tài)的對稱約束方程組有明顯優(yōu)勢。LQ 分解法和雙共軛梯度法相對本文算法而言不穩(wěn)定，精度也略有欠缺。

Fig.4 Erel（a）圖4 Erel（a）

Fig.5 Vertical slit position error（a）圖5 豎直狹縫位置誤差

圖6 是原始圖像，圖7 是圖像變換結(jié)果。把橙形變成豎直狹縫得到圖7 的變換結(jié)果。

Fig.6 Original image（a）圖6 原始圖像

Fig.7 Transform result of Fig.6圖7 圖6 的變換結(jié)果

例2：單聯(lián)通橢圓外部區(qū)域

電荷點：

約束點：

約束點分布在邊界C 上，電荷點分布在邊界C 內(nèi)部，其中i是虛數(shù)，0 ＜q＜1 是控制電荷點的位置分布。

圖8 是單聯(lián)通橢圓外部區(qū)域，構(gòu)造了40 × 40 的網(wǎng)格，圖9 是圖8 變換的結(jié)果，由保角變換把橢圓變成一條豎直狹縫，網(wǎng)格也發(fā)生一系列畸變，反映保角變換前后的角度變化，這些角度變化很好地表現(xiàn)出該算法的優(yōu)勢。

Fig.8 Problem areas（b）圖8 問題區(qū)域

圖10 和圖11 表示在公式（20）下根據(jù)兩種誤差評判指標畫出的誤差曲線圖，圖10 和圖11 畫出了本文方法和Amano 方法的誤差曲線。從圖10 可以明顯看出本文方法和Amano 方法的誤差曲線，隨著電荷數(shù)量的增加誤差都在減小，但本文方法和Amano 方法的誤差曲線還是存在區(qū)別。從實驗數(shù)據(jù)計算得到本文方法比原有方法精度平均提高約0.026 4，最大提升率為5.4%，而LQ 分解法和雙共軛梯度法在簡單邊界上的變換效果比原有Amano 方法要差一些。圖11 中雖然兩者曲線差距不太明顯，精度平均提升0.002 4，但是本文方法震動弧度不是很大，比原有方法更為穩(wěn)定。而LQ 分解法和雙共軛梯度法不穩(wěn)定，誤差一直在上下波動，所以本文方法體現(xiàn)出穩(wěn)定性好的特點。因此，本文方法不僅穩(wěn)定性好而且精度高，優(yōu)勢更強。

Fig.9 Fig.8 transformation results with the improved algorithm圖9 圖8 在改進算法下的變換結(jié)果

Fig.10 Erel（b）圖10 Erel（b）

Fig.11 Vertical slit position error（b）圖11 豎直狹縫位置誤差

圖13 是圖12 通過數(shù)值保角變換后的圖像，是一條豎直狹縫，在圖中可以看到高樓的彎曲效果，從中可以更好地看出它們之間的角度關(guān)系。

Fig.12 Original image（b）圖12 原始圖像

Fig.13 Transform result of Fig.12圖13 圖12 的變換結(jié)果

4 結(jié)語

本文提出基于改進SSOR 法的數(shù)值保角變換計算法，解決了原有方法與一些傳統(tǒng)方法在求解電荷量上的缺陷導(dǎo)致的數(shù)值保角變換結(jié)果誤差大及不穩(wěn)定等問題。從數(shù)值實驗結(jié)果可知，本文方法在單連通外部區(qū)域例證中精度和穩(wěn)定性較高，能夠應(yīng)用于圖像領(lǐng)域。

但改進的SSOR 法在雙連通及多連通上的數(shù)值效果仍有待檢驗，本文方法的一些因子常數(shù)選取是否為最佳值也需要繼續(xù)研究和改善。

圖像保角變換在醫(yī)學(xué)圖像處理中有著重要作用，未來可嘗試處理醫(yī)學(xué)圖像中的一些實際問題，如將圖像上狹小的縫隙在具有保角特性下放大觀察，將一些具有干擾性的位置變?yōu)榭珊鲆暤莫M縫。數(shù)值保角變換在流體力學(xué)中也有重要應(yīng)用，可以用變換函數(shù)繪制無壓縮流體在遇到障礙物時的曲線。