改進SSOR 迭代法的數值保角變換計算法

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

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

0 引言

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

共軛調和函數具有嚴謹的理論基礎而且容易計算。模擬電荷法是Steinberger 計算電場時提出的方法，Amano 利用模擬電荷法將數值保角變換的近似問題轉換為共軛調和函數的逼近問題。相較于Symm 提出的積分方程法和Nass?er 提出的帶有廣義Neumann kernel 的邊界積分法，基于模擬電荷法的數值保角變換計算法避免了很多數值積分微分等計算操作，可降低數值保角變換計算的復雜度和時間。

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

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

1 相關知識

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

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

其中gv和hv是一對共軛調和函數。

（1）邊界條件。

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

（3）約束條件。

（4）正則化條件。

可得：

（5）柯西條件。復變函數（4）是區域D內的解析函數，D內的任意閉曲線設為，設為只包圍Cl的封閉曲線：

故有：

可得：

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

其中ζl0為Cl內的一點。

將式（12）進行變形得：

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

2 本文方法

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

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

其中，AT A為對稱正定矩陣，ωi∈(0,2),i=1,2 為松弛因子。根據參考文獻［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 為預條件算子，一般取K 為D=diag(a1,1,a2,2,…,aN,N)或單位矩陣，i=1 時，μi取，i=2 時，μi取，x值更接近準確值，λ1，λN+1分別為Mi-Ni的最大特征值、最小特征值。

由式（18）構造兩個半步迭代公式：

根據式（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 方法收斂性如下：

證明：根據參考文獻［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相關算法。

算法流程如圖1 所示。

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

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

3 實驗結果

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

例1：當l=1 時的橙形外部區域

邊界：

邊界參數方程：

約束點：

電荷點：

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

圖2 構造40×40 網格，圖3 是橙形外部區域。約束點分布在邊界上，電荷點在區域的內部，其中i是虛數單位，0＜q＜1，是調整電荷點位置的參數。

圖3 是圖2 用改進后的算法進行數值保角變換結果。從圖中可以看出變換前后的角度關系，顯示該算法的優越性。圖中把橙形變成一條豎直的線，由圖2 和圖3 可以看到網格變換前后的變化。

Fig.2 Problem areas（a）圖2 問題區域

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

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

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

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

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

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

Fig.7 Transform result of Fig.6圖7 圖6 的變換結果

例2：單聯通橢圓外部區域

電荷點：

約束點：

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

圖8 是單聯通橢圓外部區域，構造了40 × 40 的網格，圖9 是圖8 變換的結果，由保角變換把橢圓變成一條豎直狹縫，網格也發生一系列畸變，反映保角變換前后的角度變化，這些角度變化很好地表現出該算法的優勢。

Fig.8 Problem areas（b）圖8 問題區域

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

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

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

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

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

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

Fig.13 Transform result of Fig.12圖13 圖12 的變換結果

4 結語

本文提出基于改進SSOR 法的數值保角變換計算法，解決了原有方法與一些傳統方法在求解電荷量上的缺陷導致的數值保角變換結果誤差大及不穩定等問題。從數值實驗結果可知，本文方法在單連通外部區域例證中精度和穩定性較高，能夠應用于圖像領域。

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

圖像保角變換在醫學圖像處理中有著重要作用，未來可嘗試處理醫學圖像中的一些實際問題，如將圖像上狹小的縫隙在具有保角特性下放大觀察，將一些具有干擾性的位置變為可忽視的狹縫。數值保角變換在流體力學中也有重要應用，可以用變換函數繪制無壓縮流體在遇到障礙物時的曲線。