劉吉波，王志紅

(貴州工程應(yīng)用技術(shù)學(xué)院 礦業(yè)工程學(xué)院，貴州 畢節(jié) 551700)

等值線圖是一種非常重要的專題圖型，常用的等值線圖有等高線圖、等深線圖、等溫線圖、等壓線圖等。由離散點(diǎn)生成等值線，可采用Delaunay三角剖分算法根據(jù)離散點(diǎn)構(gòu)建三角網(wǎng)，然后插值生成等值線，這種算法應(yīng)用較廣[1-5]。對(duì)于規(guī)則的格網(wǎng)狀數(shù)據(jù)，由于其數(shù)據(jù)結(jié)構(gòu)的特殊性，可不用建立三角網(wǎng)而直接生成等值線。這種方法直觀，易編程實(shí)現(xiàn)，也有很多學(xué)者研究并實(shí)現(xiàn)[6-11]。本文結(jié)合開采沉陷數(shù)據(jù)處理特殊性，充分利用對(duì)角線上的數(shù)據(jù)，提高等值線精度，實(shí)現(xiàn)規(guī)則格網(wǎng)數(shù)據(jù)等值線的快速生成。

1 數(shù)據(jù)結(jié)構(gòu)設(shè)計(jì)

設(shè)預(yù)計(jì)點(diǎn)有N行M列，分別位于網(wǎng)格點(diǎn)上，其位置由所在的行號(hào)和列號(hào)確定，如(i，j)。通過行列號(hào)的引用，每個(gè)矩形網(wǎng)格可自動(dòng)分割成兩個(gè)三角形，實(shí)際并未實(shí)施剖分，故稱之為虛擬剖分，可大幅提高算法效率，如圖1所示。

圖1 規(guī)則網(wǎng)格點(diǎn)的虛擬三角剖分

設(shè)等值距為V，則可根據(jù)網(wǎng)格角點(diǎn)數(shù)據(jù)內(nèi)插生成每條邊(包括對(duì)角線)位于等值線上的點(diǎn)。當(dāng)網(wǎng)格角點(diǎn)屬性值恰好為等值距V的整數(shù)倍時(shí)，亦即網(wǎng)格點(diǎn)也是等值線上的點(diǎn)，將使等值線構(gòu)造變復(fù)雜。為了降低算法復(fù)雜性，將該角點(diǎn)屬性值加上一個(gè)不影響等值線生成精度的極小值ε，則可避免等值線過網(wǎng)格頂點(diǎn)時(shí)需多次進(jìn)行路徑判斷問題，且不影響等值線精度。

對(duì)每個(gè)網(wǎng)格內(nèi)的數(shù)據(jù)進(jìn)行以下約定：對(duì)于網(wǎng)格(i，j)，i、j為其左下角點(diǎn)的行號(hào)和列號(hào)，上部三角形編號(hào)為Ⅰ，下部三角形編號(hào)為Ⅱ，四條邊及對(duì)角線的編號(hào)分別為①、②、③、④、⑤，每個(gè)網(wǎng)格編碼方式相同，如圖2所示。

圖2 網(wǎng)格(i，j)內(nèi)三角形和邊編號(hào)規(guī)則

2 等值線生成算法

2.1 起始三角形的確定

從編號(hào)為(0，0)的網(wǎng)格開始遍歷，找到第一個(gè)需要繪制等值線的三角形所在網(wǎng)格。首先計(jì)算三角形三條邊上當(dāng)前等值線內(nèi)插值的個(gè)數(shù)n，n可能為1，2或3。(Bb、Be)分別用來標(biāo)記一個(gè)三角形內(nèi)等值線段起點(diǎn)和終點(diǎn)所在的邊的編號(hào)。當(dāng)n=1時(shí)，Bb=Be=插值點(diǎn)所在的邊的編號(hào)；當(dāng)n=2時(shí)，第一段等值線的可能編號(hào)方式為(①，②)、(②，⑤)、(⑤，①)(對(duì)于△Ⅰ)或(③，④)、(④，⑤)、(⑤，③)(對(duì)于△Ⅱ)，每個(gè)三角形有3種情形，若為逆時(shí)針，則將Bb和Be互換；當(dāng)n=3時(shí)，則直接將這三個(gè)點(diǎn)連接構(gòu)成閉合等值線。

2.2 下一個(gè)三角形的確定

如何尋找下一個(gè)要繪制的三角形是生成等值線算法的關(guān)鍵。當(dāng)Be=①時(shí)，下一個(gè)要繪制的三角形位于網(wǎng)格(i，j-1)，編號(hào)為Ⅱ；當(dāng)Be=②時(shí)，下一個(gè)要繪制的三角形位于網(wǎng)格(i+1，j)，編號(hào)為Ⅱ；當(dāng)Be=③時(shí)，下一個(gè)要繪制的三角形位于網(wǎng)格(i，j+1)，編號(hào)為Ⅰ；當(dāng)Be=④時(shí)，下一個(gè)要繪制的三角形位于網(wǎng)格(i-1，j)，編號(hào)為Ⅰ；當(dāng)Be=⑤，位于△Ⅰ時(shí)，下一個(gè)要繪制的三角形位于網(wǎng)格(i，j)，編號(hào)為Ⅱ；當(dāng)Be=⑤，位于△Ⅱ時(shí)，下一個(gè)要繪制的三角形位于網(wǎng)格(i，j)，編號(hào)為Ⅰ。

確定了下一個(gè)三角形位置后則可計(jì)算出新三角形中另一條等值線插值點(diǎn)所在的邊，存儲(chǔ)該點(diǎn)位置信息，并將其所在的邊的編號(hào)賦值給Be。

重復(fù)上述過程，即可實(shí)現(xiàn)等值線的自動(dòng)追蹤，從而完成等值線繪制。

2.3 等值線追蹤停止條件判斷

一條等值線追蹤停止或繪制完成的條件分以下2種情況：一種是當(dāng)前等值線段的終點(diǎn)所在邊的編號(hào)Be等于整條等值線的起點(diǎn)所在邊的編號(hào)Bb，則停止追蹤，此時(shí)等值線閉合，繪制完成，如圖3中等值線A所示；另一種情況是當(dāng)前等值線段的終點(diǎn)Be所在的邊為網(wǎng)格數(shù)據(jù)的邊界線時(shí)，停止追蹤。對(duì)于后一種情況，多數(shù)情況下只是繪制了半條等值線。此時(shí)，需將已生成的等值線反向存儲(chǔ)，并令Be=Bb，重復(fù)過程(2.2)和(2.3)，生成等值線的剩余部分，從而完成整條等值線的繪制，此時(shí)等值線不閉合，如圖3中等值線B所示。

圖3 等值線跟蹤

當(dāng)生成的等值線為逆時(shí)針存儲(chǔ)時(shí)，需將等值線端點(diǎn)逆序重新存儲(chǔ)，便于統(tǒng)一管理。

2.4 特殊等值線處理

在開采沉陷數(shù)據(jù)處理中，下沉邊界等值線是由下沉值為10 mm的點(diǎn)構(gòu)成，而不是下沉值為0 mm的點(diǎn)。同時(shí)，除下沉，傾斜、水平移動(dòng)、曲率等變形都有正負(fù)之分(極少數(shù)情況下下沉也會(huì)出現(xiàn))，需要對(duì)正負(fù)零變形邊界分別處理。

對(duì)于下沉邊界，可以生成下沉值為10 mm的等值線，將其屬性賦值為0 mm即可；而對(duì)于其他變形的正負(fù)0邊界，可以用±ε(ε是個(gè)極小值)這兩條特殊等值線代替，其屬性值分別賦值為±0。

2.5 程序流程圖

根據(jù)以上算法設(shè)計(jì)，由預(yù)計(jì)數(shù)據(jù)生成等值線的程序流程圖如圖4所示。

圖4 程序流程圖

3 等值線平滑

由于預(yù)計(jì)效率要求，預(yù)計(jì)網(wǎng)格點(diǎn)密度總是有限的，導(dǎo)致直接生成的等值線會(huì)呈現(xiàn)鋸齒狀，故需進(jìn)行平滑處理。在要求不高的情況下，可采用質(zhì)心平滑算法。即對(duì)于相鄰三點(diǎn)構(gòu)成的三角形，中間點(diǎn)的坐標(biāo)用三角形質(zhì)心坐標(biāo)代替。采用平滑算法處理等值線簡單且效果較好。

4 算法應(yīng)用

在進(jìn)行煤礦開采沉陷預(yù)計(jì)時(shí)，通常將預(yù)計(jì)點(diǎn)布設(shè)成網(wǎng)格狀，根據(jù)預(yù)計(jì)網(wǎng)格點(diǎn)的變形值繪制相應(yīng)的變形曲線。根據(jù)本文闡述的等值線生成算法，采用VC++ 6.0平臺(tái)，即可實(shí)現(xiàn)規(guī)則網(wǎng)格數(shù)據(jù)等值線自動(dòng)生成，等值線生成的函數(shù)代碼如下，程序界面如圖5所示。

圖5 等值線繪制程序界面

void DrawContours(CFile ConFile)

{

int i,j,k,Flag,oneconSize;

double Value;

SP ConP;

CArraytemparray;

CString Str1,ConLable,Formats;

GetBeginRowCol(0,0);

for(i=0;i

for(j=0;j<2*(nx-1);j++)

{

GetBeginRowCol(beginrow,begincol);

i=beginrow;

j=begincol;

if(j==2*(nx-1))

break;

ConP=GetFirstP(i,j);

Formats.Format("%%.%df ",3);

ConLable.Format(Formats,ConP.v);

Value=ConP.v;

onecon.RemoveAll();

flagone=0;

flagclose=0;

if(DrawFirstTri(i,j,Value))

{

Flag=DrawNextTri(i,j,Value);

while(Flag)

Flag=DrawNextTri(nr,nc,Value);

if((IsEqual(onecon.GetAt(onecon.GetSize()-1).x,onecon.GetAt(0).x))&&(IsEqual(onecon.GetAt(onecon.GetSize()-1).y,onecon.GetAt(0).y)))

flagclose=1;

if((flagclose!=1)&&(flagone!=1))

{

be=bb;

temparray.RemoveAll();

for(k=onecon.GetSize()-1;k>=0;k--)

temparray.Add(onecon.GetAt(k));

onecon.RemoveAll();

if(be==1)

{

GetPtAndSet0(grids[i][GridCol(j)].pny,grids[i][GridCol(j)].pointsy,Value,1,&ConP);

SetGridFlag(i,j);

if(j-1>=0)

SetGridFlag(i,j-1);

}

if(be==2)

{

GetPtAndSet0(grids[i][GridCol(j)].pnx,grids[i][GridCol(j)].pointsx,Value,1,&ConP);

SetGridFlag(i,j);

if(i-1>=0)

SetGridFlag(i-1,j-1);

}

if(be==3)

{

GetPtAndSet0(grids[i][GridCol(j)+1].pny,grids[i][GridCol(j)+1].pointsy,Value,1,&ConP);

SetGridFlag(i,j);

if(j+1<2*(nx-1))

SetGridFlag(i,j+1);

}

if(be==4)

{

GetPtAndSet0(grids[i+1][GridCol(j)].pnx,grids[i+1][GridCol(j)].pointsx,Value,1,&ConP);

SetGridFlag(i,j);

if(i+1

SetGridFlag(i+1,j+1);

}

if(be==5)

{

GetPtAndSet0(grids[i][GridCol(j)].pnxy,grids[i][GridCol(j)].pointsxy,Value,1,&ConP);

SetGridFlag(i,j);

if(Odd(j))

SetGridFlag(i,j-1);

else

SetGridFlag(i,j+1);

}

Flag=DrawNextTri(i,j,Value);

while(Flag)

Flag=DrawNextTri(nr,nc,Value);

for(k=0;k

temparray.Add(onecon.GetAt(k));

onecon.RemoveAll();

for(k=0;k

onecon.Add(temparray.GetAt(k));

}

if(flagclose==1)

oneconSize=onecon.GetSize()-1;

else

oneconSize=onecon.GetSize();

if(oneconSize>=2)

{

ConFile.Write(ConLable,ConLable.GetLength());

if(flagclose==1)

{

for(k=0;k

{

Str1.Format("%.3f,%.3f ",onecon.GetAt(k).x,onecon.GetAt(k).y); ConFile.Write(Str1,Str1.GetLength());

}

Str1="C ";

ConFile.Write(Str1,Str1.GetLength());

}

else

for(k=0;k

{

Str1.Format("%.3f,%.3f ",onecon.GetAt(k).x,onecon.GetAt(k).y);

ConFile.Write(Str1,Str1.GetLength());

}

}

}

else

{

if(onecon.GetSize()>=2){

ConFile.Write(ConLable,ConLable.GetLength());

for(k=0;k

{

Str1.Format("%.3f,%.3f ",onecon.GetAt(k).x,onecon.GetAt(k).y);

ConFile.Write(Str1,Str1.GetLength());

}

Str1="C ";

ConFile.Write(Str1,Str1.GetLength());

}

}

}

}

以某煤礦開采沉陷預(yù)計(jì)為例，單一工作面開采，預(yù)計(jì)網(wǎng)格數(shù)為40×50個(gè)，可繪制下沉等值線如圖6所示。

圖6 某煤礦開采沉陷預(yù)計(jì)下沉等值線圖

5 結(jié) 論

本文闡述了根據(jù)規(guī)則格網(wǎng)數(shù)據(jù)快速生成開采沉陷預(yù)計(jì)移動(dòng)變形等值線的生成算法原理，并利用VC++ 6.0編程實(shí)現(xiàn)。研究表明，根據(jù)文中算法，程序編寫較其他算法更為簡單, 生成的等值線圖效果較好，可用于其他規(guī)則形式數(shù)據(jù)處理。該算法主要處理對(duì)象為規(guī)則數(shù)據(jù)，具有專用性，適用范圍受到限制。對(duì)于任意離散數(shù)據(jù)，可以由線性內(nèi)插法生成規(guī)則數(shù)據(jù)后再應(yīng)用本算法。