Решение обратной задачи вихретокового контроля

Содержание

Содержание..........................................................................................................................................................................	1
1.Техническоезадание......................................................................................................................................................	2
2.Анализ техническогозадания.......................................................................................................................................	3
2.1 ПрямаязадачаВТК.................................................................................................................................................	3
2.2 ОбратнаязадачаВТК.............................................................................................................................................	3
2.3 Модельзадачи.........................................................................................................................................................	4
2.4 Анализлитературы................................................................................................................................................	4
2.4.1Зарубежныеметодырешения...................................................................................................................	4
2.4.2Отечественныеметодырешения..............................................................................................................	6
3. Прямаязадача ВТКдляНВТП.....................................................................................................................................	9
3.1 УравнениеГельмгольцадля векторногопотенциала...................................................................................	10
3.2 Полевитка надмногослойнойсредой..............................................................................................................	10
3.3ВоздействиепроводящегоОК наНВТП...........................................................................................................	11
4. Обратнаязадача ВТКдляНВТП..................................................................................................................................	13
5.Некорректныезадачи....................................................................................................................................................	16
5.1 Основныеопределения.КорректностьпоАдамару......................................................................................	16
5.2КорректностьпоТихонову...................................................................................................................................	16
5.3Вариационныеметоды решениянекорректныхзадач...................................................................................	17
5.3.1 Методрегуляризации...................................................................................................................................	17
5.3.2 Методквазирешений...................................................................................................................................	17
5.3.3 Методневязки................................................................................................................................................	18
6. Нелинейноепрограммирование.................................................................................................................................	20
6.1 Метод штрафныхфункций..................................................................................................................................	20
6.2Релаксационныеметоды......................................................................................................................................	20
6.2.1 Методусловногоградиента........................................................................................................................	21
6.2.2 Методпроекцииградиента.........................................................................................................................	21
6.2.3 Методслучайногоспуска..........................................................................................................................	21
6.3 МетодмножителейЛагранжа..............................................................................................................................	21
7. Линейноепрограммирование....................................................................................................................................	23
7.1 Алгоритмсимплексногометода........................................................................................................................	23
8.Одномернаяминимизация..........................................................................................................................................	24
8.1 Алгоритмметодов..................................................................................................................................................	24
9. Результатычисленногомоделирования..................................................................................................................	25
9.1Аппроксимациипри численноммоделировании........................................................................................	25
9.2 Моделиреальныхраспределенийэлектропроводности...............................................................................	26
9.3Принципиальнаявозможностьвосстановления.............................................................................................	29
9.4Восстановлениепо зашумленнымданным......................................................................................................	29
9.5Восстановлениес учетом дополнительнойинформации.............................................................................	30
9.6Восстановлениепри различномвозбуждении................................................................................................	30
10.Заключение....................................................................................................................................................................	32
11.Литература.....................................................................................................................................................................	33
Приложение1 - Программнаяреализация...................................................................................................................	35
Приложение2 - Удельнаяэлектропроводностьматериалов....................................................................................	52
Приложение3 - Результатывосстановления................................................................................................................	53
Приложение4 -Abstract....................................................................................................................................................	78

1.Техническоезадание

Разработатьалгоритм решенияобратной задачивихретоковогоконтроля (ВТК).Объектом контроля(ОК) являютсяпроводящиенемагнитныелисты. Объектыконтроля подвергаютсятермообработке(закалка, отпуск)или насыщениювнешних слоевразличнымивеществами,что приводитк изменениюмеханических,а вследствиеэтого и электромагнитныхсвойств материалалиста по глубине.

Задачазаключаетсяв определении,в рамках допустимойпогрешности, зависимостиэлектропроводности(ЭП) от глубиныs(Н) в ОК для данногосостояния.Метод контролязаключаетсяв измеренииопределенногоколичествакомплексныхзначений вносимойЭДС на различныхчастотах спомощью накладноговихретоковогопреобразователя(НВТП).

Необходимовыбрать математическуюмодель задачи,способ аппроксимацииискомого решения,рассмотретьалгоритм решения.

Используяпрограммнуюреализацию,исследоватьповедениепогрешностиаппроксимациизависимостиs(Н)от следующихфакторов:

От величиныприборнойпогрешностиизмерения ЭДС
От видазависимостиэлектропроводности от глубиныs(Н)
От параметроваппроксимациирешения
От диапазоначастот возбужденияВТП

2.Анализ техническогозадания.

Основнаязадача вихретоковогоконтроля спомощью накладныхпреобразователейсостоит из двухподзадач:

Прямойзадачи расчетавносимой ЭДСв присутствиинемагнитногопроводящеголиста с произвольнойзависимостьюЭП по глубине.
Обратнойзадачи нахождениязависимостиЭП как функцииглубины внемагнитномпроводящемлисте по результатамизмеренийопределенногоколичествакомплексныхзначений вносимойЭДС.

2.1 Прямаязадача ВТК

ПолагаязависимостьЭП от глубиныизвестнойпроведем еекусочно-постояннуюаппроксимацию.Это позволяетсвести исходнуюзадачу к расчетуЭДС в многослойномлисте, в каждомслое которогоЭП принимаетпостоянноезначение.

Какпоказано вработе [50],подобная модельвполне адекватноописываетзадачу и даетотличное согласованиес результатамиопытов.

Рекуррентныеформулы дляпроизвольногоколичестваслоев хорошоизвестны [1-5,36,42,43,50-52].Таким образомрешение прямойзадачи в рамкахпринятой моделизатрудненийне вызывает.

2.2 Обратнаязадача ВТК

Сматематическойточки зренияобратная задачаВТК относитсяк классу некорректныхзадач[49]и ее решениенеустойчивот.е. при скольугодно малойпогрешностиисходных данных(набора измеренныхвносимых ЭДС) погрешностьрешения ( рассчитанныхлокальныхзначений ЭП) может бытьсколь угоднобольшой, а одномунабору измеренийможет отвечатьмного (формальнобесконечномного) распределенийЭП по глубине.

Припопытке расчетанекорректнойзадачи каккорректной,вычислительныйпроцесс за счетнеустойчивостисваливаетсяв заведомохудшую сторону.В нашем случаеэто означаетполучениераспределенияЭП, которое,хотя и обеспечиваеттребуемоесовпадениеизмереннойи вычисленнойЭДС, но являетсяявно нереальнымиз-за осцилляций.Следует отметить,что амплитудаи частота осцилляцийраспределенияЭП растут приувеличениичисла независимыхпараметроваппроксимацииЭП ( коэффициентовполинома вслучае полиномиальнойаппроксимации,количестваузлов присплайн-аппроксимациии т.д.).

Приналичии погрешностиизмерениявносимой ЭДС,превышающейна несколькопорядковвычислительнуюпогрешностьи на практикесоставляющейне менее (0.5-1)%от измеряемогосигнала, ситуациязначительноосложняется.

Учитываявышеизложенноедля выделенияиз множествадопустимыхраспределенийрешения, наиболееудовлетворяющегофизическойреальности,в алгоритмахрешения обратнойзадачи необходимоиспользоватьдополнительнуюаприорнуюинформацию.На практикеэто реализуетсявведениемнекоторыхкритериев,позволяющихотличить решение,отвечающеепрактике, отфизическинереального.

Длярешения обратнойзадачи ВТКпредлагалисьтри возможныестратегии[46]:

Решениебольшого числапрямых задачи табуляциярезультатовдля различныхмоделей. Измеренныеданные с помощьюнекоторыхкритериевсравниваютсяс таблицей.Подход оченьэкстенсивныйи требующийпроведенияизбыточногочисла расчетов,поэтому напрактикевстречающийсяредко.
Условнаяминимизацияневязки измеренныхи расчитанныхданных. Оченьмощный и универсальныйметод, широкораспространендля решенияобратных задачв различныхобластях техники[41,44,49].ПозволяетвосстанавливатьпроизвольноераспределениеЭП по глубине(вообще говоряпроизвольное3Dраспределение),но требуетсядовольно сложнаяпроцедурарасчета.
Аналитическоеинвертированиеядра оператораи использованиеалгоритма,зависящегоот ядра уравнения[46].Потенциальносамый малозатратныйметод, однакокак и все аналитические,применим далеконе всегда.

Внашем случаеостановимсяна втором подходе,поскольку онсочетает в себеуниверсальность,точность иотносительнуюпростоту реализации.

Вцелом процессрешения обратнойзадачи сводитсяк итерационномурешению прямойзадачи длятекущей оценкираспределенияЭП и внесениюизменений вэту оценку всоответствиис величинойневязки.

2.3 Модельзадачи

Приведемосновные положения,на основе которыхбудет построенамодель нашейзадачи:

ОКпредставляетиз себя находящуюсяв воздухе проводящуюпластину толщинойН состоящуюиз Nплоско-параллельныхслоев толщинойb_i.
Впределах каждогослоя удельнаяэлектропроводностьsимеет постоянноезначение т.е.распределениеsпо глубинеаппроксимируетсякусочно-постояннойзависимостью.
Возбуждающаяи измерительнаяобмотки ВТПзаменяютсянитевиднымимоделями. Следуетотметить, чтоэто предположениесказываетсялишь на решениипрямой задачи,а проведяинтегрированиеможно получитьвыражения длякатушек конечныхразмеров.
ДлячисленногомоделированияреальныхраспределенийЭП применимпять типоваппроксимации:сплайном,кусочно-постоянную,кусочно-линейную,экспоненциальнуюи гиперболическимтангенсом. Впроцессе решенияпрямой задачис их помощьювычисляютсязначения sв центральныхточках слоевпластины.

2.4 Анализлитературы

2.4.1Зарубежныеметоды решения

Решениюобратной задачиВТК посвященряд работ взарубежныхизданиях. Следуетотметить монографию[38],в которой рассмотреныслучаи импульсноговозбуждения,а оперируютв частотнойи временнойобластяхнапряженностьюэлектрическогополя.

Подходк решениюквазистационарныхзадач рассмотренв цикле статей[45-51].Он основан наинтегральнойпостановкезадачи с помощьюфункций Грина[31-34,39].Для иллюстрациирассмотримрешение обратнойзадачи ВТКсогласно [49].

А.Прямая задача

Определимфункцию v(r)=(s(r)- s₀)/s₀, где s(r)- произвольноераспределениепроводимости,а s₀- ее базоваявеличина. Функцияv(r)можетпредставлятьсобой как описаниепроизвольногораспределенияпроводимости(в этом случаедля удобстваполагаем s(r)=s₀вне некоторогоОК объема V,тогда v(r)отлична от нулятолько в пределахV) так инекоторогодефекта (длятрещины v(r)=-1 внутридефекта и равнанулю вне его).

Рассмотримсистему уравненийМаксвелла впредположениигармоническоговозбужденияexp(-jwt)и пренебрегаятоками смещения:

( 2.4.1)

гдеP(r)=[s(r)-s₀]ЧE(r)=s₀Чv(r)ЧE(r)- может интерпретироватьсякак плотностьдиполей эффективноготока, причинойкоторого являетсявариация s(r)-s₀.

РешениеуравненийМаксвелла можнопредставитьв виде

( 2.4.2)

гдеEⁱ(r)- возбуждающееполе, а G(r|r’)- функция Грина,удовлетворяющаяуравнениюСґСґG(r|r’)+k²ЧG(r|r’)=d(r-r’), k²=-jЧwЧm₀Чs₀, d(r-r’)- трехмернаядельта-функция.

ИмпедансВТП можно выразитькак

( 2.4.3)

гдеинтеграл беретсяпо измерительнойкатушке, J(r)- плотностьтока в возбуждающейкатушке. Применяятеорему взаимностиимпеданс можнопредставитьчерез возбуждающееполе:

( 2.4.4)

гдеинтеграл беретсяпо объему ОК.

В.Обратная задача

Пустьv(r)- оценкаистинной функцииv_true(r),Z_obs(m)- измеренныйимпеданс ВТПв точке r₀на частотевозбужденияw, m=(r₀,w)- вектор в некоторойобласти определенияM, Z[m,v]- оценка величиныZ_obs(m)на основе решенияпрямой задачи.

Определимфункционалневязки измеренныхи рассчитанныхзначений импедансаВТП как :

( 2.4.5)

Предположим,что для решенияобратной задачииспользуетсяитерационныйалгоритм типаметода спуска:vⁿ(r)=v^n-1(r)+asⁿ(r).Можнопоказать, чтов случае методанаискорейшегоспуска итерацияимеет вид:vⁿ(r)=v^n-1(r)-aЧСF[v^n-1(r)], гдеградиент функционалаСF[v] можно определитькак :

( 2.4.6)

гдеReобозначаетвещественнуючасть, * обозначаеткомплекснуюсопряженность.

Требуемыйв (2.4.6) градиентимпеданса можноопределитькак:

СZ(r)= -s₀ЧE(r)ЧE^*(r)

( 2.4.7)

гдеE^*(r)- решение уравнения

( 2.4.8)

С.Аппроксимацияпри решенииобратной задачи

Пустьэлектропроводностьмоделируетсяс помощью конечногочисла переменных(например узловыхзначений некоторойаппроксимации),а вектор рсостоит из этихпеременных.Тогда выражение(2.4.7) принимаетвид:

( 2.4.9)

где(СZ)_j- j-аякомпонентаградиентаимпеданса.

Значениеj-ойкомпонентыградиентаневязки (2.4.6) можнопредставитькак:

( 2.4.10)

Следуетобратить вниманиена то, что в случаедискретногопространстваМ (конечноечисло измерений)интеграл в(2.4.10) заменяетсясуммой.

Сучетом приведенныхпреобразованийитерация методанаискорейшегоспуска принимаетвид:

p_jⁿ=p_j^n-1-aЧ(СF^n-1)_j

( 2.4.11)

гдеn- номер итерации.

D.Пример применения

Вкачестве примерарассмотримфункцию v(r)в виде v(r)=Sc_iЧf_i(r),i=1,N, где f_i(r)- множестволинейно независимыхбазовых функцийс коэффициентамиc_i.Рассматриваякоэффициентыc_iв роли параметроваппроксимации(c_i=p_i) получимиз (2.4.9) для компонентовградиентаимпеданса:

( 2.4.12)

Вслучае проводящегоОК, состоящегоиз Nпараллельныхслоев с проводимостьюs_jраспределениеэлектропроводностипо глубинеможно представитьс помощью функцийХевисайда H(z)как s(z)=Ss_jЧ[H(z-z_j)-H(z-z_j+1)].

Подставляяв (2.4.12) базовыефункции видаf_i(z)=[H(z-z_j)-H(z-z_j+1)],получим окончательноевыражение:

( 2.4.13)

Отметимосновное преимуществотакого решения.Несмотря наопределеннуюсложностьвычисленийпри решенииинтегральныхуравнений(2.4.2-2.4.8) для расчетаградиентаимпеданса НВТПнеобходиморешить толькодве такие задачи.

2.4.2Отечественныеметоды решения

Подход,в значительноймере аналогичныйработам [45-51]был предложенв работе [41].Из-за небольшогообъема в нейуделено недостсточноевнимание вопросампрактическойреализации,объяснены невсе обозначения и не приведенырезультатычисленногомоделирования.В целом этозначительноснижает практическуюценность статьи.Приведем основныеположения этойработы.

Прямаязадача

Пустькруговой витокрадиусом ас током Iнаходится вточке P=P_s(r,j,z),jО(-p,p)вблизинемагнитногоОК, занимающегообласть V.Пусть ОК обладаетэлектрическойпроводимостьюs=s₀Чs(Р)являющейсяпроизвольнойфункцией координат.Требуется поNизмерениямвеличины э.д.с.определитьsкак функциюкоординат точекPОV.Причемi-оеизмерениеэ.д.с. будемпроводить наi-омизмерительномкруговом виткес координатамиP_i=P_i(r,j,z)i=1,Nпри неизменныхчастоте ирасположениивозбуждающеговитка.

Вобщем случаенапряженностьэлектрическогополя Еопределяетсячерез векторныймагнитныйпотенциал А,причем А= А₀+А_вн,где А₀- возбуждающий,а А_вн- вносимыйпотенциалы.

(2.4.14)

Вводяфункцию ГринаG(p,p₀)получим

(2.4.15)

Приэтом вносимаянапряженностьэлектрическогополя

E_вн=-jЧwЧA_вн

(2.4.16)

Вносимаяэ.д.с., наводимаяв i-омвитке

(2.4.17)

гдефункция ГринаG(P,P₀)имеетвид

(2.4.18)

Вдальнейшемрассмотримслучай, прикотором V-полупространство(r>0,\j\p,zс электрическойпроводимостьюs=s₀Чs(Р),гдеs(Р)-произвольнаяфункция координатыz.В этом случаевыражение(2.4.17) примет вид

(2.4.19)

гдеk²=jwm₀s₀.

Длянапряженностиэлектрическогополя Е(Р)справедливопредставление

(2.4.20)

гдеЕ⁰(р)- возбуждающееполе витка.

Послепроведениясерии из Nизмеренийвеличиныe_внвыражение(2.4.19) дает связьмежду вносимымиЭДС e_iи s(z)E(r,z).Чтобы определитьнепосредственноs=s(z),находим E(r,z)при известнойфункции s(z)E(r,z)из (2.4.20),после чегоисключаемE(r,z)из известного.

Обозначимx(p)=-k²s(z)E(r,z),а измеряемуюсовокупностьЭДС через Fi.Тогда (2.4.19)можно записатьв операторнойформе как

F= Px + d

(2.4.21)

гдеd- погрешностьизмерения.

Обратнаязадача

ПостроимфункционалФ(х)=||F-Px||²+a||x-x₀||²,где х₀- некотороеизвестноеІблизкоеІк искомомураспределение,удовлетворяющееF₀=Px₀.Образуемвариацию функционалаФ(х),используяопределениесопряженногооператора(Px,y)=(x,P*y).Для нахожденияминимума Ф(х)приравняемего вариациюdФнулю.

Вводявспомогательнуюфункцию u=x-x₀и учитываяF₀=Px₀проведем рядпреобразований.Искомое распределениеs(z)можно найтииз равенства

(2.4.22)

гденапряженностьэлектрическогополя в точкер дляизвестногораспределенияs(z)имеет вид

	(2.4.23)
	(2.4.24)

Системаалгебраическихуравнений дляопределениякоэффициентовС_iимеетвид

	(2.4.25)
, j=1,N	(2.4.26)

3. Прямаязадача ВТК дляНВТП

3.1 УравнениеГельмгольцадля векторногопотенциала

Взаимодействиепреобразователяс объектомконтроля определяетсясистемой уравненийМаксвелла вдифференциальнойформе[6]:

(3.1)

где

H	- вектор напряженностимагнитногополя
E	-вектор напряженностиэлектрическогополя
B	-вектормагнитнойиндукции
	- векторплотностиполного тока
	- векторплотноститоков проводимости
s	-удельнаяэлектрическаяпроводимость
	- векторплотноститоков смещения
D	- векторэлектрическогосмещения
	- векторплотноститоков переноса
u	- векторскорости переноса
j_стор	- векторплотностистороннихтоков

Дополнимсистему (3.1) уравнениямисвязи:

B= mЧm₀ЧH	(3.2)
B= rot A	(3.3)

где

m₀= 4ЧpЧ10^-7	-магнитнаяпостоянная
m	-относительнаямагнитнаяпроницаемость
A	-векторныйпотенциалмагнитногополя

Преобразуемсистему уравнений(3.1) с учетом следующихпредположений[4]:

ОКнеподвиженотносительноэлектромагнитногополя т.е. jпер= 0
средаизотропна иее параметрыне зависят отнапряженностейполей
воздействиясинусоидальны
последовательностьдифференцированияпо времени ипространственнымкоординатамможно изменять,а операциядифференцированиялинейна

(3.4.1)

(3.4.2)

Посколькуротор градиенталюбого скаляратождественноравен нулю,величину вскобкахвыражения(3.4.2) можно приравнятьградиентунекоторогоскаляра y, например скалярногопотенциалаэлектрическогополя :

(3.5)

Заменяявекторы напряженностимагнитногои электрическогополя в (3.4.1) черезвекторныйпотенциалмагнитногополя получаем:

graddiv A - DA= -mЧm₀Ч( s+ jЧwЧeЧe₀) Ч( grady+ jЧwЧA) + mЧm₀Чj_стор

(3.6)

Откудапосле очевидныхпреобразованийследует:

(3.7)

где

k²= w²Чm Чm₀ЧeЧe₀ - j ЧwЧm Чm₀Чs

(3.8)

Посколькувекторныйпотенциалмагнитногополя задан сточностью доградиентанекоторогоскаляра, а потенциалyс точностьюдо постояннойвеличины, имеетсявозможностьположить значениевеличины вквадратныхскобках выражения(3.7) равное нулю(так называемаякалибровкаЛоренца). Врезультатеполучаем уравнениеГельмгольцадля векторногопотенциаламагнитногополя :

(3.9)

Вдальнейшихрассужденияхиспользуемследующиепредположения:

ПолеНВТП квазистационарнов том смысле,что волновымипроцессамив воздухе можнопренебречь.Это вполнеоправдано т.к.размеры НВТПи ОК обычномного меньшедлины волныв воздухе, апотери на излучениепо сравнениюс потерями вОК малы.
Впроводящемтеле будемрассматриватьтолько волновыепроцессы,обусловленныеналичием параметровsи mт.е. токамисмещения(пропорциональнымиwЧeЧe0) как и в воздухепренебрегаем.Легко показать,что это предположениесправедливоне только дляметаллов, нои для полупроводниковыхматериаловс удельнымсопротивлениемrдо 50[ОмЧсм].В этом случаевыражение(3.8) принимаетвид :

.

3.2 Полевитка надмногослойнойсредой

Введемцилиндрическуюсистему координат( r,j,z).

Пусть:

-ток,протекающийпо нитевиднойвозбуждающейобмотке срадиусом R₁,находящейсяна расстоянииhот N-слойнойсреды
jстор=I Чd(z - h ) Чd(r- R1)

Отметим,что в силу осевойсимметриисистемы

Вцилиндрическойсистеме координатвыражение (3.9)имеет следующийвид :

(3.10)

Применяяк (3.10) преобразованиеФурье-Бесселяс ядром в видефункции Бесселяпервогопорядка имеющеевид :

получаем

(3.11)

Таккак на поверхностяхраздела слоевОК должна сохранятсянепрерывностьтангенциальныхсоставляющихвекторовнапряженностеймагнитногои электрическогополя, дополняемуравнение(3.11) граничнымиусловиями наповерхностяхслоев ОК( граничныеусловия одинаковыдля А и А* ) :

	(3.12)
	(3.13)

Решивуравнение(3.11) с учетом граничныхусловий (3.12-3.13) иприменяя крешению обратноепреобразованиеФурье-Бесселяимеющее вид:

получаем дляполупространстванад ОК

(3.14)

гдеj=j(l, m , s)- функция граничныхусловий.

3.3 ВоздействиепроводящегоОК на НВТП

Длябольшинстваинженерныхрасчетов можноиспользоватьнитевиднуюмодель обмотокНВТП использованнуюв (п 3.2). При данномупрощенииполучаем :

	- напряженностьэлектрическогополя	(3.15)
	-ЭДС наводимаяв измерительнойобмотке срадиусом R2ичислом витковw2	(3.16)

Анализируяформулу (3.14) можнозаметить, чтопервый интегралпредставляетсобой векторныйпотенциалсоздаваемыйвозбуждающейобмоткой, авторой интеграл- векторныйпотенциалвносимый ОК.В практике ВТКобычно анализируютсявносимые параметрыНВТП (напряжение,импеданс) поэтомуполучим выражениедля вычислениявносимогонапряжениекруговоготрансформаторногонакладногоВТП используя(3.15-3.16):

(3.17)

Подставляявыражение длявносимоговекторногопотенциала(3.14) в уравнение(3.17) окончательнополучаем :

(3.18)

где

w=2ЧpЧf	-круговая частотатока возбужденияI
m₀	-магнитнаяпостоянная
w_и, w_в	-числа витковизмерительнойи возбуждающейобмоток НВТП
R= Ц(R_иЧR_в)	-эквивалентныйрадиус НВТП
K_r= Ц(R_в/R_и)	-параметр НВТП
x	-переменнаяинтегрирования
h^*= (h_и+h_в)/2	-обобщенныйзазор
J₁	-функция Бесселя1 рода 1 порядка
j_m	-функция граничныхусловий

Функцияграничныхусловий дляm-слойногоОК с плоскопараллельнымислоями можетбыть вычисленапо рекуррентнойформуле[2]:

(3.19)

где

	(3.20)
	(3.21)
	(3.22)

th(z)	-гиперболическийтангенс
mm	-относительнаямагнитнаяпроницаемостьm-гослоя
bm*= 2Чt_m/ R	-относительнаятолщинаm-го слоя
t_m	-толщина m-гослоя
qm	-обобщенныйпараметр m-гослоя
j₁	-функция граничныхусловийдля нижнегополубесконечногослоя, длявоздуха ( m= 1 , e= 1 , s= 0 ) j₁= 0

Прианализе годографовдля удобстваиспользуютнормированныезависимости.Для НВТП нормировкупроизводятпо модулюмаксимальноговносимогонапряжения,которое соответствуетидеально проводящемуОК и вычисляетсяпри j_м= -1:

(3.23)

Такаянормировкаобобщает полученныерезультаты,расширяетобласть ихпримененияи делает иходнозначными.

Отметим,что для получениячасто используемогов ВТК значенияимпеданса НВТПдостаточноразделитьправую часть(3.18) на величинутока возбужденияI.

4. Обратнаязадача ВТК дляНВТП

Решениеобратной задачиВТК состоитв нахождениизависимостиs(h)распределенияэлектропроводностипо глубинепластины используянабор из Nизмеренныхс помощью НВТПвносимых напряжений.Математическиобратную задачуможно представитьинтегральнымуравнением

(4.1)

Посколькуявного методарешения уравнения (4.1) не существует,применим к немуметод квазирешения(п5.3.2). В постановкедля локальногов смысле Чебышевакритерия получимкорректнуюзадачу минимизациифункционаланевязки:

,i=1,N

(4.2)

Учетаприорнойинформациив обратнойзадачи ВТКудобно проводитьв виде интервала[s_min, s_max],которомумогут принадлежатьзначенияэлектропроводности.В этом случаеможно рассматриватьзадачу (4.2) какзадачу нелинейногопрограммированиявида:

(4.3)

Заметим,что посколькуограниченияв задаче (4.3) являютсялинейными,разумнымпредставляетсяприменениеметода условногоградиента(п6.2.1). Рассмотримпроцесс решениясистемы (4.3) впредположении,что электропроводностьаппроксимируетсяпо узловымзначениям s_j, j=1,M.

(4.4)

ЛинеаризуемфункционалФ вокрестностиисследуемойточки s⁰разложив егов ряд Тейлорас использованиемтолько первыхпроизводных.

(4.5)

Пустьy= maxФ_i’= Ф_p’і0. В этомслучае мы можемсвести задачу(4.4) к эквивалентнойзадаче линейногопрограммирования,состоящей вусловной минимизациифункции y.Рассмотримпроцесс приведениязадачи линейногопрограммированияк каноническомувиду.

Раскрываямодуль в (4.5) получаемсистему уравнений

(4.6)

Рассмотримвыражение подмодулем в (4.5) ивведем некоторыеобозначения

(4.7)

(4.8)

(4.9)

(4.10)

С учетомсистемы (4.8 - 4.10)постановказадачи (4.4) принимаетвид

(4.11)

Раскрываяскобки в (4.11) иисключая yиз первых 2Nнеравенствкроме р-гополучаем системунеравенств

(4.12)

Приведемсистему неравенств(4.12) к каноническомувиду (6.1). Для этого,в соответствиисо стандартнымподходом, запишемвсе неравенствав виде равенств,добавляя влевые частинеравенствнеотрицательныепеременныеv.

(4.13)

Вматричном видеполученнаясистема имеетвид Ax= b (4.14),где искомыйвектор-столбециз 2(N+M)+1элементовимеет вид x= ( y , z₁, ... , z_M, v₁, ... , v_2N+M)^T.В системе линейныхалгебраическихуравнений(4.13) параметрминимизацииyопределенстрокой с номеромp,которую в дальнейшембудем называтьбазовой.

Рассмотрималгоритм симплексногометода длярешения задачи(4.14):

Выборначальногобазиса - допустимогорешения (4.14). Внашем случаебазис долженсостоять из2N+Mпеременных.Удобно задатьначальныйбазис, присвоивдополнительнымпеременнымv_iзначенияправых частейb_iтех строк,в которыхкоэффициентматрицы Aпри нихравен 1.Начальноезначение параметраминимизацииyравнозначению правойчасти базовойстроки. Всеостальныекомпонентыискомого векторах принимаютсяравными нулю.
Определениепеременной,которая должнавойти в очереднойпробный базис.Для этого проводитсяанализ базовойстроки pматрицыA.Из всех положительныхэлементовстроки p,не являющихсякоэффициентамипри базисныхпеременных, выбираетсяэлемент с наибольшимзначением.Переменная,у которой этотэлемент являетсякоэффициентом,должен войтив очереднойпробный базис,т.е. за новуюбазисную переменнуюпринимаетсята, котораяимеет наибольшийвес в функцииy.Если вбазовой строкеpнет небазисныхпеременныхс положительнымикоэффициентами,то в силу неотрицательностиэлементов хследует сделатьвывод, чтооптимальномурешению, т.е.минимуму yсоответствуетвыбранныйранее базис.Вычислениязавершаютсятакже и призапрете измененияпеременныхпо ограничениям.
Определениемаксимальнойдопустимойвеличины новойбазисной переменной,не выходящейза пределыимеющихсяограничений.Вычисляютсяотношениязначений правыхчастей (4.14) ксоответствующимзначениямкоэффициентовпри новой базиснойпеременнойво всех строках,кроме базовой.При этом нерассматриваютсяотношения, вкоторых знаменательравен нулю илиотрицателен,т.е. при положительнойправой частиподобные случаисоответствуютбесконечнымзначениямпеременных.Определяетсяномер строкиq,где это отношениенаименьшее.Новой базиснойпеременнойприсваиваетсязначение отношенияв строке q.Переменная,входившая впрежний базиси определявшаясястрокой q,исключаетсяиз базиса иприравниваетсянулю. Если вовсех строках,кроме базовой,коэффициентыпри новой переменнойравны нулю илиотрицательны,то в силу неотрицательностиэлементовх и ограничениябазиса (2N+M)переменными,следует признать,что эта переменнаяне может наданном шагевычисленийвойти в базис.В этом случаенеобходимовернуться кпункту 2,не рассматриваязапрещеннуюпеременную.
Преобразованиесистемы (4.14) такимобразом, чтобыв строке qкоэффициентпри вновь введенномпараметре былравен 1,а в остальныхстроках - 0.Это достигаетсяпутем линейныхпреобразованийравенств, входящихв (4.14). Т.к. коэффициентыпри параметрах,входящих вновое пробноебазисное решение,становятсяравными 1и в каждую строкувходит толькоодин базисныйпараметр, тозначение новогобазиса определяетсяправой частьюуравнений.Далее следуетвозврат к пункту2.

Решаясистему (4.14) находимвектор s_min, соответствующийтекущему решениюзадачи(4.13). Возвращаяськ методу условногоградиентаотметим, чтонаправлениеспуска определяетсякак -s_n=s_min-s⁰, а очередноеитерационноерешение задачи(4.3) определяетсявыражениемsⁿ⁺¹=sⁿ- aЧs_n.Для полученияокончательногорезультататребуетсяопределитьоптимальнуювеличину шагаaв направлении s_n, что можноосуществитьпутем одномернойминимизациифункции s(a)=sⁿ- aЧs_nметодом золотогосечения.

5.Некорректныезадачи

5.1 Основныепонятия. Корректностьпо Адамару

В самомобщем видебольшинствообратных задачможет бытьпредставленов виде операторногоуравнения

A·x = f , xОX , fОF

( 5.1 )

гдеА - оператор,определенныйна непустоммножественекоторогометрическогопространстваХ сметрикой r_Xи действующийв метрическоепространствоFс метрикой r_F,а по заданномуэлементу fтребуетсяопределитьрешение х[10-14].

Введемв пространствеXнорму||x ||=Цеx_i²и в пространствеFнорму||f ||=Цеf_j².Заметим, чтометрики rв соответствующихпространствахбудут иметьвид r(x,y)=\x-y\.

В нашемслучае обозначенияв (5.1) имеют следующийсмысл:

Ає		-оператор,согласно которомувычисляетсявеличинаотносительногонапряжения,вносимогопластиной сэлектрическойпроводимостьюs(h)
х є	s(h )	-электрическаяпроводимостьпластины какфункция глубины
f є	U^*_вн	-величинаотносительноговносимогонапряженияНВТП

Согласноклассическогоопределениязадача (5.1) называетсякорректнойпо Адамару еслипри любойфиксированнойправой частиее решение:

существуетв Х
единственнов Х
непрерывнозависит от f

Вреальных условияхправая часть(5.1) известна всегдас некоторойпогрешностью,т.е. f= f0+df, причемобычно f0принадлежитпространствугладких функций,а погрешностьdfвыводит ее изэтого класса.Вследствиеэтого получаемпостановкузадачи, длярешения которойневозможноприменениеобычных методоврешения корректныхзадач, т.к. любойфиксированнойправой части(5.1) соответствуетбесконечноемножествонаборов исходныхданных т.е. возможныхраспределенийЭП по глубинепластины.

5.2Корректностьпо Тихонову

Задача(5.1) называетсякорректнойпо Тихоновуна множествекорректностиМ МXесли:

точноерешение задачисуществуети принадлежитМ
принадлежащееМ решениеединственнодля любой правойчасти
принадлежащееМ решениенепрерывноотносительноправой части

Вданном подходек вопросукорректностисуществованиерешения и егопринадлежностьнекоторомумножеству недоказывается,а постулируетсяв самой постановкезадачи.

Физическигипотеза опринадлежностиискомого решенияопределенномумножествукорректностиможет интерпретироватьсядля нашей задачи предположениями:

Исследуемаясреда устроенане слишкомсложно, т.е. еефизическиехарактеристики(s,m)являются достаточногладкими функциями(т.е. их можномоделироватьс помощьюаппроксимацийтипа кусочно-постоянной,кусочно-линейнойи т.п.). Предположениеосновываетсяна физическомсмысле поверхностнойобработки.
Значенияфункций находятсяво вполнеопределенныхпределах( дляs(h)истинностьданного предположениене вызываетсомнения ).

5.3Вариационныеметоды решениянекорректныхзадач

Вариационныеметоды решениянекорректныхзадач являютсянаиболееуниверсальнымииз известныхспособов решения.Практическилюбая некорректнаязадача, длякоторой разработанкакой-либометод решения,может бытьрешена такжеи вариационнымспособом[15].

Длявыбора подходящегометода решенияобратной задачирассмотримпостановкинаиболеераспространенныхвариационныхметодов в терминахвычислительнойматематикии нашей задачи.

Пустьфиксированныйнабор данныхсостоит изизмеренныхна Nчастотах Nкомплексныхзначений вносимойЭДСUiизм,текущее рассчитанноезначение которыхUi(s). Требуетсяопределитьдля выбранноготипа аппроксимацииЭП значенияМ параметроваппроксимации( обычно используютсяузловые значения).

5.3.1 Методрегуляризации

Методоснован настабилизацииневязки r(Ax,f)при помощивспомогательногонеотрицательногофункционалаW(x).Идея методасостоит в том,чтобы минимизироватьобладающийсглаживающимисвойставамифункционалФ(x,f),имеющий следующийвид:

, параметррегуляризацииa>0

(5.2)

Используяклассическийрегуляризирующийфункционалвида

втерминах нашейзадачи получаем:

(5.3)

Основноепреимуществометода состоитв регуляризациипростейшимспособом, врамкахиспользованияквадратичногофункционала.Это позволяетиспользоватьдля решениянекорректнойзадачи хорошоизвестные илегко программируемыеметоды минимизацииквадратичныхфункционалов[17].

Оборотнойстороной достоинствметода являютсяего недостатки.Требованиеминимизациинормы решенияи, как следствие,выбор гладкойреализации,в нашем случаебудет приходитьв противоречиес физикой задачии в принципене позволитнаходить решенияс выраженнымприповерхностнымизменениемЭП.

Ещеодин принципиальныйнедостатокметода состоитв постановкефункционалакак квадратичного,единого длявсех измерений.Его минимумв общем случаене гарантируетминимизациюотклонениядля произвольногоi-гоизмерения вследствиинелокальностиусловия минимизации.

Крометого, следуетучитыватьотсутствиенадежных априорныхрекомендацийпо выбору параметрарегуляризацииa.Обычно подходящиезначения aможноподобратьтолько послеряда численныхэкспериментовпо решениюоднотипныхзадач. Изменениехарактераискомого решенияприводит кнеобходимостипоиска новогозначения a.

5.3.2 Методквазирешений

Методиспользуетодну из формкритерия невязкии заключаетсяв сведенииневязки к минимумуна некоторомнепустом множествеP,содержащемподмножествоискомых решений.

Квазирешениемуравнения (5.1)на множествеPМXназываетсявсякий элементyОPдля которогосправедливоравенство r_F(Ay , f ) = inf( Ax , f ), xОP.Понятиеквазирешенияобобщает понятиерешения, а дляего существованияне требуетсяпринадлежностьрешения множествуP.

Исходяиз вышеизложенногополучаем постановкуметода в видезадачи условнойминимизациифункционалаФ(x,f):

(5.4)

Отметим,что множествоР можетиметь простойвид, напримеринтервала [x_min, x_max].

Втерминах нашейзадачи ВТКпостановказадачи (5.4) приметвид:

(5.5)

Длятого, чтобыгарантироватьминимизациюотклонениядля произвольногоi-гоизмерения,можно применитьк первому выражениюв (5.5) локальныйв смысле Чебышевакритерий, всоответствиис которым получаемокончательноевыражение :

(5.6)

Основноепреимуществометода состоитв том, что самопонятие квазирешенияснимает трудностис требованиямитихоновскойкорректности:первым (вызывающимпереопределенностьзадачи) и третьим(обычно принадлежностьприближеннойправой частиуравнения (5.1)множеству N=AMнеизвестна,а критерии этойпринадлежностичасто самибывают неустойчивы).

Кромеэтого прирассмотрениизадачи в виде(5.6) возможна постановкаминимизационнойзадачи какзадачи нелинейногопрограммированияс явно заданнымиограничениямина искомыепеременные.В этом случаенет необходимостиискажать исходныйфункционалрегуляризующимичленами какв (п5.3.1), а требованияк искомомурешению можноудовлетворить,управляяограничениямина параметрыминимизации(в нашем случае- узловые значенияЭП).

5.3.3 Методневязки

Рассмотриммножество Рформальныхрешений уравнения(5.1) Р={x: r_F(Ax , f_d) Јd},где f_d- приближеннаяправая часть(5.1), известнаяс погрешностьюd.

Вкачествеприближенногорешения (5.1) нельзябрать произвольныйэлемент множестваР, т.к.не гарантируетсяблизость Рк множествуточных решений.Для выбораприближенногорешения предлагаетсяиспользоватьстабилизирующийфункционалW(х)из (п 5.3.1) следующимобразом: W(х)=inf W(х),xОP.Этот способприводит квыбору элементовмножества Римеющих минимальнодопустимуюневязку.

Сучетом этогопостановкаметода состоитв условнойминимизациифункционалаФ(х):

(5.7)

Каки для методарегуляризацииможно использоватьстабилизирующийфункционалвида W(х)=||x||², что приводитв обозначенияхнашей задачик системе:

(5.8)

Прииспользованиилокальногов смысле Чебышевакритерия система(5.8) окончательнопримет вид:

(5.9)

6. Нелинейноепрограммирование

Содержаниенелинейногопрограммированиясоставляюттеория и методырешения задачо нахожденииэкстремумовнелинейныхфункций намножествах,определяемыхлинейными инелинейнымиограничениями(равенствамии неравенствами)[16-29].

Рассмотримнаиболеераспространенныеметоды решенияна примереосновной задачинелинейногопрограммированиявида:

(6.1)

6.1 Методштрафных функций

Идеяметода состоитв замене экстремальнойзадачи с ограничениями(6.1) на задачубезусловнойминимизацииоднопараметрическойфункции

,b>0

(6.2)

Непрерывнуюфункцию y(х)называют штрафом,если y(х)=0для хОХ и y(х)>0в противномслучае. Функцияy(х)должна бытьвыбрана такимобразом, чтобырешение задачи(6.2) сходилосьпри b®0к решению исходнойзадачи (6.1) или,по крайнеймере, стремилоськ нему.

Приведемчасто используемыевыражения дляштрафа :

, k>0	(6.3)
	(6.4)
	(6.5)

Наибольшееприменениенаходит штраф(6.3). Выражение(6.5) гарантируетконечностьметода прилюбом k>0.

Причисленнойреализацииметода штрафныхфункций возникаютпроблемы выбораначальногозначения параметраbи способа егоизменения.Сложностьсостоит в том,что выбор достаточномалого b увеличиваетвероятностьсходимостирешения (6.2) крешению (6.1), аскорость сходимостиградиентныхметодов вычисленияточек минимума(6.2), как правило,падает с убываниемвеличины b.

6.2Релаксационныеметоды

Релаксационнымметодом называютпроцесс построенияпоследовательноститочек {х_k:х_kОX , j(х_k+1) Јj(х_k) ; k=0,1... }. Основнымипредставителямиэтого классаявляются методыспуска, алгоритмкоторых состоитиз следующихшагов :

Выборначальногоприближениях₀
Выборв точке х_kнаправленияспуска -s_k
Нахождениеочередногоприближениях_k₊₁= х_k- a_kЧs_k ,где длина шагаa_k>0

Различияметодов состоятв выборе либонаправленияспуска, либоспособа движениявдоль выбранногонаправления.В последнемслучае обычноиспользуютодномернуюминимизациюфункции х_k₊₁(a)= х_k- aЧs_k(приэтом точностьвычисленияточки минимумафункции х_k₊₁(a)следует согласовыватьс точностьювычислениязначений функцииj(х))или способудвоения a(величинашага удваиваетсяпока выполняетсяусловие j(х_k+1)Јj(х_k)).

6.2.1 Методусловногоградиента

Идеяметода заключаетсяв линеаризациинелинейнойфункции j(х).В этом методевыбор направленияспуска осуществляетсяследующимобразом :

Линеаризируяфункцию j(х)в точке х_Кполучаем Ф(х)=j(х_k) + ( j(х_k)’, х- х_k)
Минимизируялинейную функциюФ(х) на множествеХ находим х_min
Направлениеспуска получаемкак -s_k= х_min- х_k

Такимобразом итерацияметода имеетвид: x_k+1=x_k+a_kЧ(s_k+1-x_k), s_k+1=argmin(Сf(x_k),x).

Основноепреимуществометода проявляетсяв случае заданиядопустимогомножества спомощью линейныхограничений.В этом случаеполучаем задачулинейногопрограммирования,решаемую стандартнымиметодами(напримерсимплексным).

6.2.2 Методпроекции градиента

Этотметод являетсяаналогом методаградиентногоспуска, используемогов задачах безограничений.Его идея состоитв проектированииточек, найденныхметодом наискорейшегоспуска, на допустимоемножество,определяемоеограничениями.Проекцией точкиyна множествоХ называетсяточка P(y)ОХтакая, что ||P(y) - y || Ј|| x - y || длявсех хОХ.Задача проектированияформализуетсякак|| x - y ||²®min, xОХ.

Выборнаправленияспуска осуществляетсяследующимобразом :

Находимточку r_k=х_k- aЧj’(х_k)
Находимпроекцию p_kточкиr_kна множествоХ
Направлениеспуска получаемкак -s_k= p_k- х_k

Такимобразом итерацияметода имеетвид:x_k+1=P_X[x_k-a_kЧСf(x_k)], где Р_X(у)- ортогональнаяпроекция точкиу намножество Х.

Дляотысканиянаправленияспуска s_kнеобходиморешить задачуминимизацииквадратичнойфункции ||r_k- х ||²на множествеХ. В общемслучае этазадача тогоже порядкасложности, чтои исходная,однако длязадач, допустимоемножествокоторых имеетпростую геометрическуюструктуру,отысканиепроекции значительноупрщается.Например, длямногомерногопараллелепипидаQ_N={xОR_N: a Јx Јb }, отысканиепроекции осуществляетсяпутем сравненияnчисели имеет видP(x)={a_i,x_ii; x_i,x_iО[a_i,b_i] ; b_i,x_i>b_i}.

6.2.3 Методслучайногоспуска

Методхарактеризуетсятем, что в качественаправленияспуска s_Kвыбираетсянекотораяреализацияn-мернойслучайнойвеличины Sс известнымзаконом распределения.Об эффективностиэтого методасудить трудно,однако благодаряиспользованиюбыстродействующихЭВМ он оказываетсяпрактическиполезным.

6.3 МетодмножителейЛагранжа

Идеяметода состоитв отысканииседловой точкифункции Лагранжазадачи (6.1). Длянахождениярешения вводитсянабор переменныхl_i, называемыхмножителиЛагранжа, исоставляетсяфункция Лагранжа,имеющая вид:

(6.6)

Алгоритмметода состоитв следующем:

Составлениефункции Лагранжа
Нахождениечастных производныхфункции Лагранжа

(6.7)

Решениесистемы из n+mуравненийвида

(6.8)

Решениямисистемы (6.8) являютсяточки, которыемогут бытьрешениямизадачи.

Выборточек, в которыхдостигаетсяэкстремум ивычислениефункции j(х)в этих точках.

7. Линейноепрограммирование

Задачалинейногопрограммированияв каноническомвиде имеетвид[15,16]:

(7.1)

Приведениек каноническомувиду любойзадачи линейногопрограммированияосуществляетсяпутем введениядополнительныхнеотрицательныхпеременных,за счет чегоограничения,имеющие виднеравенств,принимают видэквивалентныхим равенств.

Любаязадача линейногопрограммированияможет бытьрешена за конечноечисло итерацийс помощьюсимплексногометода[17,18].Следует отметить,что посколькуэтот методразработандля неотрицательныхэлементов x_j, это условиеучитываетсянеявно и в системууравнений (7.1)при численнойреализациине входит.

7.1 Алгоритмсимплексногометода

1. Приведениек каноническомувиду

2. Выборначальногобазиса

3. Проверкаоптимальностибазиса

МатрицуА можнорассматриватькак совокупностьстолбцов a^jт.е. еa^jЧx_j=bгде j=1,N.Не ограничиваяобщности можносчитать, чтобазис образуютпервые mстолбцов, тогдаостальные можнопредставитьв виде a^k=еa^jЧl_jk , j=1,mгдеl_jk.-некоторыечисла.

РассмотримкоэффициентыD_k=еc^jЧl_jk- c_kгдеj=1,mи k=1,N. Заметим,что для базовыхстолбцов D_kє0. Проверкана оптимальностьосуществляетсяследующимобразом:

D_kk=1,N	- текущийбазис оптимален
	- решениене ограниченосверху
	-существуетдругой, болееподходящийбазис

4. Составлениенового базиса

4.1 Выборэлемента длявведения вбазис.

Вбазис вводитсялюбой столбец,для которогоD_k0, обозначимего D_p

4.2 Выборэлемента исключаемогоиз базиса

Изтекущего базисаисключаетсястолбец, длякоторого минимальноотношениеb_i/A_ip,i=1,Mобозначим егоb_r/A_rp

Преобразованиевектора bи матрицыА пометоду Жордана-Гаусса

4.4Переход к пункту3

8. Одномернаяминимизация

Несмотряна кажущуюсяпростоту, дляширокого классафункций решениезадачи минимизацияфункции одногопеременногоj(х)сопряженос некоторымитрудностями.С одной стороны,в практическихзадачах частонеизвестно,является лифункция дифференцируемой.С другой стороны,задача решенияуравненияjў(х)=0может на практикеоказатьсявесьма сложной.Ввиду этогосущественноезначение приобретаютметоды минимизации,не требующиевычисленияпроизводной[15].

Посколькунас интересуетприближенноеопределениеточки минимума,то для этогоисследуютповедениефункции в конечномчисле точек,способамивыбора которыхразличаютсяметоды одномернойминимизации.

Кметодам, в которыхпри ограниченияхна количествовычисленийзначений j(х)достигаетсяв определенномсмысле наилучшаяточность, относятсяметоды Фибоначчии золотогосечения[17,18].Оба методастроятся поединой схемеи различаютсяспособом выборапараметра l.Исходнымиданными дляних являются:отрезок [a₀,b₀]содержащийточку минимумаи точностьопределенияточки минимумаe.

8.1 Алгоритмметодов

h₀= b₀- a₀, k = 1 , lО(0.5,1) , h₁= lЧh₀, h₂= h₀- h₁, c₁= a₀+ h₂, d₂= b₀- h₂
Вычислитьтекущие значенияj(c_k)и j(d_k)и действоватьв соответствиис ними:

	j(c_k) Јj(d_k)	j(c_k) >j(d_k)
a_k=	a_k-1	c_k-1
b_k=	d_k-1	b_k-1
d_k=	c_k-1
c_k=		d_k-1
h_k+2=	h_k-h_k-1	h_k-h_k-1
d_k=		b_k-h_k+2
c_k=	a_k+h_k+2

Если(h_kЈe) то x_min=min{j(c_k), j(d_k)} иначеk++и переходк шагу II

Следуетотметить, чтона каждом шагекроме первого,производитсятолько одновычислениезначения функцииj(x).

Легкопоказать, чтодля полученияоптимальнойпоследовательностиотрезков,стягивающихсяк точке минимума,необходимоположить l_k=F_k-1/F_k, где F- числоФибоначчи.

8.2 МетодФибоначчи

Решаявопрос, прикаких значенияхпараметра lза конечноечисло итерацийNмы получимотрезок минимальнойдлины, получимl=l_N=F_N-1/F_N.Иначе говоря,для поискаминимумапервоначальнонеобходимонайти числоФибоначчи Nтакое,что F_N+1(b₀-a₀)/eЈF_N+2.

8.3 Методзолотого сечения

Вреальной ситуацииначиная поискминимума мыне знаем точногочисла требуемыхитераций. Вместовычисленияlбудем выдерживатьпостоянноеотношение длининтерваловh_k-2/h_k-1= h_k-1/h_k= t.При t= (Ц5+1)/2= 1.618034 получаемметод золотогосечения.

Сравниваяприведенныеметоды прибольших значенияхNможно показать,что значениеокончательногоинтерваланеопределенностив методе золотогосечения лишьна 17% больше чемв методе Фибоначчи.

9.Результатычисленногомоделирования

9.1Аппроксимациипри численноммоделировании

Дляпостроениямоделей реальныхраспределенийЭП возможноприменениецелого рядааппроксимаций.Все они могутбыть разделенына два класса.

1.Аппроксимации,строящиесяпо набору изпроизвольногочисла узлов.Наиболеераспространенныеиз них: кусочно-постоянная,кусочно-линейнаяи сплайном. Вусловиях нашейзадачи указанныеаппроксимацииимеют несколькосущественныхнедостатков:

Результатыаппроксимацийслабо согласуютсяс реальностью.Кусочно-постояннаяи кусочно-линейнаяаппроксимациипринципиальноявляются негладкими,а аппроксимациясплайном сглаживаетвсе, в результатечего возникаютзначительныенемонотонностии всплески.
Приувеличенииколичестваузлов аппроксимациибыстро нарастаетнеустойчивостьпроцесса решенияобратной задачи,для противодействиякоторой требуетсяприменениеискусственныхприемов, негарантирующихуспеха.
Вреальных условияхмы не имеемдостовернойаприорнойинформациио величине ЭПв узлах аппроксимации,расположенныхв глубине пластины.

2.Аппроксимации,строящиесяпо значениямЭП на верхнейи нижней поверхностяхпластины инесколькимпараметрамаппроксимации.Наиболее известныеиз них: экспоненциальная,гиперболическимтангенсом игауссоидой.Аппроксимацииимеют вид:

	-аппроксимацияэкспоненциальная
	-аппроксимациягиперболическимтангенсом
	-аппроксимациягауссоидой

где

x	-координата,равна нулюна нижнейповерхностипластины иединице наверхней
s₁	- величинаэлектропроводностина верхнейповерхностипластины
s₂	- величинаэлектропроводностина нижнейповерхностипластины
a	-коэффициент,характеризующийкрутизнуэкспоненты
b	-коэффициент
g	-коэффициент,характеризующийкрутизну; g=0соответствуетслучаю слояс проводимостьюs₁и толщиной bна полупространствес проводимостьюs₂
d	-коэффициент,характеризующийкрутизну

Длянашей задачиподобныеаппроксимацииявляютсяпредпочтительными,посколькуобладают заметнымидостоинствами:

Аппроксимацииявляются монотоннымии гладкими,что хорошосогласуетсяс физическойреальностью.
Пользуясьфизическиобоснованнымирассуждениямимы можем получитьнеобходимуюаприорнуюинформациюо величинахЭП в приповерхностныхслоях пластины.
Процессрешения обратнойзадачи существенноболее устойчиви осуществляетсязначительнобыстрее

Дляиллюстрациинаших рассужденийприведем примерпримененияприведенныхвыше аппроксимацийк случаю восстановлениякусочно-линейнойфункции. Пооси абсциссотложенаотносительнаяглубина, пооси ординатэлектропроводность(МСм/м).

Награфике показаныаппроксимации:кусочнопостоянная(SIci),кусочнолинейная(SIli),сплайн(SIs),экспоненциальная(SIe),гиперболическимтангенсом(Sith),гауссоидой(SIg).

Легкозаметить, чтоаппроксимациягиперболическимтангенсомхорошо описываетприповерхностныеизменения(аналогичноэкспоненциальнойпри большомпоказателеэкспоненты).Гауссоидаможет бытьлегко воспроизведенас помощьюэкспоненциальнойаппроксимации,поэтому вдальнейшемиспользованане будет.

9.2 Моделиреальныхраспределенийэлектропроводности

Модельзадачи должнаописыватьнекоторуюпластину,подвергнутуюповерхностнойобработке. Дляопределенностизададим толщинупластины равнойдвум сантиметрам.На основе данныхиз Приложения2зададим значенияЭП вблизи нижнейи верхнейповерхностейсоответственно20(МСм/м)и 13(МСм/м).

Длярешения обратнойзадачи необходимозадать априорнуюинформациюо величине ЭПв узлах аппроксимации.В качестветаковой примеминтервал [8,25](МСм/м),полученныйвнесением 25%отклоненияот считаемыхистиннымизначений. Этоотклонениемоделируетнеточностьаприорнойинформации.

Из-заособенностейреализацииалгоритмаустойчивостьрешения сильнозависит отточности заданияЭП в узле, соответствующемнижней поверхностипластины, поэтомуограничение в нем зададиминтервалом[19,21](МСм/м).

В нашемслучае всевозможныемодели распределенийЭП могут бытьразделены надва класса.Распределенияотносящиесяк первому изних условноназовем глубинными.В них ЭП претерпеваетсущественныеизменения напротяжениивсей глубиныпластины. Второйкласс образуютраспределения,ЭП в которыхзаметно изменяетсялишь в приповерхностномслое глубинойпорядка четвертипластины., поэтомуназовем этираспределенияповерхностными.

Критериемотличия восстановленнойфункции распределенияЭП от модельнойбудем считатьвеличинуотносительнойпогрешности,посколькусравнениерезультатовс ее помощьювполне адекватноцелям нашейработы.

Следуетотметит, чтопогрешностьвосстановлениядля поверхностныхраспределенийЭП представляетпрактическийинтерес в области,примерно ограниченнойглубиной порядкачетверти пластины,что обусловленофизическимсмыслом поверхностнойобработки.Поэтому дляслучаев поверхностныхраспределенийосновное вниманиебудем уделятименно указаннымглубинам.

Дляпроверки возможностивосстановленияприповерхностныхизменений ЭПрассмотримдве базовыемодели поверхностныхраспределений.

Базоваямодель A1.

Аппроксимацияэкспонентой.

Проводимостьs₂=20[МСм/м]

Проводимостьs₁=13[МСм/м]

Показательэкспонентыa={25,38,120,200 }.

БазоваямодельA2

Аппроксимациягиперболическимтангенсом.

Проводимостьs₂=20[МСм/м]

Проводимостьs₁=6 [МСм/м]

Коэффициентb=1

Коэффициентg= {0.1, 0.05, 0.02, 0.01 }

Дляпроверки возможностивосстановленияглубинныхраспределенийЭП рассмотримдве базовыемодели глубинныхраспределений.

БазоваямодельB1

Аппроксимациякусочно-линейная.

Проводимостьзадаетсяв узлах с отсчитываемойот дна пластиныотносительнойглубиной {0,0.25,0.5,0.75,1}.

Узловыезначенияпроводимости{20,20,17.6,15.3,13},{20,20,20,16.5,13},{20,20,20,20,13} [МСм/м].

БазоваямодельB2

Аппроксимациясплайном.

Проводимостьзадаетсяв узлах с отсчитываемойот дна пластиныотносительнойглубиной {0,0.25,0.5,0.75,1}.

Узловыезначенияпроводимости{20,20,17.6,15.3,13},{20,20,20,16.5,13},{20,20,20,20,13} [МСм/м].

Заметим,что на практикеможно осуществитьдостаточноточное определениевеличины ЭПприповерхностныхслоев путемизмеренийпроводимоститрадиционнымисредствами,поэтому дополнительнорассмотриммодельныезадачи приусловии известнойЭП на верхней,а так же верхнейи нижней поверхностях.

Посколькуна практикерезультатыизмеренийвносимогонапряженияимеют определеннуюпогрешность,все моделибудем рассчитыватьэмулируя погрешностьdU=0,1,2,5%.

ДляисследованиязависимостирезультатоввосстановленияраспределенийЭП от частотывозбужденияразобьем частотныйдиапазона тричасти следующимобразом( глубиныпроникновенияприведены дляслучая постояннойЭП s=13МСм/м ):

Модели	FA1L,FB1L	Модели	FA1M,FB1M	Модели	FA1H,FB1H
f, [Гц]	h , [m]	f, [КГц]	h , [m]	f, [КГц]	h , [m]
1	0.1396	5	0.001974	55	0.0005952
10	0.04414	10	0.001396	60	0.0005699
20	0.03121	15	0.00114	80	0.0004935
50	0.01974	20	0.000987	90	0.0004653
100	0.01396	25	0.0008828	100	0.0004414
200	0.00987	30	0.0008059	200	0.0003121
500	0.006243	35	0.0007461	300	0.0002549
1000	0.004414	40	0.0006979	500	0.0001974
2000	0.003121	45	0.000658	700	0.0001668
5000	0.001974	54.1	0.0006001	1000	0.0001396

ДляисследованиязависимостирезультатоввосстановленияраспределенийЭП от числаизмеряемыхвносимых напряженийNрассмотримслучаи N={5, 10,15 }.

Низкиечастоты

f, [Гц]	1,5,10,20,35,50,100,150,200,500,750,1000,2000,3500,5000
f, [Гц]	1,10,20,50,100,200,500,1000,2000,5000
f, [Гц]	1,20,100, 500,2000

Средниечастоты

f, [КГц]	5,7.5, 10, 15,17.5, 20, 25, 27.5, 30, 35, 37.5, 40, 45, 50, 54.1
f, [КГц]	5,10, 15, 20, 25, 30, 35, 40, 45, 54.1
f, [КГц]	5,15, 25, 35, 45

Высокиечастоты

f, [КГц]	55,57.5, 60, 80, 85, 90,100, 150, 200, 300, 400, 500, 700, 850, 1000
f, [КГц]	55,60, 80, 90,100, 200, 300, 500, 700, 1000
f, [КГц]	55,80, 100, 300, 700

9.3Принципиальнаявозможностьвосстановления

ДляисследованиявозможностивосстановленияраспределенияЭП рассмотримрезультаты,полученныев предположенииналичия точныхданных (погрешностьизмеренияотсутствует).На графикахв первых четырехпунктах Приложения3 рассматриваемыезависимостипоказаны краснымцветом (исходныеданные черным).Исходя из нихможно сделатьследующиевыводы

Восстановлениес помощьюаппроксимации,использованнойпри эмуляцииизмерений(решении прямойзадачи), приводитк погрешностивосстановленияпорядка 0.1%.
Восстановлениеглубинныхраспределенийс помощьюаппроксимацийэкспоненциальнойи гиперболическимтангенсомвозможно схорошей точностью( погрешность2-5% )для приповерхностныхслоев глубинойпорядка четвертипластины.
Восстановлениеглубинныхраспределенийс помощьюаппроксимацийсплайном икусочно-линейнойвозможно схорошей точностью( погрешность2-3% ).Погрешностьвосстановленияувеличиваетсяс уменьшениемглубины.
Восстановлениеповерхностныхраспределенийс помощьюаппроксимацийсплайном икусочно-линейнойпрактическиневозможно.Имеют местоосцилляции,приводящиек погрешностям,превышающим10%.
Восстановлениеповерхностныхраспределенийс помощьюаппроксимацийэкспоненциальнойи гиперболическимтангенсомвозможно схорошей точностью(погрешность2-3%).Погрешностьвосстановленияувеличиваетсяс уменьшениемглубины, занимаемойраспределением.

9.4Восстановлениепо зашумленнымданным

Рассмотренныев данном разделерезультатыдемонстрируютвозможностьвосстановленияраспределенийЭП в реальныхусловиях. Графикипредставленыв первых четырехпунктах Приложения3.

Награфикахрассматриваемыезависимостипоказаны цветами:результатвосстановленияпри погрешностиданных равной1% - голубым,результатвосстановленияпри погрешностиданных равной2% - коричневым,результатвосстановленияпри погрешностиданных равной5% - синим.

Исходяиз них можносделать следующиевыводы:
Восстановлениеглубинныхраспределенийс помощьюаппроксимацийэкспоненциальнойи гиперболическимтангенсомвозможно схорошей точностью( погрешность2-8% )для приповерхностныхслоев глубинойпорядка четвертипластины.
Восстановлениеглубинныхраспределенийс помощьюаппроксимацийсплайном икусочно-линейнойзатруднено(погрешностьосциллируетв пределах0-10% ).Погрешностьвосстановленияувеличиваетсяс уменьшениемглубины, занимаемойраспределением.
Восстановлениеповерхностныхраспределенийс помощьюаппроксимацийсплайном икусочно-линейнойпрактическиневозможно.Имеют местоосцилляции,приводящиек погрешностям,превышающим10%.
Восстановлениеповерхностныхраспределенийс помощьюаппроксимацийэкспоненциальнойи гиперболическимтангенсомвозможно схорошей точностью(погрешность3-6%для одноименныхаппроксимацийи 7-10%в противномслучае).Погрешностьвосстановленияувеличиваетсяс уменьшениемглубины, занимаемойраспределением.

9.5Восстановлениес учетом дополнительнойинформации

С цельюулучшениярезультатоввосстановленияв реальнойобстановке,осложненнойналичием зашумленныхданных, следуетиспользоватьболее жесткиеограниченияна величинуЭП в приповерхностныхслоях.

Дляиллюстрацииприведем несколькографиков,представленныхв пятом пунктеПриложения3. Награфикахрассматриваемыезависимостипоказаны цветами:базовые ограничения- коричневым,жесткие ограниченияна верхнейповерхности- голубым, жесткиеограниченияна обоих поверхностях- малиновым.

Исходяиз полученныхрезультатовможно сделатьследующиевыводы
Восстановлениеглубинныхраспределенийс помощьюаппроксимацийэкспоненциальнойи гиперболическимтангенсомвозможно судовлетворительнойточностью дляприповерхностныхслоев глубинойпорядка четвертипластины.Дополнительныежесткие ограничениярезультатывосстановленияне улучшают.
Восстановлениеглубинныхраспределенийс помощьюаппроксимацийсплайном икусочно-линейнойзатруднено.Дополнительныежесткие ограничениярезультатывосстановленияне улучшают.
Восстановлениеповерхностныхраспределенийс помощьюаппроксимацийсплайном икусочно-линейнойпрактическиневозможно.Имеют местоосцилляции,приводящиек погрешностям,превышающим10%.
Восстановлениеповерхностныхраспределенийс помощьюаппроксимацийэкспоненциальнойи гиперболическимтангенсомвозможно судовлетворительнойточностью(погрешность6-10% ).Погрешностьвосстановленияуменьшаетсяпри использованиидополнительныеограниченийпримерно вдвое,особенно вприповерхностномслое толщинойпорядка 10%.

9.6 Восстановлениепри различномвозбуждении

Длявыбора необходимогоколичестваизмерений U_вн^*и соответствующихим частотвозбужденияВТП рассмотримтри возможныхдиапазоначастот, в каждомиз которыхисследуемслучаи пяти,десяти и пятнадцатичастот.

Награфиках вшестом пунктеПриложения3рассматриваемыезависимостипоказаны цветами:5 частот- коричневым,10частот- голубым , 15частот- малиновым.

Областьнизких частот

Исходяиз полученныхрезультатовможно сделатьследующиевыводы
Восстановлениеглубинныхраспределенийс помощьюаппроксимацийэкспоненциальнойи гиперболическимтангенсомвозможно судовлетворительнойточностью дляприповерхностныхслоев глубинойпорядка четвертипластины. Дляулучшениярезультатоввосстановленияследует использовать10частот в случаепогрешностиданных не более2% и15частот в противномслучае.
Восстановлениеглубинныхраспределенийс помощьюаппроксимацийсплайном икусочно-линейнойзатруднено.Для улучшениярезультатоввосстановленияв приповерхностномслоев глубинойпорядка четвертипластины следуетиспользовать10частот в случаепогрешностиданных не более2% и15частот в противномслучае.
Восстановлениеповерхностныхраспределенийс помощьюаппроксимацийсплайном икусочно-линейнойпрактическиневозможно.Имеют местоосцилляции,приводящиек погрешностям,превышающим10%.
Восстановлениеповерхностныхраспределенийс помощьюаппроксимацийэкспоненциальнойи гиперболическимтангенсомвозможно судовлетворительнойточностью(погрешность6-8% ).Для улучшениярезультатоввосстановленияследует использовать10частот в случаепогрешностиданных не более2% и15частот в противномслучае.

Областьсредних частот

Исходяиз полученныхрезультатовможно сделатьследующиевыводы:
Восстановлениеглубинныхраспределенийс помощьюаппроксимацийэкспоненциальнойи гиперболическимтангенсомвозможно судовлетворительнойточностью дляприповерхностныхслоев глубинойпорядка четвертипластины. Дляулучшениярезультатоввосстановленияследует использовать10частот в случаепогрешностиданных не более2% и15частот в противномслучае.
Восстановлениеглубинныхраспределенийс помощьюаппроксимацийсплайном икусочно-линейнойпрактическиневозможно.Имеют местоосцилляции,приводящиек погрешностям,превышающим10%.
Восстановлениеповерхностныхраспределенийс помощьюаппроксимацийсплайном икусочно-линейнойпрактическиневозможно.Имеют местоосцилляции,приводящиек погрешностям,превышающим10%.
Восстановлениеповерхностныхраспределенийс помощьюаппроксимацийэкспоненциальнойи гиперболическимтангенсомвозможно судовлетворительнойточностью(погрешность6-8% ).Для улучшениярезультатоввосстановленияследует использовать10частот в случаепогрешностиданных не более2% и15частот в противномслучае.

Областьвысоких частот

Исходяиз полученныхрезультатовможно сделатьследующиевыводы:
Восстановлениеглубинныхраспределенийс помощьюаппроксимацийэкспоненциальнойи гиперболическимтангенсомвозможно судовлетворительнойточностью дляприповерхностныхслоев глубинойпорядка четвертипластины. Дляулучшениярезультатоввосстановленияследует использовать15,что позволяетвосстанавливатьраспределенияс погрешностью(2-5)%.
Восстановлениеглубинныхраспределенийс помощьюаппроксимацийсплайном икусочно-линейнойпрактическиневозможно.Имеют местоосцилляции,приводящиек погрешностям,превышающим10%.
Восстановлениеповерхностныхраспределенийс помощьюаппроксимацийсплайном икусочно-линейнойпрактическиневозможно.Имеют местоосцилляции,приводящиек погрешностям,превышающим10%.
Восстановлениеповерхностныхраспределенийс помощьюаппроксимацийэкспоненциальнойи гиперболическимтангенсомвозможно судовлетворительнойточностью. Дляулучшениярезультатоввосстановленияследует использовать15частот.

10.Заключение

Порезультатампроделаннойработы можносделать следующиевыводы:
Существуетпринципиальнаявозможностьвосстановлениякак поверхностныхтак и глубинныхраспределенийЭП с погрешностью,не превышающей(2-3)%. Длявосстановленияповерхностныхраспределенийследует использоватьэкспоненциальнуюи гиперболическуюаппроксимации,а для глубинныхсплайн икусочно-постоянную(возможноиспользованиеэкспоненциальнойи гиперболическойаппрксимацийдля в приповерхностномслое глубинойпорядка четвертипластины).
Существенноеотрицательноевлияние нарезультатывосстановленияимеют погрешностьизмерения U_вн^*(не следуетиспользоватьданные с погрешностьюизмеренияболее 2%)и малая глубинараспределенияЭП (распределенияЭП сосредоточенныев приповерхностномслое глубиной менее (3-5)%восстанавливаютсяхуже).
Использованиежестких ограниченийна величинуЭП в приповерхностныхслоях оправданопри восстановленииповерхностныхраспределений,причем приналичии данныхс погрешностью,превосходящей2%,или малой глубиныраспределенияпредпочтительнеезадавать ограниченияна обеих поверхностях.При зашумленностиданных порядка(1-2)%достаточнозадать жесткиеограничениялишь на верхнейповерхности.
Внаборе частотвозбужденияВТП должныприсутствоватьнизкочастотныесоставляющие,влияние которыхособенно заметнопри работе сглубиннымираспределениямии соответствующимиаппроксимациями.Рекомендуетсяиспользоватьпорядка десятичастот, равномернораспределенныхпо частотномудиапазону(0.001ё70)КГц.В условияхвысокой погрешностиизмерений илиотчетливовыраженныхприповерхностныхизмененияхЭП заметноеположительноевлияние оказываетувеличениечисла частотвозбужденияВТП (например,до пятнадцати.).

В процессеработы надзадачей былпроведен анализлитературы,выбрана модельзадачи и способыее аппроксимации.При помощипрограммы,разработаннойсогласно предложенноймодели, былипроведенырасчеты модельныхзадач и рассмотренырезультатывосстановленияраспределенийЭП в зависимостиот основныхвлияющих факторов.

Такимобразом, цели,поставленныев техническомзадании, решеныв полном объеме.

11.Литература

Неразрушающийконтроль качестваизделий электромагнитнымиметодами,ГерасимовВГ, 1978,215
Вихретоковыйконтроль накладнымипреобразователями.,ГерасимовВГ,1985,86
Вихретоковыеметоды и приборынеразрушающегоконтроля.,РудаковВН, 1992,72
Накладныеи экранныедатчики.,СоболевВС, 1967,144
Теорияи расчет накладныхвихретоковыхпреобразователей.,ДякинВВ, 1981,135
Основыанализа физическихполей.,ПокровскийАД, 1982,89
Дефектоскопияметаллов.,ДенельАК, 1972,303
Индукционнаяструктуроскопия.,ДорофеевАЛ,1973,177
Структураи свойстваметаллов исплавов.Справочник.,ШматкоОА,1987,580
Некорректныезадачи Численныеметоды и приложения.,ГончарскийАВ,1989,198
Некорректныезадачи матфизикии анализа.,ЛаврентьевММ,1980,286
Линейныеоператоры инекорректныезадачи.,ЛаврентьевММ,1991,331
Методырешения некорректнопоставленныхзадач Алгоритмич.аспект.,МорозовВА, 1992,320
Численныеметоды решениянекорректныхзадач.,ТихоновАН,1990,230
Началатеории вычислительныхметодов,КрыловВИ,1984,260
Математическоепрограммированиев примерах изадачах.,АкуличИЛ,1993,319
Математическоепрограммирование.,КармановВГ,1986,286
Математическоепрограммирование.,ОреховаРА,1992,290
НелинейноепрограммированиеТеория и алгоритмы.,БазараМ,1982,583
Прикладноенелинейноепрограммирование.,ХиммельблауД,1975,534
Введениев методы оптимизации.,АокиМ,1977,344
Введениев оптимизацию.,ПолякБТ,1983,384
Курсметодов оптимизации.,СухаревАГ,1986,326
Практическаяоптимизация.,ГиллФ,1985,509
Численныеметоды оптимизации.,ПолакЭ,1974,367
Алгоритмырешения экстремальныхзадач.,РомановскийИВ,1977,352
Методырешения экстремальныхзадач.,ВасильевФП,1981,400
Методырешения экстремальныхзадач и ихприменениев системахоптимизации.,ЕвтушенкоЮГ, 1982,432
Численныеметоды решенияэкстремальныхзадач.,ВасильевФП,1988,549
Введениев вычислительнуюфизику.,ФедоренкоРП,1994,526
Методыматематическойфизики.,АрсенинВЯ,1984,283
Уравненияматематическойфизики.,ТихоновАН,1977
Уравненияматематическойфизики.,ВладимировВС,1988,512
Методинтегральныхуравнений втеории рассеивания.,КолтонД,1987,311
Теорияэлектромагнитногополя.,ПоливановКМ,1975,207
Eddycurrent testing. Manual on eddy current method.,Cecco VS,1981,195
Optimizationmethods with applications for PC.,Mistree F,1987,168
Electromagneticinverse profiling., TijhuisAG,1987,465
Inverseacoustic and electromagnetic scattering theory.,Colton D,1992,305
"Накладнойэлектромагнитныйпреобразовательнад объектомконтроля сизменяющимисяпо глубинеэлектрическимии магнитнымисвойствами",Касимов ГА,Кулаев ЮВ,"Дефектоскопия",1978, №6, с81-84
"Возможностипримененияметодов теориисинтеза излучающихсистем в задачахэлектромагнитногоконтроля ",Кулаев ЮВ, 1980,тематическийсборник "ТрудыМЭИ",выпуск 453, с12-18
"Analitical solutions to eddy-current probe-coil problems ", Deeds WE, Dodd CV, ІJournalof Applied PhisicsІ,1968,vol39,?3,p2829-2838
"General analysis of probe coils near stratified conductors ", Deeds WE, Dodd CV,ІInternationalJournal of Nondestructive TestingІ,1971,vol3, ?2,p109-130
"Tutorial. A review of least-squares inversion and its application togeophysical problems ",Lines LR, Treitel S,"GeophysicalProspecting ",1984, vol32, ?2,p159-186
"Eddy currentcalculations using half-space Green’s functions ", Bowler JR, ІJournalof Applied PhisicsІ,1987,vol61, ?3,p833-839
"Reconstruction of 3D conductivity variations from eddy current(electromagnetic induction ) data ", Nair SM, Rose JH, ІInverse ProblemsІ,1990, ?6,p1007-1030
"Electromagnetic induction (eddy-currents) in a conducting half-spacein the absence and presence of inhomogeneities: a new formalism ", Nair SM, Rose JH, ІJournalof Applied PhisicsІ,1990,vol68, ?12,p5995-6009
"Eddy-current probe impedance due to a volumetric flaw ", Bowler JR, ІJournalof Applied PhisicsІ,1991,vol70, ?3,p1107-1114
"Theory of eddy current inversion ", Bowler JR, Norton SJ, ІJournalof Applied PhisicsІ,1993,vol73, ?2,p501-512
"Impedance of coils over layered metals with continuously variableconductivity and permeability: Theory and experiment ", Rose JH, ІJournalof Applied PhisicsІ,1993,vol74, ?3,p2076
"Eddy-current interaction with ideal crack ", Bowler JR, ІJournalof Applied PhisicsІ,1994,vol75, ?12,p8128,8138
"Method of solution of forward problems in eddy-current testing ", Kolyshkin AA, ІJournalof Applied PhisicsІ,1995,vol77, ?10,p4903-4912

Приложение1. Программнаяреализация

Программнаяреализацияизложенногометода решенияобратной задачиВТК осуществленапри помощи компилятораBorland Pascal 7.0и состоит изшести модулей:

ErIn12.pas- исполняемыйфайл, осуществляетосновной циклпрограммы
EData.pas- содержитглобальныеданные и осуществляетчтение файлаисходных данных
EFile.pas-содержитвспомогательныефункции ииосуществляетсохранениерезультатоврасчетов
EMath.pas- осуществляетподдержкуопераций скомплекснымичислами
EDirect.pas- осуществляетрешение прямойзадачи ВТК
EMinimum.pas-осуществляетрешение обратнойзадачи ВТК

П1.1Исходные данные

Исходныеданные программыхранятся втекстовомфайле(кодировкаASCII,расширениепо умолчаниюTXT ).

HThick	- толщинапластины,[мм]
nPoints	-количествоузлов аппроксимацииэлектропроводностидля PWL,PWC,SPL аппроксимаций.В случае EXP,HTGаппроксимациивычислениезначений ЭПв них производитсяпо окончаниирасчетов
nLayers	-количество интерваловс кусочно-постояннойэлектропроводностью,на которыеразбиваетсяпластина длянепосредственногорасчета вносимойЭДС по реккурентнымформулам длямногослойнойпластины
nFreqs	-количествочастот возбуждениягармоник вносимойотносительнойЭДС
nStab	- числостабилизируемыхзначащих цифр
epsU	-погрешностьизмерения
aG	-коэффициентсжатия ограничений(aG
nApprox	- типыаппроксимациипрямой и обратнойзадач
si	- значенияпроводимостив узлах аппроксимации
siMin, siMax	-ограниченияна возможныезначенияпроводимостив узлах аппроксимациив процессерешения ОЗ
incVal	- величинаdx для численногодифференцирования
maxSteps	-максимальноечисло отрезковинтегрирования
maxX	- верхнийпредел интегрированияпри расчетеUvn*
Eps	-погрешностьинтегрированияпри расчетеUvn*
dType	- типразностнойпроизводной(=1 правая или=2 центральная)
eqlB	-толщины слоевпластиныодинаковы(b=hThick/nLayers) если eqlB>0,в противномслучае используютсякоординатыслоев из файла

П1.2Используемыеаппроксимации

Примечание.КоординатаХО[0,1]отсчитываетсяот дна пластиныдля всех аппроксимаций.

Сплайн(SPL),кусочно-линейная(PWL),кусочно-постоянная(PWC)аппроксимации.

Впроцессе расчетовищутся значенияэлектропроводностив узлах аппроксимации,причем количествоузлов увеличиваетсяот едениицыдо nPointsвцелях сохраненияустойчивостирешения.

Начальныезначения (узловыезначения ІистиннойІЭП для эмуляцииизмерений U^*_вн)задаются встолбце ІsiІфайла исходныхданных, начальныезначения ограниченийна узловыезначения ЭПв столбцахІsiMinІиІsiMaxІ(движениепо столбцусверху внизсоответствуетизменениюкоординатыот дна пластиныдо обрабатываемойповехности).

Экспоненциальнаяаппроксимация(EXP)

В случаезадания экспоненциальной аппроксимациизависимостьэлектропрводностиот толщиныпредставляетсяв виде

SIGMA = ( siE-siI)*EXP( -alfa*(1-x) ) + siI

Варьруемымипараметрамиявляютсяэектропроводностьна верхнейповерхностиsiЕ,электропроводность"на бесконечности"siIи параметрalfa.В файле исходныхданных в таблицеиз nPointsстрок с подзаголовком"si siMin siMax", информацияоб ограниченияхна параметрыsiE, siI задается впервой и nPoints-строке.Величина иограничениядля параметра alfa задаются первойстрокой в "special approximation parameters".

Аппроксимациягиперболическимтангенсом (HTG)

В случаезадания аппроксимациигиперболическимтангенсомзависимостьэлектропрводностиот толщиныпредставляетсяв виде

SIGMA = si2 + (si1-si2 )/2*{ 1 + th( ( beta-x )/gamma ) }

Величинаи ограничениядля параметровsi2,beta,gammaзадаются начинаясо второй строкив "special approximation parameters", дляsi1аналогичноsiI.

П1.3Результатырасчета

Результатырасчета помещаютсяв текстовыйфайл(кодировкаASCII,расширениепо умолчаниюLST), при этом результаткаждой итерацииотбражаетсястрокой вида:

1 0.000353 Rg= 17.003 15.639 9.697

гдепервая цифра(в данном случае1) соответствуетномеру текущейвнутреннейитерации, затемпосле текста"" идет значениетекущей абсолютнойсреднеквадратичнойневязки по всемгармоникам(в данном случае0.000353), затем после текста "Rg= ", идутискомые текущиезначения переменныхминимизации.В случае SPL,PWL,PWCаппроксимацийэто непосредственноузловые значенияэлектропроводностидля текущегоколичестваузлов, а дляEXP,HTG аппроксимацийэто параметры{ siE,siI, Alfa }или{ si1,si2, Beta, Gamma}.B качестве последнейстроки помещаютсяnPointsвычисленныхзначенийэ/проводностив равномернорасположенныхузлах пластины.

П1.4ОсновнаяпрограммаErIn

(****************************************************************************)

(* ErIn v1.42 *)

(* Eddy current inverse problem solver. *)

(* Moscow PowerEngineering Institute , Introscopy dept. *)

{****************************************************************************}

Program ErIn;{23.02.99}

Uses

DOS,CRT, EData,EMath, EDirect, EFile, EMinimum;

Var

m, mLast, i : byte; {loop counters}

procedure about; {Let me to introduce myself}

begin

clrscr;

GetTime( clk1.H,clk1.M, clk1.S, clk1.S100 ); {get start time}

writeln('***********************************************************');

writeln('* ErInv1.42 Basic * *');

writeln('***********************************************************');

end;

procedureinitParameters;

var

apDT : byte; {approximation type for direct task}

begin

apDT := nApprox SHR4; {XXXXYYYY->0000XXXX}

fHypTg:=(( apDT ANDapHypTg ) = apHypTg);

if fHypTg then

begin

si0[ 1 ]:=si[1 ]; {si1 - conductivity about bottom of slab}

si0[ 2]:=par0[ 2 ]; {si2 - conductivity about top of slab}

si0[ 3]:=par0[ 3 ]; {Beta - ratio of approx.}

si0[ 4]:=par0[ 4 ]; {Gamma- ratio of approx.}

mCur:=4;

end

else

if(( apDT AND apExp) = 0 ) then {It's not an EXP approx.}

begin

for i:=1 tonPoints do si0[ i ] :=si [ i ]; {SI data from file}

mCur:=nPoints;

end

else

begin

si0[ 1 ]:=si[1 ]; {siI - conductivity about bottom of slab}

si0[ 2 ]:=si[nPoints ]; {siE - conductivity about top of slab}

si0[ 3]:=par0[ 1 ]; {Alfa- ratio of approx.}

mCur:=3;

end;

setApproximationType(apDT ); {approx. type for direct problem}

setApproximationData(si0, mCur ); {approx. data for direct problem}

nApprox := (nApprox AND $0F ); {XXXXYYYY->0000YYYY}

fHypTg := ((nApprox AND apHypTg ) = apHypTg );

fMulti := ((nApprox AND apExp ) = 0 ) AND NOT fHypTg; {It's not an EXP approx.}

if fMulti then

begin

for i:=1 tonPoints do

begin

Gr[ 1,i]:=SiMax[ i ];

Gr[ 2,i]:=SiMin[ i ];

Rg[ i]:=( Gr[ 1,i ] + Gr[ 2,i ] )/2; {zero estimate of SI}

Rgs[ i]:=1E33; {biggest integer}

end;

mLast:=nPoints; {loop for every node of approx.}

mCur :=1; {to begin from the only node of approx}

end

else

if fHypTg then

begin

Gr[ 1,1 ]:=siMax[ 1 ]; Gr[ 2,1 ]:= siMin[ 1 ]; Rgs[ 1 ]:=1E33;

Gr[ 1,2]:=parMax[ 2 ]; Gr[ 2,2 ]:=parMin[ 2 ]; Rgs[ 2 ]:=1E33;

Gr[ 1,3]:=parMax[ 3 ]; Gr[ 2,3 ]:=parMin[ 3 ]; Rgs[ 3 ]:=1E33;

Gr[ 1,4]:=parMax[ 4 ]; Gr[ 2,4 ]:=parMin[ 4 ]; Rgs[ 4 ]:=1E33;

for i:=1 to 4do Rg[ i ]:=( Gr[ 1,i ] + Gr[ 2,i ] )/2;

mLast:=1;

mCur:=4;

end

else

begin

Gr[1,1 ]:= siMax[1]; Gr[2,1]:= siMin[1]; Rgs[ 1 ]:=1E33;

Gr[ 1,2 ]:=siMax[nPoints]; Gr[2,2]:= siMin[nPoints]; Rgs[ 2 ]:=1E33;

Gr[1,3 ]:=parMax[1]; Gr[2,3]:=parMin[1]; Rgs[ 3]:=1E33;

for i:=1 to 3do Rg[ i ]:=( Gr[ 1,i ] + Gr[ 2,i ] )/2;

mLast:=1;

mCur :=3;

end;

initConst( nLayers,parMaxH, parMaxX , parEps, parEqlB );{set probe params}

end;

procedure directTask; {emulate voltage measurements [with error]}

begin

for i:=1 to nFreqsdo

begin

getVoltage(freqs[i], Umr[ i ], Umi[ i ] ); {"measured" Uvn*}

if ( epsU >0 ) then {add measurement error}

begin

randomize;Umr[ i ]:=Umr[ i ]*( 1 + epsU*( random-0.5 ) );

randomize;Umi[ i ]:=Umi[ i ]*( 1 + epsU*( random-0.5 ) );

end;

writeln('* Voltagemeasurements have been emulated');

setApproximationType(nApprox ); {approx. type for inverse problem}

setApproximationData(Rg, mCur ); {approx. data for inverse problem}

end;

procedurereduceSILimits; {evaluate SI for m+1 points of approx. using aG}

var

x0, x1, xL, dx,Gr1, Gr2 : real;

j, k : byte;

begin

{-----------------------------get SI min/max for m+1 points of approximation}

dx:=1/(nPoints-1 );

for i:=1 to m+1do

begin

k:=1;

x1:=0;

x0:=( i-1)/m;

for j:=1 tonPoints-1 do

begin

xL:=(j-1 )/( nPoints-1 );

if( (xL

begin

k:=j;

x1:=xL;

end;

Gr[ 1,i]:=siMax[ k ] + ( siMax[ k+1 ]-siMax[ k ] )*( x0-x1 )/dx;

Gr[ 2,i]:=siMin[ k ] + ( siMin[ k+1 ]-siMin[ k ] )*( x0-x1 )/dx;

end;

{-------------------------------------get SI for m+1 points of approximation}

for i:=1 to m+1do

begin

Rg[i]:=getSiFunction((i-1)/m );

if ( Rg[i]> Gr[1,i] )then Rg[i]:=Gr[1,i];

if ( Rg[i]

if m > 1then {There're more than 1 point of approx.}

begin

Gr1:=Rg[i]+( Gr[1,i]-Rg[i] )*aG; {reduce upper bound}

Gr2:=Rg[i]-( Rg[i]-Gr[2,i] )*aG; {reduce lower bound}

if (Gr1

if (Gr2 > Gr[2,i] )then Gr[2,i]:=Gr2;

end;

setApproximationData(Rg , m+1 );

end;

procedure resultMessage; {to announce new results}

begin

if fMulti then

begin

writeln('current nodal values of conductivity');

write(' si :');for i:=1 to m do write(Rg[i] :6:3,' ');writeln;

write(' max:');for i:=1 to m do write(Gr[1,i]:6:3,' ');writeln;

write(' min:');for i:=1 to m do write(Gr[2,i]:6:3,' ');writeln;

end

else

begin

for i:=1 tonPoints do si[i]:=getSiFunction( ( i-1 )/( nPoints-1 ) );

if fHypTg then

saveHypTgResults

else

saveExpResults;

end;

procedureclockMessage; {user-friendlymessage}

begin

writeln('***********************************************************');

write( '* approximation points number :',m:3,' * Time '); clock;

writeln('***********************************************************');

end;

procedure done; {final message}

begin

Sound(222);Delay(111); Sound(444); Delay(111); NoSound; {beep}

write('* Taskprocessing time '); clock; saveTime;

writeln('* Status:Inverse problem has been successfully evaluated.');

end;

Begin

about;

loadData;

initParameters;

directTask;

for m:=1 to mLastdo

begin

if fMulti then

begin

mCur:=m;

clockMessage;

end;

doMinimization; {main part of work}

setApproximationData(Rg, mCur ); {set new approx. data}

resultMessage;

if(( fMulti)AND( m

end;

done;

End.

П1.5МодульглобальныхданныхEData

Unit EData;

Interface

Uses DOS;

Const

maxPAR = 40;{nodes of approximation max number}

maxFUN = 20;{excitation frequencies max number}

maxSPC = 4;{support approximation values number}

iterImax = 50; {maxnumber of internal iterations}

Const

apSpline = 1;{approximation type identifiers}

apHypTg = 3;

apExp = 2;

apPWCon = 4;

apPWLin = 8;

Type

Parameters =array[ 1..maxPAR ] of real; {Si,Mu data}

Functionals =array[ 1..maxFUN ] of real; {Voltage data}

SpecialPar =array[ 1..maxSPC ] of real; {Special data}

Var

hThick : real; {thickness of slab}

nPoints : integer;{nodes of approximation number within [ 0,H ]}

nLayers : integer;{number of piecewise constant SI layers within[ 0,H ]}

nFreqs : integer;{number of excitation frequencies}

nStab : integer;{required number of true digits in SI estimation}

epsU : real; {relative error of measurement Uvn*}

nApprox : byte; {approximation type identifier}

incVal : real; {step for numerical differ.}

parMaxH : integer;{max number of integration steps}

parMaxX : real; {upper bound of integration}

parEps : real; {error of integration}

derivType: byte; {1 for right; 2 for central}

Var

freqs :Functionals; {frequencies of excitment Uvn*}

Umr, Umi :Functionals; { Re(Uvn*),Im(Uvn*) for "measured" data}

Uer, Uei :Functionals; { Re(Uvn*),Im(Uvn*) for estimated data}

mu :Parameters; {relative permeability nodal values}

si, si0 :Parameters; {conductivity approximation nodal values}

siMin, siMax :Parameters; {conductivity nodal values restrictions}

par0 :SpecialPar; {alfa,si2,beta,gamma - for exp&HypTg}

parMin,parMax:SpecialPar; {-||- min/max}

zLayer :Parameters; {relative borders of slab layers [0,1]}

Var

aG :real; {scale factor for SImin/max}

Ft :real; {current discrepancy functional value}

fMulti :boolean; {TRUE if it isn't an EXP-approximation}

fHypTg :boolean; {TRUE for Hyperbolic tg approximation}

parEqlB :boolean; {TRUE if b[i]=const}

mCur :integer; {current number of approximation nodes}

inFileName :string; {data file name}

outFileName :string; {results file name}

Var

Rg : Parameters; {current SI estimation}

RgS : Parameters; {previous SI estimation}

Gr : array [ 1..2,1..maxPAR ] of real; {SI max/min}

Fh : array [1..maxPAR , 1..maxFUN ] of real; {current discrepancies}

Type

TTime=record

H, M,S, S100 : word; {hour,min,sec,sec/100}

end;

Var

clk1, clk2 : TTime; {start&finish time}

procedure loadData; {load all user-defined data from file}

Implementation

procedure loadData;

var

i,eqlB : integer;

FF : text;

begin

assign( FF,outFileName ); {clear output file}

rewrite( FF );

close( FF );

assign( FF,inFileName ); {read input file}

reset( FF );

readln( FF );

readln( FF, hThick,nPoints, nLayers, nFreqs, nStab, epsU, aG, nApprox );

readln( FF );

for i:=1 to nFreqsdo read( FF, freqs[i] );

readln( FF );

for i:=1 to nPointsdo readln( FF, si[i], siMin[i], siMax[i] );

readln( FF );

readln( FF ,incVal, parMaxH, parMaxX, parEps, derivType, eqlB );

readln( FF );

for i:=1 to maxSPCdo readln( FF, par0[i] , parMin[i] , parMax[i] );

readln( FF );

if ( eqlB=0 )then

begin

for i:=1 tonLayers+1 do read( FF, zLayer[i] );

parEqlB:=false;

end

else parEqlB:=true;

close( FF );

for i:=1 to maxPARdo mu[i]:=1;

end;

Var

str : string;

Begin

if( ParamCount = 1)then str:=ParamStr(1)

else

begin

write('EnterI/O file name, please: ');

readln( str );

end;

inFileName:=str+'.txt';

outFileName:=str+'.lst';

End.

П1.6Модульработы с файламиEFile

Unit EFile;

Interface

Uses

DOS, EData;

function isStable( ns :integer; var RG1,RG2 ) : boolean;

function saveResults(ns,iter : integer ) : boolean;

proceduresaveExpResults;

proceduresaveHypTgResults;

procedure clock;

procedure saveTime;

Implementation

Var

FF : text;

i : byte;

function decimalDegree(n:integer ) : real;{10^n}

var

s:real; i:byte;

begin

s:=1;

for i:=1 to n dos:=s*10;

decimalDegree:=s;

end;

function isStable(ns:integer ; var RG1,RG2 ) : boolean;

var

m : real;

R1 : Parametersabsolute RG1;

R2 : Parametersabsolute RG2;

begin

isStable:=TRUE;

m:=decimalDegree(ns-1 );

for i:=1 to mCur do

begin

if NOT(( ABS(R2[i]-R1[i] )*m )

RgS[i]:=Rg[i];

end;

function saveResults( ns, iter : integer ) : boolean;

var

sum : real;

begin

sum:=0;

for i:=1 to nFreqsdo sum:=sum + Fh[1,i];

sum:=SQRT(sum/nFreqs );

assign( FF ,outFileName );

append( FF );

write( iter:2,' ', sum:10:7, ' Rg=' );

write( FF , iter:2,' ', sum:10:7, ' Rg=');

for i:=1 to mCur do

begin

write( Rg[i]:6:3, ' ');

write( FF ,Rg[i]:6:3, ' ');

end;

writeln;

writeln( FF );

close( FF );

saveResults:=isStable(ns , Rgs , Rg );

end;

proceduresaveExpResults;

begin

assign( FF ,outFileName );

append( FF );

writeln( 'siE=',Rg[2]:6:3,' siI=',Rg[1]:6:3,' alfa=',Rg[3]:6:3);

writeln( FF , 'siE=',Rg[2]:6:3,' siI=',Rg[1]:6:3,' alfa=',Rg[3]:6:3);

write( ' SI:');

write( FF , ' SI:');

for i:=1 to nPointsdo

begin

write( si[i]:6:3,' ');

write( FF ,si[i]:6:3,' ');

end;

writeln;

writeln( FF );

close( FF );

end;

proceduresaveHypTgResults;

begin

assign( FF ,outFileName );

append( FF );

writeln( 'si1=',Rg[2]:6:3,' si2=',Rg[1]:6:3,' beta=',Rg[3]:6:3,'gamma=',Rg[4]:6:3);

writeln( FF , 'si1=',Rg[2]:6:3,' si2=',Rg[1]:6:3,' beta=',Rg[3]:6:3,'gamma=',Rg[4]:6:3);

write( ' SI:');

write( FF , ' SI:');

for i:=1 to nPointsdo

begin

write( si[i]:6:3,' ');

write( FF ,si[i]:6:3,' ');

end;

writeln;

writeln( FF );

close( FF );

end;

procedure clock; {t2 = t2-t1}

var

H1,M1,S1,H2,M2,S2,sec1,sec2: longint;

begin

GetTime( clk2.H,clk2.M, clk2.S, clk2.S100 ); {current time}

H2:=clk2.H;M2:=clk2.M; S2:=clk2.S; H1:=clk1.H; M1:=clk1.M; S1:=clk1.S;

sec2:= ( H2*60 + M2)*60 + S2;

sec1:= ( H1*60 + M1)*60 + S1;

if( sec2

sec2:=sec2 - sec1;

clk2.H := sec2 div3600; sec2:=sec2 - clk2.H*3600;

clk2.M := sec2 div60; sec2:=sec2 - clk2.M*60;

clk2.S := sec2;

writeln( clk2.H:2,':', clk2.M:2, ':', clk2.S:2 );

end;

procedure saveTime;

begin

assign( FF ,outFileName );

append( FF );

write( FF ,'*Processing time ',clk2.H:2, ':', clk2.M:2, ':', clk2.S:2 );

close( FF );

end;

End.

П1.7Модуль решенияпрямой задачиВТК для НВТПEDirect

{****************************************************************************}

{ Estimates Uvn* forEddy current testing of inhomogeneous multilayer slab }

{ with surface( flat )probe. }

{ It can do it usingone of five types of conductivity approximation : }

{Spline, Exponential,Piecewise constant, Piecewise linear,Hyperbolic tangent}

{****************************************************************************}

{$F+}

Unit EDirect;

Interface

Uses EData, EMath;

Type

siFunc =function( x:real ) : real;

Var

getSiFunction :siFunc; {for external getting SI estimate}

procedure initConst(par1,par2:integer; par3,par4:real; par5:boolean );

procedure getVoltage(freq : real ; var ur,ui : real ); { Uvn* = ur + j*ui }

proceduresetApproximationType( approx : byte ); { type of approx. }

proceduresetApproximationItem( SIG:real ; N : byte ); { set SIGMA[ N ]}

proceduresetApproximationData( var SIG; nVal : byte ); { SIGMA[1..nVal] }

proceduregetApproximationData( var SIG ; var N : byte ); { get SIGMA[ N ]}

Implementation

Const

PI23 = 2000*pi; {2*pi*KHz}

mu0 = 4*pi*1E-7; {magnetic const}

Var

appSigma :Parameters; {conductivity approximation data buffer}

appCount : byte; {size of conductivity approximation data buffer}

appType : byte; {conductivity approximation type identifier}

Type

commonInfo=record

w : real; {cyclical excitation frecuency}

R : real; {equivalent radius of probe}

H : real; {generalized lift-off of probe}

Kr : real; {parameter of probe}

eps : real; {error of integration}

xMax : real; {upper bound of integration}

steps : integer; {current number of integration steps}

maxsteps:integer; {max number of integration steps}

Nlay : integer; {number of layers in slab}

sigma : Parameters; {conductivity of layers}

m : Parameters; {relative permeability of layers}

b : Parameters; {thickness of layers}

zCentre: Parameters; {centre of layer}

end;

procFunc =procedure( x:real; var result:complex);

Var

siB, siC, siD :Parameters; {support for Spline approx.}

cInfo :commonInfo; {one-way access low level info}

function siSpline(x:real ) : real;{Spline approximation}

begin

if( appCount = 1)then

siSpline :=appSigma[ 1 ]

else

siSpline:=Seval(appCount, x, appSigma, siB, siC, siD);

end;

function siExp( x:real ): real;{Exponential approximation}

begin

siExp:=(appSigma[2]-appSigma[1])*EXP(-appSigma[3]*(1-x) ) + appSigma[1];

end;

function siPWConst(x:real ) : real;{Piecewise constant approximation}

var

dx, dh : real; i :byte;

begin

if( appCount = 1)then siPWConst := appSigma[ 1 ]

else

begin

dh:=1/(appCount-1 );

dx:=dh/2;

i:=1;

while( x >dx ) do

begin

i:=i + 1;

dx:=dx +dh;

end;

siPWConst:=appSigma[i ];

end;

function siPWLinear(x:real ) : real;{Piecewise linear approximation}

var

dx, dh : real;

i : byte;

begin

if( appCount = 1)then siPWLinear := appSigma[ 1 ]

else

begin

dh:=1/(appCount-1 );

dx:=0;

i:=1;

repeat

i:=i +1;

dx:=dx +dh;

until( x

siPWLinear:=appSigma[i-1]+(appSigma[i]-appSigma[i-1] )*( x/dh+2-i);

end;

function siHyperTg(x:real ) : real;{Hyperbolic tangent approximation}

begin

siHyperTg:=appSigma[2]+(appSigma[1]-appSigma[2])*(1+th((appSigma[3]-x)/appSigma[4]))/2;

end;

proceduresetApproximationType( approx : byte );

begin

appType := approx;

write('* conductivity approximation type : ');

case approx of

apSpline: begin

writeln('SPLINE');

getSiFunction:= siSpline;

end;

apExp : begin

writeln('EXP');

getSiFunction:= siExp;

end;

apPWCon : begin

writeln('PIECEWISECONST');

getSiFunction:= siPWConst;

end;

apPWLin : begin

writeln('PIECEWISELINEAR');

getSiFunction:= siPWLinear;

end;

apHypTg : begin

writeln('HYPERBOLICTANGENT');

getSiFunction:= siHyperTg;

end;

proceduresetApproximationData( var SIG ; nVal : byte );

var

Sigma : Parametersabsolute SIG; i:byte;

begin

appCount := nVal;

for i:=1 to nVal doappSigma[ i ]:=Sigma[ i ];

if( appType =apSpline )then Spline( appCount, appSigma, siB, siC, siD);

end;

proceduregetApproximationData( var SIG ; var N : byte );

var

Sigma : Parametersabsolute SIG; i:byte;

begin

N := appCount;

for i:=1 toappCount do Sigma[ i ]:=appSigma[ i ];

end;

proceduresetApproximationItem( SIG:real ; N : byte );

begin

appSigma[ N ] :=SIG;

if( appType =apSpline )then Spline( appCount, appSigma, siB, siC, siD);

end;

procedure functionFi(x:real ; var result:complex );{get boundary conditions functionvalue}

var

beta : array[1..maxPAR ]of real;

q : array[1..maxPAR ]of complex;

fi : array[0..maxPAR ]of complex;

th , z1 , z2 , z3 ,z4 , z5 , z6 , z7 : complex;

i : byte;

begin

mkComp( 0, 0, fi[0]);

with cInfo do

for i:=1 to Nlay do

begin

beta[i]:=R*sqrt(w*mu0*sigma[i] ); {calculation of beta}

mkComp(sqr(x), sqr(beta[i])*m[i], z7 ); {calculation of q, z7=q^2}

SqrtC( z7,q[i] );

mulComp( q[i],b[i], z6 ); {calculation of th,z6=q*b}

tanH( z6, th); {th=tanH(q*b)}

mkComp(sqr(m[i])*sqr(x), 0, z6 ); {z6=m2x2}

SubC( z6, z7,z5); {z5=m2x2-q2}

AddC( z6, z7,z4); {z4=m2x2+q2}

MulC( z5, th,z1); {z1=z5*th}

MulC( z4, th,z2); {z2=z4*th}

mulComp( q[i],2*x*m[i], z3 ); {z3=2xmq}

SubC( z2, z3,z4 );

MulC( z4,fi[i-1], z5 );

SubC( z1, z5,z6 ); {z6=high}

AddC( z2, z3,z4 );

MulC( z1,fi[i-1], z5 );

SubC( z4, z5,th ); {th=low}

DivC( z6, th,fi[i] );

end;

eqlComp( result,fi[ cInfo.Nlay ] );

end;

procedure funcSimple(x:real; var result:complex );{intergrand function value}

var

z : complex;

begin

with cInfo do

begin

functionFi( x,result );

mulComp(result, exp( -x*H ), z );

mulComp( z,J1( x*Kr ), result );

mulComp(result, J1( x/Kr ), z );

eqlComp(result, z );

end;

procedurefuncMax( x:real; var result:complex );{max value;When Fi[Nlay]=1}

var

z1, z2 : complex;

begin

with cInfo do

begin

mkComp( 1,0,z1);

mulComp( z1,exp(-x*H), z2 );

mulComp( z2,J1( x*Kr ), z1 );

mulComp( z1,J1( x/Kr ), result );

end;

procedure integralBS(func:procFunc ; var result:complex );{integral by Simpson}

var

z , y , tmp :complex;

hh : real;

i, iLast : word;

begin

with cInfo do

begin

hh:=xMax/steps;

iLast:=stepsdiv 2;

nullComp(tmp);

func( 0, z );

eqlComp( y, z);

for i:=1 toiLast do

begin

func( (2*i-1 )*hh , z );

deltaComp(tmp, z );

end;

mulComp( tmp,4, z );

deltaComp( y,z );

nullComp( tmp);

iLast:=iLast-1;

for i:=1 toiLast do

begin

func(2*i*hh ,z );

deltaComp(tmp, z );

end;

mulComp( tmp,2, z );

deltaComp( y,z );

func( xmax, z);

deltaComp(y,z);

mulComp( y,hh/3, z );

eqlComp(result, z );

end;

end;{I = h/3 * [ F0 +4*sum(F2k-1) + 2*sum(F2k) + Fn ]}

procedureintegral( F:procFunc; var result:complex );{integral with givenerror}

var

e , e15 : real;

flag : boolean;

delta , integ1H ,integ2H : complex;

begin

with cInfo do

begin

e15 :=eps*15;{ | I2h-I1h |/(2^k -1 )

steps:=20;

flag :=false;

integralBS( F,integ2H );

repeat

eqlComp(integ1H, integ2H );

steps:=steps*2;

integralBS(F, integ2H );

SubC(integ2H, integ1H, delta );

e:=Leng(delta );

if(e

until( ( flag) OR (steps>maxsteps) );

if( flag )then

begin

eqlComp(result, integ2H );

end

else

begin

writeln('Error:Too big number of integration steps.');

halt(1);

end;

procedure initConst(par1, par2 : integer; par3, par4 : real; par5:boolean );

var

i : byte;

bThick, dl, x :real;

const

Ri=0.02; hi=0.005; { radius and lift-off of excitation coil}

Rm=0.02; hm=0.005; { radius and lift-off of measuring coil}

begin

with cInfo do

begin

Nlay :=par1;

xMax :=par3;

maxsteps:=par2;

R :=sqrt( Ri*Rm );

H :=(hi+hm )/R;

Kr :=sqrt( Ri/Rm );

eps :=par4;

bThick :=hThick*0.002/R; {2*b/R [m]}

for i:=1 toNlay do m[i]:= mu[i];

if par5 then

begin

bThick:=bThick/NLay;

for i:=1to Nlay do b[i]:=bThick;

dl:=1/NLay;

x:=dl/2; {x grows up from bottom of slab to the top}

for i:=1to Nlay do

begin

zCentre[i]:=x;

x:=x+ dl;

end;

end

else

for i:=1to Nlay do

begin

b[i]:=(zLayer[i+1]-zLayer[i] )*bThick;

zCentre[i]:=(zLayer[i+1]+zLayer[i] )/2;

end;

procedure init( f:real);{get current approach of conductivity values}

var

i : byte;

begin

with cInfo do

begin

w:=PI23*f;

for i:=1 toNlay do sigma[i]:=getSiFunction( zCentre[i] )*1E6;

end;

procedure getVoltage(freq : real ; var ur,ui : real );

var

U, U0, Uvn, tmp :complex;

begin

init( freq );

integral(funcSimple, U ); { U =Uvn }

integral( funcMax , U0 ); { U0=Uvn max }

divComp( U,Leng(U0), Uvn ); { Uvn=U/|U0| }

mkComp( 0, 1, tmp); { tmp=( 0+j1 ) }

MulC( tmp, Uvn, U); { U= j*Uvn = Uvn* }

ur := U.re;

ui := U.im;

end;

END.

П1.8Модульрешения обратнойзадачи ВТК дляНВТПEMinimum

Unit EMinimum;

INTERFACE

Uses EData, Crt, EFile,EDirect;

proceduredoMinimization;

IMPLEMENTATION

procedure getFunctional(Reg : byte );

var

ur, ui, dur, dui,Rgt : real;

ur2, ui2:Functionals;

i, j, k : byte;

begin

getApproximationData(si , k );

setApproximationData(Rg, mCur );

case Reg of

0 : for i:=1to nFreqs do {get functional F value}

begin

getVoltage(freqs[i], ur, ui );

Uer[i ]:=ur; {we need it for dU}

Uei[i ]:=ui;

Fh[1,i]:= SQR( ur-Umr[i] ) + SQR( ui-Umi[i] );

end;

{Right:U'(i)=(U(i+1)-U(i))/h}

1 : for i:=1to mCur do {get dF/dSI[i] value}

begin

Rgt:=Rg[i]*(1+incVal ); {si[i]=si[i]+dsi[i]}

setApproximationItem(Rgt, i ); {set new si[i] value}

forj:=1 to nFreqs do

begin {get dUr/dSI,dUi/dSI}

getVoltage(freqs[ j ], ur, ui );

dur:=(ur-Uer[j] )/( Rg[i]*incVal );

dui:=(ui-Uei[j] )/( Rg[i]*incVal );

Fh[i,j]:=2*(dur*(Uer[j]-Umr[j])+dui*(Uei[j]-Umi[j]));

end;

setApproximationItem(Rg[i], i ); {restore si[i] value}

end;

{Central:U'(i)=(U(i+1)-U(i-1))/2h}

2 : for i:=1to mCur do {get dF/dSI[i] value}

begin

Rgt:=Rg[i]*(1+incVal ); {si[i]=si[i]+dsi[i]}

setApproximationItem(Rgt, i ); {set new si[i] value}

forj:=1 to nFreqs do getVoltage( freqs[j],ur2[j],ui2[j] );

Rgt:=Rg[i]*(1-incVal ); {si[i]=si[i]-dsi[i]}

setApproximationItem(Rgt, i ); {set new si[i] value}

forj:=1 to nFreqs do

begin {get dUr/dSI,dUi/dSI}

getVoltage(freqs[ j ], ur, ui );

dur:=(ur2[j]-ur )/( 2*Rg[i]*incVal );

dui:=(ui2[j]-ui )/( 2*Rg[i]*incVal );

Fh[i,j]:=2*(dur*(Uer[j]-Umr[j])+dui*(Uei[j]-Umi[j]));

end;

setApproximationItem(Rg[i], i ); {restore si[i] value}

end;

setApproximationData(si , k );

end;

proceduredoMinimization;

const

mp1Max = maxPAR +1;

mp2Max = maxPAR +2;

m2Max = 2*( maxPAR+ maxFUN );

m21Max = m2Max + 1;

n2Max = 2*maxFUN;

m1Max = maxPAR +n2Max;

n1Max = n2Max + 1;

mm1Max = maxPAR +n1Max;

minDh : real =0.001; {criterion of an exit from golden section method}

var

A : array [ 1 ..m1Max , 1 .. m21Max ] of real;

B : array [ 1 ..m1Max] of real;

Sx: array [ 1 ..m21Max] of real;

Zt: array [ 1 ..maxPAR] of real;

Nb: array [ 1 ..m1Max] of integer;

N0: array [ 1 ..m21Max] of integer;

a1, a2, dh, r, tt,tp, tl, cv, cv1, cl, cp : real;

n2, n1, mp1, mp2,mm1, m1, m2, m21 : integer;

ain : real;

apn : real;

iq : integer;

k1 : integer;

n11 : integer;

ip : integer;

iterI : integer;

i,j,k : integer;

label

102 ,103 ,104 ,105,106 ,107 ,108;

begin

n2:=2*nFreqs;n1:=n2+1; m1:=mCur+n2;

mp1:=mCur+1;mp2:=mCur+2; mm1:=mCur+n1;

m2:=2*( mCur+nFreqs); m21:=m2+1;

for k:=1 to m1Maxdo

for i:=1 tom21Max do

A[k,i]:=0;

iterI:=0;

102:

iterI:=iterI+1;

getFunctional( 0 );

for i:=1 to nFreqsdo b[i]:= -Fh[1,i];

getFunctional(derivType );

for k:=1 to mCur do

begin

Zt[k]:=Rg[k];

for i:=1 tonFreqs do

begin

A[i,k+1]:=Fh[k,i];

A[i+nFreqs,k+1]:=-A[i,k+1];

end;

for i:=1 tonFreqs do B[i]:=B[i]+Rg[k]*A[i,k+1];

end;

for i:=1 to nFreqsdo B[i+nFreqs]:=-B[i];

for i:=n1 to m1 doB[i]:=Gr[1,i-n2]-Gr[2,i-n2];

for i:=1 to m1 do

begin

if i

for k:=2 tomp1 do B[i]:=B[i]-A[i,k]*Gr[2,k-1];

A[i,1]:=-1;

if i>n2then

begin

A[i,1]:=0;

for k:=2to mp1 do

ifi-n2=k-1 then A[i,k]:=1

elseA[i,k]:=0;

end;

for k:=mp2 tom21 do

if k-mp1=ithen A[i,k]:=1

elseA[i,k]:=0;

end;

k1:=1;

for k:=1 to n2 do

if B[k1]>B[k]then k1:=k;

for k:=1 to mp1 doA[k1,k]:=-A[k1,k];

A[k1,mCur+1+k1]:=0;

B[k1]:=-B[k1];

for i:=1 to n2 do

if ik1then

begin

B[i]:=B[i]+B[k1];

for k:=1to mm1 do A[i,k]:=A[i,k]+A[k1,k];

end;

fori:=mp2 to m21 do

begin

Sx[i]:=B[i-mp1];

Nb[i-mp1]:=i;

end;

for i:=1 to mp1 doSx[i]:=0;

Sx[1]:=B[k1];

Sx[mp1+k1]:=0;

Nb[k1]:=1;

103:

for i:=2 to m21 doN0[i]:=0;

104:

fori:=m21 downto 2 do

if N0[i]=0 thenn11:=i;

for k:=2 to m21 do

if((A[k1,n11]N0[k])) then n11:=k;

if A[k1,n11]

iq:=0;

for i:=1 to m1 do

if ik1then

begin

ifA[i,n11]>0 then

begin

iq:=iq+1;

ifiq=1 then

begin

Sx[n11]:=B[i]/A[i,n11];ip:=i;

end

else

begin

ifSx[n11]>B[i]/A[i,n11] then

begin

Sx[n11]:=B[i]/A[i,n11];ip:=i;

end;

end

else

ifiq=0 then

begin

N0[n11]:=n11;

goto104;

end;

Sx[Nb[ip]]:=0;

Nb[ip]:=n11;

B[ip]:=B[ip]/A[ip,n11];

apn:=A[ip,n11];

for k:=2 to m21 doA[ip,k]:=A[ip,k]/apn;

for i:=1 to m1 do

if iipthen

begin

ain:=A[i,n11];

B[i]:=-B[ip]*ain+B[i];

for j:=1to m21 do A[i,j]:=-ain*A[ip,j]+A[i,j];

end;

for i:=1 to m1 doSx[Nb[i]]:=B[i];

goto 103;

105:

for k:=1 to mCur doSx[k+1]:=Sx[k+1]+Gr[2,k];

a1:=0;

a2:=1.;

dh:=a2-a1;

r:=0.618033;

tl:=a1+r*r*dh;

tp:=a1+r*dh;

j:=1;

108:

if j=1 then tt:=tlelse tt:=tp;

106:

for i:=1 to mCur doRg[i]:=Zt[i]+tt*(Sx[i+1]-Zt[i]);

getFunctional( 0 );

cv:=abs(Fh[1,1]);

if nFreqs>1 then

for k:=2 tonFreqs do

begin

cv1:=abs(Fh[1,k]);

if cv

end;

if (j=1) or (j=3)then cl:=cv

else cp:=cv;

if j=1 then

begin

j:=2;

goto 108;

end;

if dh

if cl>cp then

begin

a1:=tl;dh:=a2-a1; tl:=tp; tp:=a1+r*dh ; tt:=tp; cl:=cp; j:=4;

end

else

begin

a2:=tp;dh:=tp-a1; tp:=tl; tl:=a1+r*r*dh; tt:=tl; cp:=cl; j:=3;

end;

goto 106;

107:

if (iterI

end;

End.

Приложение2- Удельнаяэлектрическаяпроводимостьматериалов

Приведемсводку справочныхданных согласно[7-9].

Материал	s_min,[МСм/м]	s_max,[МСм/м]
Немагнитныестали	0.4	1.8
Бронзы(БрБ, Бр2, Бр9)	6.8	17
Латуни(ЛС59, ЛС62)	13.5	17.8
Магниевыесплавы (МЛ5-МЛ15)	5.8	18.5
Титановыесплавы (ОТ4,ВТ3-ВТ16)	0.48	2.15
Алюминиевыесплавы (В95, Д16,Д19)	15.1	26.9

Приложение4- Abstract

The inverseeddy current problem can be described as the task of reconstructingan unknown distribution of electrical conductivity from eddy-currentprobe voltage measurements recorded as function of excitationfrequency. Conductivity variation may be a result of surfaceprocessing with substances like hydrogen and carbon or surfaceheating.

Mathematical reasons andsupporting software for inverse conductivity profiling were developedby us. Inverse problem was solved for layered plane and cylindricalconductors.

Because the inverse problem isnonlinear, we propose using an iterative algorithm which can beformalized as the minimization of an error functional related to thedifference between the probe voltages theoretically predicted by thedirect problem solving and the measured probe voltages.

Numerical results were obtainedfor some models of conductivity distribution. It was shown thatinverse problem can be solved exactly in case of correctmeasurements. Good estimation of the true conductivity distributiontakes place also for measurement noise about 2 percents but in caseof 5 percent error results are worse.