黃 磊，王夢雪，李 淼

（成都理工大學 地球探測與信息技術(shù)教育部重點實驗室，四川 成都610059）

0 引 言

隕石坑是月表最明顯的形貌特征，它對人們研究月球演化、太陽系隕石撞擊歷史和空間環(huán)境等都有重要意義，在隕石坑的研究中，隕石坑影像幾何特征的信息提取是一種重要手段［1］。隨著科技的進步，很多國家發(fā)射了大量的月球探測器，如何利用月球探測器的高分辨率圖像研究月表隕石坑的幾何特征和月表年齡關(guān)系，已成為一個研究的難點與熱點問題。許多學者從不同角度進行了分析研究，如使用隕石坑大小頻度法估測月表年齡，建立表暴露年齡與單位面積中隕石坑數(shù)量關(guān)系等［2，3］。但是，這些研究只是定性描述隕石坑幾何特征，沒有對隕石坑幾何特征進行系統(tǒng)的分類以及定量描述隕石的某一幾何特征。

成都理工大學大學空間信息研究室在此研究基礎(chǔ)上，利用 “嫦娥一號”衛(wèi)星CCD 立體相機獲取的月表遙感影像，研究隕石坑的幾何特征，提出了一種 “六位數(shù)字十級分級法”，用0－9個數(shù)字量化描述隕石坑的幾何特征，為描述隕石坑的影像和幾何特征建立模型和標準。六位數(shù)字分別是直徑、邊緣清晰度、邊緣圓度、坑唇寬度、坑底結(jié)構(gòu)復(fù)雜度和高程。本文是研究邊緣圓度模型的建立以及邊緣圓度的評價算法，依據(jù)圓度誤差的定義和評定方法，分別用最小區(qū)域圓法、最小二乘圓法、最小外接圓法和最大內(nèi)接圓法對 “嫦娥一號”CCD 衛(wèi)星影像的隕石坑邊緣進行圓度誤差判定，然后利用改進的BP神經(jīng)網(wǎng)絡(luò)，組合4種評價算法得到的圓度誤差的結(jié)果作為圓度等級的評價指數(shù)。最后，把圓度評價模型加載到自主開發(fā)的數(shù)字月球平臺上做結(jié)果展示，圓度等級設(shè)為0－9級，對隕石坑的邊緣形變做10級量化。

1 基本思路

1.1 算法綜述

圓度誤差是指被測實際輪廓相對理想圓的變動量，是一項重要的形位精度指標，是評定許多其它形位誤差的基礎(chǔ)，關(guān)于圓度誤差的評定方法研究很多，本文所用的圓度考量因子是圓度誤差。根據(jù)國際標準化組織新頒發(fā)的標準規(guī)定了有5種圓度誤差的評定方法，分別是最小二乘圓法（the least square circle method，LSCM）、最小區(qū)域圓法（the minimum zone circle method，MZCM）、最小外接圓法（the minimum circumscribed circle method，MCCM）、最大內(nèi)接圓法（the maximum inscribed circle method，MICM）和切比雪夫擬合圓法［4］。本文提出一種基于圓度誤差的隕石坑邊緣圓度評定算法，該算法對圖像增強和降噪平滑后的圖像做邊緣提取，然后用圓度誤差計算方法進行判定，得出邊緣圓度量化評價結(jié)果。對圓度誤差計算方法的選擇，本文選用的評定方法分別是最小二乘法、最小區(qū)域法、最小外接圓法和最大內(nèi)接圓法。最后，用標準數(shù)據(jù)對神經(jīng)網(wǎng)絡(luò)進行訓(xùn)練，將4種評價結(jié)果作為網(wǎng)絡(luò)的輸出，得到量化的規(guī)則度評價等級，算法流程如圖1所示。

圖1 邊緣圓度評價模型算法流程

1.2 邊緣圓度的定義

邊緣圓度的判定是對月球隕石坑的CCD 衛(wèi)星圖像的邊緣進行圓度誤差的判定，圓度誤差作為隕石坑形變的定性判定，誤差越大，邊緣圓度就越差，誤差越小，邊緣圓度就越好。以 “六位數(shù)字十級分級法”的角度出發(fā)，從定量來看，分為0 到9 級，圓度誤差逐漸減小，邊緣圓度逐漸提高。

1.3 圖像處理方法

1.3.1 圖像預(yù)處理“嫦娥一號”CCD 圖像的預(yù)處理對邊緣圓度的判定有影響，本文圖像的預(yù)處理主要采用圖像增強和雙邊濾波。圖像增強是根據(jù)局部灰度值的統(tǒng)計數(shù)據(jù)對選取原始圖像進行增強處理，增大圖像的對比度，有選擇地圖像中隕石坑的幾何特征，同時抑制圖像中的其它信息。

為消除圖像噪聲對月表隕石坑的影響，采用非線性、非迭代的雙邊濾波方法處理圖像。雙邊濾波在處理相鄰像素時，既要考慮到距離上的鄰近關(guān)系，還要考慮灰度上的相似度，通過對這兩者的非線性組合，可以自適應(yīng)地進行濾波。圖像在濾除噪聲的同時，能很好地保留圖像的邊緣信息。預(yù)處理結(jié)果如圖2所示。

圖2 原始圖像和處理后效果圖像

1.3.2 邊緣提取

邊緣提取的精度對邊緣圓度的判定有影響，本文分析多種邊緣檢測的方法后，采用Canny算法提取月球隕石坑的邊緣。由于Canny算子的計算量大，實時性較差，檢測邊緣的精度沒有達到單像素級，因而采用改進型Canny邊緣檢測算法提取邊緣及細化［5］。值得注意的是，因為太陽高度角、CCD 鏡頭拍攝角度等原因，隕石坑內(nèi)經(jīng)常會出現(xiàn)高亮區(qū)和陰影區(qū)，二者同時出現(xiàn)且它們之間的連線與太陽光入射矢量基本保持平行，如圖2（a）所示。特別是在小型隕石坑內(nèi)，由于半徑小，當太陽光高度角接近90°時，會出現(xiàn)大面積的高亮區(qū)和陰影區(qū)，幾乎占到坑底的一般面積，并覆蓋大部分坑唇。提取的隕石坑邊緣中既包括了真實的邊緣信息，又包括了隕石坑內(nèi)部高亮區(qū)和陰影區(qū)的偽邊緣信息以及其它地形特征的邊緣信息，如圖3（a）所示。通過檢測隕石坑圖像的梯度方向，發(fā)現(xiàn)在隕石坑內(nèi)部的偽邊緣與隕石坑的真實邊緣的梯度方向正好相反，真實隕石坑的邊緣梯度方向與光源方向的夾角小于90°，但是偽邊緣梯度方向與光源方向的夾角大于90°，即真實邊緣滿足式（1）

圖3 偽邊緣去除效果

2 邊緣圓度評價算法

邊緣圓度的評價實質(zhì)上是根據(jù)圓度誤差的定義構(gòu)造函數(shù)模型，然后進行函數(shù)的優(yōu)化求解［6，7］。最小區(qū)域圓法（MZCM）是以最小區(qū)域圓（MZC）作為評定的基準圓的方法，按此方法求得的圓度誤差值ΔZZ

式中： （Ri）max、 （Ri）min——隕石坑邊緣上取得的樣點相對最小區(qū)域圓（MZC）的最大、最小偏離值；Ri——隕石坑邊緣上取得的樣點相對最小區(qū)域圓（MZC）的偏離值。同樣將最小二乘圓（LSC）、最小外接圓（MCC）、最大內(nèi)接圓（MIC）作為評定的基準圓求得的圓度誤差分別為ΔZs、ΔZc、ΔZi。

2.1 圓度誤差的數(shù)學模型

根據(jù)圓的標準方程，最小區(qū)域圓、最小二乘圓、最小外接圓和最大內(nèi)接圓的一般方程見式（3）

（1）最小二乘圓的標準方程見式（3），在隕石坑的邊緣上取得n （n ＞3） 個樣點，各個點的坐標記為xi、yi（i＝1，2，…，n） ，則各個點到最小二乘圓的圓心的距離為

根據(jù)最小二乘圓法構(gòu)造的目標函數(shù)見式（5），該函數(shù)的實質(zhì)是求解圓心的坐標xc、yc和半徑Rc的最優(yōu)化問題，求解結(jié)果、、即為最小二乘圓的圓心坐標值和半徑值

（2）最小區(qū)域圓的標準方程見式（3），隕石坑的邊緣上各個樣點到最小區(qū)域圓的圓心的距離為

根據(jù)最小區(qū)域圓法構(gòu)造的目標函數(shù)見式（7），該函數(shù)的實質(zhì)是求解圓心的坐標xc、yc的最優(yōu)化問題，求解結(jié)果、即為最小區(qū)域圓的圓心坐標值

式中： ［Ri］max、 ［Ri］min——隕石坑邊緣上取得的樣點與圓心距離的最大值和最小值。

（3）最小外接圓的標準方程見式（3），隕石坑的邊緣上各個樣點到最小外接圓的圓心的距離為Ri見式（6）。根據(jù)最小外接圓法構(gòu)造的目標函數(shù)見式（8），該函數(shù)的實質(zhì)是求解圓心的坐標xc、yc的最優(yōu)化問題，求解結(jié)果、即為最小外接圓的圓心坐標值

（4）最大內(nèi)接圓的標準方程見式（3），隕石坑的邊緣上各個樣點到最大內(nèi)接圓的圓心的距離為Ri見式（6）。根據(jù)最大內(nèi)接圓法構(gòu)造的目標函數(shù)見式（9），該函數(shù)的實質(zhì)是求解圓心的坐標xc、yc的最優(yōu)化問題，求解結(jié)果、即為最大內(nèi)接圓的圓心坐標值

2.2 圓度誤差的求解

圓度誤差的求解過程為：

（1）在隕石坑的邊緣上取得的n （n ＞3） 個樣點，各個點的坐標記為xi、yi（i＝1，2，…，n） ，選取其中3個點，分別記為P1（x1，y1）、P1（x2，y21）、P1（x3，y3），將它們帶入方程組 （10），可以得到圓心坐標初值x0、y0和半徑初值R0

（2）將求到的圓心坐標初值x0、y0和半徑初值R0代入構(gòu)造函數(shù)式（5），同時將圓心坐標初值x0、y0分別代入構(gòu)造函數(shù)式 （7）、式 （8）、式 （9），然后進行迭代求解，分別得到滿足迭代求解精度的函數(shù)解、、和、，即最小二乘圓的圓心坐標值、半徑值和最小區(qū)域圓、最小外接圓、最大內(nèi)接圓的圓心坐標值。

（4）分別將求得的最小二乘圓法、最小區(qū)域圓法、最小外接圓法、最大內(nèi)接圓法相對應(yīng)的Ri（i＝1，2，…，n） 的值代入式（2），即可分別得到圓度誤差ΔZs、ΔZZ、ΔZc、ΔZi。

3 基于BP網(wǎng)絡(luò)的邊緣圓度等級判定

上述的4種圓度評定方法由于使用不同的理想圓，評定結(jié)果會有所差別，考慮到月表的復(fù)雜性和評價結(jié)果的局部誤差，只對4種結(jié)果進行加權(quán)求和做歸一化處理得到的圓度的等級缺乏準確性。如何建立一種標準的圓度等級劃分模型對檢測樣本進行判定，本文借助于BP神經(jīng)網(wǎng)絡(luò)良好的復(fù)雜非線性問題處理和糾錯能力，設(shè)計一種適用于月表情況BP神經(jīng)網(wǎng)絡(luò)判定算法［8，9］。

3.1 建立學習樣本

“嫦娥一號”CCD 相機拍攝的遙感影像在月球的中低緯度地區(qū)的成像質(zhì)量更高，本文的實驗對象選擇這個區(qū)域。首先選取10個標準隕石坑為樣本，樣本的標準為隕石坑周圍地形起伏不大，坑唇較窄，坑底較為平緩，然后對這些樣本進行人為的模糊處理，劃分成圓度好、一般、差3 個圓度等級。通過4種評價方法分別對樣本進行評價，再將不同圓度等級的樣本評價結(jié)果通過線性內(nèi)插為10 個評價值，對應(yīng)相應(yīng)的0－9級。然后將得到的100個評價值用以建立學習樣本S＝s i，（ ）j4×100，神經(jīng)網(wǎng)絡(luò)的輸入數(shù)據(jù)應(yīng)進行歸一化，再把神經(jīng)網(wǎng)絡(luò)的期望輸出設(shè)置為圓度Y ＝［0，0.1…，0.9，…，0，0.1，…，0.9］1×100。傳統(tǒng)的最大最小值法易使數(shù)據(jù)歸一化為傳遞函數(shù)Sigmoid的極大值0和極小值1，導(dǎo)致連接權(quán)值應(yīng)足夠大才能使神經(jīng)網(wǎng)絡(luò)的輸出值與之相匹配，進而使訓(xùn)練次數(shù)增多和訓(xùn)練時間變長。為了避免這個現(xiàn)象，可以把區(qū)間數(shù)據(jù)規(guī)范到［0.15，0.85］區(qū)間，本文使用改進型歸一化公式，見式（11），最小二乘圓法評價作為逆指標，歸一化公式見式（12）

3.2 神經(jīng)網(wǎng)絡(luò)訓(xùn)練

BP神經(jīng)網(wǎng)絡(luò)分為三層，分別為輸入層、中間層（或隱含層）、輸出層［10］。其中，輸入層X1、X2、X3、X4分別是4種評價算法的值，中間層（或隱含層）的各個神經(jīng)元通過Wij（i＝1，2，…，n；j＝1，2，…，p） 和輸入層建立權(quán)重連接，通過Vjk（j＝1，2，…，p；k＝1，2，…，q）和輸出層建立連接，輸出層Yi（i＝0，1，…，9）是最終的評級結(jié)果，如圖4所示。

從圖4可知，神經(jīng)網(wǎng)絡(luò)中各層單元都與上一層或下一層建立連接，但是各層之間并無連接。神經(jīng)網(wǎng)絡(luò)學習的基本步驟為：

（1）正向傳播，將輸入層輸入的學習樣本通過權(quán)值連接送入中間層，經(jīng)中間層的傳遞函數(shù)激活，再通過連接權(quán)值送入到輸出層，經(jīng)輸出層的傳遞函數(shù)激活，的得到輸出結(jié)果；

（2）反向調(diào)整，BP 神經(jīng)網(wǎng)絡(luò)是一種前饋型神經(jīng)網(wǎng)絡(luò)，定義了一種誤差函數(shù)E，若E 不滿足預(yù)定精度，則需要不斷調(diào)整連接權(quán)值和閥值減少誤差；

（3）循環(huán)調(diào)整，調(diào)整連接權(quán)值和閥值，直到誤差函數(shù)E滿足要求且需要設(shè)置迭代上限避免學習陷入死循環(huán)。

圖4 神經(jīng)網(wǎng)絡(luò)

如圖5所示，其中，圖5（a）是神經(jīng)網(wǎng)絡(luò)的誤差變化曲線，誤差曲線表明訓(xùn)練了972步后，輸出誤差小于誤差上限0.001，圖5（b）是神經(jīng)誤差的誤差檢驗結(jié)果圖，是為了檢驗神經(jīng)網(wǎng)絡(luò)的準確性，再次把學習樣本輸入網(wǎng)絡(luò)進行檢驗，得到的散點圖，散點圖表明圓度的評價等級檢驗的絕對誤差都小于0.1，因而神經(jīng)網(wǎng)絡(luò)的準確性得到了保證。

圖5 神經(jīng)網(wǎng)絡(luò)的訓(xùn)練結(jié)果

3.3 評價結(jié)果展示

本文的實驗選取了30個隕石坑，為了使圓度判定的結(jié)果與主觀的判定相互印證，從中又選取了10幅能夠主觀上判斷圓度等級的影像在本文中展示，相關(guān)圖像和判定數(shù)據(jù)見圖6和表1。表1中的神經(jīng)網(wǎng)絡(luò)測試結(jié)果對應(yīng)圓度判定的等級，根據(jù)學習樣本和測試集的神經(jīng)網(wǎng)絡(luò)測試結(jié)果，神經(jīng)網(wǎng)絡(luò)測試結(jié)果的值與邊緣圓度判定等級劃分的對應(yīng)表見表2。

圖6 測試圖片集 （按降序圓度減?。?/p>

最后，將本文設(shè)計的邊緣圓度的評價模型加載到數(shù)字月球平臺，評價模型可以對月表的隕石坑進行邊緣圓度判定，結(jié)果展示如圖7所示。

表2 圓度等級劃分的數(shù)值區(qū)間

圖7 邊緣圓度評價示例

4 結(jié)束語

本文基于圓度誤差的理論和求解方法，提出了一種月表隕石坑邊緣圓度評價算法，通過BP神經(jīng)網(wǎng)絡(luò)綜合了最小二乘圓法、最小區(qū)域圓法、最小外接圓法和最大內(nèi)接圓法4種算法的評價結(jié)果，進而得到邊緣圓度等級。通過邊緣圓度算法建立評價模型，并且對實驗樣本進行了評價，實驗結(jié)果和人工對比結(jié)果基本一致。本文的邊緣圓度評價結(jié)果最終可以和其它5個判定因子一起，建立了全月性的月表的隕石坑幾何影像特征庫，為研究月表地質(zhì)和探月工程的后續(xù)工作提供了一個研究平臺。由于運用的 “嫦娥一號”全月球影像數(shù)據(jù)分辨率為120m，不能精確的對中小型和小型的隕石坑的進行邊緣圓度等級判定，未來的研究工作是利用 “嫦娥二號”7m 分辨率影像數(shù)據(jù)，根據(jù)隕石坑大小的級別，分別建立隕石坑的邊緣圓度等級判定模型。

［1］LI Kun，LIU Jianjun，MOU Lingli，et al.Crater retention age calculation of lunar geological unit using small－scale craters［J］.Journal of Jilin University（Earth Science Edition），2012，42 （2）：452－459 （in Chinese）.［李坤，劉建軍，牟伶俐，等.利用小型撞擊坑測算月球地質(zhì)單元撞擊年齡 ［J］.吉林大學學報地球科學版，2012，42 （2）：452－459.］

［2］Morota T，Haruyama J，miyamoto H，et al.Formation age of the lunar crater Giordano Bruno［J］.Meteorities ＆Planetary Science，2009，44 （8）：1115－1120.

［3］Melosh H J.Planetary surface processes ［M］.Cambridge University Press，2011.

［4］FANG Qinlin.Research on roundness evaluation algorithm and programming ［D］.Wuhan：Huazhong University of Science and Technology，2007 （in Chinese）.［方沁林.圓度誤差評定的算法研究與軟件設(shè)計 ［D］.武漢：華中科技大學，2007.］

［5］FENG Junhua，CUI Hutao，CUI Pingyuan，et al.Autonomous crater detection and matching on planetary surface ［J］.Chinese Journal of Aeronautics，2010，31 （9）：1858－1863 （in Chinese）.［馮軍華，崔祜濤，崔平遠，等.行星表面隕石坑檢測與匹配方法［J］.航空學報，2010，31 （9）：1858－1863.］

［6］ZHANG Yumei.Research on roundness on－machine measurement based on multiscale analysis［D］.Changchun：Jilin University，2009 （in Chinese）.［張玉梅.基于多尺度分析的圓度誤差在線檢測研究 ［D］.長春：吉林大學，2009.］

［7］ZHANG Chunyang，LEI Xianqing，LI Jishun，et al.Method for roundness error evaluation based on geometry optimization［J］.Journal of Mechanical Engineering，2010，46 （12）：8－12（in Chinese）.［張春陽，雷賢卿，李濟順，等.基于幾何優(yōu)化的圓度誤差評定算法 ［J］.機械工程學報，2010，46 （12）：8－12.］

［8］WEN Chao.Research on image－based cultural relic retrieval and ontology－based annotation ［D］.Xi’an：Northwest University，2012 （in Chinese）.［溫超.基于圖像的文物檢索與本體標注方法研究 ［D］.西安：西北大學，2012.］

［9］WAN Dingsheng，HU Yuting，REN Xiang.Applied research of BP network with feedback input［J］.Computer Engineering and Design，2010，31 （2）：398－400 （in Chinese）.［萬定生，胡玉婷，任翔.帶反饋輸入BP 神經(jīng)網(wǎng)絡(luò)的應(yīng)用研究 ［J］.計算機工程與設(shè)計，2010，31 （2）：398－400.］

［10］WANG Qianqian，HUANG Zhiwen，LIU Kai，et al.Classification of plastics with laser－induced breakdown spectroscopy based on principal component analysis and artificial neural network model［J］.Spectroscopy and Spectral Analysis，2012，32 （12）：3179－3182 （in Chinese）.［王茜倩，黃志文，劉凱，等.基于主成分和人工神經(jīng)網(wǎng)絡(luò)的激光誘導(dǎo)擊穿光譜塑料分類識別方法研究 ［J］.光譜學與光譜分析，2012，32 （12）：3179－3182.］