汪強(qiáng)強(qiáng), 劉海飛, 柳建新, 李 星, 施昕祎
(1.中南大學(xué) 地球科學(xué)與信息物理學(xué)院,長(zhǎng)沙 410083;2.有色資源與地質(zhì)災(zāi)害探查湖南省重點(diǎn)實(shí)驗(yàn)室,長(zhǎng)沙 410083)
地球物理反演問題大多數(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)證了本文算法的有效性。
地球反演中常因線性反演方程的嚴(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
所謂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)。
前面提出的奇異值分解求解病態(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)。 高階希爾伯特(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.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)的效果較好。 通過對(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)空間。3 數(shù)值實(shí)驗(yàn)分析
3.1 病態(tài)問題分析
3.2 正則化方法的測(cè)試與分析
4 結(jié)論