鄭史雄， 朱進波， 唐 煜， 郭俊峰

(1.西南交通大學土木工程學院，四川成都610031;2.西南石油大學土木工程與建筑學院，四川成都610500)

顫振理論認為橋梁顫振是一種氣動彈性不穩(wěn)定現(xiàn)象，當氣流經(jīng)過不規(guī)則橋梁斷面時，流體與結構產生相互激勵作用，隨著來流風速的增大，氣流對結構的輸送能量趨近于結構振動所耗散的能量，表現(xiàn)為運動耦合方程中系統(tǒng)振動阻尼由正值逼近0，而當風速超過一定值后，系統(tǒng)振動阻尼變?yōu)樨撝?，橋梁振動幅度加大并趨于運動發(fā)散，線性顫振發(fā)生.針對顫振現(xiàn)象的研究，Scanlan等［1-2］首先基于風洞模型試驗，引入顫振導數(shù)的概念，描述非定常氣動自激力，應用于振動微分方方程中，建立了經(jīng)典的自激力框架模型［3-4］.在此模型基礎上，許多研究者通過求解該框架模型，研究分析結構的顫振機理，先后提出了不同的二維顫振分析方法.

總體上，二維顫振分析方法可以分成兩類:一類求解方法是二自由度復模態(tài)特征值解法［5-6］，其基本思想是將二維顫振問題轉變?yōu)榍蠼鈴吞卣髦祮栴}，以某階模態(tài)下的阻尼比為0時，作為系統(tǒng)發(fā)散的依據(jù)，認為顫振發(fā)生，對應的風速為顫振臨界風速.該方法通過設置風速級數(shù)量和控制復頻率的收斂性誤差，可以得到方程組的精確解.為深入了解顫振機理，分析顫振導數(shù)在顫振過程中的貢獻，Matsumot等［7-8］對一系列簡單二維斷面進行了系統(tǒng)研究，提出了另一種求解方法，即顫振分步分析法(step by step analysis);項海帆和楊詠昕等［9-10］在此基礎上提出了描述顫振自由度參與程度的合理方法，能夠定量描述顫振形態(tài).該方法以某一自由度牽連運動阻尼比為0作為發(fā)散的根據(jù)，基于不同自由度間的激勵-反饋原理來解耦系統(tǒng)方程，能夠研究橋梁斷面的顫振驅動機理.

上述兩類方法都是頻率迭代方法，即在某級風速下，選取初始的頻率值代入頻率迭代方程，求得該級風速下的振動頻率與系統(tǒng)的阻尼，然后又逐級變化風速搜索直到出現(xiàn)阻尼比為0時的風速與頻率值.丁泉順等［11］已經(jīng)指出上述兩類方法結果的一致性.近年來，研究者分別對兩類方法進行了優(yōu)化發(fā)展，如葛耀君等［12］提出的二維顫振直接分析法就較好地提高了傳統(tǒng)復模態(tài)特征值解法的計算效率.

本文將遺傳混合算法應用于二維二自由度耦合顫振分步分析中，提出了基于遺傳混合算法的分析方法.該方法將二維顫振分析轉變?yōu)榍蠼夥蔷€性方程組的問題，引入具有自組織、自適應和自學性能的遺傳算法［13］，不需要人為選取初始頻率，自動搜索每級風速下全局最優(yōu)的各自由度振動頻率解，避免了迭代算法陷入局部收斂甚至不收斂的情況.同時，為了避免在搜索最優(yōu)解時的誤差，再引入良好的最優(yōu)算法——L-M 算法［14］——進行局部收斂修正.算例的計算結果驗證了該方法的可靠度和適用性.

1 二維二自由度耦合顫振分析方法

基于Scanlan的經(jīng)典自激力框架模型，首先給出自激力作用下二維橋梁節(jié)段僅做豎向和扭轉兩個方向的顫振運動方程:

式中:h、α分別為橋梁斷面在豎向和扭轉方向的位移;mh和I分別為結構豎向廣義質量和扭轉廣義質量慣矩;ξh0和ξα0分別為結構豎向和扭轉兩個自由度運動的結構阻尼比;ωh0和ωα0分別為結構豎向和扭轉兩個自由度運動的結構固有頻率，其中，廣義質量、廣義質量慣矩、結構固有頻率都為常數(shù);ρ為空氣密度;v為來流平均風速，一般可以視來流為均勻流，風速大小根據(jù)實際情況確定;B為橋梁實際斷面寬度;i=1，2，…，4)分別為無量綱氣動升力導數(shù)和氣動力矩導數(shù)，被稱之為顫振導數(shù);K為無量綱的折減頻率，K=Bω/v，ω為系統(tǒng)振動圓頻率.

求解上述方程時，二維二自由度耦合顫振分步分析方法將方程解耦為扭轉牽連運動方程和豎向牽連運動方程，其振動頻率分別對應為系統(tǒng)扭轉頻率和系統(tǒng)豎向運動頻率.此處扭轉牽連運動不僅包含扭轉自由度運動，還包括豎向自由度運動，豎向牽連運動也是如此.

當二維橋梁節(jié)段運動系統(tǒng)做單自由度扭轉運動時，運動方程形式變得格外簡單，如式(2).

式中:Mse(α，α)為扭轉運動自身所產生的氣動力矩.

式中:Lse(h，h)為豎向運動自身所產生的氣動力;Lse(h，α)為由豎向和扭轉運動間的激勵-反饋效應引起的氣動力.

因為顫振為自激發(fā)散振動，對于每個牽連運動方程，可以轉換為自由運動的形式，故將式(3)和式(4)牽連運動方程等式右端移到等式左端，結合為自由運動的形式，如式(5).

至此，成功地將一個耦合的振動系統(tǒng)解耦，系統(tǒng)是否發(fā)散取決于某一牽連運動方程的阻尼比.通過頻率與阻尼比迭代求解方程組(5).

系統(tǒng)豎向運動頻率為

系統(tǒng)扭轉牽連阻尼比為

系統(tǒng)豎向牽連阻尼比為

式(6)～(9)清楚地用顫振導數(shù)表達出了系統(tǒng)阻尼與系統(tǒng)剛度，且能夠反映不同狀態(tài)下系統(tǒng)阻尼所包含的運動系統(tǒng)自身阻尼和各分項氣動阻尼的影響貢獻.

無量綱系數(shù)為

不同耦合運動間的相位差角為

式(11)可描述為將系統(tǒng)運動分解為兩種牽連主運動，當系統(tǒng)以某一單自由度運動為結構主運動時，其運動方程稱為主運動方程，此時自激力框架下的另一方程稱為約束耦合方程，系統(tǒng)的某一單自由度運動會產生耦合氣動力，并以外荷載的形式通過約束耦合方程強制激起橋梁節(jié)段的另一自由度的運動，形成的耦合氣動力再反饋作用給主運動系統(tǒng)，構成了各自由度運動之間的激勵-反饋機制［15-16］.

在式(1)～(11)的基礎上，通過逐級搜索的方式找尋顫振臨界點，具體的做法如下:風速從初始風速開始搜索，每級風速下求解式(6)～(7)中的系統(tǒng)扭轉牽連運動頻率和豎向牽連運動頻率，計算并記錄系統(tǒng)的各自由度牽連阻尼比.隨著風速的增加，系統(tǒng)某一自由度阻尼比由正轉負時，即認為系統(tǒng)發(fā)生顫振，當前風速的大小即該系統(tǒng)發(fā)生顫振的臨界風速.值得一提的是，上述方法計算每級風速下的系統(tǒng)運動頻率是通過給定圓頻率初值代入式(6)和式(7)中，以不動點的迭代形式求得.根據(jù)上述原理，將二維二自由度耦合顫振分析法程序化，其計算流程見圖1.

圖1 傳統(tǒng)的分步分析法流程Fig.1 Flowchart of the traditional method

2 傳統(tǒng)數(shù)值解法的收斂性問題

2.1 構建非線性方程組

采用上述方法分析時，每級風速下的各自由度振動頻率都需要經(jīng)過式(6)和式(7)來迭代求得，以此得到啟發(fā)，將顫振問題的分析轉變?yōu)榍蠼庀铝卸蔷€性方程組的問題:

方程組(12)中顫振導數(shù)、各自由度相位差角和無量綱系數(shù)都是關于系統(tǒng)振動圓頻率的函數(shù)，故可以表示為

式中:X=(ωh，ωα).

求解上述形式的非線性方程組時，大量數(shù)值迭代法被提出，其主要思想是通過給定的初值代入迭代公式中逐次逼近準確值.文獻［17］中給出了很多求解迭代的方法，有基于不動點迭代的簡單迭代法，即原始耦合分析方法中使用的迭代法，本文中稱為原始方法，如式(6)和式(7)所示;有為了增加方程組求解收斂速度而提出的基于泰勒展開級數(shù)思想的牛頓法;同時也涌現(xiàn)出大量的改善方法，如修正牛頓法、擬牛頓法等.牛頓法及擬牛頓法等具有比簡單迭代法更高的收斂效率，完全適用于本文的求解過程，此類方法統(tǒng)一稱為傳統(tǒng)方法.

2.2 傳統(tǒng)數(shù)值解法的收斂性問題

使用傳統(tǒng)方法解非線性方程組時，首先要弄清楚方程組的性質.如線性方程組Ax=b，在理論上只要A－1存在，則解存在且唯一.而對于非線性方程組而言，情況變得尤為復雜，其可以有多解、唯一解或無解，最為重要的是迭代的收斂性與所取的閉區(qū)域有關，即傳統(tǒng)方法僅具有局部收斂性.

通過一個較簡單的非線性方程組來說明上述問題.

方程組(14)的真解為 x*=(0.7，0.7)，若取初值 x0=(0.1，0.1)作為初始點代入牛頓法［18］中收斂到點(1.550 7，－0.737 8)，迭代結果如圖 2 所示;若初始點取 x0=(0.9，0.9)或 x0=x*+0.5，則能夠收斂到精確解，迭代結果如圖3所示.

圖2 初值不在局部收斂域中的牛頓迭代法結果Fig.2 Results of the Newton iteration method with initial value not in the local convergence domain

上述類似方法都可以解決二維耦合顫振分析方法中每級風速下頻率的計算問題，不同迭代法的優(yōu)劣性也會因Jacobian矩陣的計算和函數(shù)求值的復雜性而不同，但都僅具有局部收斂性，初值選取不利會導致最終的迭代結果遠離真解的情況.而用選取不動點迭代方法時，將式(14)迭代形式變?yōu)?/p>

圖3 初值在局部收斂域中的牛頓迭代法結果Fig.3 Results of the Newton iteration method with initial value in the local convergence domain

方程組(15)為Jacoobi迭代格式，選取多個初值都做發(fā)散變化不曾收斂，其收斂過程如圖4所示.

圖4 Jacoobi迭代格式迭代結果Fig.4 Iteration results of Jacoobi iterative format

方程組(14)也可以轉變?yōu)?/p>

方程組(16)為Gauss-Seidel格式，選取真解附近的初始值(0.71，0.70)，迭代過程出現(xiàn)周期震蕩性的不收斂，如圖5所示.可見簡單迭代法的收斂性受迭代形式與格式的制約.

圖5 Gauss-Seidel格式迭代結果Fig.5 Iteration results of Gauss-Seidel iterative format

對于不復雜的非線性方程組，傳統(tǒng)的數(shù)值解法計算局限性都已經(jīng)很明顯了，而本文中要求解的方程組(12)很不規(guī)則，并且參數(shù)隨著風速變化，想要從理論上研究其在某個取值區(qū)域的收斂性質是十分困難的.嚴格來說，在方程性質尚不明確的條件下，采用局部收斂的迭代格式(傳統(tǒng)方法)進行求解是存在一定風險的，可能并不能收斂到真解.為規(guī)避這一風險，研究全局收斂的數(shù)值求解方法顯得很有必要.

3 遺傳混合算法

結合遺傳算法和傳統(tǒng)方法的優(yōu)點提出了遺傳混合算法，規(guī)避了傳統(tǒng)方法避免局部收斂甚至發(fā)散的情況.遺傳算法具有群體搜索和全局收斂性，不需要人為選取初始點，具有自行適應性，克服了傳統(tǒng)算法對初始點的敏感性問題.同時引入L-M算法對遺傳算法計算結果下的優(yōu)良個體進行局部搜索，克服遺傳算法局部搜索緩慢的問題.

3.1 遺傳算法

結合求解的非線性方程組，介紹遺傳算法［19］的基本步驟.首先結合方程組(13)將求解方程組等價于求解如下的最小值優(yōu)化問題:

式中:D為方程組解的區(qū)間.當f(ω)=0時，對應的ω即為方程組的解.通常利用目標函數(shù)來定義適應度函數(shù)，而適應度為非負值，適應度值越大對應的解越合理.所以對于求極小值問題，將目標函數(shù)變換為適應度函數(shù)，如式(18).

式中:C為正常數(shù).

對上述優(yōu)化問題實現(xiàn)遺傳算法求解時，需要先進行解的編碼與解碼.由于實數(shù)編碼容易過早收斂，本文采用具有穩(wěn)定性高、種群多樣性等優(yōu)點的二進制編碼，即將問題的解即系統(tǒng)扭轉頻率值分別空間映射到位串空間Bl上，其中，下標l為位數(shù)，Bl取0或1，每串編碼上的每一位稱為基因，所在位置稱為基因位.在此基礎上進行遺傳算法操作，得到的解經(jīng)過解碼還原成實數(shù)解并進行適應度評價.

遺傳算法的基本步驟及解釋:

(1)隨機產生N個個體構成二進制編碼的初始群體，表示為 P(0)={ω1，ω2，…，ωN}，其中ωi={ωhi，ωαi}，i=1，2，…，N，ωi被形象地稱為染色體，其上基因的組合構成染色體不同的表現(xiàn)型.

(2)根據(jù)目標函數(shù)f(ω)計算當前群體中各個體的適應度，即ffit(ωi).

(3)判斷算法終止條件是否滿足，若滿足則轉(8).

(4)根據(jù)各個體的適應度執(zhí)行選擇操作，本文采用基于各個體適應度概率決定它們繼續(xù)繁衍還是消亡的輪盤式選擇方法.

個體的選擇概率為pi=ffit(ωi) ∑ffit(ωi)，，生成一個［0，1］內的隨機數(shù) r，若 p1+p2+…+pi－1＜r＜p1+p2+…+pi，則選擇個體 i.

(5)按交叉概率pc執(zhí)行交叉操作.將選擇操作后的某兩個父代個體的部分片段加以替換重組形成新的子代，在二進制編碼中選擇交叉點，體現(xiàn)為部分字符串的互換.

(6)按變異概率pm執(zhí)行變異操作.以一定的概率選擇父代個體的某一基因座，通過改變該基因座的基因變化為新個體，在二進制編碼中體現(xiàn)為某段字符串的取反，即若原來為0，則變?yōu)?，若為1則取0.

(7)若已得到由N個新個體構成的新一代群體，則轉(2)，否則轉(4).

(8)輸出搜索結果，終止.

綜上，遺傳算法［19］可以形式化描述為

式中:l為二進制編碼長度，取決于解的取值范圍和精度要求;

s為選擇、交叉和變異的策略方法;

g為遺傳算子，通常包括交叉算子和變異算子;

p為遺傳算子的執(zhí)行概率，通常包括交叉概率和變異概率;

f為適應度函數(shù);

t為迭代次數(shù)即最大進化代數(shù)或者終止準則.

依照上述流程，根據(jù)相關經(jīng)驗選取種群規(guī)模、染色體長度、最大進化代數(shù)、交叉概率和變異概率等遺傳參數(shù)，能夠有效地找出非線性方程組(12)的全局最優(yōu)解.

3.2 L-M 算法

遺傳算法雖有很好的全局收斂性，但局部搜索能力較差，在接近準確解時，搜索進程變得緩慢甚至停止.本文引入基于約束最小二乘思想的L-M算法對遺傳算法求出的近似解進行局部修正.首先與方程組(17)相同，確定求解目標函數(shù):式中:F(ω)=(f1(ω)，f2(ω))T，ω=(ωh，ωα)T.

設已知第 k次的迭代點ωk，則由泰勒公式可知:Δ

fu(ω)≈fu(ωk)+(fu(ωk))T(ω－ωk)，u=1，2.所以 F(ω)≈F(ωk)+A(ωk)(ω－ωk).其中:

記 Ak=A(ωk)，則有

式中:dk=ω－ωk.

記 φ(ω)=［Akdk+F(ωk)］T［Akdk+F(ωk)］，那么求解φ(ω)的極小點則近似原問題的極小點.而為克服雅克比矩陣的奇異致線性搜索得不到進一步下降，故將目標函數(shù)(20)轉為如式(23)的信賴域模型.

式中:hk為信賴域半徑.

這個方程的解可由解如式(24)得到.

適當調整 λk，保持 AT(ωk)A(ωk)+λkI正定，從而

若 (AT(ωk)A(ωk))－1AT(ωk)F(ωk) ≤ hk，則 λk=0;否則 λk＞0.由于 AT(ωk)A(ωk)+λkI正定，從而式(25)產生的方向一直是下降方向，最終能夠搜索到精確解.

該算法結合傳統(tǒng)高斯-牛頓法和梯度下降法的優(yōu)點，同時增加對角元改善雅克比矩陣奇異性，避免了傳統(tǒng)的Gauss-Newton算法不迭代情況，其快速的局部搜索能力得以充分發(fā)揮，具有很強的適用性.

3.3 遺傳混合算法的應用

遺傳算法擁有強大的全局搜索能力，其理論比較成熟且得到了數(shù)學上的證明，已然逐步應用到工程實踐領域，而L-M算法克服了雅克比矩陣奇異的問題且擁有較快的收斂速度，兩種方法組合的遺傳混合算法成為一種可行高效的方法.基于遺傳混合算法的二維二自由度耦合分步分析方法的思路如下:

(1)以系統(tǒng)振動頻率為未知數(shù)，根據(jù)式(10)和式(11)分別表示各個相位角與無量綱量.

(2)以折算頻率表示各顫振導數(shù).若風速已知，則顫振導數(shù)僅與系統(tǒng)振動頻率有關.

(3)輸入風速v.

(4)求解關于系統(tǒng)頻率的二元非線性方程組(12).首先應用遺傳算法全局搜索最優(yōu)近似解，其次將近似解作為L-M算法迭代的初始點做局部精細搜索.

(5)將求解的系統(tǒng)頻率精確解代入式(8)和式(9)中計算系統(tǒng)阻尼比.若阻尼比都大于0，則風速增加一個定值Δv，返回(3)，重新計算;若某一阻尼比等于0，則程序終止，對應風速為顫振臨界風速.

4 數(shù)值算例

為驗證上述顫振直接分析方法的可靠性和適用性，用該方法對一平板的耦合顫振問題進行分析，并將分析結果與傳統(tǒng)不動點迭代分析方法的結果進行比較.

二自由度平板斷面基本參數(shù)如下:豎向振動的固有頻率為0.120 1 Hz，扭轉振動的的固有頻率為0.280 1 Hz，斷面寬度 B=38.8 m，每米長度廣義質量m=22 600 kg/m，廣義質量慣性矩 Im=3 667 100 kg·m2/m，空氣密度 ρ=1.225 kg/m3，結構豎向和扭轉振動阻尼比均取為0.005，平板的顫振導數(shù)采用Theodorsen解［20］.分別用上述兩種方法進行運算，顫振分析結果列于表1，隨風速變化的系統(tǒng)振動頻率與阻尼比如圖6～7所示.

從圖6～7和表1中數(shù)據(jù)可知，在取定的多個檢測風速節(jié)點處，本文方法與傳統(tǒng)方法的計算結果誤差都低于0.1‰，結果幾乎一致，對應的顫振臨界風速完全一致.而本文建立的遺傳混合算法，由于具有全局收斂的先天優(yōu)勢，故更適宜應用于大跨度橋梁的顫振分析中.

表1 橋梁斷面顫振分析結果對比Tab.1 Comparison results of the bridge flutter analysis

圖6 隨風速變化的系統(tǒng)振動圓頻率比較Fig.6 System vibration frequency with the changing wind speed

圖7 隨風速變化的系統(tǒng)牽連阻尼比比較Fig.7 System damping ratio with the changing wind speed

5 結 論

通過對傳統(tǒng)的二維二自由度耦合分步分析方法中數(shù)學問題的研討，發(fā)現(xiàn)應用其中的數(shù)值迭代方法存在局部收斂性的問題.

將二維二自由度耦合顫振分析轉變?yōu)榍蠼怅P于系統(tǒng)振動頻率的非線性方程組問題，使得求解的方式多元化.引入了不需要自行選取初始值的遺傳混合算法進行全局搜索和局部精細收斂運算，建立基于遺傳混合算法的二維耦合顫振分步分析方法.

理論推導與算例分析均表明，本文所提出的基于遺傳混合算法的二維耦合顫振分步分析方法具有精度高、無條件收斂的優(yōu)點.

［1］ SCANLAN R H，TOMKO J J.Airfoil and bridge deck flutter derivatives［J］. Journal of Engineering Mechanics，1971，97(6):1717-1737.

［2］ HUSTON D R.Flutter derivatives from 14 generic deck sections［C］∥ Proc. ofASCE Structure Congress.Orlando:ASCE，1987:281-292.

［3］ 陳政清，于向東.大跨橋梁顫振自激力的強迫振動法研究［J］.土木工程學報，2002，35(5):34-41.CHEN Zhenqing，YU Xiangdong.A new method for measuring flutter self-excited forces of long-span bridges［J］.China Civil Engineering Journal， 2002，35(5):34-41.

［4］ 李志國，王騎，伍波，等.不同攻角下扁平箱梁顫振機理研究［J/OL］.西南交通大學學報，優(yōu)先出版.(2017-09-27) ［2017-10-28］.http://kns.cnki.net/kcms/detail/51.1277.U.20170927.2111.034.html.［5］ WILDE K，F(xiàn)UJINO Y，MASUKAWA J.Time domain modeling of brige deck flutter［J］. Structural Engineering Earthquake Engineering，1996，13(2):93-104.

［6］ XIANG H，ZHANG R.On mechanism of flutter and unified flutter theory of bridges［M］∥Wind Engineering Into the 21st Century.Rotterdam: ［s.n.］，1999:1069-1074.

［7］ MATSUMOTO M，HAMASAKI H，YOSHIZUMI F.On flutter stability of decks for super long-span bridge［J］.Structural Engineering EarthquakeEngneering， 1997，14:185-200.

［8］ MATSUMOTO M. Flutter instability and recent development in stabilization of structures［J］.Journal of Wind Engineering and Industrial Aecodynamics，2007，95(9/10/11):888-907.

［9］ 楊詠昕.大跨度橋梁二維顫振機理及其應用研究［D］.上海:同濟大學，2002.

［10］ 楊詠昕，葛耀君，項海帆.大跨度橋梁典型斷面顫振機理［J］.同濟大學學報，2006，34(4):455-460.YANG Yongxin，GE Yaojun，XIANG Haifan.Flutter mechanism of representative sections for long-span bridges［J］.JournalofTongjiUniversity， 2006，34(4):455-460.

［11］ 丁泉順，朱樂東.橋梁主梁斷面氣動耦合顫振分析與顫振機理研究［J］.土木工程學報，2007，40(3):69-74.DING Quanshun， ZHU Ledong. Aerodynamically coupling flutter analysis and flutter mechanism for bridge deck sections［J］.ChinaCivilEngineering Journal，2007，40(3):69-74.

［12］ 邵亞會，葛耀君，柯世堂.超大跨度懸索橋二維顫振頻域直接分析方法［J］.哈爾濱工業(yè)大學學報，2011，43(8):119-123.SHAO Yahui，GE Yaojun，KE Shitang.Straight forward method for two-dimensional flutter analysis of superlong-span suspension bridge in frequency domain［J］.Journal of Harbin Institute of Technology University，2011，43(8):119-123.

［13］ 羅亞中，袁端才，唐國金.求解非線性方程組的混合遺傳算法［J］.計算力學學報，2005，22(1):109-114.LUO Yazhong，YUAN Duancai，TANG Guojin.Hybrid genetic algorithm for solving systems of nonlinear equations［J］. Chinese Journal of Computational Mechanics，2005，22(1):109-114.

［14］ 張鴻燕，耿征.Levenberg-Marquard算法的一種新解釋［J］.計算機工程與應用，2009，45(19):5-9.ZHANG Hongyan，GENG Zheng.Novel interpretation for Levenberg-Marquardtalgorithm［J］.Computer Engineering and Applications，2009，45(19):5-9.

［15］ 項海帆，葛耀君，朱樂東，等.現(xiàn)代橋梁抗風理論與實踐［M］.北京:人民交通出版社，2005:206-216.

［16］ 葛耀君.大跨度懸索橋抗風［M］.北京:人民交通出版社，2011:259-269.

［17］ 朱小飛.迭代法解非線性方程組［D］.合肥:合肥工業(yè)大學，2014.

［18］ 晁玉翠.求解非線性方程組的修正牛頓法［D］.哈爾濱:哈爾濱工業(yè)大學，2007.

［19］ 屈穎爽.基于擬牛頓法和遺傳算法求解非線性方程組的混合算法［D］.長春:吉林大學，2005.

［20］ 陳政清.橋梁風工程［M］.北京:人民交通出版社，2005:67-73.