景區(qū)滅火的數(shù)學(xué)模型(共27頁(yè))

1、精選優(yōu)質(zhì)文檔-傾情為你奉上西南交通大學(xué)峨眉校區(qū)2017年全國(guó)大學(xué)生數(shù)學(xué)建模競(jìng)賽第二次預(yù)選賽試題 A題：景區(qū)滅火的數(shù)學(xué)模型姓名劉琦王楊學(xué)號(hào)專業(yè)交運(yùn)城軌電氣工程聯(lián)系電話QQ摘要本文以景區(qū)防火問題為背景，在通過多種方式插值方法修補(bǔ)殘缺等高線的基礎(chǔ)上，結(jié)合等高線高程信息建立景區(qū)三維地形模型并巧妙求出其地表總面積，最后通過圖論模型制定最優(yōu)滅火路線，達(dá)到及時(shí)控制和消滅火情的目的。問題一中，面對(duì)有5處需要修補(bǔ)的殘缺等高線圖，我們先根據(jù)各處殘缺是否垂直于X軸的不同特點(diǎn)，通過人機(jī)交互分式將5處殘缺線分離出來，并對(duì)原圖進(jìn)行去噪處理，進(jìn)而將問題轉(zhuǎn)換為一元插值問題。根據(jù)問題特點(diǎn)，我們選擇了三次樣條插值和三次多項(xiàng)式插值

2、，并對(duì)二者插值結(jié)果進(jìn)行比較，發(fā)現(xiàn)三次多項(xiàng)式插值效果更好，因此我們選擇三次多項(xiàng)式插值結(jié)果作為最終修補(bǔ)結(jié)果。最后我們將分區(qū)域修補(bǔ)好的等高線與去噪后的等高線重疊并保存，便可得到修補(bǔ)完全的等高線圖。問題二中，我們?cè)趩栴}一得到完整等高線圖的基礎(chǔ)上，首先將8種高程信息不同的等高線從完整等高線種分離出來，并通過坐標(biāo)變換，將地圖上已知等高線的各點(diǎn)坐標(biāo)轉(zhuǎn)換為實(shí)際點(diǎn)坐標(biāo)。這里考慮DTM（Digital Terrain Model）數(shù)字地面模型，由于高程值的集合是已知的，每一條等高線對(duì)應(yīng)一個(gè)已知的高程值，可以得出其等高線模型，再使用MATLAB工具箱對(duì)多種不同插值方法得到的三維圖進(jìn)行比較，最后通過立方插值和最近點(diǎn)插

3、值的混合插值，得到了更多地表曲面的點(diǎn)再進(jìn)行擬合處理，進(jìn)而得到了較為理想的景區(qū)三維地形模型。對(duì)于景區(qū)地表面積的求解，我們可以通過對(duì)景區(qū)地表曲面進(jìn)行曲面積分進(jìn)行求解，但這種方式需要景區(qū)地表曲面表達(dá)式，在多項(xiàng)式擬合效果不理想的情況下，我們采用了三角網(wǎng)模型（TIN），以直代取，巧妙合理而又不失準(zhǔn)確性的求出了景區(qū)地表面積。問題三我們將確定最優(yōu)的滅火路線問題轉(zhuǎn)換為求解景區(qū)地表曲面指定兩點(diǎn)之間的最短曲線段問題。我們引入圖論中賦權(quán)圖的概念，將曲面上任意兩點(diǎn)間曲線段的長(zhǎng)度作為權(quán)值，賦給其對(duì)應(yīng)投影所對(duì)應(yīng)的網(wǎng)格的邊，進(jìn)而將求三維立體空間內(nèi)最短曲線段的問題轉(zhuǎn)換為求平面賦權(quán)圖網(wǎng)絡(luò)中指定頂點(diǎn)間具有最小權(quán)的路問題，采用迪克

4、斯特拉（Dijkstra）算法，按距A點(diǎn)從近到遠(yuǎn)的順序，依次求得A點(diǎn)到圖中各頂點(diǎn)的最短路和距離，直至B點(diǎn)。但計(jì)算量是遠(yuǎn)超計(jì)算機(jī)所能承受的，為了減少計(jì)算機(jī)計(jì)算量，我們采用了逐步優(yōu)化的方法，對(duì)512*512個(gè)網(wǎng)格，先以8個(gè)網(wǎng)格為邊長(zhǎng)，64個(gè)小網(wǎng)格為一個(gè)大網(wǎng)格，進(jìn)行求解，進(jìn)而對(duì)每個(gè)大網(wǎng)格的路線進(jìn)行二次優(yōu)化，以減少誤差。首先利用Excel處理數(shù)據(jù)，通過C+編寫程序，最后得到了平面賦權(quán)圖中A點(diǎn)到B點(diǎn)具有最小權(quán)的路，并利用MATLAB將其對(duì)應(yīng)的曲線在三維地表曲面中表示了出來。 關(guān)鍵詞： 等高線 DTM 圖論 Dijkstra算法 三角網(wǎng) C+ MATLAB專心-專注-專業(yè)一、問題提出1. 知識(shí)背景隨著時(shí)代

5、和經(jīng)濟(jì)的發(fā)展，森林公園等自然景區(qū)逐漸成為了人們出游時(shí)一個(gè)不可或缺的選擇，防火等安全問題也顯得日益重要。本題基于景區(qū)防火的問題，要求在補(bǔ)全等高線的基礎(chǔ)上建立其三維地形圖，并且通過構(gòu)建數(shù)學(xué)模型制定最優(yōu)滅火路線，達(dá)到及時(shí)控制和消滅火情的目的。2. 需要解決的問題現(xiàn)旅游管理部門想在景區(qū)發(fā)生火災(zāi)時(shí)能及時(shí)控制和消滅火情，請(qǐng)你利用附件提供的數(shù)據(jù)通過建立數(shù)學(xué)模型解決下面三個(gè)問題：（1） 本題中給出的等高圖存在局部破損，利用數(shù)學(xué)模型修補(bǔ)好該等高圖。（2） 在修補(bǔ)好等高圖的基礎(chǔ)上，通過數(shù)學(xué)模型構(gòu)建出景區(qū)的三維地形圖，同時(shí)計(jì)算其地表面積。（3） 等高圖中B點(diǎn)發(fā)生火災(zāi)，從固定的消防站A出發(fā)去滅火，通過建立模型確定最優(yōu)

6、的滅火路線。附件：景區(qū)地形等高圖說明：1.該圖水平及豎直方向以10m每像素為單位，山高以50m為單位。2.實(shí)際圖形見附件，為512×512像素。 二.基本假設(shè)1. 不考慮曲面坡度等其他因素對(duì)于救火員速度的影響，救火員在景區(qū)地表行走時(shí)速度恒定，即在空間內(nèi)任意兩點(diǎn)內(nèi)救火員速度均保持不變。2.救火員接收到報(bào)警電話后就立即出發(fā)，景區(qū)表面每一處均不會(huì)阻礙救火員行進(jìn)。3.計(jì)算面積時(shí)以50米海拔處為景區(qū)地平線。3. 符號(hào)說明符號(hào)意義單位備注網(wǎng)格的x坐標(biāo)無全局變量網(wǎng)格的y坐標(biāo)無全局變量網(wǎng)格的z坐標(biāo)無全局變量S景區(qū)地表總面積全局變量G賦權(quán)圖無全局變量V頂點(diǎn)集無全局變量E邊的集合無全局變量W鄰接矩陣無全

7、局變量到這條邊是否在最短路徑上無0-1矩陣結(jié)點(diǎn)到對(duì)應(yīng)的空間曲線長(zhǎng)度m全局變量結(jié)點(diǎn)到對(duì)應(yīng)的空間曲線d對(duì)應(yīng)的權(quán)值無全局變量D兩個(gè)頂點(diǎn)之間具有的權(quán)無全局變量4 問題分析4.1問題背景隨著我國(guó)的旅游行業(yè)進(jìn)一步發(fā)展，選擇觀賞自然景區(qū)的人越來越多，在對(duì)一些旅游景區(qū)的消防安全以及滅火救援的要求上也有所提高。制定合理的滅火路線，及時(shí)控制和消滅火情，從而達(dá)到保障游客人身安全和保護(hù)自然景觀的目的，直接影響了我國(guó)社會(huì)的安定和諧，同時(shí)對(duì)旅游經(jīng)濟(jì)的快速發(fā)展也具有重大意義。4.2分析問題在問題一中，將補(bǔ)全等高線轉(zhuǎn)化為一元插值問題，利用插值法進(jìn)行曲線修補(bǔ)，最終得到完整的等高圖。 在問題二中，需要在第一問得到的完整等高線基礎(chǔ)

8、上得到景區(qū)地形三維模型，我們可以采用二維插值的方式進(jìn)行求解，對(duì)于表面積的求解，我們可以采用對(duì)曲面進(jìn)行曲面積分的方法或采用三角網(wǎng)模型，以曲代直，進(jìn)而得到景區(qū)地表表面積。在問題三中，需要我們制定最優(yōu)的救火路線，我們可以轉(zhuǎn)換為求解地表曲面上A點(diǎn)到B點(diǎn)的之間最短曲線段的問題，進(jìn)而可以將曲面進(jìn)行投影，引入賦權(quán)圖概念，將最短曲線段的求解進(jìn)一步轉(zhuǎn)換為投影平面內(nèi)求解對(duì)應(yīng)兩點(diǎn)的最小權(quán)路的問題。一元插值缺失的等高圖二維插值問題一補(bǔ)全的等高圖問題二景區(qū)三維地形圖最短路徑問題問題二景區(qū)地表表面積問題三的最優(yōu)救火路徑4.3數(shù)據(jù)的前期處理 由于5處缺口分別分布在4條等高線上，所以要對(duì)圖像進(jìn)行預(yù)處理，即將五處缺口分離為4張

9、新的圖像。對(duì)于圖像分割，我們可以采取灰度閾值分割法、基于區(qū)域的分割方法以及基于邊緣的分割方法等等，但本題給出的圖像較為簡(jiǎn)單，同時(shí)考慮到效率問題，所以采用人機(jī)交互的方法，利用畫圖板的橡皮擦手動(dòng)分割圖像。分割后可以得到以下4張圖片。圖4.1 圖4.2 圖4.3 圖4.45 模型的建立與求解5.1問題一模型建立與求解5.1.1問題一的分析 由數(shù)據(jù)的前期處理，我們將五處缺口分別分散在四張圖片中，圖中只保留了缺口部分的黑點(diǎn)。然后利用MATLAB中的imread函數(shù)可以讀出其灰度矩陣，則對(duì)等高線的補(bǔ)全就轉(zhuǎn)化為了一元插值問題，利用插值法可以分別補(bǔ)全每個(gè)缺口，最后將得到的四張修補(bǔ)后的圖片與最初的缺失等高圖進(jìn)行

10、圖層重疊，重疊后可得到完整的等高圖。5.1.2問題一模型的建立 由問題知我們需要讀出每張圖片的灰度矩陣，然后選擇當(dāng)中部分點(diǎn)進(jìn)行插值法，首先采用人為選點(diǎn)的方法，但在同一圖片選取不同組的點(diǎn)進(jìn)行插值，對(duì)結(jié)果進(jìn)行比較后發(fā)現(xiàn)，人為選點(diǎn)的方法具有很大的主觀性，會(huì)產(chǎn)生較大的誤差，所以該方法不可取。于是利用圖片中保留的所有點(diǎn)，即插值的對(duì)象為整個(gè)矩陣。我們分別采用對(duì)線性插值、三次多項(xiàng)式插值、三次樣條插值等多種方法進(jìn)行插值運(yùn)算，對(duì)曲線進(jìn)行修補(bǔ)，比較其插值結(jié)果，發(fā)現(xiàn)三次多項(xiàng)式插值的結(jié)果更加合理，所以本問題全部采用三次多項(xiàng)式插值的方法。 5.13問題一的模型求解 對(duì)等高線的補(bǔ)全可以轉(zhuǎn)化為一元插值問題，利用MATLAB

11、中imread函數(shù)讀出每張圖片的灰度矩陣，選取合適的插值方法即可解決。1. 計(jì)算思路：經(jīng)過插值結(jié)果的比較，發(fā)現(xiàn)三次多項(xiàng)式插值更加適合本題，用MATLAB中的imread函數(shù)讀出每張圖片的灰度矩陣，然后運(yùn)用三次多項(xiàng)式插值，可以補(bǔ)全圖片中的缺失部分，其余剩下的圖片采取同樣的方法可以補(bǔ)全缺失部分，最后將處理后的圖片與原等高圖在同一窗口中顯示，即得補(bǔ)全的等高圖。2. 計(jì)算步驟：Step1：利用imread函數(shù)讀出圖1的灰度矩陣；Step2：圖4.1缺失部分可以由灰度矩陣得到的坐標(biāo)點(diǎn)數(shù)據(jù)進(jìn)行三次多項(xiàng)式插值得到，圖一的補(bǔ)全完成；Step3：對(duì)圖4.2的操作重復(fù)前面兩步即可完成；Step4：圖4.3、圖4.

12、4首先重復(fù)第一步，但由于圖中缺失部分附近的等高線彎曲程度太大，導(dǎo)致灰度矩陣在讀取時(shí)無法一一對(duì)應(yīng)，所以在讀出其灰度矩陣后要先對(duì)矩陣進(jìn)行轉(zhuǎn)置，然后重復(fù)第二步，圖三圖四的補(bǔ)全完成；Step5：圖4.1-4.4的補(bǔ)全均完成，將其與原等高圖在同一窗口顯示，得到最終補(bǔ)全的等高圖。3. 計(jì)算結(jié)果：經(jīng)過上述的處理后，將最終補(bǔ)全得到的等高圖保存為bmp，單色格式，圖片如下（程序見附錄）： 圖5.1.1 三次樣條插值結(jié)果 圖5.1.2三次多項(xiàng)式插值結(jié)果（效果更好） 由多種方式插值得到的圖像可知，相比于三次樣條插值，三次多項(xiàng)式插值得到的完整等高線更為光滑，效果也更好。因此最終我們選擇三次多項(xiàng)式插值的結(jié)果作為等高線的

13、補(bǔ)全圖。5.1.3問題一結(jié)果的分析與驗(yàn)證對(duì)于補(bǔ)全后的等高圖，我們可以采用分析點(diǎn)的差異度來判斷曲線的補(bǔ)全是否合理。對(duì)于點(diǎn)的分析包括該點(diǎn)的位置和其對(duì)應(yīng)的切向量，然后對(duì)兩個(gè)點(diǎn)進(jìn)行切向量和差異函數(shù)的計(jì)算，差異函數(shù)值越低則表示圖像中的兩個(gè)點(diǎn)越接近事實(shí)。5.2.1問題二模型的分析問題要求我們?cè)趩栴}一的基礎(chǔ)上，根據(jù)補(bǔ)全的等高線圖建立景區(qū)地表的三維地形圖，并求出景區(qū)地表面積。通過第一問對(duì)于等高線的修補(bǔ)，結(jié)合等高線的高程信息，我們可以得到景區(qū)地表的離散坐標(biāo)，即我們可以用一組有序數(shù)值陣列形式表示地面高程信息。但我們所得到的點(diǎn)比較稀疏，數(shù)量還遠(yuǎn)遠(yuǎn)不能擬合出地表曲面。因此，我們需要進(jìn)行二維插值，進(jìn)行二維插值模型。根據(jù)

14、已有點(diǎn)信息插值得到更多地表曲面上的點(diǎn)，然后MATLAB進(jìn)行擬合曲面。得到擬合后的地表曲面后，意味著我們得到了地表曲面上所有點(diǎn)的坐標(biāo)信息以及曲面的表達(dá)式，緊接著我們便建立三角形網(wǎng)模型，采用以直代取，將曲面面積用足夠多的小三角形平面面積去表示，最后得到的三角形面積總和即為景區(qū)地表曲面面積。而我們欲得到景區(qū)地表三維模型，首先應(yīng)根據(jù)像素與比例尺求得地圖中各等高線的實(shí)際高度，我們首先將修補(bǔ)好的等高線圖根據(jù)各等高線高程不同進(jìn)行分割，進(jìn)而采用MATLAB循環(huán)賦值的方法進(jìn)行數(shù)據(jù)的前期處理（程序見附錄）。等高線1-8如下圖所示： 5.2.2問題二模型建立首先我們先采用MATLAB內(nèi)部二維插值工具Curve Fi

15、tting Tool中的多種插值方式Cubic、v4，以及可以求出曲面表達(dá)式的polynomial等插值方法進(jìn)行求解，求得的三維地形效果差不多，但是根據(jù)可以求出曲面表達(dá)式的polynoimal插值方法得到的曲面效果卻不理想，因?yàn)閜olynoimal插值法x,y最高次為二元五次表示式，故我們用Custom Equation自擬表達(dá)式來提高x,y次數(shù)，發(fā)現(xiàn)隨著次數(shù)的增高，效果并沒有變得更理想。 圖5.2.1 Cubic插值 圖5.2.2 V4插值圖5.2.3 二元五次表達(dá)式擬合曲面 圖5.2.4 二元六次表達(dá)式擬合曲面 因此我們需要重新選取方法，考慮到求解表面積需要插值點(diǎn)坐標(biāo)以及本問題特點(diǎn)，我們根

16、據(jù)部分等高線高程進(jìn)行二維插值，而二維插值分為結(jié)點(diǎn)插值與散亂結(jié)點(diǎn)插值，因本問題中曲面上的點(diǎn)并不是全部落在網(wǎng)絡(luò)節(jié)點(diǎn)上，因此我們采用散亂結(jié)點(diǎn)插值，即本問題可簡(jiǎn)化為： 已知個(gè)結(jié)點(diǎn)，求點(diǎn)處的插值。 對(duì)于上述問題。MATLAB中提供了插值函數(shù)gridata，其格式為：其中，x,y,z均為n維向量，指明所給數(shù)據(jù)點(diǎn)的橫坐標(biāo)，縱坐標(biāo)和豎坐標(biāo)；向量是給定的網(wǎng)格點(diǎn)的橫坐標(biāo)和縱坐標(biāo)；返回值為網(wǎng)格處的函數(shù)值。為方向不同的向量，即一個(gè)是行向量，另一個(gè)是列向量。為保證準(zhǔn)確性以及減小誤差，我們采用混合插值的方式，進(jìn)行立方差值和最近點(diǎn)插值的混合插值，通過編寫程序可得到三維地形圖（程序見附錄）：我們首先得到了以50米為“地平線”

17、的景區(qū)地形三維圖，即得到了海拔高度50米以上的景區(qū)地形圖。圖5.2.5 50米為地平線的景區(qū)三維圖但題目要求我們建立整個(gè)景區(qū)的地形三維模型。對(duì)于海拔小于50米的部分，我們首先想到的是將其他部分都假設(shè)海拔高度為0，但這種想法會(huì)產(chǎn)生大量斷崖，顯然不符合實(shí)際。而實(shí)際情況應(yīng)該為其余部分海拔高度逐漸上升，進(jìn)而與50米處平緩連接。因此，我們考慮了其余區(qū)域，將其也列入插值區(qū)域，重新進(jìn)行插值，并得到了重新插值后的三維地形圖及其俯視圖：圖5.2.6 景區(qū)地形三維圖 圖5.2.7 景區(qū)地形三維圖（俯視角） 由最終插值所得到的景區(qū)地形三維地形圖及其俯視圖可知，最終插值將海拔50米以下的區(qū)域考慮進(jìn)去，彌補(bǔ)了第一次插值

18、的缺陷，建立了（）的整個(gè)景區(qū)的地形三維圖。成功建立了景區(qū)的地形三維圖后，我們利用已建立的模型以及數(shù)據(jù)進(jìn)行求解景區(qū)地表表面積，我們首先想到的方法就是通過多項(xiàng)式擬合的方式求出曲面的表達(dá)式，進(jìn)而利用曲面積分便可求得地表曲面的面積，即：便可求得地表曲面的面積。但是從表達(dá)式進(jìn)行擬合曲面效果圖來看，可以發(fā)現(xiàn)擬合效果并不理想，誤差太大，因此我們用表達(dá)式進(jìn)行曲面積分而求出景區(qū)表面積的方法是不可取的，所以我們建立三角形網(wǎng)模型來近似求得景區(qū)地表表面積。在本模型中我們先將景區(qū)地表劃分為n個(gè)小三角形，以直代曲，求景區(qū)面積就近似等于這n個(gè)小三角形的面積之和。因?yàn)槿c(diǎn)可以確定一個(gè)平面，而且等高圖為512*512像素的圖像

19、，容易想到將整個(gè)景區(qū)表面積劃分為512*512個(gè)網(wǎng)格，選取其中一個(gè)網(wǎng)格示意圖如下： 其中A,B,C,D分別代表上述網(wǎng)格的四個(gè)頂點(diǎn)。Matlab中的cross函數(shù)，用于返回兩個(gè)向量的數(shù)量積，為了方便編程，所以采用向量的數(shù)量積來對(duì)三角形的面積進(jìn)行求解，則面積為同理，可以求得面積為：那么一個(gè)網(wǎng)格的面積為最終得到景區(qū)地表面積的為：由此，對(duì)所有的小三角形進(jìn)行遍歷求和，便可近到整個(gè)景區(qū)地表面積S 5.2.3問題二模型的求解1. 計(jì)算思路：需要在第一問得到的完整等高線基礎(chǔ)上得到景區(qū)地形三維模型，我們可以采用二維插值的方式進(jìn)行求解，對(duì)于表面積的求解，我們可以采用對(duì)曲面進(jìn)行曲面積分的方法或采用三角網(wǎng)模型，以曲代

20、直，進(jìn)而得到景區(qū)地表表面積。2. 計(jì)算步驟：Step1：將第一問得到的完整等高線按照高程分離，得到8張分離后的等高線。Step2：利用混合插值的方法還原其三維地形圖。Step3：將景區(qū)表面積劃分三角網(wǎng)，以直代曲，通過matlab編程對(duì)所有三角形面積求和最終得到景區(qū)地表表面積。3. 計(jì)算結(jié)果： 由我們?cè)谥暗牡降亩S插值結(jié)果，編寫遍歷循環(huán)程序（程序見附錄）對(duì)所有小三角形面積進(jìn)行求和。最后得到整個(gè)景區(qū)地表表面積近似為： 5.2.4問題二結(jié)果的分析與驗(yàn)證1.三維地形圖的結(jié)果分析：在構(gòu)建景區(qū)的三維地形圖時(shí)，采用了多種插值方法對(duì)圖形進(jìn)行繪制，采用MATLAB內(nèi)部二維插值工具Curve Fitting T

21、ool中的多種插值方式Cubic、v4，以及可以求出曲面表達(dá)式的polynomial等插值方法進(jìn)行求解，求得的三維地形大致符合要求，但是在光滑度等方面較差。通過比較我們采用混合插值，將多種方法結(jié)合，最終得到滿意的三維地形圖。2.景區(qū)地表表面積的結(jié)果分析：在問題二中計(jì)算面積時(shí)，我們并沒有對(duì)等高線區(qū)域以外的部分進(jìn)行條件限制，若假設(shè)等高線區(qū)域以外的點(diǎn)高度都是50米，利用相同的算法，計(jì)算可得其面積為將兩種處理方式進(jìn)行對(duì)比，計(jì)算其差異度由此可見，兩種處理計(jì)算計(jì)算出的表面積的差異度很小。假設(shè)該景區(qū)全為平地，那么其表面積，因?yàn)闃?gòu)建出的景區(qū)地形三維圖存在坡面，所以其地表表面積一定大于全為平地情況下的表面積，通

22、過比較：。結(jié)合上述兩種方法的分析驗(yàn)證，我們可以確定計(jì)算出的地表表面積是正確的。5.3.1問題三問題分析 問題三在問題二的基礎(chǔ)上，要求建立一條救火的最佳路徑。我們假設(shè)救火員在空間任意兩點(diǎn)之間的速度是相同的，因此，我們可以將問題轉(zhuǎn)換為最短路徑的問題。因?yàn)榇司然鹇窂綖榭臻g中位于地表曲面上的一條曲線，我們可先將位于固定曲面上兩點(diǎn)的最短曲線轉(zhuǎn)換為二維頂點(diǎn)之間最短路問題的數(shù)學(xué)規(guī)劃模型。我們引入賦權(quán)網(wǎng)絡(luò)圖的概念，通過景區(qū)地表曲面在平面內(nèi)的投影區(qū)域構(gòu)造賦權(quán)有向圖。5.3.2問題三模型建立構(gòu)造賦權(quán)圖:（1）以景區(qū)地表曲面在平面內(nèi)的投影區(qū)域上劃分網(wǎng)格中的結(jié)點(diǎn)為網(wǎng)絡(luò)圖的結(jié)點(diǎn)，記做。頂點(diǎn)集，這里表示各個(gè)結(jié)點(diǎn)。（2）以

23、景區(qū)地表曲面在平面內(nèi)的投影區(qū)域上劃分網(wǎng)格中的結(jié)點(diǎn)為網(wǎng)絡(luò)圖的邊。記做。為邊的集合，（3）以各條邊在景區(qū)地表曲面上對(duì)應(yīng)的曲線段的長(zhǎng)度作為其相應(yīng)邊的權(quán)值，即為網(wǎng)絡(luò)圖中結(jié)點(diǎn)到的權(quán)值。我們?cè)O(shè)投影到頂點(diǎn)的空間坐標(biāo)為，投影到頂點(diǎn)的空間坐標(biāo)為，根據(jù)假設(shè)，我們近似以曲面上兩點(diǎn)之間的線段長(zhǎng)度來代替兩點(diǎn)之間的曲線段長(zhǎng)度，即：（4） 鄰接矩陣是表示頂點(diǎn)之間相鄰關(guān)系的矩陣，以鄰邊矩陣表示網(wǎng)絡(luò)圖，記做，因本模型建立的為賦權(quán)圖，有： 我們?cè)O(shè)救火員起點(diǎn)A為,終點(diǎn)B為,我們引入0-1變量：我們的目標(biāo)函數(shù)是求賦D權(quán)圖上指定的兩個(gè)頂點(diǎn)之間具有最小權(quán)的路。即最短路徑上各條邊的長(zhǎng)度權(quán)值最小，即： 故最終我們建立的數(shù)學(xué)規(guī)劃表達(dá)式為： 5

24、.3.3問題三模型的求解1.計(jì)算思路：?jiǎn)栴}三在問題二的基礎(chǔ)上，要求建立一條救火的最佳路徑。假設(shè)在空間任意兩點(diǎn)之間的速度是相同的，因此將問題轉(zhuǎn)換為最短路徑的問題。因?yàn)榇司然鹇窂綖榭臻g中位于地表曲面上的一條曲線，我們可先將位于固定曲面上兩點(diǎn)的最短曲線轉(zhuǎn)換為二維頂點(diǎn)之間最短路問題的數(shù)學(xué)規(guī)劃模型。我們引入賦權(quán)網(wǎng)絡(luò)圖的概念，通過景區(qū)地表曲面在平面內(nèi)的投影區(qū)域構(gòu)造賦權(quán)有向圖。確定起點(diǎn)A和終點(diǎn)B利用Dijkstra算法確定行走路徑。迪克斯拉特（Dijkstra）算法，作為求解最短路徑一種成熟的算法，其基本思想是按距離從近到遠(yuǎn)為順序，依次求得到G的各項(xiàng)點(diǎn)的最短路和距離，直至，算法結(jié)束。為避免重復(fù)并保留每一步的

25、計(jì)算信息，采用了標(biāo)號(hào)算法。下面為該算法：Setp1:令=0,對(duì)，令。Setp2:對(duì)每個(gè)，用代替，這里用來表示頂點(diǎn)之間邊的權(quán)值。計(jì)算，把達(dá)到這個(gè)最小值的一個(gè)頂點(diǎn)記為，令。Setp3:若，則停止；若，則用代替，轉(zhuǎn)Setp2。算法結(jié)束時(shí)，從到各頂點(diǎn)的距離由的最后一次標(biāo)號(hào)給出。在進(jìn)入之前的標(biāo)號(hào)叫T標(biāo)號(hào)，進(jìn)入之后的標(biāo)號(hào)叫P標(biāo)號(hào)。算法就是不斷修改各項(xiàng)頂點(diǎn)的T標(biāo)號(hào)，直至獲得P標(biāo)號(hào)。2.計(jì)算步驟： 問題三我們采用Dijkstra算法，但由于數(shù)據(jù)量過大，如果對(duì)于全局采用Dijkstra算法，運(yùn)算量會(huì)很大，導(dǎo)致算法效率不高。因此，我們采用逐步優(yōu)化的搜索方法，從而達(dá)到減少運(yùn)算量的目的。逐步優(yōu)化的步驟如下：Step1

26、：輸入景區(qū)地形的數(shù)值模型，降低網(wǎng)格密度，取80m為單位網(wǎng)格長(zhǎng)度Step2；在step1得出的稀疏網(wǎng)格中，利用Dijkstra算法求出最優(yōu)路徑，將此路徑上的點(diǎn)作為控制點(diǎn)，轉(zhuǎn)Step3。Step3：利用step2中的控制點(diǎn)，將每個(gè)80m*80m的控制節(jié)點(diǎn)分為64個(gè)10m*10m的網(wǎng)格，局部?jī)?yōu)化，即二次優(yōu)化。Step4：將每個(gè)控制節(jié)點(diǎn)與其局部?jī)?yōu)化的節(jié)點(diǎn)結(jié)合，最終得出最優(yōu)路徑。3.計(jì)算結(jié)果：我們首先利用C+軟件編程我們得到了在稀疏網(wǎng)格中節(jié)點(diǎn)的距離矩陣（程序見附錄），再通過MATLAB軟件編程我們得到了最佳救火路線在地形三維圖上的表示，如下圖所示（程序見附錄）： 5.3.4問題三模型的改進(jìn)在求解問題三時(shí)

27、，我們不考慮曲面坡度等其他因素對(duì)于救生員速度的影響，假定救火員在景區(qū)地表行走時(shí)速度恒定。但在實(shí)際情況中，由于受到坡度、地面平坦度和樹木分布等因素，這種假設(shè)是不合理的。所以問題三轉(zhuǎn)化為一個(gè)多目標(biāo)決策問題，我們根據(jù)層次分析法的原理，可以建立篩選最優(yōu)路線的層次結(jié)構(gòu)模型。模型兼顧距離、坡度、地面平坦度、樹木分布、火勢(shì)蔓延等因素，確定各方案的優(yōu)劣順序，進(jìn)行最優(yōu)路線分析，構(gòu)建AHP判斷矩陣并計(jì)算權(quán)重，最終得到最優(yōu)路線。結(jié)果表明，改進(jìn)后的模型對(duì)于救火路線方案進(jìn)行優(yōu)選更加科學(xué)、方便。6、 模型的評(píng)價(jià)與推廣6.1 模型的評(píng)價(jià)優(yōu)點(diǎn)：1. 問題一、問題二繪制圖形時(shí)采用了多種插值方法進(jìn)行比較，最終得出最理的結(jié)果。2.

28、 問題三的模型簡(jiǎn)單易懂，并且考慮到效率問題，采取了二次優(yōu)化，大大減少了計(jì)算量。缺點(diǎn)：1. 在問題三中沒有考慮坡度等因素，在實(shí)際情況中其實(shí)是不允許的。2. 由于給出的等高圖較為簡(jiǎn)單，所以在處理圖像時(shí)很多都是人工進(jìn)行操作，對(duì)于那些等高線密集的圖像，這種處理方式是很不方便的。6.2 模型的推廣 對(duì)于問題一中使用的三次多項(xiàng)式插值，通過利用一些特征節(jié)點(diǎn)，用多項(xiàng)式擬合的方法來產(chǎn)生平滑的插值曲線，其易操作且計(jì)算量不大，對(duì)于溫度、高程、地下水位高度和污染濃度等問題同樣適用。對(duì)于問題二，在使用插值的基礎(chǔ)上，利用軟件重建其三維地形圖，由平面圖得到其三維圖，在機(jī)械三維、建筑三維、美術(shù)三維等方面都是極其重要的。問題三

29、對(duì)于生活生產(chǎn)中最短路徑的問題中具有參考意義。7、 參考文獻(xiàn)1數(shù)學(xué)建模算法與應(yīng)用 司守奎，孫璽菁編著 國(guó)防工業(yè)出版社 2016.72MATLAB在數(shù)學(xué)建模中的應(yīng)用（第二版）卓金武著 北京航空航天大學(xué)出版社 2014.93MATLAB教程 張志涌，楊祖櫻等編著 北京航空航天大學(xué)出版社 2015.14基于等高線的三角形建模及真實(shí)感地形重建 翁巧琳,姜昱明 西安電子科技大學(xué)計(jì)算機(jī)學(xué)院，陜西 西安 5基于MSDAE的等髙線地圖斷線重連技術(shù)研究汪梁 西安電子科技大學(xué)計(jì)算機(jī)學(xué)院 TP3116DEM數(shù)據(jù)提取關(guān)鍵技術(shù)研究 喬利軍 國(guó)防科技大學(xué)研究生院 2005.117彩色地形圖中等高線斷裂重建算法的研究 李曉清

30、 西安電子科技大學(xué)碩士學(xué)位論文 2014.12八、 附錄附錄清單：附錄1：求解問題一的MATLAB程序附錄2：求解問題二的MATLAB程序附錄3：求解問題三的C+程序附錄1 問題一中補(bǔ)全等高線MATLAB程序8.1.1三次樣條插值程序%3圖4圖 無需轉(zhuǎn)置 a=ones(512,512);imshow(a);hold onfor n=3:4imageName=strcat(num2str(n),'.bmp');A=imread(imageName);k=1;for j=1:512for i=1:512if A(i,j)=0B(1,k)=j;B(2,k)=i;k=k+1;break

31、endendendB;x=B(1,:);y=B(2,:);xi=min(B(1,:):1:max(B(1,:);yi=interp1(x,y,xi,'spline');plot(xi,yi,'k')axis(1,512,1,512)hold onclearendhold offF1=getframe;imwrite(F1.cdata,'test1.bmp')c=imread('test1.bmp')d=im2bw(c,0.8);imshow(d)imwrite(d,'test3.bmp')%1圖2圖 要轉(zhuǎn)置a=on

32、es(512,512);imshow(a);hold onfor n=1:2imageName=strcat(num2str(n),'.bmp');A=imread(imageName);A=A'k=1;for j=1:512for i=1:512if A(i,j)=0B(1,k)=j;B(2,k)=i;k=k+1;breakendendendB;x=B(1,:);y=B(2,:);xi=min(B(1,:):1:max(B(1,:);yi=interp1(x,y,xi,'spline');plot(xi,yi,'k')axis(1,51

33、2,1,512)hold onclearendhold offF1=getframe;imwrite(F1.cdata,'test2.bmp')c=imread('test2.bmp')d=im2bw(c,0.8);imshow(d)imwrite(d,'test4.bmp')%3張圖合并到一張圖里并保存a=imread('quzaoyuantu.bmp')b=imread('test3.bmp')c=imread('test4.bmp')c=c'for i=1:512for j=1:512

34、if a(i,j)&b(i,j)&c(i,j)=1d(i,j)=1;elsed(i,j)=0;endendenddimshow(d)imwrite(d,'zuizhongtu.bmp')8.1.2三次多次式插值程序%3圖4圖 無需轉(zhuǎn)置 a=ones(512,512);imshow(a);hold onfor n=3:4imageName=strcat(num2str(n),'.bmp');A=imread(imageName);k=1;for j=1:512for i=1:512if A(i,j)=0B(1,k)=j;B(2,k)=i;k=k+1

35、;breakendendendB;x=B(1,:);y=B(2,:);xi=min(B(1,:):1:max(B(1,:);yi=interp1(x,y,xi,'cubic');plot(xi,yi,'k')axis(1,512,1,512)hold onclearendhold offF1=getframe;imwrite(F1.cdata,'test1.bmp')c=imread('test1.bmp')d=im2bw(c,0.8);imshow(d)imwrite(d,'test3.bmp')%1圖2圖 要轉(zhuǎn)

36、置a=ones(512,512);imshow(a);hold onfor n=1:2imageName=strcat(num2str(n),'.bmp');A=imread(imageName);A=A'k=1;for j=1:512for i=1:512if A(i,j)=0B(1,k)=j;B(2,k)=i;k=k+1;breakendendendB;x=B(1,:);y=B(2,:);xi=min(B(1,:):1:max(B(1,:);yi=interp1(x,y,xi,'cubic');plot(xi,yi,'k')axis(

37、1,512,1,512)hold onclearendhold offF1=getframe;imwrite(F1.cdata,'test2.bmp')c=imread('test2.bmp')d=im2bw(c,0.8);imshow(d)imwrite(d,'test4.bmp')%3張圖合并到一張圖里并保存a=imread('quzaoyuantu.bmp')b=imread('test3.bmp')c=imread('test4.bmp')c=c'for i=1:512for j=1

38、:512if a(i,j)&b(i,j)&c(i,j)=1d(i,j)=1;elsed(i,j)=0;endendenddimshow(d)imwrite(d,'zuizhongtu.bmp') 附錄2 問題二中建立三維模型及求得地表表面積MATLAB程序clear C=; a=0; z=; for n=1:8imageName=strcat(num2str(n),'.bmp');A=imread(imageName);k=1;for i=1:512for j=1:512if A(i,j)=0B(1,k)=i;B(2,k)=j;k=k+1;end

39、endendB;C=C B;B=;x=C(1,:);y=C(2,:);D=size(C);b=D(:,2)-a;a=D(:,2);E=50*n*ones(1,b);z=z E;endxmm=minmax(x)ymm=minmax(y)xi=xmm(1):xmm(2);yi=ymm(1):ymm(2);zi1=griddata(x,y,z,xi,yi','cubic');zi2=griddata(x,y,z,xi,yi','nearest');zi=zi1;zi(isnan(zi1)=zi2(isnan(zi1)subplot(1,2,1),plo

40、t(x,y,'*')subplot(1,2,2),mesh(xi,yi,zi)s=0;for i=1:403 for j=1:367 a=i*10,j*10,zi(i,j); b=i*10,j*10+10,zi(i,j+1); c=i*10+10,j*10,zi(i+1,j); d=i*10+10,j*10+10,zi(i+1,j+1); s1=0.5*norm(cross(b-a,c-a); s2=0.5*norm(cross(b-d,c-d); s=s+s1+s2; endendsst=5120*5120-3680*4040+s求解結(jié)果：s = 1.5174e+07st =

41、 2.6521e+07附錄3 求解問題三中最短路徑的MATLAB程序8.3.1求解距離矩陣的C+程序：#include<iostream>#include<math.h>#include<iomanip>using namespace std;void main(void)int qd,zd;long double qz,qz1,qz2,qz3,qz4,qz5,qz6,qz7,qz8,qz9,qz10;int i,n,m;cout<<"請(qǐng)輸入長(zhǎng)度為10邊的個(gè)數(shù)n&qu

42、ot;cin>>n;for(i=0;i<n;i+)        cin>>qdi>>zdi>>qz1i>>qz2i;qzi=pow(double(10*10+pow(abs(qz1i-qz2i),2),(double)1/2)+pow(double(10*10+pow(abs(qz2i-qz3i),2),(double)1/2)+pow(double(10*10+pow(abs(qz3i-qz4i),2),(double)1/2)+pow(doub

43、le(10*10+pow(abs(qz4i-qz5i),2),(double)1/2)+pow(double(10*10+pow(abs(qz5i-qz6i),2),(double)1/2)+pow(double(10*10+pow(abs(qz6i-qz7i),2),(double)1/2)+pow(double(10*10+pow(abs(qz7i-qz8i),2),(double)1/2)+pow(double(10*10+pow(abs(qz8i-qz9i),2),(double)1/2)+pow(double(10*10+pow(abs(qz9i-qz10i),2),(double)

44、1/2);cout<<"nnnnnnnnnnnnnnn計(jì)算完成nnnnnnnnnnnnnnn"for(i=0;i<n;i+)cout<<setw(10)<<qdi<<setw(10)<<zdi<<setw(10)<<qzi<<endl; 8.3.2求解救火路線的C+程序：#include <iostream>#include <limits>using namespace std;struct Nod

45、e  /定義表結(jié)點(diǎn)  int adjvex; /該邊所指向的頂點(diǎn)的位置  int weight;/ 邊的權(quán)值  Node *next; /下一條邊的指針;struct HeadNode / 定義頭結(jié)點(diǎn)    int nodeName; / 頂點(diǎn)信息    int inDegree; / 入

46、度    int d; /表示當(dāng)前情況下起始頂點(diǎn)至該頂點(diǎn)的最短路徑,初始化為無窮大    bool isKnown; /表示起始頂點(diǎn)至該頂點(diǎn)的最短路徑是否已知,true表示已知，false表示未知    int parent; /表示最短路徑的上一個(gè)頂點(diǎn)    Node *link; /指向第一條依附該頂點(diǎn)的邊的指針;/G表示指向頭結(jié)點(diǎn)數(shù)組的第一個(gè)結(jié)點(diǎn)的指針

47、/deNum表示結(jié)點(diǎn)個(gè)數(shù)/arcNum表示邊的個(gè)數(shù)void createGraph(HeadNode *G, int nodeNum, int arcNum)   cout << "開始創(chuàng)建圖(" << nodeNum << ", " << arcNum  << ")"

48、 << endl;  /初始化頭結(jié)點(diǎn)  for (int i = 0; i < nodeNum; i+)     Gi.nodeName = i+1; /位置0上面存儲(chǔ)的是結(jié)點(diǎn)v1,依次類推    Gi.inDegree = 0; /入度為0    G

49、i.link = NULL;    for (int j = 0; j < arcNum; j+)     int begin, end, weight;    cin >> begin >> end >> weight; &#

50、160;  / 創(chuàng)建新的結(jié)點(diǎn)插入鏈接表    Node *node = new Node;    node->adjvex = end - 1;    node->weight = weight;    +Gend-1.inDegree; /入度加1   &

51、#160;/插入鏈接表的第一個(gè)位置    node->next = Gbegin-1.link;    Gbegin-1.link = node;  void printGraph(HeadNode *G, int nodeNum)   for (int i = 0; i < nodeNum; i

52、+)     cout << "結(jié)點(diǎn)v" << Gi.nodeName << "的入度為"    cout << Gi.inDegree << ", 以它為起始頂點(diǎn)的邊為: "    Node *node =&

53、#160;Gi.link;    while (node != NULL)       cout << "v" << Gnode->adjvex.nodeName << "(權(quán):" << node->weight << ")"&

54、#160;<< "  "      node = node->next;        cout << endl;  /得到begin->end權(quán)重int getWeight(HeadNode *G, int begin, int end) &

55、#160; Node *node = Gbegin-1.link;  while (node)     if (node->adjvex = end - 1)       return node->weight;        node =

56、0;node->next;  /從start開始，計(jì)算其到每一個(gè)頂點(diǎn)的最短路徑void Dijkstra(HeadNode *G, int nodeNum, int start)   /初始化所有結(jié)點(diǎn)  for (int i = 0; i < nodeNum; i+)     Gi.d = INT_MAX;

57、 /到每一個(gè)頂點(diǎn)的距離初始化為無窮大    Gi.isKnown = false; / 到每一個(gè)頂點(diǎn)的距離為未知數(shù)    Gstart-1.d = 0; /到其本身的距離為0  Gstart-1.parent = -1; /表示該結(jié)點(diǎn)是起始結(jié)點(diǎn)  while(true)    /= 如果所有的結(jié)點(diǎn)的最短距離都已知,&#

58、160;那么就跳出循環(huán)   int k;   bool ok = true; /表示是否全部ok   for (k = 0; k < nodeNum; k+)      /只要有一個(gè)頂點(diǎn)的最短路徑未知,ok就設(shè)置為false     if (!Gk.isKnown

59、)        ok = false;       break;            if (ok) return;   /=   /= 搜索未知結(jié)點(diǎn)中d最小的,將其變?yōu)閗nown   

60、/= 這里其實(shí)可以用最小堆來實(shí)現(xiàn)   int i;   int minIndex = -1;   for (i = 0; i < nodeNum; i+)      if (!Gi.isKnown)        if (m

61、inIndex = -1)         minIndex = i;       else if (GminIndex.d > Gi.d)         minIndex = i;     

62、      /=   cout << "當(dāng)前選中的結(jié)點(diǎn)為: v" << (minIndex+1) << endl;     GminIndex.isKnown = true; /將其加入最短路徑已知的頂點(diǎn)集     / 將以minIndex為起

63、始頂點(diǎn)的所有的d更新     Node *node = GminIndex.link;     while (node != NULL)        int begin = minIndex + 1;       int end

64、 = node->adjvex + 1;       int weight = getWeight(G, begin, end);       if (GminIndex.d + weight < Gend-1.d)