基于分形理論的雷電先導(dǎo)建模與仿真畢設(shè)論文終極版

1、<p>　　基于分形理論的雷電先導(dǎo)建模與仿真</p><p><b>　　摘要</b></p><p>　　雷電是一種自然放電現(xiàn)象，它是由雷云(帶電的云層)對地面建筑物及大地的自然放電引起的。雷電放電的隨機(jī)性以及雷電先導(dǎo)在其產(chǎn)生和發(fā)展過程中易受周圍環(huán)境影響的特點導(dǎo)致雷電活動對電力系統(tǒng)產(chǎn)生巨大危害。由于雷電放電的隨機(jī)性，需要用統(tǒng)計學(xué)的方法來進(jìn)行防雷性能的評估。

2、為此采用分形模型模擬了先導(dǎo)的發(fā)展過程，在二維和三維中描述雷電放電通道的分支和扭曲。通過比較自然界中雷電放電和仿真結(jié)果的分形維數(shù)，可以找到最接近實際的分形系數(shù)。與傳統(tǒng)的電氣幾何模型、改進(jìn)電氣幾何模型以及先導(dǎo)發(fā)展模型對比，分形模型不僅能夠給出逼真的先導(dǎo)通道圖形，而且能夠從統(tǒng)計學(xué)意義上對一些不合理的側(cè)面繞擊現(xiàn)象作出解釋。</p><p>　　本課題基于對雷電發(fā)展過程準(zhǔn)靜態(tài)場的假設(shè)，采用有限差分法作為電場計算的方法，將偏

3、微分方程的問題化為線性方程組問題，并用迭代法進(jìn)行線性方程的數(shù)值解。結(jié)合以上計算出的空間電場，采用分形方法完成了二維和三維雷電先導(dǎo)放電通道的建模與仿真，并實現(xiàn)了雷電先導(dǎo)發(fā)展過程及相應(yīng)的空間電勢圖形靜態(tài)和動態(tài)顯示。在此基礎(chǔ)上，本課題還討論了避雷針(線)對閃電軌跡的影響。</p><p>　　關(guān)鍵詞：分形模型；三維建模；雷電先導(dǎo)；有限差分法；MATLAB</p><p>　　Modeling a

4、nd Simulation of Lightning Leader based on the Fractal Theory</p><p><b>　　Abstract</b></p><p>　　Lightning is a natural discharge phenomenon, which is caused by thundercloud’s (electr

5、ically charged clouds) natural discharge to ground-level structures and earth. Lightning causes great damage to power system due to the stochastic nature of lightning discharges and the vulnerability of lightning leaders

6、 to various environments in the process of inception and propagation. As a result of the randomness of lightning discharge, it’s necessary to introduce statistical methods into the evaluation of</p><p>　　Thi

7、s topic is based on the assumption of quasi-static field in lightning development process, using the finite difference method to calculate the electric field, changing partial differential equations into linear equations

8、, using the iterative method to get the numerical solution of linear equations. Combining with the above electric field, this topic uses the fractal theory to model and simulate the two-dimensional and three-dimensional

9、lightning leader discharge channel, achieving development</p><p>　　Keywords: Fractal model; Three-dimensional modeling; Lightning leader; Finite difference method; MATLAB</p><p><b>　　目 錄&l

10、t;/b></p><p><b>　　摘要I</b></p><p>　　AbstractII</p><p><b>　　第1章 緒論1</b></p><p>　　1.1 課題背景1</p><p>　　1.2 國內(nèi)外研究動態(tài)2</p>

11、<p>　　1.3 本文主要內(nèi)容和難點3</p><p>　　1.3.1 工作內(nèi)容4</p><p>　　1.3.2 課題難點4</p><p>　　第2章 有限差分法和格點式隨機(jī)分形模型5</p><p>　　2.1 有限差分法概述5</p><p>　　2.2 DBM概述7</p>

12、;<p>　　2.1.1 先導(dǎo)發(fā)展規(guī)則7</p><p>　　2.1.2 分形先導(dǎo)發(fā)展模型關(guān)鍵參數(shù)的選取8</p><p>　　第3章 雷電先導(dǎo)放電通道的二維模型10</p><p>　　3.1 雷電先導(dǎo)放電的物理過程及參數(shù)10</p><p>　　3.1.1雷電放電概述10</p><p>

13、　　3.3 程序流程13</p><p>　　3.3.1 模塊114</p><p>　　3.3.2 模塊214</p><p>　　3.3.3 模塊316</p><p>　　3.4 動畫效果和電勢的變化19</p><p>　　3.4.1 動畫效果19</p><p>　　3.4

14、.2 電勢的變化21</p><p>　　第4章 雷電先導(dǎo)放電通道的三維模型22</p><p>　　4.1 模型的建立與簡化22</p><p>　　4.2 程序流程22</p><p>　　4.2.1 模塊1 初始情況22</p><p>　　4.2.2 模塊2 計算電勢23</p>&

15、lt;p>　　4.2.3 模塊3 計算概率，決定下一個擊穿點24</p><p>　　4.3 動畫效果和電勢的變化24</p><p>　　4.3.1 動畫效果24</p><p>　　4.3.2 電勢的變化25</p><p>　　第5章 避雷針(線)對先導(dǎo)通道發(fā)展的影響28</p><p>　　5

16、.1 避雷針對先導(dǎo)通道發(fā)展的影響28</p><p>　　5.1.1 二維模擬避雷針對雷電先導(dǎo)的影響28</p><p>　　5.1.2 三維模擬避雷針對雷電先導(dǎo)的影響29</p><p>　　5.2 避雷線對架空線路的保護(hù)作用32</p><p>　　第6章 結(jié)論與展望36</p><p><b&g

17、t;　　6.1 結(jié)論36</b></p><p><b>　　6.2 展望36</b></p><p><b>　　參考文獻(xiàn)37</b></p><p><b>　　致謝39</b></p><p><b>　　第1章 緒論</b>&

18、lt;/p><p><b>　　1.1 課題背景</b></p><p>　　雷電是伴有閃電和雷鳴的一種雄偉壯觀而又有點令人生畏的放電現(xiàn)象[]，是發(fā)生在大氣層中的聲、光、電物理現(xiàn)象，常見的雷電現(xiàn)象是一部分帶電云層與另一部分帶異種電荷的云層與大地之間的迅猛的放電過程[]。</p><p>　　空中的塵埃、冰晶等物質(zhì)在云層中翻滾運動的時候，經(jīng)過一些復(fù)雜

19、過程，使這些物質(zhì)分別帶上了正電荷與負(fù)電荷。經(jīng)過運動，帶上相同電荷的質(zhì)量較重的物質(zhì)會到達(dá)云層的下部(一般為負(fù)電荷)，帶上相同電荷的質(zhì)量較輕的物質(zhì)會到達(dá)云層的上部(一般為正電荷)。這樣，同性電荷的匯集就形成了一些帶電中心，當(dāng)異性帶電中心之間的空氣被其強大的電場擊穿時，就形成“云間放電”(即閃電)。</p><p>　　帶負(fù)電荷的云層向下靠近地面時，地面的凸出物、金屬等會被感應(yīng)出正電荷，隨著電場的逐步增強，雷云向下形成

20、下行先導(dǎo)，地面的物體形成向上閃流，二者相遇即形成對地放電。這就容易造成雷電災(zāi)害。雷電災(zāi)害是“聯(lián)合國國際減災(zāi)十年”公布的最嚴(yán)重的十種自然災(zāi)害之一。最新統(tǒng)計資料表明，雷電造成的損失已經(jīng)上升到自然災(zāi)害的第三位。全球每年因雷擊造成人員傷亡、財產(chǎn)損失不計其數(shù)。據(jù)不完全統(tǒng)計，我國每年因雷擊以及雷擊負(fù)效應(yīng)造成的人員傷亡達(dá)3000～4000人，財產(chǎn)損失在50億元到100億元人民幣。</p><p>　　雷電災(zāi)害所涉及的范圍幾乎遍

21、布各行各業(yè)?，F(xiàn)代電子技術(shù)的高速發(fā)展，帶來的負(fù)效應(yīng)之一就是其抗雷擊浪涌能力的降低。以大規(guī)模集成電路為核心組件的測量、監(jiān)控、保護(hù)、通信、計算機(jī)網(wǎng)絡(luò)等先進(jìn)電子設(shè)備廣泛運用于電力、航空、國防、通信、廣電、金融、交通、石化、醫(yī)療以及其它現(xiàn)代生活的各個領(lǐng)域，以大型CMOS集成元件組成的這些電子設(shè)備普遍存在著對暫態(tài)過電壓、過電流耐受能力較弱的缺點，暫態(tài)過電壓不僅會造成電子設(shè)備產(chǎn)生誤操作，也會造成更大的直接經(jīng)濟(jì)損失和廣泛的社會影響。</p>

22、<p>　　雷電的發(fā)生是不能避免，不能阻止的，因為這是自然現(xiàn)象。但是，我們可以設(shè)法避雷，設(shè)法使雷電天氣對我們不造成災(zāi)害。在雷雨天氣，我們有時會—看到一些些參天大樹被雷電擊倒，而周圍的一些高塔、高樓等高層建筑卻安然無恙。這是因為這些參天大樹由于受到帶大量電荷的云層的感應(yīng)，也帶上大量電荷，積累的電荷過多時，樹就被擊倒。而同樣倩況下．高層建筑安然無恙恐怕要歸功于避雷針了。在許多高塔上都有于個金屬做的、狀如針頭向上直立的東西，這就

23、是避雷針。原來，避雷針應(yīng)叫“引雷針”。在雷雨天氣，高樓上空出現(xiàn)帶電云層時，迅雷針和高樓頂部都被感應(yīng)上大量電荷，由于避雷針針頭是尖的，而靜電感應(yīng)時，導(dǎo)體尖端總是聚集了最多的電荷。這樣，避雷針就聚集了大部分電荷。避雷針又與這些帶電云層形成了一個電容器，由于它較尖，即這個電容器的兩極板正對面積很小，電容也就很小，也就是說它所能容納的電荷很少。而它又聚集了大部分電荷，所以，當(dāng)云層上電荷較多時，避雷針與云層之間的空氣就很容易被擊穿，成為導(dǎo)體。這樣

24、，帶電云層與避雷針形成通路，而避雷針又是接地的。避雷針就可以把云層上的電荷導(dǎo)人大地，使其不對高層建筑構(gòu)成危險，保證了它的安全。</p><p>　　國內(nèi)外超特高壓輸電線路的運行統(tǒng)計結(jié)果表明，雷擊跳閘已經(jīng)成為特高壓線路故障的主要原因，而繞擊則是造成雷擊跳閘的主要因素。另外，超高壓線路雷擊跳閘故障中很大部分也是由繞擊造成的[]。隨著特高壓電網(wǎng)的建設(shè)，如何防護(hù)輸電線路的繞擊逐漸成為線路防雷分析中需要考慮的主要問題，而如

25、何正確的描述雷電的發(fā)展過程又決定了防雷分析的準(zhǔn)確性和合理性。自然界的雷電放電是一種超長間隙的放電現(xiàn)象。由于其范圍巨大，其發(fā)展過程中會受到各種不確定因素的影響，如自然界的背景噪聲、空間電場、氣象條件等，從而呈現(xiàn)出彎曲分叉的隨機(jī)放電路徑，給研究者和工程師研究雷擊過程帶來了很大的困難。目前我國輸電線路的雷電防護(hù)設(shè)計及其性能估算以我國電力行業(yè)標(biāo)準(zhǔn)DL/T620-1997《交流電氣裝置的過電壓保護(hù)和絕緣配合》為依據(jù)。規(guī)程法的繞擊計算未考慮雷電發(fā)展

26、過程、雷電流大小和地面傾角對屏蔽效果的影響，而是根據(jù)經(jīng)驗和小電流試驗?zāi)Ｐ驮囼灲Y(jié)果提出的綜合平均法，常不能反映線路的具體特性。因此，非常需要找到一種最合理的雷電分析模型。</p><p>　　1.2 國內(nèi)外研究動態(tài)</p><p>　　目前，輸電線路雷擊屏蔽模型主要有：電氣幾何模型(EGM)[、]、改進(jìn)電氣幾何模型[、]、先導(dǎo)發(fā)展模型(LPM)[、、]。</p><p&g

27、t;　　對于EGM模型及改進(jìn)EGM模型，將雷電流幅值與雷擊距離聯(lián)系起來，將雷擊輸電線路用幾何作圖的方法來描述，繞擊率與雷電流幅值相關(guān)，能夠考慮線路結(jié)構(gòu)和雷電參數(shù)等對繞擊率的影響，其通過計算雷電擊距，以繞擊弧的方式對輸電線路繞擊性能進(jìn)行評估，但該方法沒有考慮雷電放電的分散性及線路運行電壓對線路繞擊耐雷性能的影響。雖然在不同線路參數(shù)及環(huán)境下的應(yīng)用中提出了一些修正的方法，但其防雷設(shè)計與超特高壓輸電線路實際運行結(jié)果仍存在較大的誤差[]。但該方法

28、的理論是由保護(hù)角較大、桿塔高度較低的線路運行經(jīng)驗總結(jié)而來的，在超高壓線路設(shè)計時，理論計算的繞擊率總是遠(yuǎn)小于實際運行值。如20世紀(jì)90年代中后期，華東地區(qū)連年屢次發(fā)生雷害事故，而遭受雷害事故的桿塔多為ZM1型桿塔，其邊相導(dǎo)線保護(hù)角為7.2°，滿足EGM的完全屏蔽條件。根據(jù)現(xiàn)場調(diào)查情況，經(jīng)分析認(rèn)為：由于放電的分散性，超出EGM繞擊區(qū)定位的雷閃仍可能發(fā)生繞擊[3]，這一現(xiàn)象顯然無法用EGM來解釋。</p><p&

29、gt;　　LPM模型基于實驗室試驗的基礎(chǔ)上，闡述了雷電發(fā)展的整個過程，是較EGM模型更為先進(jìn)的雷電屏蔽模型，LPM模型提出輸電線路繞擊側(cè)面距離及屏蔽失效寬度這兩個基本參數(shù)使得模型在對輸電線路繞擊耐雷性能的評估有較大的適用性，該模型考慮了上下行先導(dǎo)相對發(fā)展的實際過程，是對雷電發(fā)展過程的近似模擬，但對雷電先導(dǎo)的發(fā)展僅考慮了雷電放電的確定性因素，認(rèn)為雷電先導(dǎo)一直沿先導(dǎo)發(fā)展空間最大場強方向發(fā)展，忽略了隨機(jī)性因素，其模型中用于表征輸電線路繞擊耐雷

30、水平的屏蔽失效寬度這一參數(shù)也無法體現(xiàn)出雷云放電的分散性。因此傳統(tǒng)的LPM模型及其線路繞擊耐雷性能評估方法具有一定的局限性。國內(nèi)學(xué)者基于雷電放電和長空氣間隙放電物理過程的相似性，建立了實驗室的等比例縮小輸電線路雷擊模型[、]，通過輸電線路繞擊概率曲線的形式對輸電線路繞擊耐雷性能進(jìn)行研究。但實驗室試驗中長間隙放電的先導(dǎo)發(fā)展并不充分，如何將實驗室數(shù)據(jù)應(yīng)用到實際的防雷設(shè)計還有待深入地研究。</p><p>　　因此，尋求

31、一種能夠充分考慮雷電放電確定性及隨機(jī)性因素的輸電線路繞擊屏蔽模型及其繞擊耐雷性能評估方法對特高壓的防雷設(shè)計及運行至關(guān)重要?？紤]到雷電先導(dǎo)的發(fā)展呈現(xiàn)樹狀和狹窄的隨機(jī)放電通道結(jié)構(gòu)[]，該結(jié)構(gòu)有明顯的自相似特征，而分形正是描述該結(jié)構(gòu)的最好數(shù)學(xué)方法。分形理論是當(dāng)今世界十分風(fēng)靡和活躍的新理論、新學(xué)科。分形概念是美籍?dāng)?shù)學(xué)家曼德布羅特(B.B.Mandelbort)首先提出的。分形理論既是非線性科學(xué)的前沿和重要分支，又是一門新興的橫斷學(xué)科。作為一種方

32、法論和認(rèn)識論，其啟示是多方面的：一是分形整體與局部形態(tài)的相似，啟發(fā)人們通過認(rèn)識部分來認(rèn)識整體，從有限中認(rèn)識無限；二是分形揭示了介于整體與部分、有序與無序、復(fù)雜與簡單之間的新形態(tài)、新秩序；三是分形從一特定層面揭示了世界普遍聯(lián)系和統(tǒng)一的圖景。</p><p>　　分形先導(dǎo)發(fā)展模型是一種基于分形理論的雷電先導(dǎo)模型，它在傳統(tǒng)雷電先導(dǎo)模型的基礎(chǔ)上，充分考慮雷電上下行先導(dǎo)路徑的分形特性以及先導(dǎo)發(fā)展過程的隨機(jī)性和確定性因素，將

33、LPM模型的理論基礎(chǔ)和分形理論相結(jié)合而建立的一種基于分形理論的超特高壓輸電線路耐雷性能評估方法。Niemeyer,Pientronero和Wielsmann于1983年研究了平行玻璃之間SF6中的電擊現(xiàn)象，他們得到了多枝杈的徑向放電圖形(1µs的30kV電脈沖)，他們對圖形進(jìn)行了計算機(jī)模擬，提出了介電擊穿模型(dielectric breakdown model，DBM)。繼Niemeyer等人提出DBM模型后[]，不少研究者

34、又將其用在對雷電閃道的模擬上[、、、、、]。分形模擬不僅其結(jié)果能夠體現(xiàn)通道的分岔和彎曲，而且其過程也與雷電放電過程基本類似。閃電現(xiàn)象屬于氣體電介質(zhì)擊穿現(xiàn)象，用類似DBM模型的方法進(jìn)行計算機(jī)模擬是比較自然的想法，因此一些學(xué)者以此為依據(jù)來模擬自然界中的閃電現(xiàn)象。金秀儒的《用計算機(jī)模擬閃電形成的嘗試》即為一次不錯的嘗試[]，但作者最終沒有得出一個完善的模擬。中國科學(xué)技術(shù)大學(xué)物理2班的王寧遠(yuǎn)在他的基礎(chǔ)上更近一步[]，用</p>&

35、lt;p>　　1.3 本文主要內(nèi)容和難點</p><p>　　本課題的主要任務(wù)是基于MATLAB軟件編程，采用分形方法實現(xiàn)二維和三維雷電先導(dǎo)放電通道的建模與仿真，充分了解介質(zhì)電擊穿的具體原理以及性質(zhì)，在此基礎(chǔ)上利用其核心公式進(jìn)行條件的確立與簡化，找到合適的邊界條件并進(jìn)行建模。模型首先通過在二維空間和三維空間劃分網(wǎng)格把實際問題離散化，然后通過適當(dāng)?shù)乃惴▽﹂W電形成過程進(jìn)行模擬。MATLAB編程時采用有限差分法不

36、斷更新電位，并結(jié)合隨機(jī)函數(shù)使程序不停地尋找生長點，直到閃電擊穿到地面。</p><p>　　1.3.1 工作內(nèi)容</p><p>　　(1)充分了解介質(zhì)電擊穿的具體原理以及性質(zhì)[]；</p><p>　　(2)理解有限差分法的具體原理[23]；</p><p>　　(3)學(xué)會使用MATLAB軟件，尤其是MATLAB程序設(shè)計與圖形處理[]；&l

37、t;/p><p>　　(4)編寫二維程序，實現(xiàn)雷電先導(dǎo)打到地面或邊界后程序結(jié)束運行；</p><p>　　(5)編寫三維程序，實現(xiàn)雷電先導(dǎo)打到地面或邊界后程序結(jié)束運行；</p><p>　　(6)程序中添加雷電擊中地面的時間；</p><p>　　(7)實現(xiàn)雷電的動態(tài)效果；</p><p>　　(8)通過圖形動態(tài)顯示電勢的

38、變化；</p><p>　　(9)研究避雷針（線）對雷電先導(dǎo)的影響；</p><p>　　(10)研究一根輸電線上場強的變化。</p><p>　　1.3.2 課題難點</p><p>　　(1)對介質(zhì)電擊穿的核心公式進(jìn)行條件的確立與簡化，找到合適的邊界條件；</p><p>　　(2)如何使程序不停地尋找生長點；&l

39、t;/p><p>　　(3)編程實現(xiàn)動態(tài)效果；</p><p>　　(4)編程實現(xiàn)三維電勢的動態(tài)變化。</p><p>　　第2章 有限差分法和格點式隨機(jī)分形模型</p><p>　　2.1 有限差分法概述</p><p>　　經(jīng)典場的邊值問題在數(shù)學(xué)上表達(dá)為泊松方程和拉普拉斯方程，但解偏微分方程往往是比較困難的[23]。

40、但是很多時候?qū)τ诰唧w問題我們需要的不是解析解，而是數(shù)值解，所以可以考慮用連續(xù)變量離散化的方法求出數(shù)值解，在足夠的精度上進(jìn)行逼近，這就引出了有限差分法。</p><p><b>　　有限差分法：</b></p><p>　　微分： (2-1)</p><p>　　用有限的h

41、代替，使得</p><p><b>　　(2-2)</b></p><p><b>　　差分的種類：</b></p><p><b>　　一階差分：</b></p><p><b>　　或者或者</b></p><p><b&

42、gt;　　二階差分：</b></p><p><b>　　(2-3)</b></p><p>　　設(shè)U(x,y,z)為空間電勢的函數(shù)。</p><p><b>　　泊松方程：</b></p><p><b>　　(2-4)</b></p><p&

43、gt;　　使用二階差分代替泊松方程中的拉普拉斯算符，有：</p><p><b>　　(2-5)</b></p><p>　　Σ表示分別對三個變元求差分之和，以下相同。</p><p>　　矩陣(數(shù)組)是計算機(jī)中重要的數(shù)據(jù)結(jié)構(gòu)，為了方便用矩陣去存儲數(shù)據(jù)，我們網(wǎng)格去劃分空間，從而不僅使方程化為簡單的有限差分形式，而且這樣的數(shù)據(jù)類型在計算機(jī)中易于儲

44、存和運算。那么h=k=l=1，并且令f(x,y,z)=u(x,y,z)就有：</p><p><b>　　(2-6)</b></p><p>　　這就是泊松方程的有限差分形式，以下估計該方程的精度：</p><p>　　由泰勒公式，易知有以下結(jié)果：</p><p><b>　　(2-7)</b>&l

45、t;/p><p><b>　　若考慮離散的點：</b></p><p><b>　　(2-8)</b></p><p><b>　　上述六式相加</b></p><p><b>　　(2-9)</b></p><p><b>

46、　　代入</b></p><p><b>　　得：</b></p><p><b>　　(2-10)</b></p><p>　　(2-1)式的精度為三階，四階及更高階項被略去。</p><p><b>　　若滿足</b></p><p>&

47、lt;b>　　有</b></p><p><b>　　(2-11)</b></p><p>　　此即拉普拉斯方程的有限差分形式。</p><p>　　這里，我們通過有限差分的方法，把偏微分方程在三階精度下簡化為形式易于計算的代數(shù)方程，從而使之易于在計算機(jī)上實現(xiàn)。</p><p>　　注：有時我們需要解二

48、維靜電場，則方程退化為：</p><p><b>　　(2-12)</b></p><p><b>　　2.2 DBM概述</b></p><p>　　2.1.1 先導(dǎo)發(fā)展規(guī)則</p><p>　　為簡潔起見，以二維平板放電為例，介紹DBM的主要算法和具體步驟[24]。圖2-1中給出了運行若干步的模

49、擬放電圖樣，其中黑點表示放電通道中的點，黑點以實線相連形成通道。</p><p>　　圖2-1 平板放電過程的DBM模擬方法示意圖</p><p>　　模擬區(qū)域內(nèi)的電場在準(zhǔn)靜態(tài)近似下，電位滿足Poisson方程。將上下邊界和放電通道都作為第一類邊界條件，即其中格點的電位是固定的。于是無需考慮空間電荷，因為其效應(yīng)已包含在邊界條件中。圖2-1中每一格點電位可由離散的Laplace方程求取，即

50、：</p><p><b>　　(2-13)</b></p><p>　　每一步的放電過程是鄰近已生成的放電通道的未放電格點(i′,j′)，如果與放電通道中某點(i,j)之間的平均場強超過放電臨界場強Ec，即：</p><p><b>　　(2-14)</b></p><p>　　該點就是可能放電點

51、，(2)式中di′j′,ij為兩點之間的距離。可能放電點在圖2-1中用白圈表示。一個可能放電點成為新的通道點的概率p與局部的場強有關(guān)，其表達(dá)式如(3)式所示：</p><p><b>　　(2-15)</b></p><p>　　其中(i,j)和(i′,j′)分別為一個已放電點和一個可能放電點，η稱為發(fā)展概率指數(shù)，放電通道由(i,j)向(i′,j′)發(fā)展的概率與兩點之

52、間平均場強的η次方成正比。(3)式右端分母是對所有(i,j)和(i′,j′)的組合求和，需要說明的是，每一個通道點和每一個與之鄰近的可能發(fā)展點都構(gòu)成一個這樣的組合。按計算出的概率隨機(jī)抽取下一步的放電發(fā)展方向，如果新的放電通道是由(i,j)向(i′,j′)發(fā)展，假設(shè)沿著放電通道電場強度為Ech，則(i′,j′)處更新后的電位為：</p><p><b>　　(2-16)</b></p&g

53、t;<p>　　然后開始下一步的模擬，如此循環(huán)，直到通道聯(lián)通上下邊界，表示此時間隙擊穿。</p><p>　　2.1.2 分形先導(dǎo)發(fā)展模型關(guān)鍵參數(shù)的選取</p><p>　　由(2-15)式可以看出，η是影響上式的關(guān)鍵參數(shù)[]。通過計算不同取值下模擬出的雷電先導(dǎo)通道的分形維數(shù)，對比實際觀測到的雷電先導(dǎo)通道的分形維數(shù)，可以確定η的合理取值。按照(2-15)式計算出的發(fā)展概率隨機(jī)

54、抽取下一步的發(fā)展方向，如此循環(huán)，直到間隙擊穿。</p><p>　　分形維數(shù)反映復(fù)雜圖形占有空間的有效性，是分形圖形有序性與無序性，確定性與隨機(jī)性的量度。由式(2-15)式可見，分形參數(shù)η其值大小影響了分形先導(dǎo)通道周圍各處發(fā)展概率，決定了雷電先導(dǎo)發(fā)展的隨機(jī)性及確定性。選取合適的分形參數(shù)η值對整個雷電路徑模擬的精確性有非常重要的影響。因此，本文對各種分形參數(shù)η值下雷電放電路徑進(jìn)行模擬，以期得出分形參數(shù)η值對雷電模擬

55、路徑分形維數(shù)的影響，確定分形先導(dǎo)發(fā)展模型合適的分形參數(shù)η值。</p><p>　　基于盒維法的雷電發(fā)展分形維數(shù)為：</p><p><b>　　(2-17)</b></p><p>　　式中，Db為分形維數(shù)，ε為測量尺度，N(ε)為分形圖形能夠被邊長為ε的正方形完全覆蓋的最小個數(shù)。</p><p>　　基于盒維法，文[2

56、6]對不同分形參數(shù)η值下雷電模擬路徑的分形維數(shù)進(jìn)行了計算，如表2-1所示。并選取了3個不同η值下雷電路徑模擬圖(見圖2-2)。</p><p>　　表2-1 不同η值下的分形維數(shù)Db</p><p>　　當(dāng)η取值較小時，通過式(2-15)可以看出，先導(dǎo)頭部各方向發(fā)展概率趨于相同。因此，導(dǎo)致先導(dǎo)通道分支過多，如圖2-2(a)所示。此時，隨機(jī)性在雷電先導(dǎo)的發(fā)展中占主要作用，相應(yīng)的分形維數(shù)較大

57、。反之當(dāng)η值取值較大時，先導(dǎo)頭部最大電場方向獲得了很大的發(fā)展概率，而其它方向的先導(dǎo)發(fā)展概率較小，先導(dǎo)主要沿最大場強方向發(fā)展，分支發(fā)展的概率也相對較小，如圖2-2(c)所示。此時，確定性在雷電先導(dǎo)的發(fā)展中占主要地位，相應(yīng)的分形維數(shù)較小。對于實際放電分形維數(shù)，Charles I Richman 運用盒維法對雷云放電圖片進(jìn)行分形維數(shù)的測量得出，實際雷云放電路徑的分形維數(shù)在1.05~1.4之間，平均值為1.213。而由表2-1可以看出，當(dāng)η=1

58、時，雷電發(fā)展分形維數(shù)Db=1.20±0.04，與實際放電路徑的分形維數(shù)最接近。</p><p>　　圖2-2 雷電先導(dǎo)模擬圖</p><p>　　但是對于η值得選取，其他文獻(xiàn)還有不同的見解，例如文[23]提出根據(jù)生活經(jīng)驗，η大約取2或3之間時應(yīng)該更符合客觀情況。由于這不是本文研究的重點，所以在下文的編程中均采取文[23]的觀點，取η=3。</p><p>

59、　　第3章 雷電先導(dǎo)放電通道的二維模型</p><p>　　3.1 雷電先導(dǎo)放電的物理過程及參數(shù)</p><p>　　3.1.1雷電放電概述</p><p>　　雷電放電包括雷云對大地、雷云對雷云和雷云內(nèi)部的放電現(xiàn)象[]。大多數(shù)雷電放電是在雷云與雷云之間進(jìn)行的，只有少數(shù)是對地進(jìn)行的。我們主要研究雷云對大地的放電，因為這是造成雷害事故的主要因素。</p>

60、<p>　　雷云帶有大量電荷，由于靜電感應(yīng)作用，在雷云下方的地面或地面上的物體將感應(yīng)聚集與雷云相反的電荷，雷云與大地間就形成了電場。當(dāng)雷云中的電荷逐漸積累，達(dá)到一定的電荷密度，使其表面空間的電場強度足夠大(25~30kV∕cm)時，就發(fā)展局部放電。如果此地最大場強方向主要是對地的，就會發(fā)展對地的放電，形成下行雷。</p><p>　　雷電的極性是按照從雷云流入大地的電荷極性決定的。廣泛實測表明，90

61、%左右是負(fù)極性雷。</p><p>　　正雷出現(xiàn)的機(jī)會較少，故對正雷的研究也較少。正雷的下行先導(dǎo)階段沒有明顯的逐級發(fā)展的特征，正雷通常只有一個分量，十分難得有兩個分量的。正雷一般有很長的波頭(長達(dá)幾百微秒)和很長的波尾(長達(dá)上千微秒)，這樣，它所傳遞的電荷可能比多分量的負(fù)雷還多得多，而其電流陡度則比相應(yīng)的負(fù)雷小得多。</p><p>　　若低面上存在特別高聳的導(dǎo)電性能良好的物體時，也可能首

62、先從該物體頂端出發(fā)發(fā)展向上的先導(dǎo)(這種雷稱為上行雷)。但上行先導(dǎo)到達(dá)雷云時，一般不會發(fā)生主放電過程，這是因為雷云的放電性能比大地差得多，難以在極短時間內(nèi)供應(yīng)為中和先導(dǎo)通道中電荷所需的極大的主放電電流，而只能向雷云深處發(fā)展多分支的云中先導(dǎo)，通過寬廣區(qū)域的電暈流注，從分散的水性質(zhì)點上卸下電荷，匯集起來，以中和上行先導(dǎo)中的部分電荷。這樣的放電過程顯然只能是較緩和的，而不可能具有大沖擊電流的特性，而放電電流一般不足千安，而延續(xù)時間則較長，有的可

63、能長達(dá)10-1s。此外，上行先導(dǎo)從一開始就出現(xiàn)分支的概率較大。</p><p>　　3.1.2 雷電的先導(dǎo)過程及參數(shù)</p><p>　　關(guān)于雷電先導(dǎo)的情況[27]，主要由展開照相得來，它與長間隙擊穿的先導(dǎo)有定性的相似，也是由先導(dǎo)通道、先導(dǎo)頭部和流注區(qū)三部分組成。先導(dǎo)過程延續(xù)約幾毫秒，以逐級發(fā)展的、高電導(dǎo)的、高溫的、具有高電位的先導(dǎo)通道將雷云到大地的氣隙擊穿。</p><

64、;p>　　整塊雷云可以有若干個電荷中心[]，負(fù)電荷中心位于雷云的下部，距地面500~10000m的范圍內(nèi)，直接擊向地面的放電通常從負(fù)電荷中心邊緣開始。</p><p>　　下行負(fù)先導(dǎo)具有明顯的分級發(fā)展的特點[27]。每級長度在10~200m范圍內(nèi)，平均為40m.。相鄰兩級間歇時間為30~90μs，平均為60μs。分級前伸速度約為(1~5)×107m/s，延續(xù)約為1μs。由于存在間歇，其總的平均速度

65、約為(1~8)×105m/s，先導(dǎo)過程歷時約0.005~ 0.01s[]。下行負(fù)先導(dǎo)向地面推進(jìn)時還可能出現(xiàn)一些分支，但通常只有其中的一支到達(dá)地面。</p><p>　　梯級先導(dǎo)形成的原因一般解釋為[29]：在先導(dǎo)通道形成之前有一看不見的引路先導(dǎo)(流注)。先導(dǎo)頭部有很強的電場，引路先導(dǎo)在電場的驅(qū)動下向下發(fā)展，由于空間電荷區(qū)的消耗，又因先導(dǎo)頂端聚集了大量正電荷，使得先導(dǎo)通道頂端處的電場大為減弱甚至出現(xiàn)反向。

66、此時引路先導(dǎo)便停止發(fā)展，而電子會退到先導(dǎo)頭部的正電荷區(qū)，產(chǎn)生強烈的復(fù)合，形成很強的光電離，先導(dǎo)就沿著引路先導(dǎo)形成的通道向前發(fā)展，完成一次梯級過程，其頭部強電場又將恢復(fù)。這個過程重復(fù)下去，就是梯級先導(dǎo)過程經(jīng)典的物理圖象。從另一個角度看，下行雷的放電是一種超長間隙放電，放電過程也由先導(dǎo)-流注系來完成。通常的下行負(fù)雷擊從雷云下部負(fù)電荷中心的某處開始，以梯級先導(dǎo)放電的形式向地面發(fā)展，各梯級的方向有明顯的垂直分量，但一般來說，每個梯級在空間的取向

67、是隨機(jī)的，最終使放電通道成為大家習(xí)見的曲折形狀。各個梯級的方向很可能由放電通道頂端周圍的局部電場分布所確定，所以雷電的第一次閃擊幾乎總是有重重分枝，而這些分枝也一樣按梯級形式發(fā)展，以幾十米的梯級向下延伸。因此，雷電的發(fā)展具有明顯的分形特性。</p><p>　　很難直接測出雷云電荷中心對地的電勢[27]，但可以間接推算出：先導(dǎo)根部(在云中)對地的電勢約為50~100MV；先導(dǎo)頭部對地電勢為20~80MV。<

68、/p><p>　　當(dāng)下行雷電導(dǎo)從雷云向建筑物方向發(fā)展時，從接地的建筑物上可能產(chǎn)生向上的迎面先導(dǎo)。產(chǎn)生初始迎面先導(dǎo)的條件與上述產(chǎn)生上行雷的條件相似。區(qū)別在于：前者不僅考慮雷云電場的作用，還應(yīng)考慮下行先導(dǎo)電荷造成局部電場的作用。</p><p>　　迎面先導(dǎo)在相當(dāng)大的程度上影響著下行先導(dǎo)的發(fā)展路線并決定雷擊點的所在，所以，它在對地雷擊的發(fā)展中具有很重要的意義。</p><p&g

69、t;　　直擊到平地的雷電，實際上不存在從地面向上發(fā)展的迎面先導(dǎo)。實際觀測統(tǒng)計指出，平原地區(qū)高度在200m以上的建筑物上，觀測到相當(dāng)多的上行雷電。</p><p><b>　　3.2 模型的建立</b></p><p>　　3.2.1 理想的模擬場景、步驟和參數(shù)</p><p>　　在文[24]中只考慮比較常見的下行負(fù)極性雷的情形，給出了一個理想

70、的雷電先導(dǎo)模型。雷電通道的形成分為下行先導(dǎo)和上行先導(dǎo)兩部分，因而模擬時在一般氣體放電模擬的基礎(chǔ)上考慮多條上行先導(dǎo)起始、發(fā)展以及上、下行先導(dǎo)頭部相接(末躍)等的實現(xiàn)。下行先導(dǎo)從雷云底部起始，當(dāng)其逼近大地時，上行先導(dǎo)會從地面及輸電線路的某位置起始并向上發(fā)展。模擬中，大地以及地面上構(gòu)筑物(如桿塔、避雷線等)設(shè)為零電位，雷云底部等效成一個平板電極，電位為?U0。另外要在雷云底部隨機(jī)選擇一點，產(chǎn)生一段初始的下行先導(dǎo)通道。在每一時間步長，檢查下行先

71、導(dǎo)頭部與地面之間的平均場強，如果超過擊穿臨界值Eg，則雷電將直擊大地。檢查地面目標(biāo)上可能產(chǎn)生上行先導(dǎo)的點處的場強，如果其中某點與一鄰近點之間的平均場強大于上行先導(dǎo)起始的臨界場強Eu，則此處出現(xiàn)擊穿，上行先導(dǎo)起始。對已經(jīng)起始的每條上行先導(dǎo)，按照類似于下行先導(dǎo)的方式發(fā)展。計算每條上行先導(dǎo)的頭部與下行先導(dǎo)的頭部之間的平均場強，設(shè)其最大值E(k)對應(yīng)于第k條上行先導(dǎo)。如果E(k)不小于Ef(Ef為末躍場強判據(jù))，則發(fā)生末躍。單次模擬結(jié)果顯然具有

72、很大的隨機(jī)性。為使得模擬具有實際價值，須多次模擬得到統(tǒng)計</p><p>　　3.2.2 模型的簡化</p><p>　　雖然現(xiàn)有雷電放電路徑發(fā)展觀察資料表明，雷電先導(dǎo)發(fā)展過程表現(xiàn)出梯級向下發(fā)展的特點。但由于在計算機(jī)先導(dǎo)發(fā)展模擬中，忽略先導(dǎo)發(fā)展梯級間隔對先導(dǎo)路徑的選擇及最后雷擊目標(biāo)物選擇等關(guān)鍵問題并沒有明顯影響。因此，本模型將雷電先導(dǎo)梯級發(fā)展簡化為連續(xù)向下發(fā)展。上文指出，通過對實際雷云向地

73、面放電過程的觀察得出，雷電先導(dǎo)梯級長度為10~ 200m。對于計算機(jī)模擬計算，將計算步長設(shè)置小于1個梯級發(fā)展長度便可滿足模擬先導(dǎo)路徑精度要求，進(jìn)一步減小計算步長可以提高模擬精度，但過小的計算步長會造成較大的計算量。因此,本模型綜合考慮模擬精度及計算機(jī)計算量,將計算步長取為雷電梯級長度的最小值，即為10m。</p><p>　　圖3-1 二維格點電勢圖</p><p>　　本模型是一個從概率

74、角度去描述介質(zhì)電擊穿的模型。如圖3-1所示，每一個點都有可能擊穿離它最近的八個點，每擊穿一個點，擊穿點的電勢就可以認(rèn)為和擊穿源相同。由于雷雨天氣中云層頂部帶正電，底部帶負(fù)電，所以會使得地面成為高電勢，云層成為低電勢，所以在計算中我們不妨假設(shè)云層的電勢為0，地面的電勢為50MV(中等雷電)。天地間簡化為一個平面。閃電的路徑為擊穿的點，由于閃電從云端一直到地面，第一個被擊穿的點在云端，可以認(rèn)為所有擊穿的點的電勢都為0，其概率滿足(2-15)

75、式。但是因為我們處理的都是電勢，所以改寫(2-15)式可以得到：</p><p><b>　　(3-1)</b></p><p>　　從點(i,j)出發(fā)，如果僅僅考慮對(i+1,j)，(i-1,j)，(i,j-1)，(i,j+1)點的擊穿，很明顯，當(dāng)網(wǎng)格均勻時，所有di都相等，(3-1)式就蛻變?yōu)橄率剑?lt;/p><p><b>　　(

76、3-2)</b></p><p>　　再考慮加入(i+1,j-1)，(i-1,j-1)，(i+1,j+1)，(i-,j+1)四個點，考慮到(3-1)式中距離di 的影響，我們對這四個點只要把φ改寫為φ∕就可以直接應(yīng)用式(3-2)求解。</p><p>　　本課題模擬的二維雷電先導(dǎo)從1000米高空開始，分級前伸速度300000m/s。我們假設(shè)第一個擊穿的點在云端的中間，每擊穿一個

77、點，我們把已被擊穿的點計入邊界，這滿足拉普拉斯方程，所以我們再用有限差分法迭代重新計算平面的電勢分布，然后依照新算出來的電勢再配合隨機(jī)函數(shù)得到下一個擊穿的點。每個已被擊穿的點都可能會作為生長點，使最近的八個點之一被擊穿。最終擊穿點到達(dá)邊界則計算結(jié)束。因為計算電勢需要完整的邊界條件，所以我們做一個近似的假定，我們認(rèn)為在整個過程中平面左右兩個邊界上的電勢一直不會改變，即滿足最初的線性變化的分布。</p><p>&l

78、t;b>　　3.3 程序流程</b></p><p>　　根據(jù)以上的理論和假設(shè)可以做出程序流程圖，如圖3-2。</p><p>　　顯然，該算法包括以下3大模塊：</p><p><b>　　模塊1：初始情況</b></p><p><b>　　模塊2：計算電勢</b></p

79、><p>　　模塊3：計算概率，決定下一個擊穿點</p><p>　　圖3-2 程序流程圖</p><p><b>　　3.3.1 模塊1</b></p><p>　　clc %清除命令窗口中顯示的所有輸入和輸出</p><p>　　tic;

80、 %啟動秒表計時器，輸出tic使用后的流</p><p>　　時間，與toc結(jié)合記錄程序運行時間</p><p>　　lx=100;ly=100; %定義網(wǎng)格的最大行、列數(shù) </p><p>　　w=2/(1+sin(pi/(lx-1))); %加速收斂

81、因子</p><p>　　v1=zeros(ly,lx); %v1初始化，使其成為全零矩陣</p><p>　　v=v1; %將v1賦給v</p><p>　　r=3; %發(fā)展概率指數(shù)</p><p>

82、;　　c=3e5; %分級前伸速度</p><p>　　T=0; %時間矩陣</p><p>　　tt=0; %返回雷電先導(dǎo)歷時</p><p><b>　　temp=0;</b><

83、/p><p>　　e1=zeros(1,97); %求取一條線路上相鄰兩點之間的場強，放在e1中</p><p>　　%-----------------------------------------使電勢從上到下線性遞增----------------------------------------</p><p>　　for i

84、=1:ly</p><p>　　for j=1:lx</p><p>　　v1(i,j)=(ly-i)*50/(ly-1); </p><p><b>　　end </b></p><p><b>　　end </b></p><p>　　o=sqrt(2);<

85、/p><p>　　v(1,:)=1;v(:,1)=1;v(ly,:)=1;v(:,lx)=1; %邊界</p><p>　　k=0; %迭代次數(shù)</p><p>　　i=99; %起始點的縱坐標(biāo)</p><p>　　

86、j=fix(lx/2); %起始點的橫坐標(biāo)</p><p>　　第3行中天地間開辟了一個100×100的空間，由于上節(jié)中指出模擬的雷電從1000米高空開始，時間步長設(shè)置為10米，所以每一電位坐標(biāo)代表10米。</p><p><b>　　3.3.2 模塊2</b></p><p>　　

87、命xi為第i點電勢，Laplace方程的離散形式可以寫成：</p><p><b>　　(3-3)</b></p><p>　　則對全部n個點，可以列n個方程：</p><p><b>　　(3-4)</b></p><p>　　方程組中，第i個方程中xj前的系數(shù)aj如此給定：當(dāng)j=i時，aj=-4

88、；當(dāng)j為i鄰近的四個點之一時，aj=1；其他情況，aj=0。即有n個方程，每個方程中的n個元素中，除了這個點和與它最鄰近的4個點外，都為0元；把給定的m個邊界條件(包括擊穿點)的電勢值代入方程組，剩余的n-m獨立方程，恰能給出剩余的n-m個未知網(wǎng)格的電勢；這樣就把求各點電勢問題轉(zhuǎn)化為求解線性方程組的問題了。</p><p>　　由于我們用到的差分公式本身就是近似公式，求出該方程的精確解沒有實際意義，我們只需要一個

89、符合精度要求的近似解即可，我們所以我們應(yīng)該考慮一個效率更高的，更容易實現(xiàn)的算法，那么可以考慮使用迭代法[23]。</p><p><b>　　迭代程序如下：</b></p><p>　　while(maxt>1e-5) %精度要求，達(dá)到精度跳出循環(huán)</p><p><b>　　k=k+1; </b></p

90、><p><b>　　maxt=0; </b></p><p>　　for m=2:ly-1 </p><p>　　for n=2:lx-1; </p><p>　　if v(m,n)==0 </p><p>　　v2(m,n)=(v1(m,n+1)+v1(m+1,n)+v1(m-1,n)+v1(m,

91、n-1))/4;</p><p>　　t=abs(v2(m,n)-v1(m,n)); </p><p>　　if(t>maxt) maxt=t; </p><p><b>　　end </b></p><p><b>　　end</b></p><p><b

92、>　　end</b></p><p><b>　　end</b></p><p><b>　　v1=v2;</b></p><p><b>　　end</b></p><p>　　為了加快迭代速度，可以提前使用新值進(jìn)行迭代，因為在上述程序中，我們的掃描賦值方式

93、是一行一行掃描，在對v2(i，j)賦值時，它周圍四個點其中有2個已經(jīng)被掃描過了，即已經(jīng)獲得了新的數(shù)值，這個數(shù)值應(yīng)該更優(yōu)，所以我們盡量使用新算出的數(shù)值進(jìn)行迭代，因此，只要把上述代碼中：</p><p>　　v2(i,j)=(v1(i,j+1)+v1(i+1,j)+v1(i-1,j)+v1(i,j-1))/4</p><p><b>　　改為：</b></p>

94、<p>　　v2(i,j)=(v1(i,j+1)+v1(i+1,j)+v2(i-1,j)+v2(i,j-1))/4</p><p>　　就可以將迭代次數(shù)減少。</p><p>　　為了能使迭代次數(shù)進(jìn)一步減少，考慮讓每次迭代獲得更多的增量，可以采用逐次超松弛法[]：</p><p><b>　　那么將：</b></p>

95、<p>　　v2(i,j)=(v1(i,j+1)+v1(i+1,j)+v2(i-1,j)+v2(i,j-1))/4</p><p>　　變形為：v2(i,j)=v1(i,j)+(v1(i,j+1)+v1(i+1,j)+v2(i-1,j)+v2(i,j-1)-4*v1(i,j))/4</p><p>　　再考慮：v2(i,j)=v1(i,j)+w*(v1(i,j+1)+v1(i+

96、1,j)+v2(i-1,j)+v2(i,j-1)-4*v1(i,j))/4</p><p>　　其中w稱為“加速收斂因子”，是一個供選擇的參數(shù)，其值在1~2之間，該方法的快慢與w有著明顯的關(guān)系。實踐表明，如果w選得好，可以較快地加快迭代速度。經(jīng)驗表明正方形場域由正方形網(wǎng)格劃分，每邊節(jié)點數(shù)為(p+1)時，最佳收斂因子w為</p><p><b>　　(3-5)</b>&

97、lt;/p><p>　　程序最終的電位迭代程序為：</p><p>　　while(maxt>1e-5) </p><p><b>　　k=k+1; </b></p><p><b>　　maxt=0; </b></p><p>　　for m=2:ly-1 </p&

98、gt;<p>　　for n=2:lx-1; </p><p>　　if v(m,n)==0 </p><p>　　v2(m,n)=v2(m,n)+w*(v1(m,n+1)+v1(m+1,n)+v2(m-1,n)+v2(m,n-1)-4*v2(m,n))/4;</p><p>　　t=abs(v2(m,n)-v1(m,n)); </p>

99、<p>　　if(t>maxt) </p><p><b>　　maxt=t; </b></p><p><b>　　end </b></p><p><b>　　end</b></p><p><b>　　end</b></p

100、><p><b>　　end</b></p><p><b>　　v1=v2;</b></p><p><b>　　end</b></p><p><b>　　3.3.3 模塊3</b></p><p>　　結(jié)合式(3-2)，設(shè)變量X

101、表示可能擊穿的n個點，則這n個點的概率可以寫成：</p><p><b>　　(3-6)</b></p><p>　　我們以列表的形式表示出：</p><p>　　表3-1 可能擊穿點的概率分布表</p><p>　　從中顯而易見， 作為概率函數(shù)的pi滿足：</p><p><b>　　

102、(3-7)</b></p><p>　　于是，我們命F(x)作為概率分布函數(shù)有：</p><p><b>　　(3-8)</b></p><p>　　由rand(1,1)生成一個隨機(jī)數(shù)，這個數(shù)必然在[0,1]區(qū)間中隨機(jī)選取，而且這個數(shù)必然對應(yīng)于唯一一個可能擊穿的點，由此即可判定x所在的區(qū)間，進(jìn)而判定pi以及所對應(yīng)的ai，這樣，我們就

103、確定了下一個擊穿點，命這個擊穿點Φ=0。</p><p><b>　　相關(guān)程序如下：</b></p><p><b>　　p=0;</b></p><p>　　for i=2:ly-1 </p><p>　　for j=2:lx-1 </p><p>　　if (v(i,j)

104、==2) </p><p>　　if ( v(i-1,j)==2) </p><p><b>　　p1=0; </b></p><p>　　else p1=v1(i-1,j)^r; </p><p><b>　　end; </b></p><p>　　if (v(i+1,j)

105、==2) </p><p><b>　　p2=0; </b></p><p>　　else p2=v1(i+1,j)^r;</p><p><b>　　end; </b></p><p>　　if (v(i,j-1)==2) </p><p><b>　　p3=0;

106、 </b></p><p><b>　　else</b></p><p>　　p3=v1(i,j-1)^r;</p><p><b>　　end; </b></p><p>　　if (v(i,j+1)==2) </p><p><b>　　p4=0;

107、</b></p><p><b>　　else</b></p><p>　　p4=v1(i,j+1)^r;</p><p><b>　　end; </b></p><p>　　if ( v(i-1,j-1)==2) </p><p><b>　　p5=0

108、; </b></p><p><b>　　else</b></p><p>　　p5=(v1(i-1,j-1)/o)^r; </p><p><b>　　end; </b></p><p>　　if ( v(i-1,j+1)==2) p6=0; </p><p>

109、　　else p6=(v1(i-1,j+1)/o)^r; </p><p>　　end; if ( v(i+1,j-1)==2) </p><p><b>　　p7=0; </b></p><p>　　else p7=(v1(i+1,j-1)/o)^r; </p><p><b>　　end; </b&g

110、t;</p><p>　　if ( v(i+1,j+1)==2) p8=0; </p><p>　　else p8=(v1(i+1,j+1)/o)^r; </p><p>　　end; p=p1+p2+p3+p4+p5+p6+p7+p8+p; </p><p><b>　　end</b></p><p

111、><b>　　end;</b></p><p><b>　　end; </b></p><p>　　q=rand(1,1); %取0~1之間的隨機(jī)數(shù)</p><p><b>　　s=0; </b></p><p>　　%---------

112、-------------------------------計算每一個可能擊穿點的概率-----------------------------------</p><p>　　for i=2:ly-1 </p><p>　　for j=2:lx-1 </p><p>　　if (v(i,j)==2) </p><p>　　if ( v(i

113、-1,j)==2) p1=0; </p><p>　　else p1=(v1(i-1,j)^r/p); end; </p><p>　　if (v(i+1,j)==2) p2=0;</p><p>　　else p2=(v1(i+1,j)^r/p);end; </p><p>　　if (v(i,j-1)==2) p3=0;</p>

114、;<p>　　else p3=(v1(i,j-1)^r/p);end; </p><p>　　if (v(i,j+1)==2) p4=0; </p><p>　　else p4=(v1(i,j+1)^r/p);end; </p><p>　　if ( v(i-1,j-1)==2) p5=0; </p><p>　　else p5

115、=((v1(i-1,j-1)/o)^r)/p; end; </p><p>　　if ( v(i-1,j+1)==2) p6=0; </p><p>　　else p6=((v1(i-1,j+1)/o)^r)/p; end; </p><p>　　if ( v(i+1,j-1)==2) p7=0; </p><p>　　else p7=((v

116、1(i+1,j-1)/o)^r)/p; end; </p><p>　　if ( v(i+1,j+1)==2) p8=0; </p><p>　　else p8=((v1(i+1,j+1)/o)^r)/p; end; </p><p>　　%--------------------------------------結(jié)合隨機(jī)函數(shù)選擇下一個擊穿點-----------

117、-----------------------</p><p>　　s=s+p1+p2+p3+p4+p5+p6+p7+p8; </p><p>　　if (q>=s-p1-p2-p3-p4-p5-p6-p7-p8) && (q<s-p2-p3-p4-p5-p6-p7-p8) </p><p>　　m=i-1; n=j; x0=j;y0=i