汪強(qiáng)強(qiáng)， 劉海飛， 柳建新， 李 星， 施昕祎

(1.中南大學(xué) 地球科學(xué)與信息物理學(xué)院，長(zhǎng)沙 410083；2.有色資源與地質(zhì)災(zāi)害探查湖南省重點(diǎn)實(shí)驗(yàn)室，長(zhǎng)沙 410083)

0 引言

地球物理反演問題大多數(shù)是非線性多元函數(shù)的極值問題，由于涉及高階偏導(dǎo)數(shù)矩陣的計(jì)算和存儲(chǔ)，采用直接法對(duì)其進(jìn)行求解尤為困難，必須通過線性化手段將非線性問題轉(zhuǎn)換為解方程組的問題[1]。將非線性問題進(jìn)行泰勒級(jí)數(shù)展開并略去二次以上的項(xiàng)可以得到線性方程組，但得到的線性方程組往往是病態(tài)的[2]。利用奇異值分解法先將病態(tài)線性方程組進(jìn)行分解，然后采用正則化技術(shù)削弱小奇異值的影響，可以使求得的解更好地逼近真實(shí)解，能夠達(dá)到有效求解病態(tài)方程組的目的[3]。目前削弱小奇異值影響的正則化技術(shù)主要有:Wiggins零化小奇異值法和Tikhonov濾波正則化方法[4]，基于Tikhonov濾波正則化方法求取病態(tài)線性方程組的正則解，尋找最佳正則化因子成為奇異值分解法成功求解病態(tài)線性方程組的關(guān)鍵。Hansen[5]提出了L-曲線的概念，并證明了點(diǎn)(ln‖b-Axλ‖2,ln‖xλ‖2)所構(gòu)成的L-曲線上曲率最大的點(diǎn)對(duì)應(yīng)最佳正則化因子的位置；O’ Leary等[6]利用三次樣條平滑的L-曲線計(jì)算曲率從而確定最佳角點(diǎn)；J.Longina Castellanos等[7]采用三角形法確定L-曲線的最佳角點(diǎn)；Hansen等[8]提出了一種自適應(yīng)修剪算法，確定L-曲線最佳角點(diǎn)的效果較好。筆者在前人工作的基礎(chǔ)上，給出了一種病態(tài)線性方程組正則化因子的分步選取方法，即先根據(jù)L-曲線的局部特征找到可能的候選角點(diǎn)，再?gòu)暮蜻x角點(diǎn)中找到最佳角點(diǎn)，進(jìn)而確定最佳正則化因子，再將其代入濾波正則化方程即得病態(tài)線性方程組的全局最優(yōu)解。最后，通過希爾伯特方程組對(duì)算法進(jìn)行測(cè)試，驗(yàn)證了本文算法的有效性。

1 奇異值分解法

地球反演中常因線性反演方程的嚴(yán)重病態(tài)特征導(dǎo)致計(jì)算結(jié)果存在較大誤差，奇異值分解法是求解病態(tài)線性方程組主要方法之一。對(duì)于線性方程組

Ax=b

(1)

(2)

當(dāng)λ大于某數(shù)之后，若A+λI正定，則可以得到原問題的直接解而無需奇異值分解。但是λ的取值難以確定，而奇異值分解無需判定矩陣A是否對(duì)稱正定，對(duì)線性方程組的系數(shù)矩陣進(jìn)行奇異值分解后，可以通過分解的結(jié)果得到不同λ的正則解。A的奇異值分解為

(3)

其中：左矩陣U=(u1,u2,…,un)，且U∈Rm×n，UTU=In,In為單位矩陣；右矩陣V=(v1,v2,…,vn)，且V∈Rn×n，VTV=In；對(duì)角陣S=diag(σ1,σ2,…,σn)，且奇異值σi以σ1≥σ2≥…≥σn≥0依次遞減。則方程(1)的解x為

(4)

其中：A+為系數(shù)矩陣A的廣義逆；S+為對(duì)角陣S的廣義逆。從式(4)可以看出，隨著i的增加奇異值σi逐漸趨于零時(shí)或者說矩陣A的條件數(shù)σ1/σn很大時(shí)，直接求解Ax=b很難求得精確解，采用正則化削弱奇異值的影響，可以使所得解更接近精確解。經(jīng)正則化后，線性方程組的解為[9]式(5)。

(5)

選擇合適的正則化因子λ，可以壓制較小奇異值的影響。

圖1 L-曲線示意圖Fig.1 Schematic diagram of L-curve

2 基于L-曲線尋找最佳正則化因子

2.1 L-曲線

所謂L-曲線是指點(diǎn)(ln‖b-Axλ‖2,ln‖xλ‖2)在直角坐標(biāo)系中構(gòu)成的曲線圖，其中‖xλ‖2為正則化解的范數(shù)，‖b-Axλ‖2為正則化解的殘差范數(shù)。根據(jù)Tikhonov正則化方法，當(dāng)正則化因子λ增大時(shí)，正則化解的范數(shù)‖xλ‖2減小，相應(yīng)的殘差范數(shù)‖b-Axλ‖2增大，反之正則化解的范數(shù)‖xλ‖2增大，殘差范數(shù)‖b-Axλ‖2減小，如圖1所示。當(dāng)λ取到較為合理的值時(shí)，正則化解的范數(shù)‖xλ‖2和殘差范數(shù)‖b-Axλ‖2能夠達(dá)到較好的平衡，這個(gè)正則化因子就是最佳正則化因子，它對(duì)應(yīng)L-曲線的最佳角點(diǎn)。

基于L-曲線選取最佳角點(diǎn)的方法主要有最大曲率法、三角形法等，O’Leary等[6]提出，如果平滑連續(xù)的L-曲線由(ρ,η)表示，ρ和η分別表示解的范數(shù)和對(duì)應(yīng)的殘差范數(shù)，并且ρ和η為二階可導(dǎo)的函數(shù)，那么可以通過式(6)直接計(jì)算曲率k(λ)；對(duì)于由許多離散點(diǎn)組成的L-曲線，先擬合一條關(guān)于離散點(diǎn)的三次樣條曲線，再通過曲率計(jì)算公式計(jì)算曲率，從而確定最大曲率點(diǎn)，最大曲率點(diǎn)對(duì)應(yīng)L-曲線的最佳角點(diǎn)。

(6)

(7)

最佳角點(diǎn)對(duì)應(yīng)定向面積為負(fù)且角度小于7π/8的頂點(diǎn)。

盡管最大曲率法和三角形法都可以基于L-曲線選定最佳角點(diǎn)，但它們各自都存在不足之處。對(duì)于由離散點(diǎn)構(gòu)成的L-曲線，最大曲率法通過擬合一條關(guān)于離散點(diǎn)的三次樣條曲線來計(jì)算曲率，但這條樣條曲線對(duì)離散點(diǎn)的分布很敏感，并且有時(shí)不能完全考慮所有的離散點(diǎn)，導(dǎo)致計(jì)算得到的最大曲率點(diǎn)并非對(duì)應(yīng)最佳正則化因子。三角形方法在選定最佳角點(diǎn)時(shí)，位于L-曲線垂直分支上的離散點(diǎn)都參與了計(jì)算，導(dǎo)致計(jì)算量大，復(fù)雜性高。針對(duì)以上問題，我們提出了一種正則化因子的分步選取方法，考慮L-曲線的斜率突變區(qū)可能存在局部角點(diǎn)這一問題，基于一系列子L-曲線選擇候選角點(diǎn)并在這些候選角點(diǎn)中選擇L-曲線的最佳角點(diǎn)。

2.2 分步選取法

前面提出的奇異值分解求解病態(tài)方程組的正則化因子的分步選取方法，主要分為兩個(gè)階段，①在L-曲線上抽取不同數(shù)目的點(diǎn)構(gòu)建新的L-曲線，尋找每條新L-曲線的最佳角點(diǎn)作為候選角點(diǎn)，可構(gòu)成一個(gè)候選角點(diǎn)列表；②從第一階段得到的候選角點(diǎn)列表中選出最佳角點(diǎn)，在候選角點(diǎn)列表中，最佳角點(diǎn)應(yīng)該是到達(dá)垂直分支之前的最后一個(gè)候選角點(diǎn)。如果候選角點(diǎn)分布于L-曲線的垂直分支和水平分支上，那么最佳角點(diǎn)為L(zhǎng)-曲線到達(dá)垂直分支前的一個(gè)候選角點(diǎn)；如果曲線的垂直分支上沒有任何候選角點(diǎn)，那么最后一個(gè)候選點(diǎn)為最佳角點(diǎn)。

2.2.1 搜索候選角點(diǎn)

在算法的第一階段中,基于新的L-曲線搜索候選角點(diǎn)主要包括兩個(gè)準(zhǔn)則:①基于角度搜索候選角點(diǎn);②基于最小距離搜索候選角點(diǎn)。

1)基于角度搜索候選角點(diǎn)：對(duì)于L-曲線上n個(gè)離散點(diǎn)按照λ減小的順序依次編號(hào)為Pi(i=1,2,…n)，用向量vi,i+1(i=1,2,…,n-1)表示點(diǎn)Pi與點(diǎn)Pi+1構(gòu)成的向量。首先計(jì)算向量的模長(zhǎng)并對(duì)向量作歸一化處理使得‖vi,i+1‖2=1，再選擇其中模長(zhǎng)較長(zhǎng)的p(p=min(5,n-1))個(gè)向量，這些向量構(gòu)成新的L-曲線，計(jì)算相鄰向量之間的叉積

delta=(vx,x+1)1(vy,y+1)2-(vx,x+1)2(vy,y+1)1

(8)

其中：x

2)基于最小距離搜索候選角點(diǎn)：定義水平向量v=(-1,0)T，對(duì)于上述子L-曲線上的p個(gè)向量，計(jì)算向量與水平向量v之間的叉積，公式下角標(biāo)的1和2表示向量的第1個(gè)和第2個(gè)坐標(biāo)

-(vj,j+1)2

(9)

|Δ|的大小反映歸一化向量在垂直方向上分量的大小，在這些向量中選擇垂直分量最小的向量vlh,lh+1作為具有水平特征的向量，選擇垂直分量最大的向量vlv,lv+1作為具有垂直特征的向量，但具有垂直特征的向量必須位于水平特征向量左側(cè)。那么“原點(diǎn)o”可以定義為Plh和所處的水平線與向量vlv,lv+1所在的直線的交點(diǎn)，然后在整個(gè)L-曲線上的所有離散點(diǎn)中選擇與該原點(diǎn)o具有最小歐幾里得距離的點(diǎn)作為候選角點(diǎn)[10]。

2.2.2 選定最佳角點(diǎn)

判斷第一階段得到的候選角點(diǎn)列表，如果列表為空，那么認(rèn)為該L-曲線不存在最佳角點(diǎn)，否則將點(diǎn)P1添加到候選角點(diǎn)列表中，并對(duì)列表按照索引值由小到大排序，假設(shè)原來列表中候選角點(diǎn)個(gè)數(shù)為k，那么它們構(gòu)成k個(gè)向量vj,j+1(j=1,2,…,k)，對(duì)向量作歸一化處理，判斷每個(gè)向量的垂直變化量與水平變化量,即正則化解范數(shù)的變化量與對(duì)應(yīng)的殘差范數(shù)的變化量，如果正則化解范數(shù)變化量不小于相應(yīng)的殘差范數(shù)變化量，那么認(rèn)為該向量具有垂直特征，否則認(rèn)為該向量具有水平特征。如果所有的向量都是具有水平特征，那么選擇候選角點(diǎn)列表中的最后一個(gè)候選角點(diǎn)作為L(zhǎng)-曲線的最佳角點(diǎn)，如果存在具有垂直特征的向量，那么選擇到達(dá)垂直特征向量前的最后一個(gè)候選角點(diǎn)作為最佳角點(diǎn)。

3 數(shù)值實(shí)驗(yàn)分析

3.1 病態(tài)問題分析

高階希爾伯特(Hilbert)矩陣方程組具有嚴(yán)重病態(tài)特征，為測(cè)試上述方法求解病態(tài)線性方程組的能力，采用100階希爾伯特方程組為例進(jìn)行測(cè)算。希爾伯特矩陣方程組如式(10)所示。

Ax=b

采用Oligo 7.0軟件設(shè)計(jì)引物。針對(duì)SSR位點(diǎn)查找的結(jié)果，在二、三、四、五、六重復(fù)單元的位點(diǎn)中分別設(shè)計(jì)出10對(duì)引物，共設(shè)計(jì)出50對(duì)引物，委托上海生工生物工程股份有限公司進(jìn)行合成。通過瓊脂糖凝膠電泳和變性聚丙烯酰胺凝膠電泳篩選出具多態(tài)性的引物。篩選獲得的多態(tài)引物將在其5’端添加羥基熒光素標(biāo)記（FAM），進(jìn)行熒光引物的合成。

(10)

圖2 100階希爾伯特矩陣的奇異值分布圖Fig.2 Singular value distribution of the 100-order Hilbert matrix

100階希爾伯特方程組經(jīng)奇異值分解后的主對(duì)角元素如圖2所示，從圖2中可以看出，奇異值元素非負(fù)且降序排列，曲線首支下降較快，尾支曲線近于水平且奇異值很小，病態(tài)問題非常嚴(yán)重，如果不對(duì)奇異值進(jìn)行處理，其求解效果幾乎無異于其它常規(guī)方法。通過給定不同的正則化因子可以削弱方程的病態(tài)程度，正則化因子λ通過λ(0)=1,λ(i)=λ(i-1)/100,i=1、2、…、9給出，方程組的求解誤差以最大絕對(duì)誤差ξmax=max(|1-xi|),i=1、2、…、100為評(píng)判標(biāo)準(zhǔn)。

圖3 不同正則化因子對(duì)應(yīng)的最大絕對(duì)誤差曲線Fig.3 Maximum absolute error curve corresponding to different regularization factors

圖4 不同正則化因子的100階希爾伯特方程組的解Fig.4 Solution of the 100-order Hilbert equations with different regularization factors

圖5 L-曲線圖Fig.5 L-curve

圖6 誤差曲線與 L-曲線的曲率Fig.6 Error curve and curvature of the L-curve

圖7 曲線的拐點(diǎn)區(qū)域局部放大曲線Fig.7 Partial magnification curve of the inflection point of the L-curve

給定不同的正則化因子并逐一求解希爾伯特方程組，解的誤差表現(xiàn)出明顯變化，如圖3所示。從圖3中可以看出,隨著正則化因子減小，解的最大誤差逐漸減小并達(dá)到最??；當(dāng)正則化因子繼續(xù)減小時(shí)，解的最大誤差又逐漸增大；當(dāng)正則化因子減小到很小如1E-18時(shí)，相當(dāng)于沒有做奇異值處理。給定不同的正則化因子并逐一求解希爾伯特方程組，其對(duì)應(yīng)的解如圖4所示。從圖中可以看出，當(dāng)正則化因子取為1E-4～1E-6時(shí)幾乎接近精確解，解的最大誤差基本在1E-2～1E-3內(nèi)變化。

3.2 正則化方法的測(cè)試與分析

3.2.1 曲率法計(jì)算結(jié)果

利用文獻(xiàn)[6]給出的最大曲率法對(duì)圖5的L-曲線計(jì)算曲率，繪制誤差曲線與曲率曲線疊合如6所示，可以看出,平均均方誤差最小時(shí)對(duì)應(yīng)的正則化因子在1E-12～1E-11之間，而最大曲率的對(duì)應(yīng)的正則化因子約為1E-15，兩者不重合，說明曲率法計(jì)算的正則化因子不是“最佳”正則化因子，為此需進(jìn)一步分析原因。

取圖5拐點(diǎn)區(qū)的數(shù)據(jù)重新繪制曲線如圖7所示，數(shù)據(jù)區(qū)的正則化因子在1E-10～1E-16范圍內(nèi)變化?？梢钥闯?，在L-曲線的斜率突變區(qū)，曲線并不平滑，導(dǎo)致樣條曲線擬合可能存在誤差，使得計(jì)算得到的最大曲率點(diǎn)出現(xiàn)偏差，沒有對(duì)應(yīng)“最佳”正則化因子。

表1 候選角點(diǎn)的坐標(biāo)(ρ,η)和正則化因子λ

圖8 逐步尋優(yōu)法搜索最佳角點(diǎn)的結(jié)果Fig.8 Results of searching the optimal corner point by stepwise optimization method

圖9 L-曲線拐點(diǎn)處放大的結(jié)果Fig.9 The result of amplification at the inflection point of the l-curve

3.2.2 分步選取法計(jì)算結(jié)果

第一步：在L-曲線上選取候選角點(diǎn)。L-曲線上的30個(gè)離散點(diǎn)構(gòu)成29個(gè)向量，選取模長(zhǎng)最長(zhǎng)的5個(gè)向量，根據(jù)相鄰兩向量間的角度大小關(guān)系判斷是否存在候選角點(diǎn)，并將搜索得到的候選角點(diǎn)添加至列表；增加選取的向量的個(gè)數(shù)，重復(fù)判斷角度大小關(guān)系來搜索候選角點(diǎn)。候選角點(diǎn)的坐標(biāo)如表1所示，候選角點(diǎn)的位置如圖9中紅色圓點(diǎn)標(biāo)識(shí)所示。

第二步：在候選角點(diǎn)中選取最優(yōu)點(diǎn)。表中所示5個(gè)候選角點(diǎn)與第1個(gè)離散點(diǎn)構(gòu)成5個(gè)向量，對(duì)向量作歸一化處理，判斷每個(gè)向量的垂直變化量與水平變化量的大小關(guān)系，如果垂直變化量大于水平變化量，那么認(rèn)為該向量具有垂直特征，否則認(rèn)為該向量具有水平特征。如果所有的向量都是具有水平特征，那么選擇候選角點(diǎn)列表中的最后一個(gè)候選角點(diǎn)作為L(zhǎng)-曲線的最佳角點(diǎn)，如果存在具有垂直特征的向量，那么選擇到達(dá)垂直特征向量前的候選角點(diǎn)作為最佳角點(diǎn)。從5個(gè)候選角點(diǎn)中搜索最佳角點(diǎn)為(ρ,η) = (-63.670 2, 4.605 17)，即對(duì)應(yīng)著最佳正則化因子，圖8中紫色方塊為最佳正則化因子的位置。由圖8和圖9可以看出，分步選取法搜索最佳角點(diǎn)的效果較好。

4 結(jié)論

通過對(duì)基于奇異值分解法求解病態(tài)線性方程組的正則因子分步選取方法的研究，得到以下幾點(diǎn)結(jié)論：

1)奇異值分解法是求解病態(tài)線性方程組的最有效的方法之一，通過給定最佳正則化因子，可以使病態(tài)線性方程組的求解結(jié)果接近真實(shí)解。

2)基于L-曲線給出了一種給定最佳正則化因子的分步選取方法，從測(cè)試結(jié)果來看，結(jié)果較好，并且優(yōu)于最大曲率法。

3)可以將分步選取方法與奇異值分解法相結(jié)合，用于求解小規(guī)模的地球物理反演問題。對(duì)于大規(guī)模地球物理反問題，由于奇異值分解法的分解過程的耗費(fèi)時(shí)間與方程階數(shù)呈指數(shù)遞增，在加快計(jì)算效率方面,仍有較大的改進(jìn)空間。