1、 關于中國南海島礁巡航路線方案的探討 姓名：劉發強 2016年11月8日 摘要南海島礁巡航路線方案的確定實質上是一個求取平面有障礙區域兩點之間最短路徑問題和TSP問題的混合。針對前者，本文借助ESPO算法的思想，將外國非法侵占島礁所確定的12海里敏感區域抽象為正多邊形，在對所有頂點及起始點之間任一條路徑的合法性判定（借助計算幾何中點、線段、多邊形的位置關系）之后，利用Dijkstra算法生成障礙密集區（南沙群島）我方任意兩島礁之間的最短路徑，后對我方全體（西沙、南沙、中沙）島礁利用Floyd算法生成最短路徑完全圖的鄰接矩陣，為后續處理奠定基礎。針對后者，利用自適應的蟻群算法求取最優解，經檢驗效

2、果較好。同時借助openCV提供的圖像處理和畫圖函數對算法的運行和結果進行了可視化呈現，易于理解調試。理論上，本文的解決方案對于這類問題有很廣泛的適用性。關鍵詞：TSP問題；ESPO算法；Dijkstra算法；Floyd算法；openCV；蟻群算法；算法可視化 緒論 隨著中國海權意識的逐步覺醒和軍事力量的日益強大以及美國重返亞太戰略的實施，南海態勢在多方的博弈下波譎云詭，不安定不穩定性極大。為了維護我國的主權和領土完整，對我國實際控制的島礁進行周期性的補給和巡航勢在必行。因此，在周圍情況極其復雜的南海確定一個最優的巡航路徑非常關鍵。本文主要借助ESPO算法和蟻群算法對我國南海實際控制的包括西沙

3、、中沙、部分南沙在內的23個島礁進行了路徑規劃，比較圓滿的解決了上述問題，且具有良好的可擴展性和適用性。本文的基本假設：1、 地球表面近似為球面，在經線上，每緯度約111千米，在緯線上，每經度約為111*cos千米，其中為緯度。由于所考查區域（緯度從東經111°到東經117°，北緯8°到北緯18°）面積相對較小且在低緯度，故近似為平面即每經度111千米。利用百度地圖的測距工具和坐標拾取系統得出的結果與對平面上兩點歐幾里得距離相比最大相對誤差低于1%，這對于本問題的解決足夠準確。2、 島礁抽象為其幾何中心處一點。3、 本文考察西沙群島、中沙群島的所有島礁和

4、南沙群島的7個島礁。由于中沙群島大部分在中部呈一緊湊環裝分布，為了簡化且不失準確性，將大環礁抽象為環北一點和環南一點。4、 本文所規劃的路徑除不可避免之外不進入外國侵占島礁所形成的12海里敏感區域，由于東門礁、南薰礁距離鴻庥島和九章群礁過近，在此兩處將12海里放寬為6海里。5、 由于部分島礁之間平面歐幾里得距離遠大于其面積尺度，故在無障礙區即西沙、中沙我方島礁的最短距離可以由平面歐幾里得距離代替。并且最優解中，可以預見路徑穿越島礁的可能性極小，再次印證上述假定的合理性。6、 在構造最短路徑完全圖的過程中，對于非南沙與非南沙之間、非南沙和南沙之間其距離認定為其平面歐幾里得距離，南沙群島各島礁距離

5、利用ESPO算法生成，由于南沙距離非南沙相對于其尺度很大，故近似比較準確。方案形成過程南海島礁經緯度采集和部分島礁之間的測距借助百度地圖及其提供的坐標拾取系統和測距工具，獲取了我方23個島礁、他方18個島礁的經緯度信息和用于連通南沙我方部分島礁（障礙外部）與中南暗沙、大環礁、中建島、盤石嶼等島礁的距離。并利用MATLAB計算了我方非南沙各島礁之間的歐氏距離，從而生成了一部分最短路徑完全圖。由于表格較大，不在此呈現，詳見表格。圖 1南海部分島礁南沙群島有障礙區域我方島礁最短路徑的形成由假設知，各敏感區域是由以外國侵占島礁幾何中心為圓心半徑為12海里的區域，為凸集。若可行區域的邊界是光滑曲面。則最

6、短路徑必由下列弧組成，它們或者是空間中的自然最短曲線，或者是可行區域的邊界弧。而且，組成最短路徑的各段弧在連接點處必定相切。上述定理由J.W.Craggs在1973年證明。由于光滑曲線在計算機中不便于存儲和處理，本文將敏感區域近似為其內接正多邊形（本文中最終使用以12邊形近似的數據），并基于上述定理近似推斷在障礙區兩點之間的最短路徑必有以下線段組成，他們或是障礙多邊形的邊，或是障礙多邊形各頂點及起始點之間的合法線段組成，其中合法線段為不與任何障礙多邊形相交或在其外部的線段。因此，在上述分析的支持下，試圖用計算機解決這個平面簡單凸多邊形障礙區域最短路徑問題首先要解決如何判斷平面上點、線段與多邊形

7、的位置關系。這是一個基本的計算幾何問題，對于線段和多邊形位置關系的判定，首先要解決線段與線段之間、點與多邊形之間的位置關系判定方法。對于點與多邊形，可以利用射線法，即由該點為端點任意引一條射線，若射線與多邊形各邊相交奇數次，則在內部，其他在外部（在邊上也認為在外部）；對于線段與線段之間，可利用正交變換（便于處理）將其中一條變換到X軸且一段與原點重合，進而判斷另一條線段與X軸的交點與第一條線段的關系，若在其上則相交，否則不想交（在本問題中如此判斷合適）。所以對于線段與多邊形，若線段與多邊形任一條邊相交，則線段與多邊形相交，否則進而判斷線段中點與多邊形的位置關系（由于多邊形凸，在這種情況下，中點與

8、多邊形的位置關系和線段與多邊形關系等價），若在外，則在外，反之易得。上述算法可以判定某條邊的合法性，若不合法，則距離為無窮大，否則為歐氏距離。在上述算法的支撐下，可以求得所有障礙多邊形頂點及起始點互連而成的所有距離。再借助Dijkstra算法生成我方7個島礁的最短路徑。Dijkstra算法的思想是若從某點A到連通圖上其他點最短路徑生成是遞增順序的的，則到下一點X待生成的最短路徑或者是A-X，或者經過已生成最短路徑的點。下面是結果呈現，數據在表格中給出。圖 2 6邊形近似圖 3 12邊形近似圖 4南沙所有最短路徑 圖 5永暑礁到東門礁三、測量中沙、西沙部分島礁到南沙7個島礁的距離由于中沙、西沙基

9、本被我方控制，且距離南沙島礁較遠，因此，中沙、西沙到南沙部分島礁的距離可以用歐氏距離代替，下面分別測量了從中建島、盤石嶼、中南暗沙、中沙大環礁、黃巖島到美濟礁、渚碧礁、永暑礁的距離，用于連通整個賦值圖。然后利用Floyd算法生成最短路徑完全圖Dist2323和路徑中繼前驅Path2323,其中，Distij為島礁i與島礁j的最短距離，它或是歐氏距離，或是ESPO算法生成的距離，或是經過中繼獲得的最短距離。Pathij表是從島礁i到島礁j最短路徑中的中級前驅，如果沒有中繼則為i。借助Path矩陣可以在生成TSP路徑后提供詳細路徑信息。圖 5用于連通南沙與西沙、中沙的路徑1、 利用蟻群算法求解di

10、st矩陣的最優TSP解蟻群算法是求解TSP問題的經典算法，是一種元啟發式搜索算法，借助群體智能來實現組合優化。但是參數的選取并不容易，算法陷入局部最優的可能性很大，影響全局搜索能力。幾個重要的參數包括信息素重要程度指標alpha、可見度（距離倒數）重要程度指標beta、信息素揮發因子dispersefactor、螞蟻數量antnum和螞蟻圈最大信息素Q。其中alpha=0時，算法退化為貪心算法，beta=0時，算法退化為隨機搜索。在實驗過程中試驗了大量參數，發現螞蟻數量過多越容易陷入局部最優。最終選取alpha=0.8，beta=2，antnum=34，dispersefactor=0.5，并

11、引入dispersefactor的陷入局部最優的自適應調節機制即dispersefactor<=(dispersefactor*limit)0.5，該迭代修正將最終使dispersefactor趨于limit不至于修正過度，取limit為0.4。下面是最優解，為3235.49千米。所有圖片都在附件中給出。 結論 本文完滿的解決了既定問題，具有良好的適用性。 參考文獻 王萬良 人工智能及其應用（第三版）群智能算法及其應用 Blog：線與多邊形位置關系 附錄(源程序）1、 Floyd算法 path=zeros(23,23); dist=zeros(23,23); for i=1:23 for

12、 j=1:23 dist(i,j)=original(i,j); if dist(i,j)<inf path(i,j)=i; else path(i,j)=-1; end if i=j dist(i,j)=0; end end end for k=1:23 for i=1:23 for j=1:23 if dist(i,j)>dist(i,k) + dist(k,j) dist(i,j) = dist(i,k) + dist(k,j); path(i,j) = path(k,j); end end end End二:ESPO算法含Dijkstra算法、幾何判定算法、可視化（運行平臺

13、需要支持opencv）#include<iostream>#include<opencv2opencv.hpp>#define pi 3.1415926#define inf 9999999999using namespace std;cv:Mat mapaa;class pointpublic:double x, y;string name;point()x = y = 0;point(double x, double y)this->x = x;this->y = y;void setname(string name)this->name = na

14、me;class Linepublic:point start, end;double length;Line()Line(point start, point end)this->start = start;this->end = end;length = sqrt(start.x - end.x)*(start.x - end.x) + (start.y - end.y) *(start.y - end.y);void setlength()length = sqrt(start.x - end.x)*(start.x - end.x) + (start.y - end.y)

15、*(start.y - end.y);template <int edgenum>class Polygon/whose methods appropriate for sample convex polygonspublic:int edgen;point vertaxedgenum;Line edgeedgenum;Polygon()void setPolygon(point vertax)edgen = edgenum;for (int i = 0; i < edgenum; i+)this->vertaxi = vertaxi;for (int i = 0; i

16、 < edgenum; i+)edgei.start = vertaxi;edgei.end = vertax(i + 1) % edgen;edgei.setlength();int intersect(Line line1, Line line2)subtract(line1.end, line1.start);subtract(line2.end, line1.start);subtract(line2.start, line1.start);subtract(line1.start, line1.start);double theta = -arctan(line1.end.x,

17、 line1.end.y);rotate(line1.end, theta);rotate(line2.start, theta);rotate(line2.end, theta);if (line2.start.y*line2.end.y < 0)double x0 = line2.start.x - line2.start.y*(line2.end.x - line2.start.x) / (line2.end.y - line2.start.y);if (x0 - line1.end.x)*x0 < 0)return 1;if (fabs(x0) < 0.000001

18、| fabs(x0 - line1.end.x) < 0.000001)return 2;else if (fabs(line2.start.y*line2.end.y) < 0.000001)return 3;elsereturn 0;point intermediate(point p1, point p2)point temp;temp.x = (p1.x + p2.x) / 2;temp.y = (p1.y + p2.y) / 2;return temp;bool interconnect(Line line)for (int i = 0; i < edgen; i+

19、)if (line.length < 0.00001)return false;if (intersect(line, edgei) = 1)return true;if (inside(intermediate(line.start, line.end) = false)return false;elsereturn true;bool inside(point p)int cnt = 0;int lock = 0;point auxi;auxi.x = -500;auxi.y = -300;Line ray(p, auxi);for (int i = 0; i < edgen;

20、 i+)Line temp = edgei;int value = intersect(ray, temp);if (fabs(temp.length) < 0.00001)return false;if (fabs(p.y - temp.start.y)*(temp.end.x - temp.start.x) - (temp.end.y - temp.start.y)*(p.x - temp.start.x) <= 0.00001)/on th edgereturn false;if (value = 1|value=3)cnt+;else if (value = 2)if (l

21、ock = 0)cnt+;lock = 1;elselock = 0;if (cnt % 2 = 1)return true;elsereturn false;double arctan(double x, double y)if (x = 0 && y = 0)cout << "error" << endl;elseif (x > 0)return atan(y / x);else if (x = 0)if (y > 0)return pi / 2;else if (y < 0)return -pi / 2;els

22、e if (y >= 0)return atan(y / x) + pi;else if (y < 0)return atan(y / x) - pi;void rotate(point &p1, double theta)point temp = p1;p1.x = temp.x*cos(theta) - temp.y*sin(theta);p1.y = temp.y*cos(theta) + temp.x*sin(theta);void subtract(point& p1, point p2)p1.x -= p2.x;p1.y -= p2.y;template

23、<int edgenum>class enemypublic:point center;double radius;int edgen;string name;point vertaxedgenum;void setenemy(point enemy,double radius,string name)this->center = enemy;this->radius = radius;this->name = name;edgen = edgenum;for (int i = 0; i < edgen; i+)vertaxi.x = center.x +

24、radius*cos(2*pi / edgenum* i);vertaxi.y = center.y + radius*sin(2*pi / edgenum * i);enemy()void setradius(double radius)this->radius = radius;for (int i = 0; i < edgen; i+)vertaxi.x = center.x + radius*cos(2*pi / edgenum * i);vertaxi.y = center.y + radius*sin(2*pi / edgenum * i);void drawarea(

25、)for (int i = 0; i < edgen; i+)line(mapaa, cv:Point(vertaxi.x, vertaxi.y), cv:Point(vertax(i + 1) % edgen.x, vertax(i + 1) % edgen.y), cv:Scalar(052, 25, 65);double euclidean_metric(point p1, point p2)return sqrt(p1.x - p2.x)*(p1.x - p2.x) +(p1.y - p2.y)*(p1.y - p2.y);voidespo(enemy<12>*ene

26、, int enenum, point china,int chinanum)Polygon<12> *poly = new Polygon<12>enenum;for (int i = 0; i < enenum; i+)polyi.setPolygon(enei.vertax);enei.drawarea();int vertaxnum = enenum * ene0.edgen + chinanum;point*vertax = new pointvertaxnum;for (int i = 0; i < chinanum; i+)vertaxi =

27、chinai;for (int i = 0; i < enenum * ene0.edgen; i+)vertaxi + chinanum = enei / ene0.edgen.vertaxi % ene0.edgen;double *dist = new double*vertaxnum;for (int i = 0; i < vertaxnum; i+)disti = new doublevertaxnum;for (int i = 0; i < vertaxnum; system("cls"), cout<<"已完成"

28、;<<i<<"/"<<vertaxnum,i+)for (int j = 0; j < vertaxnum;j+)if (j = i)distij = inf;continue;Line temp(vertaxi, vertaxj);int lock = 0;for (int t = 0; t < enenum; t+)if (erconnect(temp)=true)lock = 1;break;if (lock = 1)distij = distji = inf;else distij = distji

29、= euclidean_metric(vertaxi, vertaxj);for (int cnc = 0; cnc < chinanum; cnc+)int *generated = new intvertaxnum;for (int i = 0; i < vertaxnum; i+)generatedi = 0;int *pre = new intvertaxnum;double *mindist = new doublevertaxnum;for (int i = 0; i < vertaxnum; i+)generatedi = 0;if (distcnci <

30、 inf)prei = cnc;elseprei = -1;mindisti = distcnci;mindistcnc = 0;for (int i = 0; i < vertaxnum; i+)double min = inf;int temp = -1;for (int j = 0; j < vertaxnum; j+)if (generatedj = 0 && mindistj <= min)min = mindistj;temp = j;generatedtemp = 1;for (int j = 0; j <vertaxnum; j+)if

31、(!generatedj && mindistj>mindisttemp + disttempj)mindistj = mindisttemp + disttempj;prej = temp;for (int i = cnc+1; i < chinanum; i+)/cv:Mat mappp = mapaa.clone();int pathnum = 0;int *path = new intvertaxnum;cout << "從 " << << " 到 " &

32、lt;< << "距離是" << mindisti/800*397.7 <<"千米"<< endl;int index = prei;pathpathnum = i;while (index != -1)pathnum+;pathpathnum = index;index = preindex;for (int ij = 0; ij < pathnum; ij+)circle(mapaa, cv:Point(int)vertaxpathij.x, (int)vertaxp

33、athij.y), 5, cv:Scalar(i*10, i, i*12);line(mapaa, cv:Point(int)vertaxpathij.x, (int)vertaxpathij.y), cv:Point(int)vertaxpath(ij + 1) % (pathnum + 1).x, (int)vertaxpath(ij + 1) % (pathnum + 1).y), cv:Scalar(i*10, 45, 255);ostringstream os;os << mindisti / 800 * 397.7;string str= os.str();/cv:im

34、show("wall1", mapaa);cv:imwrite("完整.bmp", mapaa);/getchar();int main()double length = 394.8, width = 267.7, baselong = 112.28, baselati = 11.20, radius = 44.11, mapl = 800, maph = 536;string obstacle = "中業島","大現礁","小現礁","福祿礁","鴻庥島"

35、;,"九章群礁","安達礁","舶蘭礁","三角礁","畢生礁","雙黃沙洲","庫歸礁","仁愛礁","仙娥礁","馬歡島","五方礁","火艾礁","蒙自礁" ;double enelong = 114.29,113.85,113.99,113.62,114.36,114.43,114.71,114.58,115.34,113.67,

36、114.42,114.61,115.86,115.45,115.86,115.72,114.90,114.78 ;double enelati = 11.03,10.05,9.94,10.18,10.15,9.82,10.30,10.37,10.18,8.94,10.64,10.81,9.7,9.35,10.66,10.47,10.83,11.09 ;string china = "赤瓜礁","美濟礁","永暑礁","華陽礁","南薰礁","渚碧礁","東門礁&qu

37、ot; ;double cnlong = 114.28,115.54,112.96,112.83,114.22,114.08,114.57 ;double cnlati = 9.68,9.86,9.59,8.85,10.17,10.88,9.92 ;point*chin = new point7;enemy<12> enearea18;mapaa = cv:imread("南沙諸島.png");for (int i = 0; i < 18; i+)eneareai.setenemy(point(enelongi - baselong) * 111 / le

38、ngth*mapl, (baselati - enelatii) * 111 / width*maph),radius,obstaclei);enearea4.setradius(20);enearea5.setradius(20);for (int i = 0; i < 7; i+)chini.x = (cnlongi - baselong) * 111 / length*mapl;chini.y = (baselati - cnlatii) * 111 / width*maph;chini.setname(chinai);espo(enearea, 18, chin, 7);三、蟻群

39、算法含可視化#include<iostream>#include<stdlib.h>#include<time.h>#include<opencv2opencv.hpp>using namespace std;#define size 23int Q = 1000;double Q0 = 100;#define antnum 34double alpha = 0.8;double beta = 2;#define cir 10000double dispersefactor = 0.5;#define original 100cv:Mat bac

40、kground;double distsizesize;double pheromonesizesize;int presizesize;string island = "永興島","西沙洲","銀礫灘", "東島", "甘泉島", "濱湄灘", "湛涵灘", "玉琢礁", "華光礁", "盤石嶼", "中建島", "浪花礁", "黃巖島"

41、;, "中沙環礁北", "中沙環礁南","中南暗沙", "永暑礁", "華陽礁", "赤瓜礁", "美濟礁", "東門礁", "南薰礁","渚碧礁" ;double longitude = 112.35,112.25,112.24,112.73,111.68,112.47,112.56,112.03,111.71,111.79,111.2,112.52,117.76,114.71,114.06,1

42、15.44,112.96,112.83,114.28,115.54,114.57,114.22,114.08 ;double latitude = 16.84,16.98,16.76,16.63,16.51,16.33,16.4,16.33,16.22,16.03,15.77,16.03,15.15,16.07,15.51,13.94,9.59,8.85,9.68,9.86,9.92,10.17,10.88 ;class pointpublic:double x, y;string name;point()void setinfor(string name, double x, double

43、y)this->name = name;this->x = x;this->y = y;point islansize;class Antpublic:double possibilitysize;int currentpos;int time;double pathlength;int pathsize;Ant()void init()currentpos = rand() % size;time = 0;pathlength = 0;path0 = currentpos;for (int i = 0; i < size; i+)possibilityi = 0;bo

44、ol haspassed(int nextpos)for (int i = 0; i < time; i+)if (nextpos = pathi)return true;return false;int nextpos()double pm = 0;for (int i = 0; i < size; i+)if (haspassed(i) = false)pm += possibilityi=pow(pheromonecurrentposi, alpha)*pow(distcurrentposi, -beta);else possibilityi = 0;int posssize

45、;int sum = 0;for (int i = 0; i < size; i+)possibilityi /= pm;sum+=possi = (int)(possibilityi*cir);int x;int randnum = rand() % sum;sum = 0;for (int i = 0; i < size; i+)sum += possi;if (randnum < sum)return i;int i = 0;while (haspassed(i) = true)i+;return i;void move()time+;int next= nextpos

46、();pathtime = next;pathlength += distcurrentposnext;currentpos = next;if (time = size - 1)pathlength += distnextpath0;time = 0;bool hasgonethrough(int x, int y)for (int i = 0; i < time; i+)if (pathi = x&&path(i + 1) % size = y)return true;return false;void getpath(int x, int y, int&nu

47、m, int path)num = 0;int index = y;path0 = y;donum+;index = prexindex;pathnum =index; while (index!= x);for (int i = 0; i <=num/ 2; i+)int temp = pathi;pathi = pathnum - i;pathnum - i = temp;class antColonypublic:Ant antantnum;double optimal;int x;int n0;int opathsize;int cnt;int flag;antColony()v

48、oid init()flag = 0;cnt = 0;optimal = inf;n0 = 5;void findoptimal()double min = optimal;for (int i = 0; i < antnum; i+)if (anti.pathlength < min)min = anti.pathlength;x = i;for (int j = 0; j < size; j+)opathj = antx.pathj;flag = 0;if (fabs(min - optimal) < 0.000001)flag+;optimal = min;voi

49、d run(int n)for (int i = 0; i < n; i+)round(1000);show();draw();void show()cout << optimal << "千米" << " "for (int i = 0; i < size; i+)cout << << " "cout << endl;cout << alpha << " " << beta << " " << dispersefactor << endl;void update()for (int i = 0; i < size; i+)for (int j = 0; j < size; j+)double delta = 0;for (int k = 0; k < antnum; k+)if (antk.hasgone